/* wxyz.gp */
/* Four polynomial sequences w,x,y,z are discrete
versions of the four Jacobi theta functions or
the four Weierstrass sigma functions. */
/* 17 May 2016 Michael Somos */
/* 24 Dec 2018 Michael Somos */
/* This file is a PARI/GP program in the GP language. See the
PARI/GP home page at URL http://pari.math.u-bordeaux.fr/
This file is available at URL http://grail.eecs.csuohio.edu/~somos/wxyz.gp
Mathematica version is at URL http://grail.eecs.csuohio.edu/~somos/wxyz.wl
The inspiration of my research was a solution of the
congruent number 5 problem on pp. 419-427 in the book:
Uspensky and Heaslet, Elementary Number Theory, 1939.
I placed the sequences A129206,..,A129209 in the OEIS.
My original code dates back to about 2002 partly based
on Kerawala's 1947 article on Poncelet's porism.
The sequences are homogenous polynomials in the
variables X,Y,Z with scale factors of W,x0,y0,z0
for the four sequences respectively. Thus, the w
sequence terms all have a factor of W, while the x
sequence terms all have a factor of x0, and so on.
The n-th sequence terms all have a factor of Q^n^2.
Note that w is an elliptic divisibility sequence.
The w sequence is an odd sequence while the x,y,z
sequences are all even sequences. Thus, w is an
analog of theta_1, while x,y,z are analogs of the
other three Jacobi theta functions. Also, w is an
analog of the Weierstrass sigma function, while
x,y,z are analogs of the other sigma functions.
More precisely, given numbers t and |q|<1, while
x0 = theta_2(0, q), y0 = theta_3(0, q),
z0 = theta_4(0, q), W = theta_1(t, q),
X = theta_2(t, q)/x0, Y = theta_3(t, q)/y0,
Z = theta_4(t, q)/z0, and Q = 1, then we get
wn(n) = theta_1(n*t, q), xn(n) = theta_2(n*t, q),
yn(n) = theta_3(n*t, q), zn(n) = theta_4(n*t, q).
There is a similar result for the four Weierstrass
sigma functions. I have versions of the Weierstrass
zeta function and its first few derivatives which I
notate beginning with the letter "w", but they are
rational functions. I have also polynomial versions
beginning with the letter "W" which are more closely
related to the the sigma function polynomial sequences.
Note that I have introduced several constants with
more or less arbitrary names. The g2,g3 are just the
Weierstrass invariants. The ex,ey,ez correspond to
the Weierstrass e1,e2,e3. The DD is the discriminant
Delta. The j,J correspond to Klein's modular function.
The p1,p2,p3,p4,p5 are my invariants of generalized
Somos-4 sequences.
Note that it is possible to take the "derivative" of
any expression using a derivative function that I
named Du, and use it to verify identities for the
derivatives of the Weierstrass elliptic functions
and the Jacobi elliptic functions including analogs
of the Jacobi Zeta and Epsilon functions.
Note that not all possible identities between the
elliptic functions are satisfied by these polynomial
sequences. For exmaple, the important identity for
the 4th power of theta null functions does not hold.
That is, the identity 0 = x0^4 +y0^4 -z0^4 is clearly
not valid since the variables x0, y0, z0 are free.
Note that this is a work in progress and there may
be slight changes in detail but almost everything
in here is in a workable state which is unlikely to
change in future. Please inform me of any errors you
find. I may write fuller documentation if there is
any real interest in my work detailed here.
Note the following special case: if
W = u^2-v^2, X = (u^2+v^2)/2, Y = Z = u*v,
x0 = 2, y0 = z0 = Q = 1, then we have
wn(n) = (u*v)^(n^2-n) * (u^(2*n)-v^(2*n)),
xn(n) = (u*v)^(n^2-n) * (u^(2*n)+v^(2*n)),
and yn(n) = zn(n) = (u*v)^(n^2).
My results are comparable to those given in 1948
and 1950 articles by Morgan Ward in the American
Journal of Mathematics, but different in origin
and scale. More directly comparable are formulas
of the Chudnovsky's from an article in Advances
in Applied Mathematics from 1986 on page 418.
In 1878 Edouard Lucas published a memoir on the
arithmetization of elliptic functions. He may have
been working towards the generalization of this
special case without success. He refers to a
memoir of Moutard that was included in the 1862
volume 1 of Poncelet's book _Applications ..._
pp. 535-560, but apparently he did not really
understand its implications for his quest.
More details about this work of Lucas is contained
in Hugh C. Williams, _E'douard Lucas and Primality
Testing_, 1998, Wiley, pp. 451-457. On page 453 is
h{m+n}h{m-n}=h{m+1}h{m-1}h(n)^2-h{n+1}h{n-1}h(m)^2
which is similar to one of our own equations and
this quote "Lucas mentioned that these formulas
belonged to the theory of elliptic functions,...".
The memoir by Moutard on page 542 defined three
sequences of homogeneous polynomials using two
sets duplication equations. Equation (6) for even
and (7) for odd cases. For example, b_{2p}+c_{2p} =
2b^2_pc^2_p. The three sequences are slightly
special cases of x_n, y_n, z_n. Thus, a_n = x_n
with Q = x_0 = y_0 = z_0 = 1, x_1 = a, y_1 = b,
z_1 = c.
Note that some chk() computations are commented
out because they take too much time to complete.
If that is not a concern for you, you can remove
the \\ comment chars at the beginning of the line.
*/
/* weights
0: Q
1: W,X,Y,Z,Wz
*/
{addhelp(degrees,"degrees(): prints a list of values ordered by degree.")};
{degrees() = print1(
"0: W,X,Y,Z,Wz,JZ,x0,y0,z0,ex,ey,ez,j,J,g2,g3, wz1,wz2,wz3,wp1,wp2,wp3,...\n"
"1: w1,x1,y1,z1,X1,Y1,Z1,x1/x0,y1/y0,z1/z0,i1, Wz1\n"
"2: XY,XZ,YZ,f1,p5,i2,h1,hx,hy,hz,Ex,Ey,Ez, Wp1\n"
"3: b2,i3,iq, Wpp1\n"
"4: G2,f2,r2,b1,p4,h2,w2,x2,y2,z2,l1,m1, W2p1\n"
"5: W3p1\n"
"6: G3,Dr,f3,fq,h3,hq,b5,r3,p1,j1,jx,jy,jz, W4p1\n"
"7: W5p1\n"
"8: b3,c1,p2,l2,m2,w3/w1,x3/x1,y3/y1,z3/z1, Wp2\n"
"9: n1,w3/w1,x3/x0,y3/y0,z3/z0\n"
"12: DD,d1,b4,p3,l3,lq,m3,w4/w2,k1,kx,ky,kz, Wpp2\n"
"16: l4,l5,w4,x4,y4,z4, W2p2\n"
"18: n2,nx,ny,nz, Wp3\n"
"20: W3p2\n"
"24: c2,c3,c4,c5,c6,w5/w1,x5/x1,y5/y1,z5/z1, W4p2\n"
"25: w5\n"
"27: nq,n3,w6/w3, Wpp3\n"
"28: W5p2\n"
"32: x6/x2,y6/y2,z6/z2, Wp4\n"
"36: w6,x6,y6,z6, W2p3\n"
"45: W3p3\n"
"48: c8,w8/w4,c7,w7/w1,x7/x1,y7/y1,z7/z1, Wpp4\n"
"64: x8/x0, W2p4\n"
"72: c9,c10\n"
"75: Wpp5\n"
"80: W3p4\n"
"96: c12\n"
"108: Wpp6\n"
"120: c11\n"
)};
/* The minimum runtime stack required. You may want more */
NN = 90*10^6;
/* This is to avoid stack overflow while computing */
if(default(parisize)