module wstrass;
% Author: James H. Davenport.
fluid '(!*exp
!*gcd
!*mcd
!*structure
!*uncached
!*keepsqrts % Forced SIMP to keep square roots around
!*tra
!*trmin
intvar
listofallsqrts
listofnewsqrts
magiclist
previousbasis
sqrt!-intvar
sqrtflag
sqrts!-in!-integrand
taylorasslist
taylorvariable
thisplace
zlist);
global '(coates!-fdi);
exports simpwstrass,weierstrass!-form;
imports gcdn,sqrtsinplaces,
makeinitialbasis,mkvec,completeplaces,integralbasis,
normalbasis,mksp,multsq,xsubstitutesq,taylorform,taylorevaluate,
coatessolve,checkpoles,substitutesq,removecmsq,printsq,interr,
terpri!*,printplace,finitise,fractional!-degree!-at!-infinity,
!*multsq,fdi!-print,fdi!-upgrade,fdi!-revertsq,simp,newplace,
xsubstitutep,sqrtsinsq,removeduplicates,!*exptf,!*multf,
!*multsq,!*q2f,mapvec,upbv,coates!-lineq,addsq,!*addsq;
symbolic procedure simpwstrass u;
begin
scalar intvar,sqrt!-intvar,taylorvariable,taylorasslist;
scalar listofallsqrts,listofnewsqrts;
scalar sqrtflag,sqrts!-in!-integrand,tt,u;
scalar !*keepsqrts,!*exp,!*gcd,!*mcd,!*structure,!*uncached;
!*keepsqrts:=t; % Else nothing will work
!*exp := !*gcd := !*mcd := !*uncached := t;
!*structure := nil; % Algebraic code certainly wants this off:
% keeping it on inhibits sqrt(z)^2 -> z
tt:=readclock();
sqrtflag:=t;
taylorvariable:=intvar:=car u;
sqrt!-intvar:=mvar !*q2f simpsqrti intvar;
u:=for each v in cdr u
collect int!-simp v;
sqrts!-in!-integrand:=sqrtsinsql(u,intvar);
u:=errorset!*('(weierstrass!-form sqrts!-in!-integrand),t);
if atom u
then return u
else u:=car u;
printc list('time,'taken,readclock()-tt,'milliseconds);
printc "New x value is:";
printsq car u;
u:=cdr u;
printc "New y value is:";
printsq car u;
u:=cdr u;
printc "Related by the equation";
printsq car u;
return car u
end;
put('wstrass,'simpfn,'simpwstrass);
symbolic procedure weierstrass!-form sqrtl;
begin
scalar sqrtl2,u,x2,x1,vect,a,b,c,d,lhs,rhs;
if !*tra or !*trmin then <<
printc "Find Weierstrass form for elliptic curve defined by:";
for each u in sqrtl do
printsq simp u >>;
sqrtl2:=sqrts!-at!-infinity sqrtl;
sqrtl2:=append(car sqrtl2,
for each u in cdr sqrtl2 collect
u.u);
% one of the places lying over infinity
% (after deramification as necessary).
x2:=coates!-wstrass(list sqrtl2,list(-3),intvar);
% Note that we do not multiply by the MULTIPLICITY!-FACTOR
% since we genuinely want a pole of order -3 irrespective
% of any ramification problems.
if !*tra then <<
printc "Function with pole of order 3 (x2) is:";
printsq x2 >>;
x1:=coates!-wstrass(list sqrtl2,list(-2),intvar);
if !*tra then <<
printc "Function with pole of order 2 (x1) is:";
printsq x1 >>;
vect := mkvec list(1 ./ 1,
x1,
x2,
!*multsq(x1,x1),
!*multsq(x2,x2),
!*multsq(x1,!*multsq(x1,x1)),
!*multsq(x1,x2));
u:=!*lcm!*(!*exptf(denr x1,3),!*multf(denr x2,denr x2)) ./ 1;
for i:=0:6 do
putv(vect,i,!*q2f !*multsq(u,getv(vect,i)));
if !*tra then <<
printc "List of seven functions in weierstrass!-form:";
mapvec(vect,function printsf) >>;
vect := wstrass!-lineq vect;
if vect eq 'failed then
interr "Linear equation solving failed in Weierstrass";
% printsq(addsq(getv(vect,0),addsq(!*multsq(getv(vect,1),x1),
% addsq(!*multsq(getv(vect,2),x2),
% addsq(!*multsq(getv(vect,3),!*multsq(x1,x1)),
% addsq(!*multsq(getv(vect,4),!*multsq(x2,x2)),
% addsq(!*multsq(getv(vect,5),exptsq(x1,3)),
% !*multsq(getv(vect,6),
% !*multsq(x1,x2)))))))));
x2:=!*addsq(!*multsq(!*multsq(2 ./ 1,getv(vect,4)),x2),
!*addsq(!*multsq(x1,getv(vect,6)),
getv(vect,2)));
putv(vect,4,!*multsq(-4 ./ 1,getv(vect,4)));
a:=!*multsq(getv(vect,4),getv(vect,5));
b:=!*addsq(!*multsq(getv(vect,6),getv(vect,6)),
!*multsq(getv(vect,3),getv(vect,4)));
c:=!*addsq(!*multsq(2 ./ 1,!*multsq(getv(vect,2),getv(vect,6))),
!*multsq(getv(vect,1),getv(vect,4)));
d:=!*addsq(!*multsq(getv(vect,2),getv(vect,2)),
!*multsq(getv(vect,0),getv(vect,4)));
lhs:=!*multsq(x2,x2);
rhs:=!*addsq(d,!*multsq(x1,
!*addsq(c,!*multsq(x1,
!*addsq(b,!*multsq(x1,a))))));
if lhs neq rhs then <<
printsq lhs;
printsq rhs;
interr "Previous two unequal - consistency failure 1" >>;
u:=!*lcm!*(!*lcm!*(denr a,denr b),!*lcm!*(denr c,denr d));
if u neq 1 then <<
% for now use u**2 whereas we should be using the least
% square greater than u**2 (does it really matter).
x2:=!*multsq(x2,u ./ 1);
u:=!*multf(u,u) ./ 1;
a:=!*multsq(a,u);
b:=!*multsq(b,u);
c:=!*multsq(c,u);
d:=!*multsq(d,u) >>;
if (numr b) and not (quotf(numr b,3)) then <<
% multiply all through by 9 for the hell of it.
x2:=multsq(3 ./ 1,x2);
u:=9 ./ 1;
a:=multsq(a,u);
b:=multsq(b,u);
c:=multsq(c,u);
d:=multsq(d,u) >>;
x2:=!*multsq(x2,a);
x1:=!*multsq(x1,a);
c:=!*multsq(a,c);
d:=!*multsq(!*multsq(a,a),d);
lhs:=!*multsq(x2,x2);
rhs:=!*addsq(d,!*multsq(x1,!*addsq(c,!*multsq(x1,!*addsq(b,x1)))));
if lhs neq rhs then <<
printsq lhs;
printsq rhs;
interr "Previous two unequal - consistency failure 2" >>;
b:=quotf(numr b,3) ./ 1;
x1:=!*addsq(x1,b);
d:=!*addsq(d,!*addsq(multsq(2 ./ 1,!*multsq(b,!*multsq(b,b))),
negsq !*multsq(c,b)));
c:=!*addsq(c,multsq((-3) ./ 1,!*multsq(b,b)) );
% b:=nil ./ 1; % not used again.
if !*tra then <<
printsq x2;
printsq x1;
printc "with coefficients";
printsq c;
printsq d;
rhs:=!*addsq(d,
!*addsq(!*multsq(c,x1),
!*multsq(x1,!*multsq(x1,x1)) ));
lhs:=!*multsq(x2,x2);
if lhs neq rhs then <<
printsq lhs;
printsq rhs;
interr "Previous two unequal - consistency failure 3" >> >>;
return weierstrass!-form1(c,d,x1,x2)
end;
symbolic procedure !*lcm!*(u,v); !*multf(u,quotf(v,gcdf(u,v)));
symbolic procedure weierstrass!-form1(c,d,x1,x2);
begin scalar b,u;
u:=gcdf(numr c,numr d);
% We will reduce by anything whose square divides C
% and whose cube divides D.
if not numberp u then begin
scalar cc,dd;
u:=jsqfree(u,mvar u);
u:=cdr u;
if null u
then return;
% We found no repeated factors.
for each v in u do
for each w in v do
while (cc:=quotf(numr c,multf(w,w))) and
(dd:=quotf(numr d,exptf(w,3)))
do <<
c:=cc ./ 1;
d:=dd ./ 1;
x1:=!*multsq(x1,1 ./ w);
x2:=!*multsq(x2,1 ./ multf(w,simpsqrt2 w)) >>;
u:=gcdn(algint!-numeric!-content numr c,
algint!-numeric!-content numr d)
end;
b:=2;
while not(b*b > u) do begin
scalar nc,nd,uu;
nc:=0;
while cdr(uu:=divide(u,b))=0 do <<
nc:=nc+1;
u:=car uu >>;
if nc < 2
then go to next;
uu:=algint!-numeric!-content numr d;
nd:=0;
while cdr(uu:=divide(uu,b))=0 do <<
nd:=nd+1;
uu:=car uu >>;
if nd < 3
then go to next;
nc:=min(nc/2,nd/3);
% re-normalise by b**nc.
uu:=b**nc;
c:=multsq(c,1 ./ (uu**2));
d:=multsq(d,1 ./ (uu**3));
x1:=multsq(x1,1 ./ uu);
x2:=multsq(x2,1 ./ (uu*b**(nc/2)) );
if not evenp nc
then x2:=!*multsq(x2,!*invsq simpsqrti b);
next:
b:=nextprime(b)
end;
u:=!*kk2q intvar;
u:=!*addsq(!*addsq(d,multsq(c,u)),exptsq(u,3));
if !*tra or !*trmin then <<
printc "Standard form is y**2 = ";
printsq u >>;
return list(x1,x2,simpsqrtsq u)
end;
symbolic procedure sqrts!-at!-infinity sqrtl;
begin
scalar inf,hack,sqrtl2,repeating;
hack:=list list(intvar,'expt,intvar,2);
inf:=list list(intvar,'quotient,1,intvar);
sqrtl2:=list sqrt!-intvar;
while sqrt!-intvar member sqrtl2 do <<
if repeating
then inf:=append(inf,hack);
newplace inf;
sqrtl2:=for each v in sqrtl conc
sqrtsinsq(xsubstitutep(v,inf),intvar);
repeating:=t >>;
sqrtl2:=removeduplicates sqrtl2;
return inf.sqrtl2
end;
symbolic procedure coates!-wstrass(places,mults,x);
begin
scalar thisplace,u,finite!-hack,save,places2,mults2;
if !*tra or !*trmin then <<
princ "Find function with zeros of order:";
printc mults;
if !*tra then
princ " at ";
terpri!*(t);
if !*tra then
mapc(places,function printplace) >>;
% finite!-hack:=placesindiv places;
% FINITE!-HACK is a list of all the substitutors in PLACES;
% u:=removeduplicates sqrtsintree(finite!-hack,x,nil);
% if !*tra then <<
% princ "Sqrts on this curve:";
% terpri!*(t);
% superprint u >>;
% algnos:=removeduplicates for each j in places collect basicplace j;
% if !*tra then <<
% printc "Algebraic numbers where residues occur:";
% superprint algnos >>;
finite!-hack:= finitise(places,mults);
% returns list (places,mults,power of x to remove.
places2:=car finite!-hack;
mults2:=cadr finite!-hack;
finite!-hack:=list(places,mults,caddr finite!-hack);
coates!-fdi:=fractional!-degree!-at!-infinity u;
if coates!-fdi iequal 1
then return !*multsq(wstrassmodule(places2,mults2,x,finite!-hack),
!*p2q mksp(x,caddr finite!-hack));
if !*tra
then fdi!-print();
newplace list(intvar . thisplace,
list(intvar,'expt,intvar,coates!-fdi));
places2 := for each j in places2 collect fdi!-upgrade j;
save:=taylorasslist;
u:=wstrassmodule(places2,
for each u in mults2 collect u*coates!-fdi,
x,finite!-hack);
taylorasslist:=save;
u:=fdi!-revertsq u;
return !*multsq(u,!*p2q mksp(x,caddr finite!-hack))
end;
symbolic procedure wstrassmodule(places,mults,x,finite!-hack);
begin
scalar pzero,mzero,u,v,basis,sqrts,magiclist,mpole,ppole;
% MAGICLIST holds the list of extra unknowns created in JHDSOLVE
% which must be found in CHECKPOLES (calling FINDMAGIC).
sqrts:=sqrtsinplaces places;
if !*tra then <<
princ "Sqrts on this curve:";
superprint sqrts >>;
u:=places;
v:=mults;
while u do <<
if 0<car v
then <<
mzero:=(car v).mzero;
pzero:=(car u).pzero >>
else <<
mpole:=(car v).mpole;
ppole:=(car u).ppole >>;
u:=cdr u;
v:=cdr v >>;
basis:=mkvec makeinitialbasis ppole;
u:=completeplaces(ppole,mpole);
basis:=integralbasis(basis,car u,cdr u,x);
basis:=normalbasis(basis,x,0);
u:=coatessolve(mzero,pzero,basis,force!-pole(basis,finite!-hack));
% This is the list of special constraints needed
% to force certain poles to occur in the answer.
previousbasis:=nil;
if atom u
then return u;
v:= checkpoles(list u,places,mults);
if null v
then return 'failed;
if not magiclist
then return u;
u:=removecmsq substitutesq(u,v);
% Apply the values from FINDMAGIC.
if !*tra or !*trmin then <<
printc "Function is";
printsq u >>;
magiclist:=nil;
if checkpoles(list u,places,mults)
then return u
else interr "Inconsistent checkpoles"
end;
symbolic procedure force!-pole(basis,finite!-hack);
begin
scalar places,mults,u,ans;
places:=car finite!-hack;
mults:=cadr finite!-hack;
finite!-hack:=caddr finite!-hack;
u:=!*p2q mksp(intvar,finite!-hack);
basis:=for each v in basis collect multsq(u,v);
while places do <<
u:=for each v in basis collect
taylorevaluate(taylorform xsubstitutesq(v,car places),
car mults);
mults:=cdr mults;
places:=cdr places;
ans:=u.ans >>;
return ans
end;
symbolic procedure wstrass!-lineq vect;
begin
scalar zlist,powlist,m,rightside,v;
scalar zero,one;
zero:=nil ./ 1;
one:=1 ./ 1;
for i:=0:6 do
zlist:=varsinsf(getv(vect,i),zlist);
zlist:=intvar . findzvars(zlist,nil,intvar,nil);
for i:=0:6 do
putv(vect,i,f2df getv(vect,i));
for i:=0:6 do
for each u in getv(vect,i) do
if not ((tpow u) member powlist)
then powlist:=(tpow u).powlist;
m:=for each u in powlist collect begin
scalar v;
v:=mkvect 6;
for i:=0:6 do
putv(v,i,(lambda u;
if null u
then zero
else tc u)
assoc(u,getv(vect,i)));
return v
end;
v:=mkvect 6;
for i:=0:6 do
putv(v,i,zero);
putv(v,4,one);
% we know that coefficient e is non-zero.
m:=mkvec (v.m);
v:=upbv m;
rightside:=mkvect v;
putv(rightside,0,one);
for i:=1:v do
putv(rightside,i,zero);
return coates!-lineq(m,rightside)
end;
% This is same as NUMERIC!-CONTENT in the EZGCD module, but is included
% here so that that module doesn't need to be loaded.
symbolic procedure algint!-numeric!-content form;
%Find numeric content of non-zero polynomial.
if domainp form then abs form
else if null red form then algint!-numeric!-content lc form
else begin
scalar g1;
g1 := algint!-numeric!-content lc form;
if not (g1=1)
then g1 := gcddd(g1,algint!-numeric!-content red form);
return g1
end;
endmodule;
end;