Part A: the map f, the pull-back, curvature, flux and holonomy (quadrature)
  ok   f commutes with the Z^2 translations (eps=0.3, c=1.7)
  ok   det df = 1 (orientation and area preserving)
  ok   f^*alpha = eps cos(2pi y) dx + eps c cos^2(2pi y) dy
  ok   harmonic parts: E_eps (0,0), f^*E_eps (0, eps c/2) = (0, 0.255)
  ok   sup |F| = 2 pi |eps| = 1.88496
  ok   integral of alpha over f(y-circle) = eps c/2; over the y-circle = 0
  ok   curvature flux through the homotopy = eps c/2
  ok   f commutes with the Z^2 translations (eps=0.05, c=-2.0)
  ok   det df = 1 (orientation and area preserving)
  ok   f^*alpha = eps cos(2pi y) dx + eps c cos^2(2pi y) dy
  ok   harmonic parts: E_eps (0,0), f^*E_eps (0, eps c/2) = (0, -0.05)
  ok   sup |F| = 2 pi |eps| = 0.314159
  ok   integral of alpha over f(y-circle) = eps c/2; over the y-circle = 0
  ok   curvature flux through the homotopy = eps c/2
Part B: twisted signature operator on T^2 (Fourier-Galerkin, x-modes |m|<=3, y-modes |j|<=NY)
  ok   tau^2 = 1 on Lambda^*(R^2) (x) C
  -- eps=0.05, c=1.0: sup|F| = 0.3142, r = |eps| sqrt(1+c^2) = 0.07071, d = 0.025
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [7.867312e-18 1.009734e-15 1.800853e-15 6.650457e-15]; 5th = 6.2828; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [2.425294e-15 4.390804e-15 1.267822e-14 2.603299e-14]; 5th = 6.2828; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.024998 0.024998 0.024998 0.024998]; 5th = 6.2583; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.024998 0.024998 0.024998 0.024998]; 5th = 6.2583; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.02499842 in [e^(-|eps|/pi) d, d] = [0.02460526, 0.02500000]
  ok   both operators: exactly four eigenvalues in [-r, r]
  -- eps=0.01, c=1.0: sup|F| = 0.06283, r = |eps| sqrt(1+c^2) = 0.01414, d = 0.005
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [1.908350e-15 2.069930e-15 3.418511e-15 6.755441e-15]; 5th = 6.2832; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [9.851450e-18 3.590487e-15 7.441698e-15 1.881201e-14]; 5th = 6.2832; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.005 0.005 0.005 0.005]; 5th = 6.2782; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.005 0.005 0.005 0.005]; 5th = 6.2782; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.00499999 in [e^(-|eps|/pi) d, d] = [0.00498411, 0.00500000]
  ok   both operators: exactly four eigenvalues in [-r, r]
  -- eps=0.2, c=0.7: sup|F| = 1.257, r = |eps| sqrt(1+c^2) = 0.2441, d = 0.07
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [3.130538e-29 1.503852e-16 4.209791e-15 1.437528e-14]; 5th = 6.2768; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [1.055128e-25 1.530774e-15 4.735360e-15 1.821370e-14]; 5th = 6.2768; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.069929 0.069929 0.069929 0.069929]; 5th = 6.2153; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.069929 0.069929 0.069929 0.069929]; 5th = 6.2153; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.06992909 in [e^(-|eps|/pi) d, d] = [0.06568255, 0.07000000]
  ok   both operators: exactly four eigenvalues in [-r, r]
  -- eps=0.5, c=2.0: sup|F| = 3.142, r = |eps| sqrt(1+c^2) = 1.118, d = 0.5
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [1.552022e-18 6.896181e-16 9.862194e-15 1.063109e-14]; 5th = 6.2440; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [2.388839e-18 9.182568e-16 1.252587e-14 1.902626e-14]; 5th = 6.2440; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.496768 0.496768 0.496768 0.496768]; 5th = 5.7985; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.496768 0.496768 0.496768 0.496768]; 5th = 5.7985; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.49676772 in [e^(-|eps|/pi) d, d] = [0.42643210, 0.50000000]
  ok   both operators: exactly four eigenvalues in [-r, r]
  -- eps=1.0, c=1.0: sup|F| = 6.283, r = |eps| sqrt(1+c^2) = 1.414, d = 0.5
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [2.217783e-20 4.867301e-18 2.112189e-16 2.020598e-14]; 5th = 6.1330; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [4.509335e-17 8.482388e-16 3.166441e-15 2.863513e-14]; 5th = 6.1330; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.487263 0.487263 0.487263 0.487263]; 5th = 5.8441; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.487263 0.487263 0.487263 0.487263]; 5th = 5.8441; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.48726265 in [e^(-|eps|/pi) d, d] = [0.36368867, 0.50000000]
  ok   both operators: exactly four eigenvalues in [-r, r]
  -- eps=0.001, c=2.0: sup|F| = 0.006283, r = |eps| sqrt(1+c^2) = 0.002236, d = 0.001
     E_eps   NY=12: dim ker = 4; 4 smallest |lambda| = [2.700911e-22 3.336382e-17 3.418609e-16 1.378904e-14]; 5th = 6.2832; #eigs in [-r,r] = 4
     E_eps   NY=20: dim ker = 4; 4 smallest |lambda| = [1.369640e-15 1.724137e-15 1.742322e-14 1.814423e-14]; 5th = 6.2832; #eigs in [-r,r] = 4
  ok   E_eps: dim ker D = 4
  ok   E_eps: kernel in x-mode 0, tau-split 2+2, parity split even 2 + odd 2
  ok   E_eps: dim(ker ∩ Omega^k) = [0, 2, 0], dim(projection of ker to Omega^k) = [2, 2, 2]
  ok   E_eps: all nonzero |lambda| >= 4
     f^*E_eps NY=12: dim ker = 0; 4 smallest |lambda| = [0.001 0.001 0.001 0.001]; 5th = 6.2822; #eigs in [-r,r] = 4
     f^*E_eps NY=20: dim ker = 0; 4 smallest |lambda| = [0.001 0.001 0.001 0.001]; 5th = 6.2822; #eigs in [-r,r] = 4
  ok   f^*E_eps: dim ker D = 0
  ok   f^*E_eps: four eigenvalues +-sigma, sigma = 0.00100000 in [e^(-|eps|/pi) d, d] = [0.00099968, 0.00100000]
  ok   both operators: exactly four eigenvalues in [-r, r]
Part C: kernel criterion (Lemma 2.2) for connections depending on x and y (2D Fourier-Galerkin)
  ok   (A,B)=(0,0): smallest singular values of 2dbar, 2d: 1.24e-16, 1.12e-17; (A,B) in 2piZ^2: True
  ok   (A,B)=(6.283,0): smallest singular values of 2dbar, 2d: 1.74e-16, 9.02e-17; (A,B) in 2piZ^2: True
  ok   (A,B)=(6.283,-6.283): smallest singular values of 2dbar, 2d: 9.40e-17, 2.01e-16; (A,B) in 2piZ^2: True
  ok   (A,B)=(12.57,6.283): smallest singular values of 2dbar, 2d: 9.46e-17, 2.73e-16; (A,B) in 2piZ^2: True
  ok   (A,B)=(0.3,0): smallest singular values of 2dbar, 2d: 2.99e-01, 2.99e-01; (A,B) in 2piZ^2: False
  ok   (A,B)=(0,0.3): smallest singular values of 2dbar, 2d: 2.99e-01, 2.99e-01; (A,B) in 2piZ^2: False
  ok   (A,B)=(3.142,3.142): smallest singular values of 2dbar, 2d: 4.22e+00, 4.22e+00; (A,B) in 2piZ^2: False
  ok   (A,B)=(6.283,0.05): smallest singular values of 2dbar, 2d: 4.98e-02, 4.98e-02; (A,B) in 2piZ^2: False
Part D: parameter range: dim ker D(f^*E_eps) = 4 iff eps c in 4 pi Z (eps = 0.5)
  ok   c = 25.13274: eps c/(4 pi) = 1.00000, dim ker = 4 (expected 4)
  ok   c = 25.03274: eps c/(4 pi) = 0.99602, dim ker = 0 (expected 0)
  ok   c = 25.23274: eps c/(4 pi) = 1.00398, dim ker = 0 (expected 0)
  ok   c = 50.26548: eps c/(4 pi) = 2.00000, dim ker = 4 (expected 4)
  ok   c = 12.56637: eps c/(4 pi) = 0.50000, dim ker = 0 (expected 0)
  ok   c = 0.00000: eps c/(4 pi) = 0.00000, dim ker = 4 (expected 4)
Part E: degree-wise quantities on T^2 (twisted Laplacians Delta_k = d d^* + d^* d)
     E_eps   eps=0.01, c=1.0: dim ker Delta_k = [0, 2, 0]; #eigenvalues of Delta_k in [0,1] = [1, 2, 1]
  ok   E_eps: dim ker Delta_k as in the paper
  ok   E_eps: small-eigenvalue counts equal the Betti numbers of T^2 (computation)
     f^*E_eps eps=0.01, c=1.0: dim ker Delta_k = [0, 0, 0]; #eigenvalues of Delta_k in [0,1] = [1, 2, 1]
  ok   f^*E_eps: dim ker Delta_k as in the paper
  ok   f^*E_eps: small-eigenvalue counts equal the Betti numbers of T^2 (computation)
     E_eps   eps=0.2, c=0.7: dim ker Delta_k = [0, 2, 0]; #eigenvalues of Delta_k in [0,1] = [1, 2, 1]
  ok   E_eps: dim ker Delta_k as in the paper
  ok   E_eps: small-eigenvalue counts equal the Betti numbers of T^2 (computation)
     f^*E_eps eps=0.2, c=0.7: dim ker Delta_k = [0, 0, 0]; #eigenvalues of Delta_k in [0,1] = [1, 2, 1]
  ok   f^*E_eps: dim ker Delta_k as in the paper
  ok   f^*E_eps: small-eigenvalue counts equal the Betti numbers of T^2 (computation)
Part F: variants (real rank-2 bundle, spin Dirac operator, T^4 = T^2 x T^2)
  ok   real rank-2 bundle, E_eps: dim ker D = 8 (expected 8)
  ok   real rank-2 bundle, f^*E_eps: dim ker D = 0 (expected 0)
  ok   Clifford relations for the spin representation of Cl(R^2)
  ok   spin Dirac operator (periodic spin structure), E_eps: dim ker = 2 (expected 2)
  ok   spin Dirac operator (periodic spin structure), f^*E_eps: dim ker = 0 (expected 0)
  ok   tau^2 = 1 on Lambda^*(R^4) (x) C
  ok   T^4, pr1^*E_eps: dim ker D = 16 (tau-split 8+8); dim(ker ∩ Omega^k) = [0, 2, 4, 2, 0]
  ok   T^4: tau-split 8+8 and degree-wise (0,2,4,2,0)
  ok   T^4, (f x id)^*pr1^*E_eps: dim ker D = 0 (tau-split 0+0); dim(ker ∩ Omega^k) = [0, 0, 0, 0, 0]
Part G: explicit kernel elements of D(E_eps) on a grid (spectral differentiation)
  ok   eps=0.05: dbar(s+) = 0 and d(s-) = 0 (residuals 1.8e-13, 1.8e-13)
  ok   eps=0.05: D annihilates the four explicit forms (max residual 1.8e-13)
  ok   eps=0.7: dbar(s+) = 0 and d(s-) = 0 (residuals 1.9e-13, 1.4e-13)
  ok   eps=0.7: D annihilates the four explicit forms (max residual 1.9e-13)

ALL CHECKS PASSED
