Learning with GP: Verification of identities

[Jacques Gélinas <jacquesg00@hotmail.com>]

How can GP detect typos and sign errors in published identities ?
A quick test is also needed after transforming formulas or using
special cases of general identities. However, "===" and "==" are
not adequate for floating-point, summations, infinite series,
limits, etc., when equality conditions must be relaxed somehow.

We first store in "verif.gp" some routines for such comparisons.

\\-------- File: verif.gp --------------------- Start cut

          \\ relative proportion of equal bits
eqr(x,y) = if(x==y,1,exponent(if(!x,y,1-y/x))/exponent(0.));
          \\ floating-point equality at r*precision
fleq(x,y,r=7/8) = if(x==y, oo, eqr(0,\
         normlp(Pol(x-y))/normlp(Pol(x+y)))) >= r;

IS(e,msg="") = print(if(e,"Pass: ","Fail: "),msg);
ISF(x,y,msg,r=7/8) = IS( fleq(x,y,r),msg);
          \\ f(n)=g(n) on [1,N]
ISN(f,g,msg,r=7/8,N=8) = IS( N == 
             sum(n=1,N, fleq(f(n),g(n),r)), msg);
          \\ f(m,n)=g(m,n) on [1,M]x[1,N]
ISMN(f,g,msg,r=7/8,M=8,N=8) = IS( M*N == 
  sum(m=1,M, sum(n=1,N, fleq(f(m,n),g(m,n),r))), msg);

\\-------- File: verif.gp --------------------- End paste

Next, we define "ckverif()" to perform mininal tests
of the previous routines, using identities published
by John Blissard, the inventor in the 1860s of the 
representative notation (known as umbral calculus).

\\-------------------------------------------- Start cut

ckverif() = {
                          \\ floating-point
  IS(1.===1,"1.===1"); IS(1.==1,"1.==1");
  ISF( Pi, 4*atan(1), "Pi=4*atan(1) (r=oo)",oo);
  print( "eqr(Pi,4*atan(1)) = ", eqr(Pi,4*atan(1)) );
  ISF( Pi, 4*atan(1), "Pi=4*atan(1) (r=1)",1);

                          \\ vectors
  ISF( vector(8,n, -bernfrac(2*n)/2/n),\
       vector(8,n, zeta(1-2*n)), "bernfrac (r=oo)",oo);
  ISF( vector(8,n, -bernfrac(2*n)/2/n),\
       vector(8,n, zeta(1-2*n)), "bernfrac (r=1)",1);
       
                           \\ finite sums [1]
  ISN( n->sum(k=1,n,k), n->n*(n+1)/2, "1+...+n (r=oo)",oo);
  ISMN((m,n)->sum(k=0,n, (-1)^k*binomial(n,k)/(m+k)),
       (m,n)->n!/prod(k=0,n,m+k), "Blissard 124 (r=oo)",oo);

                           \\ Maclaurin series [1]
  my( A(k,n)=sum(j=1,k,(-1)^(k-j)*binomial(n,j-1)/(k+1-j)));
  ISN( n->(1+x)^n*log(1+x), n->sum(k=1,n, A(k,n)*x^k )
     + n!*sum(k=1, default(seriesprecision),
     (-1)^(k-1)*(k-1)!/(n+k)!*x^(n+k) ), "Blissard 125 (r=oo)",oo);

                           \\ infinite series [1]
  my( L130(n,N) = sum(k=1,N, (-1.)^(k-1)*(k-1)!*n!/(n+k)!) );
  my( R130(n) = 2^n*log(2) - sum(k=1,n, A(k,n)) );
  ISF( L130(16,4000), R130(16), "Blissard 130 (r=3/4)", 3/4);
}

ckverif();

\\-------------------------------------------- End paste

Output for GP 2.12.0 on Windows/Cygwin-64 (38 digits) :

Fail: 1.===1
Pass: 1.==1
Fail: Pi=4*atan(1) (r=oo)
eqr(Pi,4*atan(1)) = 127/128
Pass: Pi=4*atan(1) (r=1)
Fail: bernfrac (r=oo)
Pass: bernfrac (r=1)
Pass: 1+...+n (r=oo)
Pass: Blissard 124 (r=oo)
Pass: Blissard 125 (r=oo)
Pass: Blissard 130 (r=3/4)

------------------------------ Features used
!,==,===,Pi,atan,exp,log,binomial,default,
seriesprecision,exponent,normlp,sum,prod,vector,
bernfrac,print,if,my,anonymous closure

------------------------------ References
[1] John Blissard 1862 On the generalization of certain
    trigonometrical formulae, Messenger of Mathematics I, 124-136
    https://archive.org/details/oxfordcambridge00unkngoog/page/n136
[2] Tutorial: http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.12.0/tutorial.pdf
[3] Reference card: http://pari.math.u-bordeaux.fr/pub/pari/manuals/2.12.0/refcard.pdf

Jacques Gélinas         Next: Verification of identities - Romik numbers