reduce-historical/r33/algint.red

module afactor;

% Author: James H. Davenport.

fluid '(!*galois !*noextend !*sqfree afactorvar listofnewsqrts
        monicpart);

global '(!*trfield);

exports afactor;
imports exptf,ordop,!*multf,addf,makemainvar,algebraicsf,divsf,contents;
imports quotf!*,negf,sqfr!-norm2,prepf,gcdinonevar,algint!-subf,!*q2f;
imports jfactor,printsf;

% internal!-fluid '(monicpart);

symbolic procedure afactor(u,v);
  % Factorises U over the algebraics as a polynomial in V (=afactorvar).
begin
  scalar afactorvar,!*noextend,!*sqfree;
  % !*sqfree is known to be square free (from sqfr-norm).
  !*noextend:=t; % else we get recursion.
  afactorvar:=v;
  if !*trfield
    then <<
    princ "We must factorise the following over: ";
    for each u in listofnewsqrts do <<princ u; princ " " >>;
    terpri();
    printsf u >>;
  v:=algfactor u;
  if !*trfield then <<
    printc "factorises as ";
    mapc(v,function printsf) >>;
  return v
  end;


symbolic procedure algfactor2(f,a);
if null a
  then for each u in jfactor(f,mvar f) collect numr u
  else if algebraicsf f
    then algfactor3(f,a)
    else begin
      scalar w;
      if !*trfield then <<
        princ "to be factorized over ";
        for each u in a do << princ u; princ " " >>;
        terpri();
        printsf f >>;
      if  (!*galois neq 2) and
          (numberp red f) and
          (not numberp argof car a)
        then return algfactor2(f,cdr a);
        % assumes we need never express a root of a number in terms of
        % non-numbers.
      w:=algfactor2(f,nil);
      if null cdr w
        then return algfactor3(f,a)
        else return 'partial.w
      end;


symbolic procedure algfactor3(f,a);
begin
  scalar ff,w,gg,h,p;
  w:=sqfr!-norm2(f,mvar f,car a);
  !*sqfree:=car w;
  w:=cdr w;
  ff:=algfactor2(!*sqfree,cdr a);
  if car ff eq 'partial then <<
    p:='partial;
    ff:=cdr ff >>;
  if null cdr ff
    then return list f;
    %does not factor.
  a:=car a;
  gg:=cadr w;
  w:=list list(afactorvar,'plus,afactorvar,prepf car w);
  h:=for each u in ff
       collect (!*q2f algint!-subf(gcdinonevar(u,gg,afactorvar),w));
  if p eq 'partial
    then h:=p.h;
  return h
  end;

symbolic procedure algfactor u;
begin
  scalar a,aa,z,w,monicpart;
  z:= makemainvar(u,afactorvar);
  if ldeg z iequal 1
    then return list u;
  z:=lc z;
  if z iequal 1
    then go to monic;
  if algebraicsf z
    then u:=!*multf(u,numr divsf(1,z));
    % this de-algebraicises the top coefficient.
  u:=quotf!*(u,contents(u,afactorvar));
  z:=makemainvar(u,afactorvar);
  if lc z neq 1
    then if lc z iequal -1
      then u:=negf u
      else <<
        w:=lc z;
        u:=makemonic z >>;
monic:
  aa:=listofnewsqrts;
  if algebraicsf u
    then go to normal;
  a:=cdr aa;
  % we need not try for the first one, since algfactor2
  % will do this for us.
  z:=t;
  while a and z do begin
    scalar alg,v;
    alg:=car a;
    a:=cdr a;
    v:=algfactor3(u,list alg);
    if null cdr v
      then return;
    if car v eq 'partial
      then v:=cdr v;
      % we do not mind if the factorisation is only partial.
    a:=mapcan(v,function algfactor);
    z:=nil
    end;
  monicpart:=w;
  if null z
    then if null w
      then return a
      else return mapcar(a,function demonise);
normal:
  z:=algfactor2(u,aa);
  monicpart:=w;
  if null cdr z or (car z neq 'partial)
    then if null w
      then return z
      else return mapcar(z,function demonise);
  % does not factor.
  if null w
    then return mapcan(cdr z,function algfactor)
    else return for each u in z conc
                  algfactor demonise u;
  end;


symbolic procedure demonise u;
% Replaces afactorvar by afactorvar*monicpart in u.
if atom u
  then u
  else if afactorvar eq mvar u
    then addf(demonise red u,
              !*multf(lt u .+ nil,exptf(monicpart,ldeg u)))
    else if ordop(afactorvar,mvar u)
      then u
      else addf(demonise red u,
                !*multf(!*p2f lpow u,demonise lc u));


symbolic procedure makemonic u;
% U is a makemainvar'd polynomial.
begin
  scalar v,w,x,xx;
  v:=mvar u;
  x:=lc u;
  xx:=1;
  w:=!*p2f lpow u;% the monic term.
  u:=red u;
  for i:=(isub1 ldeg w) step -1 until 1 do begin
    if atom u
      then go to next;
    if mvar u neq v
      then go to next;
    if ldeg u iequal i
      then w:=addf(w,!*multf(lc u,
                     !*multf(!*p2f lpow u,xx)));
    u:=red u;
  next:
    xx:=!*multf(x,xx)
    end;
  w:=addf(w,!*multf(u,xx));
  return w
  end;

% unfluid '(monicpart);

endmodule;


module algfn;

% Author: James H. Davenport.

% Check if an expression is in a pure algebraic extension of
% Q(all "constants")(var).


exports algfnpl,algebraicsf;

imports simp,interr,dependsp,dependspl;

symbolic procedure algfnp(pf,var);
   if atom pf then t
    else if not atom car pf then interr "Not prefix form"
    else if car pf eq '!*sq then algfnsq(cadr pf,var)
      else if car pf eq 'expt
       then if not algint!-ratnump caddr pf
              then (not dependsp(cadr pf,var))
                and (not dependsp(caddr pf,var))
             else algfnp(cadr pf,var)
    else if not memq(car pf,'(minus plus times quotient sqrt))
           % JPff fiddle
     then not dependspl(cdr pf,var)
    else algfnpl(cdr pf,var);

symbolic procedure algfnpl(p!-list,var);
   null p!-list or algfnp(car p!-list,var) and algfnpl(cdr p!-list,var);

symbolic procedure algfnsq(sq,var);
   algfnsf(numr sq,var) and algfnsf(denr sq,var);

symbolic procedure algfnsf(sf,var);
   atom sf
 or algfnp(mvar sf,var) and algfnsf(lc sf,var) and algfnsf(red sf,var);

symbolic procedure algint!-ratnump q;
   if atom q then numberp q
    else car q eq 'quotient and (numberp cadr q) and (numberp caddr q);

symbolic procedure algebraicsf u;
   if atom u then nil
    else algebraicp mvar u or algebraicsf lc u or algebraicsf red u;

symbolic procedure algebraicp u;
   if atom u then nil
    else if car u eq 'expt then 1 neq denr simp caddr u
    else car u eq 'sqrt;

endmodule;


module algnums;

% Author: James H. Davenport.

exports denr!-algno;


symbolic procedure denr!-algno u;
% Returns the true denominator of the algebraic number u.
begin
  scalar sqlist,n,m,u!*!*j,d,isub1n;
  u!*!*j:=1 ./ 1;
  sqlist:=sqrtsinsq(u,nil);
  sqlist:=multbyallcombinations(list(1 ./ 1),
                               for each u in sqlist
                                 collect !*kk2q u);
  n:=0;
  sqlist:=for each u in sqlist collect
    (numr u) . (n:=iadd1 n);
    % format is of an associtaion list.
  n:=length sqlist;
  m:=mkvect n;
  isub1n:=isub1 n;
  for i:=0:n do
    putv(m,i,mkvect2(n,nil ./ 1));
  putv(getv(m,0),cdr assoc(1,sqlist),1 ./ 1);
  % initial matrix is now set up.
  for j:=1:n do begin
    scalar v,w;
    u!*!*j:=!*multsq(u!*!*j,u);
    dump!-sqrts!-coeffs(u!*!*j,sqlist,getv(m,j));
    v:=firstlinearrelation(m,n);
    if null v
      then return;
    if last!-non!-zero v > j
      then return;
    if (w:=getv(v,j)) neq (1 ./ 1)
      then <<
        w:=!*invsq w;
        for i:=0:j do
          putv(v,i,!*multsq(w,getv(v,i))) >>;
    m:=v;
    n:=j;
    return
    end;
  % Now m is a monic polynomial, minimal for u, of degree n.
  d:=1;
  for i:=0:isub1 n do begin
    scalar v,prime;
    v:=denr getv(m,i);
    prime:=2;
loop:
    if v = 1
      then return;
    if not zerop cdr divide(v,prime)
      then prime:=nextprime(prime)
      else <<
        d:=d*prime;
        for i:=0:n do
          putv(v,i,multsq(getv(v,i),1 ./ (prime ** (n-i)) )) >>;
    go to loop;
    end;
  return d;
  end;


symbolic procedure dump!-sqrts!-coeffs(u,sqlist,vec);
begin
  scalar w;
  dump!-sqrts!-coeffs2(numr u,sqlist,vec,1);
  u:=1 ./ denr u;
  if denr u neq 1
    then for i:=0:upbv vec do
      if numr(w:=getv(vec,i))
        then putv(vec,i,!*multsq(u,w));
  end;


symbolic procedure dump!-sqrts!-coeffs2(u,sqlist,vec,sqrtssofar);
if null u
  then nil
  else if numberp u
    then putv(vec,cdr assoc(sqrtssofar,sqlist),u)
    else <<
      dump!-sqrts!-coeffs2(red u,sqlist,vec,sqrtssofar);
      dump!-sqrts!-coeffs2(lc u,sqlist,vec,!*multf(sqrtssofar,
                                                   !*k2f mvar u)) >>;


symbolic procedure last!-non!-zero vec;
begin
  scalar n;
  for i:=0:upbv vec do
    if numr getv(vec,i)
      then n:=i;
  return n
  end;

endmodule;


module antisubs;

% Author: James H. Davenport.

exports antisubs;

imports purge,interr,dependsp;


symbolic procedure antisubs(place,x);
% Produces the inverse substitution to a substitution list.
begin
  scalar answer,w;
  while place and
        (x=caar place) do<<
    w:=cdar place;
    % w is the substitution rule.
    if atom w
      then if w neq x
        then interr "False atomic substitution"
        else nil
      else answer:=(x.anti2(w,x)).answer;
    place:=cdr place>>;
  if null answer
    then answer:=(x.x).answer;
  return answer
  end;


symbolic procedure anti2(eexpr,x);
%Produces the function inverse to the eexpr provided.
if atom eexpr
  then if eexpr eq x
    then x
    else interr "False atom"
  else if car eexpr eq 'plus
    then deplus(cdr eexpr,x)
    else if car eexpr eq 'minus
      then subst(list('minus,x),x,anti2(cadr eexpr,x))
      else if car eexpr eq 'quotient
        then if dependsp(cadr eexpr,x)
          then if dependsp(caddr eexpr,x)
            then interr "Complicated division"
            else subst(list('times,caddr eexpr,x),x,anti2(cadr eexpr,x))
          else if dependsp(caddr eexpr,x)
            then subst(list('quotient,cadr eexpr,x),x,
                       anti2(caddr eexpr,x))
            else interr "No division"
        else if car eexpr eq 'expt
          then if caddr eexpr iequal 2
            then subst(list('sqrt,x),x,anti2(cadr eexpr,x))
            else interr "Unknown root"
          else if car eexpr eq 'times
            then detimes(cdr eexpr,x)
            else if car eexpr eq 'difference
              then deplus(list(cadr eexpr,list('minus,caddr eexpr)),x)
              else interr "Unrecognised form in antisubs";



symbolic procedure detimes(p!-list,var);
% Copes with lists 'times.
begin
  scalar u,v;
  u:=deplist(p!-list,var);
  v:=purge(u,p!-list);
  if null v
    then v:=var
    else if null cdr v
      then v:=list('quotient,var,car v)
      else v:=list('quotient,var,'times.v);
  if (null u) or
     (cdr u)
    then interr "Weird multiplication";
  return subst(v,var,anti2(car u,var))
  end;


symbolic procedure deplist(p!-list,var);
% Returns a list of those elements of p!-list which depend on var.
if null p!-list
  then nil
  else if dependsp(car p!-list,var)
    then (car p!-list).deplist(cdr p!-list,var)
    else deplist(cdr p!-list,var);


symbolic procedure deplus(p!-list,var);
% Copes with lists 'plus.
begin
  scalar u,v;
  u:=deplist(p!-list,var);
  v:=purge(u,p!-list);
  if null v
    then v=var
    else if null cdr v
      then v:=list('plus,var,list('minus,car v))
      else v:=list('plus,var,list('minus,'plus.v));
  if (null u) or
     (cdr u)
    then interr "Weird addition";
  return subst(v,var,anti2(car u,var))
  end;

endmodule;


module coates;
 
% Author: James H. Davenport.
 
fluid '(intvar magiclist nestedsqrts previousbasis sqrt!-intvar
        taylorasslist thisplace);
 
global '(!*tra !*trmin coates!-fdi);
 
exports coates,makeinitialbasis,checkpoles,multbyallcombinations;
 
 
 
 
 
symbolic procedure coates(places,mults,x);
begin
  scalar u,tt;
  tt:=readclock();
  u:=coates!-hpfsd(places,mults);
  if !*tra or !*trmin then
    printc  list ('coates,'time,readclock()-tt,'milliseconds);
  return u
  end;
 
 
 
symbolic procedure coates!-real(places,mults);
begin
  scalar thisplace,u,v,save;
  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) >>;
%  v:=placesindiv places;
    % V is a list of all the substitutors in PLACES;
%  u:=mkunique sqrtsintree(v,intvar,nil);
%  if !*tra then <<
%    princ "Sqrts on this curve:";
%    terpri!*(t);
%    superprint u >>;
%  algnos:=mkunique mapcar(places,function basicplace);
%  if !*tra then <<
%    printc "Algebraic numbers where residues occur:";
%    superprint algnos >>;
  v:=mults;
  for each uu in places do <<
    if (car v) < 0
      then u:=(rfirstsubs uu).u;
    v:=cdr v >>;
  thisplace:=list('quotient,1,intvar);
  if member(thisplace,u)
    then <<
      v:= finitise(places,mults);
      % returns list (places,mults,power of intvar to remove.
      u:=coates!-real(car v,cadr v);
      if atom u
        then return u;
      return multsq(u,!*p2q mksp(intvar,caddr v)) >>;
% It is not sufficient to check the current value of U in FRACTIONAL...
% as we could have zeros over infinity JHD 18/8/86;
  for each uu in places do
    if rfirstsubs uu = thisplace
      then u:=append(u,mapcar(cdr uu,function car));
  coates!-fdi:=fractional!-degree!-at!-infinity u;
% Do we need to blow everything up by a factor of two (or more)
% to avoid fractional powers at infinity?
  if coates!-fdi iequal 1
    then return coatesmodule(places,mults,intvar);
  if !*tra
    then fdi!-print();
  places:=mapcar(places,function fdi!-upgrade);
  save:=taylorasslist;
  u:=coatesmodule(places,
    mapcar(mults,function (lambda u;u*coates!-fdi)),
                  intvar);
  taylorasslist:=save;
% u:=fdi!-revertsq u;
% That previous line is junk, I think (JHD 22.8.86)
% just because we blew up the places doesn't mean that
% we should deflate the function, because that has already been done.
  return u
  end;
 
 
 
symbolic procedure coatesmodule(places,mults,x);
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 >>;
  % ***time-hack-2***;
  if previousbasis then basis:=previousbasis
    else 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,nil);
    % The NIL is the list of special constraints needed
    % to force certain poles to occur in the answer.
  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 "These values give the function";
    printsq u >>;
  magiclist:=nil;
  if checkpoles(list u,places,mults)
    then return u
    else interr "Inconsistent checkpoles"
  end;
 
 
 
symbolic procedure makeinitialbasis places;
begin
  scalar u;
  u:=multbyallcombinations(list(1 ./ 1),
                           for each u in getsqrtsfromplaces places
         collect !*kk2q u);
  if !*tra then <<
    printc "Initial basis for the space m(x)";
    mapc(u,function printsq) >>;
  return u
  end;
 
 
 
symbolic procedure multbyallcombinations(u,l);
% Produces a list of all elements of u,
% each multiplied by every combination of elements of l.
if null l
  then u
  else multbyallcombinations(nconc(multsql(car l,u),u),cdr l);
 
 
 
symbolic procedure checkpoles(basis,places,mults);
% Checks that the BASIS really does have all the
%  poles in (PLACES.MULTS).
begin
  scalar u,v,l;
  go to outer2;
outer:
  places:=cdr places;
  mults:=cdr mults;
outer2:
  if null places
    then return if magiclist
                  then findmagic l
                  else t;
  if 0 leq car mults
    then go to outer;
  u:=basis;
inner:
  if null u
    then <<
      if !*tra
        then <<
          princ "The answer from the linear equations did";
          printc " not have the poles at:";
          printplace car places >>;
      return nil >>;
  v:=taylorform xsubstitutesq(car u,car places);
  if taylorfirst v=car mults
    then <<
      if magiclist
        then l:=taylorevaluate(v,car mults) . l;
      go to outer >>;
  if taylorfirst v < car mults
    then interr "Extraneous pole introduced";
  u:=cdr u;
  go to inner
  end;
 
 
 
symbolic procedure coates!-hpfsd(oplaces,omults);
begin
  scalar mzero,pzero,mpole,ppole,fun,summzero,answer,places,mults;
  places:=oplaces;
  mults:=omults;
  % Keep originals in case need to use COATES!-REAL directly.
  summzero:=0;
    % holds the sum of all the mzero's.
  while places do <<
    if 0<car mults
      then <<
        summzero:=summzero + car mults;
        mzero:=(car mults).mzero;
        pzero:=(car places).pzero >>
      else <<
        mpole:=(car mults).mpole;
        ppole:=(car places).ppole >>;
    places:=cdr places;
    mults:=cdr mults >>;
  if summzero > 2 then begin
    % We want to combine a zero/pole pair
    % so as to reduce the total index before calling coates!-real
    % on the remaining zeros/poles.
    scalar nplaces,nmults,f,multiplicity,newpole,sqrts,fz,zfound,mult1;
    sqrts:=getsqrtsfromplaces ppole;
    if !*tra or !*trmin then <<
      princ "Operate on divisor:";
      printc append(mzero,mpole);
      printc "at";
      mapc(pzero,function printplace);
      mapc(ppole,function printplace) >>;
iterate:
    nplaces:=list car pzero;
    multiplicity:=car mzero;
    nmults:=list 1;
    if cdr ppole
      then <<
        nplaces:=(car ppole) . ( (cadr ppole) . nplaces);
        multiplicity:=min(multiplicity,- car mpole,- cadr mpole);
        nmults:=(-1) .((-1) . nmults) >>
      else <<
        nplaces:=(car ppole) . nplaces;
        multiplicity:=min(multiplicity,(- car mpole)/2);
        nmults:=(-2) . nmults >>;
    previousbasis:=nil;
    f:=coates!-real(nplaces,nmults);
    if atom f
      then <<
 if !*tra or !*trmin then
          printc "Failure: must try whole divisor";
 return coates!-real(oplaces,omults) >>;
%    newpole:=removezero(findzeros(f,sqrts),car pzero).
    fz:=findzeros(f,sqrts);
    zfound:=assoc(car pzero,fz);
    if not zfound
       then interr "Didn't seem to find the zeros we looked for";
    if cdr zfound > car mzero
       then interr "We found too many zeros";
    fz:=delete(zfound,fz);
    if !*tra or !*trmin then <<
      printc "Replaced by the pole";
      if fz then prettyprint fz
       else <<terpri(); prin2t "The zero we were already looking for">>;
      princ multiplicity;
      printc " times" >>;
    mult1:=car mzero - multiplicity * cdr zfound;
    if mult1 < 0
 then << printc "A zero has turned into a pole";
         multiplicity:= car mzero / cdr zfound ;
  mult1:=remainder(car mzero, cdr zfound); >>;
    if zerop mult1
      then <<
        mzero:=cdr mzero;
        pzero:=cdr pzero >>
      else rplaca(mzero,mult1);
    if null cdr ppole
      then <<
        if zerop (car mpole + 2*multiplicity)
          then <<
            ppole:=cdr ppole;
     mpole:=cdr mpole >>
          else rplaca(mpole,car mpole + 2 * multiplicity) >>
      else <<
        if zerop (cadr mpole + multiplicity)
          then <<
            ppole:=(car ppole) . (cddr ppole);
            mpole:=(car mpole) . (cddr mpole) >>
          else rplaca(cdr mpole,cadr mpole + multiplicity);
        if zerop (car mpole + multiplicity)
          then <<
            ppole:=cdr ppole;
            mpole:=cdr mpole >>
          else rplaca(mpole,car mpole + multiplicity) >>;
    while fz do <<
      newpole:=caar fz;
      mult1:=multiplicity*(cdar fz);
      if newpole member pzero
        then begin
          scalar m,p;
          while newpole neq car pzero do <<
            m:=(car mzero).m;
            mzero:=cdr mzero;
            p:=(car pzero).p;
            pzero:=cdr pzero >>;
          if mult1 < car mzero then <<
            mzero:=(car mzero - mult1) . cdr mzero;
            mzero:=nconc(m,mzero);
            pzero:=nconc(p,pzero);
            return >>;
          if mult1 > car mzero then <<
            ppole:=newpole.ppole;
            mpole:=(car mzero - mult1) . mpole >>;
          mzero:=nconc(m,cdr mzero);
          pzero:=nconc(p,cdr pzero)
          end
        else if newpole member ppole then begin
          scalar m,p;
          m:=mpole;
          p:=ppole;
          while newpole neq car p do <<
            p:=cdr p;
            m:=cdr m >>;
          rplaca(m,car m - mult1)
          end
        else <<
          mpole:=nconc(mpole,list(-mult1));
          ppole:=nconc(ppole,list newpole) >>;
      fz:=cdr fz >>;
    f:=mk!*sq f;
    if multiplicity > 1
      then answer:=list('expt,f,multiplicity).answer
      else answer:=f.answer;
    summzero:=0;
    for each x in mzero do summzero:=summzero+x;
    if !*tra then <<
      princ "Function is now: ";
      printc append(mzero,mpole);
      printc "at";
      mapc(pzero,function printplace);
      mapc(ppole,function printplace) >>;
    if summzero > 2
      then go to iterate;
    end;
  fun:=coates!-real(nconc(pzero,ppole),
                    nconc(mzero,mpole));
  if null answer
    then return fun
    else answer:=(mk!*sq fun).answer;
  return !*k2q('times.answer);
    % This is not valid, but we hope that it will be unpicked;
    % (e.g. by SIMPLOG) before too much damage is caused.
  end;
 
 
 
symbolic procedure removezero(l,place);
if place member l
  then (lambda u; if null cdr u
                    then car u
                    else interr "Removezero") delete(place,l)
  else interr "Error in removezeros";
 
 
 
symbolic procedure findzeros(sq,sqrts);
begin
  scalar u,potentials,answer,n;
  u:=denr sqrt2top invsq sq;
  potentials:=for each v in jfactor(u,intvar) collect begin
    scalar w,place;
    w:=makemainvar(numr v,intvar);
    if ldeg w neq 1
      then interr "Can't cope";
    if red w
      then place:=list(intvar,'plus,intvar,prepsq(negf red w ./ lc w))
      else place:=intvar . intvar;
      % This IF .. ELSE .. added JHD 3 Sept 1980.
    return place
    end;
  potentials:=list(intvar,'quotient,1,intvar).potentials;
  for each place in potentials do begin
    scalar slist,nestedsqrts;
    place:=list place;
    newplace place;
    u:=substitutesq(sq,place);
    while involvesq(u,sqrt!-intvar) do begin
      scalar z;
      z:=list list(intvar,'expt,intvar,2);
      place:=nconc(place,z);
      newplace place;
      u:=substitutesq(u,z);
      end;
    slist:=sqrtsinsq(u,intvar);
    for each v in sqrts do
      slist:=union(slist,sqrtsinsq(xsubstitutesq(!*kk2q v,place),
       intvar));
    slist:=sqrtsign(slist,intvar);
    for each s in slist do
      if (n:=taylorfirst taylorform substitutesq(u,s)) > 0
        then answer:=(append(place,s).n).answer;
    return answer;
    end;
  if null answer
    then interr "No zero found";
  return answer
  end;
 
endmodule;


module coatesid;
 
% Author: James H. Davenport.
 
fluid '(intvar magiclist nnn taylorasslist taylorvariable);
 
global '(!*tra);
 
exports coatessolve,vecprod,coates!-lineq;
 
imports !*invsq,!*multsq,negsq,!*addsq,swap,check!-lineq,non!-null!-vec,
 printsq,sqrt2top,mapvec,mksp,vecsort,addsq,mkilist,mkvec,mapply,
 taylorformp,xsubstitutesq,taylorform,taylorevaluate,multsq,
 invsq,removecmsq;
 
symbolic procedure coatessolve(mzero,pzero,basis,normals);
begin
  scalar m,n,rightside,nnn;
% if null normals
%   then normals:=list mkilist(basis,1 ./ 1);
%     This provides the default normalisation,
%     viz merely a de-homogenising constraint;
% No it doesn't - JHD May 1983 and August 1986.
% This may be precisely the wrong constraint, as can be seen from
% the example of SQRT(X**2-1).  Fixed 19/8/86 to amend COATESMATRIX
% to insert a normalising constraint if none is provided.
  nnn:=max(length normals,1);
  basis:=mkvec basis;
  m:=coatesmatrix(mzero,pzero,basis,normals);
  n:=upbv m;
  rightside:=mkvect n;
  for i:=0:n do
    putv(rightside,n-i,(if i < nnn
                       then 1
                       else nil) ./ 1);
  n:=coates!-lineq(m,rightside);
  if n eq 'failed
    then return 'failed;
  n:=removecmsq vecprod(n,basis);
  if !*tra then <<
    printc "Answer from linear equation solving is ";
    printsq n >>;
  return n
  end;
 
 
 
symbolic procedure coatesmatrix(mzero,pzero,basis,normals);
% NORMALS is a list of the normalising constraints
% that we must apply.  Thypically, this is NIL, and we have to
% invent one - see the code IF NULL NORMALS ...
begin
  scalar ans,n1,n2,j,w,save,nextflag,save!-taylors,x!-factors,
         normals!-ok,temp;
  save!-taylors:=mkvect isub1 length pzero;
  save:=taylorasslist;
  normals!-ok:=nil;
  n1:=upbv basis;
  n2:=isub1 mapply(function plus2,mzero) + max(length normals,1);
    % the number of constaints in all (counting from 0).
  taylorvariable:=intvar;
  if !*tra then <<
    printc "Basis for the functions with precisely the correct poles";
    mapvec(basis,function printsq) >>;
  ans:=mkvect n2;
  for i:=0:n2 do
    putv(ans,i,mkvect n1);
  for i:=0:n1 do begin
    scalar xmz,xpz,k;
    xmz:=mzero;
    k:=j:=0;
    xpz:=pzero;
    while xpz do <<
      newplace basicplace car xpz;
      if nextflag
        then w:=taylorformp list('binarytimes,
     getv(save!-taylors,k),
     getv(x!-factors,k))
 else if not !*tra
          then w:=taylorform xsubstitutesq(getv(basis,i),car xpz)
          else begin
            scalar flg,u,slists;
            u:=xsubstitutesq(getv(basis,i),basicplace car xpz);
            slists:=extenplace car xpz;
            for each w in sqrtsinsq(u,intvar) do
              if not assoc(w,slists)
                then flg:=w.flg;
            if flg then <<
              printc "The following square roots were not expected";
              mapc(flg,function superprint);
              printc "in the substitution";
              superprint car xpz;
              printsq getv(basis,i) >>;
            w:=taylorform xsubstitutesq(u,slists)
            end;
      putv(save!-taylors,k,w);
      k:=iadd1 k;
      for l:=0 step 1 until isub1 car xmz do <<
        astore(ans,j,i,taylorevaluate(w,l));
        j:=iadd1 j >>;
      if null normals and j=n2 then <<
 temp:=taylorevaluate(w,car xmz);
 astore(ans,j,i,temp);
 % The defaults normalising condition is that the coefficient
 % after the last zero be a non-zero.
 % Unfortunately this too may fail (JHD 21.3.87) - check for it later
 normals!-ok:=normals!-ok or numr temp >>;
      xpz:=cdr xpz;
      xmz:=cdr xmz  >>;
    nextflag:=(i < n1) and
       (getv(basis,i) = multsq(!*kk2q intvar,getv(basis,i+1)));
    if nextflag and null x!-factors then <<
      x!-factors:=mkvect upbv save!-taylors;
      xpz:=pzero;
      k:=0;
      xmz:=invsq !*kk2q intvar;
      while xpz do <<
 putv(x!-factors,k,taylorform xsubstitutesq(xmz,car xpz));
        xpz:=cdr xpz;
        k:=iadd1 k >> >>
    end;
  if null normals and null normals!-ok then <<
     if !*tra
       then printc "Our default normalisation condition was vacuous";
     astore(ans,n2,n1,1 ./ 1)>>;
  while normals do <<
    w:=car normals;
    for k:=0:n1 do <<
      astore(ans,j,k,car w);
      w:=cdr w >>;
    j:=iadd1 j;
    normals:=cdr normals >>;
  tayshorten save;
  return ans
  end;
 
 
symbolic procedure printmatrix(ans,n2,n1);
if !*tra
  then <<
    printc "Equations to be solved:";
    for i:=0:n2 do begin
      if null getv(ans,i)
        then return;
      princ "Row number ";
      princ i;
      for j:=0:n1 do
        printsq getv(getv(ans,i),j)
      end >>;
 
 
 
symbolic procedure vecprod(u,v);
begin
  scalar w,n;
  w:=nil ./ 1;
  n:=upbv u;
  for i:=0:n do
    w:=addsq(w,!*multsq(getv(u,i),getv(v,i)));
  return w
  end;
 
 
 
symbolic procedure coates!-lineq(m,rightside);
begin
  scalar nnn,n;
  nnn:=desparse(m,rightside);
  if nnn eq 'failed
    then return 'failed;
  m:=car nnn;
  if null m
    then <<
      n:=cddr nnn;
      goto vecprod >>;
  rightside:=cadr nnn;
  nnn:=cddr nnn;
  n:=check!-lineq(m,rightside);
  if n eq 'failed
    then return n;
  n:=jhdsolve(m,rightside,non!-null!-vec nnn);
  if n eq 'failed
    then return n;
  for i:=0:upbv n do
    if (m:=getv(nnn,i))
      then putv(n,i,m);
vecprod:
  for i:=0:upbv n do
    if null getv(n,i) then putv(n,i,nil ./ 1);
  return n
  end;
 
 
 
symbolic procedure jhdsolve(m,rightside,ignore);
% Returns answer to m.answer=rightside.
% Matrix m not necessarily square.
begin
  scalar n1,n2,ans,u,row,swapflg,swaps;
  % The SWAPFLG is true if we have changed the order of the
  % columns and need later to invert this via SWAPS.
  n1:=upbv m;
  for i:=0:n1 do
    if (u:=getv(m,i))
      then (n2:=upbv u);
  printmatrix(m,n1,n2);
  swaps:=mkvect n2;
  for i:=0:n2 do
    putv(swaps,i,n2-i);
    % We have the SWAPS vector, which should be a vector of indices,
    % arranged like this because VECSORT sorts in decreasing order.
  for i:=0:isub1 n1 do begin
    scalar k,v,pivot;
  tryagain:
    row:=getv(m,i);
    if null row
      then go to interchange;
    % look for a pivot in row.
    k:=-1;
    for j:=0:n2 do
      if numr (pivot:=getv(row,j))
        then <<
          k:=j;
          j:=n2 >>;
    if k neq -1
      then goto newrow;
    if numr getv(rightside,i)
      then <<
        m:='failed;
        i:=sub1 n1; %Force end of loop.
        go to finished >>;
    % now interchange i and last element.
interchange:
    swap(m,i,n1);
    swap(rightside,i,n1);
    n1:=isub1 n1;
    if i iequal n1
      then goto finished
      else goto tryagain;
  newrow:
    if i neq k
      then <<
        swapflg:=t;
        swap(swaps,i,k);
          % record what we have done.
        for l:=0:n1 do
          swap(getv(m,l),i,k) >>;
        % place pivot on diagonal.
    pivot:=sqrt2top negsq !*invsq pivot;
    for j:=iadd1 i:n1 do begin
      u:=getv(m,j);
      if null u
        then return;
      v:=!*multsq(getv(u,i),pivot);
      if numr v then <<
        putv(rightside,j,
     !*addsq(getv(rightside,j),!*multsq(v,getv(rightside,i))));
        for l:=0:n2 do
   putv(u,l,!*addsq(getv(u,l),!*multsq(v,getv(row,l)))) >>
      end;
  finished:
    end;
  if m eq 'failed
    then go to failed;
    % Equations were inconsistent.
  while null (row:=getv(m,n1)) do
    n1:=isub1 n1;
  u:=nil;
  for i:=0:n2 do
    if numr getv(row,i)
      then u:='t;
  if null u
    then if numr getv(rightside,n1)
      then go to failed
      else n1:=isub1 n1;
      % Deals with a last equation which is all zero.
  if n1 > n2
    then go to failed;
    % Too many equations to satisfy.
  ans:=mkvect n2;
  n2:=n2 - ignore;
  if n1 < n2 then <<
    if !*tra then <<
      printc "The equations do not completely determine the functions";
      printc "Matrix:";
      mapvec(m,function superprint);
      printc "Right-hand side:";
      superprint rightside >>;
    for i:=iadd1 n1:n2 do <<
      u:=gensym();
      magiclist:=u.magiclist;
      putv(ans,i,!*kk2q u) >>;
    if !*tra then printc "If in doubt consult an expert">>;
  % now to do the back-substitution.
  for i:=n1 step -1 until 0 do begin
    row:=getv(m,i);
    if null row
      then return;
    u:=getv(rightside,i);
    for j:=iadd1 i:n2 do
      u:=!*addsq(u,!*multsq(getv(row,j),negsq getv(ans,j)));
    putv(ans,i,!*multsq(u,sqrt2top !*invsq getv(row,i)))
    end;
  if swapflg
    then vecsort(swaps,list ans);
  return ans;
failed:
  if !*tra then printc "Unable to force correct zeroes";
  return 'failed
  end;
 
 
 
symbolic procedure desparse(matrx,rightside);
begin
  scalar vec,changed,n,m,zero,failed;
  zero := nil ./ 1;
  n:=upbv matrx;
  m:=upbv getv(matrx,0);
  vec:=mkvect m;
  % for i:=0:m do putv(vec,i,zero);   %%% initialize - ach
  changed:=t;
  while changed and not failed do begin
    changed:=nil;
    for i:=0:n do
      if changed or failed
        then i:=n   % and hence quit the loop.
        else begin
          scalar nzcount,row,pivot;
          row:=getv(matrx,i);
          if null row
            then return;
          nzcount:=0;
          for j:=0:m do
            if numr getv(row,j)
              then <<
                nzcount:=iadd1 nzcount;
                pivot:=j >>;
          if nzcount = 0
            then if null numr getv(rightside,i)
              then return putv(matrx,i,nil)
              else return (failed:='failed);
          if nzcount > 1
            then return nil;
          nzcount:=getv(rightside,i);
          if null numr nzcount
            then <<
              putv(vec,pivot,zero);
       go to was!-zero >>;
   nzcount:=!*multsq(nzcount,!*invsq getv(row,pivot));
          putv(vec,pivot,nzcount);
          nzcount:=negsq nzcount;
          for i:=0:n do
            if (row:=getv(matrx,i))
              then if numr (row:=getv(row,pivot))
  then putv(rightside,i,!*addsq(getv(rightside,i),
      !*multsq(row,nzcount)));
was!-zero:
          for i:=0:n do
            if (row:=getv(matrx,i))
              then putv(row,pivot,zero);
          changed:=t;
          putv(matrx,i,nil);
          swap(matrx,i,n);
          swap(rightside,i,n);
          end;
    end;
  if failed
    then return 'failed;
  changed:=t;
  for i:=0:n do
    if getv(matrx,i)
      then changed:=nil;
  if changed
    then matrx:=nil;
    % We have completely solved the equations by these machinations.
  return matrx.(rightside.vec)
  end;
 
 
symbolic procedure astore(a,i,j,val);
   putv(getv(a,i),j,val);
 
endmodule;


module findmagc;

% Author: James H. Davenport.

fluid '(magiclist);

global '(!*tra);

symbolic procedure findmagic l;
begin
  scalar p,n,pvec,m,intvec,mcount,temp;
  % L is a list of things which must be made non-zero by means of
%   a suitable choice of values for the variables in MAGICLIST;
  l:=for each u in l collect
     << mapc(magiclist,function (lambda v;
                                 if involvesf(denr u,v)
                                   then interr "Hard findmagic"));
        numr u >>;
  if !*tra then <<
    printc "We must make the following non-zero:";
    mapc(l,function printsf);
    princ "by suitable choice of ";
    printc magiclist >>;
  % Strategy is random choice in a space which has only finitely
%   many singular points;
  p:=0;
  n:=isub1 length magiclist;
  pvec:=mkvect n;
  putv(pvec,0,2);
  for i:=1:n do
    putv(pvec,i,nextprime getv(pvec,isub1 i));
  % Tactics are based on Godel (is this a mistake ??) and let P run
%   through numbers and take the prime factorization of them;
  intvec:=mkvect n;
loop:
  p:=iadd1 p;
  if !*tra then <<
    princ "We try the number ";
    printc p >>;
  m:=p;
  for i:=0:n do <<
    mcount:=0;
    while zerop cdr (temp:=divide(m,getv(pvec,i)) ) do <<
      mcount:=iadd1 mcount;
      m:=car temp >>;
    putv(intvec,i,mcount) >>;
  if m neq 1
    then go to loop;
  if !*tra then <<
    printc "which corresponds to ";
    superprint intvec >>;
  m:=nil;
  temp:=magiclist;
  for i:=0:n do <<
    m:=((car temp).getv(intvec,i)).m;
    temp:=cdr temp >>;
  % M is the list of substitutions corresponding to this value of P;
  temp:=l;
loop2:
  if null numr algint!-subf(car temp,m)
    then go to loop;
  temp:=cdr temp;
  if temp
    then go to loop2;
  if !*tra then <<
    printc "which corresponds to the values:";
    superprint m >>;
  return m
  end;

endmodule;


module findres;

% Author: James H. Davenport.

fluid '(!*gcd
        basic!-listofallsqrts
        basic!-listofnewsqrts
        intvar
        listofallsqrts
        listofnewsqrts
        nestedsqrts
        sqrt!-intvar
        taylorvariable);

global '(!*tra !*trmin);

exports find!-residue,findpoles;
imports sqrt2top,jfactor,prepsq,printplace,simpdf,involvesf,simp;
imports stt,interr,mksp,negf,multf,addf,let2,substitutesq,subs2q,quotf;
imports printsq,clear,taylorform,taylorevaluate,involvesf,!*multsq;
imports sqrtsave,sqrtsinsq,sqrtsign;

symbolic procedure find!-residue(simpdl,x,place);
  % evaluates residue of simpdl*dx at place given by x=place(y).
begin
  scalar deriv,nsd,poss,slist;
  listofallsqrts:=basic!-listofallsqrts;
  listofnewsqrts:=basic!-listofnewsqrts;
  deriv:=simpdf(list(place,x));
  if involvesf(numr deriv,intvar)
    then return residues!-at!-new!-point(simpdl,x,place);
  if eqcar(place,'quotient) and (cadr place iequal 1)
    then goto place!-correct;
  place:=simp list('difference,intvar,place);
  if involvesq(place,intvar)
    then interr "Place wrongly formatted";
  place:=list('plus,intvar,prepsq place);
place!-correct:
  if car place eq 'plus and caddr place = 0
    then place:=list(x.x)
    else place:=list(x.place);
  % the substitution required.
  nsd:=substitutesq(simpdl,place);
  deriv:=!*multsq(nsd,deriv);
  % differential is deriv * dy, where x=place(y).
  if !*tra then <<
    printc "Differential after first substitution is ";
    printsq deriv >>;
  while involvesq(deriv,sqrt!-intvar)
    do <<
      sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,place);
      nsd:=list(list(x,'expt,x,2));
      deriv:=!*multsq(substitutesq(deriv,nsd),!*kk2q x);
      % derivative of x**2 is 2x, but there's a jacobian of 2 to
      % consider.
      place:=nconc(place,nsd) >>;
  % require coeff x**-1 in deriv.
  nestedsqrts:=nil;
  slist:=sqrtsinsq(deriv,x);
  if !*tra and nestedsqrts
    then printc "We have nested square roots";
  slist:=sqrtsign(slist,intvar);
  % The reversewoc is to ensure that the simpler sqrts are at
  % the front of the list.
  % Slist is a list of all combinations of signs of sqrts.
  taylorvariable:=x;
  for each branch in slist do <<
    nsd:=taylorevaluate(taylorform substitutesq(deriv,branch),-1);
    if numr nsd
      then poss:=(append(place,branch).nsd).poss >>;
  poss:=reversewoc poss;
  if null poss
    then go to finished;
  % poss is a list of all possible residues at this place.
  if !*tra
    then <<
      princ "Residues at ";
      printplace place;
      printc " are ";
      mapc(poss, function (lambda u; <<
                       printplace car u;
                       printsq cdr u >>)) >>;
finished:
  sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,place);
  return poss
  end;


symbolic procedure residues!-at!-new!-point(func,x,place);
begin
  scalar place2,tempvar,topterm,a,b,xx;
  if !*tra then <<
    printc "Find residues at all roots of";
    superprint place >>;
  place2:=numr simp place;
  topterm:=stt(place2,x);
  if car topterm = 0
    then interr "Why are we here?";
  tempvar:=gensym();
  place2:=addf(place2,
               multf(!*p2f mksp(x,car topterm),negf cdr topterm));
  % The remainder of PLACE2.
  let2(list('expt,tempvar,car topterm),
       subst(tempvar,x,prepsq(place2 ./ cdr topterm)),
       nil,t);
  place2:=list list(x,'plus,x,tempvar);
  !*gcd:=nil;
    % No unnecessary work: only factors of X worry us.
  func:=subs2q substitutesq(func,place2);
  !*gcd:=t;
  xx:=!*k2f x;
  while (a:=quotf(numr func,xx)) and (b:=quotf(denr func,xx))
    do func:=a ./ b;
  if !*tra then <<
    printc "which gives rise to ";
    printsq func >>;
  if null a
    then b:=quotf(denr func,xx);
    % because B goes back to the last time round that WHILE loop.
  if b then go to hard;
  if !*tra then printc "There were no residues";
  clear tempvar;
  return nil;
  % *** thesis remark ***
%   This test for having an X in the denominator only works
%   because we are at a new place, and hence (remark of Trager)
%   if we have a residue at one place over this point, we must have one
%   at them all, since the places are indistinguishable;
hard:
  taylorvariable:=x;
  func:=taylorevaluate(taylorform func,-1);
  printsq func;
  interr "so far"
  end;


symbolic procedure findpoles(simpdl,x);
begin
  scalar simpdl2,poles;
  % finds possible poles of simpdl * dx.
  simpdl2:=sqrt2top simpdl;
  poles:=jfactor(denr simpdl2,x);
  poles:=mapcar(poles,function prepsq);
  % what about the place at infinity.
  poles:=list('quotient,1,x).poles;
  if !*tra or !*trmin
    then <<
      printc "Places at which poles could occur ";
      for each u in poles do
        printplace list(intvar.u) >>;
  return poles
  end;

endmodule;


module finitise;

% Author: James H. Davenport.

fluid '(intvar);

global '(!*tra);

exports finitise;
imports newplace,getsqrtsfromplaces,interr,completeplaces2,sqrtsign;
imports mkilist,extenplace;


symbolic procedure finitise(places,mults);
begin
  scalar placesmisc,multsmisc,m,n,sqrts;
  scalar places0,mults0,placesinf,multsinf;
  newplace list (intvar.intvar);
    % fix the disaster with 1/sqrt(x**2-1)
    % (but with no other 1/sqrt(x**2-k).
  sqrts:=getsqrtsfromplaces places;
  placesmisc:=places;
  multsmisc:=mults;
  n:=0;
  while placesmisc do <<
    if eqcar(rfirstsubs car placesmisc,'quotient)
        and (n > car multsmisc)
      then <<
        n:=car multsmisc;
        m:=multiplicity!-factor car placesmisc >>;
    placesmisc:=cdr placesmisc;
    multsmisc:=cdr multsmisc >>;
  if n = 0
    then interr "Why did we call finitise ??";
  % N must be corrected to allow for our representation of
  % multiplicities at places where X is not the local parameter.
  n:=divide(n,m);
  if not zerop cdr n and !*tra
    then printc
     "Cannot get the poles moved precisely because of ramification";
   if (cdr n) < 0
     then n:=(-1) + car n
     else n:=car n;
        % The above 3 lines (as a replacement for the line below)
        % inserted JHD 06 SEPT 80.
%  n:=car n;
% ***** not true jhd 06 sept 80 *****;
    % This works because, e.g., DIVIDE(-1,2) is -1 remainder 1.
    % Note that N is actually negative.
  % We now wish to divide by X**N, thus increasing
  % the degrees of all infinite places by N and
  % decreasing the degrees of all places lying over 0.
  while places do <<
    if atom rfirstsubs car places
      then <<
        places0:=(car places).places0;
        mults0:=(car mults).mults0 >>
      else if car rfirstsubs car places eq 'quotient
        then <<
          placesinf:=(car places).placesinf;
          multsinf:=(car mults).multsinf >>
        else <<
          placesmisc:=(car places).placesmisc;
          multsmisc:=(car mults).multsmisc >>;
    places:=cdr places;
    mults:=cdr mults >>;
  if places0
    then <<
      places0:=completeplaces2(places0,mults0,sqrts);
      mults0:=cdr places0;
      places0:=car places0;
      m:=multiplicity!-factor car places0;
      mults0:=for each u in mults0 collect u+n*m >>
    else <<
      places0:=for each u in sqrtsign(sqrts,intvar)
                 collect (intvar.intvar).u;
      mults0:=mkilist(places0,n * (multiplicity!-factor car places0))>>;
  placesinf:=completeplaces2(placesinf,
                             multsinf,
                             for each u in extenplace car placesinf
                               collect lsubs u);
  multsinf:=cdr placesinf;
  placesinf:=car placesinf;
  while placesinf do <<
    m:=multiplicity!-factor car placesinf;
    if (car multsinf) neq n*m
      then <<
        placesmisc:=(car placesinf).placesmisc;
        multsmisc:=(car multsinf -n*m).multsmisc >>;
      % This test ensures that we do not add places
      % with a multiplicity of zero.
    placesinf:=cdr placesinf;
    multsinf:=cdr multsinf >>;
  return list(nconc(places0,placesmisc),
              nconc(mults0,multsmisc),
              -n)
  end;


symbolic procedure multiplicity!-factor place;
begin
  scalar n;
  n:=1;
  for each u in place do
    if (lsubs u eq intvar) and
        eqcar(rsubs u,'expt)
      then n:=n*(caddr rsubs u);
  return n
  end;

endmodule;


module fixes;

% Author: James H. Davenport.

fluid '(!*nosubs asymplis!* dmode!*);

global '(ncmp!*);

% The standard version of SUBF messes with the order of variables before
% calling SUBF1, something we can't afford, so we define a new version.

symbolic procedure algint!-subf(a,b); algint!-subf1(a,b);

symbolic procedure algint!-subsq(u,v); 
   quotsq(algint!-subf(numr u,v),algint!-subf(denr u,v));

symbolic procedure algint!-subf1(u,l);
   %U is a standard form,
   %L an association list of substitutions of the form
   %(<kernel> . <substitution>).
   %Value is the standard quotient for substituted expression.
   %Algorithm used is essentially the straight method.
   %Procedure depends on explicit data structure for standard form;
   if domainp u
     then if atom u then if null dmode!* then u ./ 1 else simpatom u
          else if dmode!* eq car u then !*d2q u
          else simp prepf u
    else begin integer n; scalar kern,m,w,x,xexp,y,y1,z;
        z := nil ./ 1;
    a0: kern := mvar u;
        if m := assoc(kern,asymplis!*) then m := cdr m;
    a:  if null u or (n := degr(u,kern))=0 then go to b
         else if null m or n<m then y := lt u . y;
        u := red u;
        go to a;
    b:  if not atom kern and not atom car kern then kern := prepf kern;
        if null l then xexp := if kern eq 'k!* then 1 else kern
         else if (xexp := algint!-subsublis(l,kern)) = kern
                   and not assoc(kern,asymplis!*)
          then go to f;
    c:  w := 1 ./ 1;
        n := 0;
        if y and cdaar y<0 then go to h;
        if (x := getrtype xexp) then typerr(x,"substituted expression");
        x := simp xexp;
        % SIMP!* here causes problem with HE package;
        x := reorder numr x ./ reorder denr x;
        % needed in case substitution variable is in XEXP;
        if null l and kernp x and mvar numr x eq kern then go to f
         else if null numr x then go to e;   %Substitution of 0;
        for each j in y do
         <<m := cdar j;
           w := multsq(exptsq(x,m-n),w);
           n := m;
           z := addsq(multsq(w,algint!-subf1(cdr j,l)),z)>>;
    e:  y := nil;
        if null u then return z
         else if domainp u then return addsq(algint!-subf1(u,l),z);
        go to a0;
    f:  sub2chk kern;
        for each j in y do 
           z := addsq(multpq(car j,algint!-subf1(cdr j,l)),z);
        go to e;
    h:  %Substitution for negative powers;
        x := simprecip list xexp;
    j:  y1 := car y . y1;
        y := cdr y;
        if y and cdaar y<0 then go to j;
    k:  m := -cdaar y1;
        w := multsq(exptsq(x,m-n),w);
        n := m;
        z := addsq(multsq(w,algint!-subf1(cdar y1,l)),z);
        y1 := cdr y1;
        if y1 then go to k else if y then go to c else go to e
     end;

symbolic procedure algint!-subsublis(u,v);
   begin scalar x;
      return if x := assoc(v,u) then cdr x
              else if atom v then v
              else if car v eq '!*sq then
                      list('!*sq,algint!-subsq(cadr v,u),caddr v)
%    Previous two lines added by JHD 7 July 1982.
%    without them, CDRs in SQ expressions buried inside;
%    !*SQ forms are lost;
              else if flagp!*!*(car v,'subfn) 
               then algint!-subsubf(u,v)
              else for each j in v collect algint!-subsublis(u,j)
   end;

symbolic procedure algint!-subsubf(l,expn);
   %Sets up a formal SUB expression when necessary;
   begin scalar x,y;
      for each j in cddr expn do
         if (x := assoc(j,l)) then <<y := x . y; l := delete(x,l)>>;
      expn := sublis(l,car expn)
                 . for each j in cdr expn 
                       collect algint!-subsublis(l,j);
        %to ensure only opr and individual args are transformed;
      if null y then return expn;
      expn := aconc!*(for each j in reversip!* y
                     collect list('equal,car j,aeval cdr j),expn);
      return mk!*sq if l then algint!-simpsub expn
                     else !*p2q mksp('sub . expn,1)
   end;

symbolic procedure algint!-simpsub u;
   begin scalar !*nosubs,w,x,z;
    a:  if null cdr u
          then <<if getrtype car u or eqcar(car u,'equal) 
                   then typerr(car u,"scalar");
                 u := simp!* car u;
                 z := reversip!* z;   % to put replacements in same
                                      % order as input.
                 return quotsq(algint!-subf(numr u,z),
                               algint!-subf(denr u,z))>>;
        !*nosubs := t;  % We don't want left side of eqns to change.
        w := reval car u;
        !*nosubs := nil;
        if getrtype w eq 'list
          then <<u := append(cdr w,cdr u); go to a>>
         else if not eqexpr w then errpri2(car u,t); 
        x := cadr w;
        if null getrtype x then x := !*a2k x;
        z := (x . caddr w) . z;
        u := cdr u;
        go to a;
   end;

endmodule;


module fracdi;

% Author: James H. Davenport.

fluid '(basic!-listofallsqrts basic!-listofnewsqrts expsub intvar
        sqrt!-intvar);

global '(coates!-fdi);

exports fdi!-print,fdi!-revertsq,fdi!-upgrade,
   fractional!-degree!-at!-infinity;

% internal!-fluid '(expsub);

symbolic procedure fdi!-print();
<< princ "We substitute";
   princ intvar;
   princ "**";
   princ coates!-fdi;
   princ " for ";
   princ intvar;
   printc " in order to avoid fractional degrees at infinity" >>;


symbolic procedure fdi!-revertsq u;
if coates!-fdi iequal 1
  then u
  else (fdi!-revert numr u) ./ (fdi!-revert denr u);


symbolic procedure fdi!-revert u;
if not involvesf(u,intvar)
  then u
  else addf(fdi!-revert red u,
            !*multf(fdi!-revertpow lpow u,
                    fdi!-revert lc u));


symbolic procedure fdi!-revertpow pow;
if not dependsp(car pow,intvar)
  then (pow .* 1) .+ nil
  else if car pow eq intvar
    then begin
      scalar v;
      v:=divide(cdr pow,coates!-fdi);
      if zerop cdr pow
        then return (mksp(intvar,car pow) .* 1) .+ nil
        else interr "Unable to revert fdi";
      end
    else if eq(car pow,'sqrt)
      then simpsqrt2 fdi!-revert !*q2f simp argof car pow
      else interr "Unrecognised term to revert";


symbolic procedure fdi!-upgrade place;
begin
  scalar ans,u,expsub,n;
  n:=coates!-fdi;
  for each u in place do
    if eqcar(u:=rsubs u,'expt)
      then n:=n / caddr u;
      % if already upgraded, we must take account of this.
  if n = 1
    then return place;
  expsub:=list(intvar,'expt,intvar,n);
  ans:=nconc(basicplace place,list expsub);
  expsub:=list expsub; % this prevents later nconc from causing trouble.
  u:=extenplace place;
  while u do begin
    scalar v,w,rfu;
    v:=fdi!-upgr2 lfirstsubs u;
    if v iequal 1
      then return (u:=cdr u);
    if eqcar(rfu:=rfirstsubs u,'minus)
      then w:=argof rfu
      else if eqcar(rfu,'sqrt)
        then w:=rfu
        else interr "Unknown place format";
    w:=fdi!-upgr2 w;
    if w iequal 1
      then interr "Place collapses under rewriting";
    if eqcar(rfu,'minus)
      then ans:=nconc(ans,list list(v,'minus,w))
      else ans:=nconc(ans,list(v.w));
    u:=cdr u;
    return
    end;
  sqrtsave(basic!-listofallsqrts,
           basic!-listofnewsqrts,
           basicplace ans);
  return ans
  end;


symbolic procedure fdi!-upgr2 u;
begin
  scalar v,mv;
  v:=substitutesq(simp u,expsub);
  if denr v neq 1
    then goto error;
  v:=numr v;
loop:
  if atom v
    then return v;
  if red v
    then go to error;
  mv:=mvar v;
  if (not dependsp(mv,intvar)) or (mv eq intvar)
    then <<
      v:=lc v;
      goto loop >>;
  if eqcar(mv,'sqrt)
    then if sqrtsinsf(lc v,nil,intvar)
      then go to error
      else return mv
    else go to error;
error:
  printc "*** Format error ***";
  princ "unable to go x:=x**";
  printc coates!-fdi;
  superprint u;
  rederr "Failure to make integral at infinity"
  end;


symbolic procedure fractional!-degree!-at!-infinity sqrts;
if sqrts
  then lcmn(fdi2 car sqrts,fractional!-degree!-at!-infinity cdr sqrts)
  else 1;


symbolic procedure fdi2 u;
   % Returns the denominator of the degree of x at infinity
   % in the sqrt expression u.
begin
  scalar n;
  u:=substitutesq(simp u,list list(intvar,'quotient,1,intvar));
  n:=0;
  while involvesq(u,sqrt!-intvar) do <<
    n:=iadd1 n;
    u:=substitutesq(u,list list(intvar,'expt,intvar,2)) >>;
  return (2**n)
  end;


symbolic procedure lcmn(i,j);
  i*j/gcdn(i,j);

% unfluid '(expsub);

endmodule;


module genus;

% Author: James H. Davenport.

fluid '(!*galois
        gaussiani
        intvar
        listofallsqrts
        listofnewsqrts
        nestedsqrts
        previousbasis
        sqrt!-intvar
        sqrt!-places!-alist
        sqrtflag
        sqrts!-in!-integrand
        taylorasslist
        taylorvariable);

global '(!*tra !*trmin);

symbolic procedure simpgenus u;
begin
  scalar intvar,sqrt!-intvar,taylorvariable,taylorasslist;
  scalar listofnewsqrts,listofallsqrts,sqrt!-places!-alist;
  scalar list!-of!-all!-sqrts,list!-of!-new!-sqrts;
  scalar sqrtflag,sqrts!-in!-integrand,tt,u,simpfn;
  tt:=readclock();
  sqrtflag:=t;
  taylorvariable:=intvar:=car u;
  simpfn:=get('sqrt,'simpfn);
  put('sqrt,'simpfn,'proper!-simpsqrt);
  sqrt!-intvar:=mvar !*q2f simpsqrti intvar;
  listofnewsqrts:= list mvar gaussiani; % Initialise the SQRT world.
  listofallsqrts:= list (argof mvar gaussiani . gaussiani);
  u:=for each v in cdr u
            collect simp!* v;
  sqrts!-in!-integrand:=sqrtsinsql(u,intvar);
  u:=!*n2sq length differentials!-1 sqrts!-in!-integrand;
  put('sqrt,'simpfn,simpfn);
  printc list('time,'taken,readclock()-tt,'milliseconds);
  return u
  end;
put('genus,'simpfn,'simpgenus);

symbolic procedure differentials!-1 sqrtl;
begin
  scalar asqrtl,faclist,places,v,nestedsqrts,basis,
         u,n,hard!-ones,sqrts!-in!-problem;
    % HARD!-ONES  A list of all the factors of our equations which do
    % not factor, and therefore such that we can divide the whole of
    % our INTBASIS by their product in order to get the true INTBASIS,
    % since these ones can cause no complications.
  asqrtl:=for each u in sqrtl
            collect !*q2f simp argof u;
  if !*tra or !*trmin then <<
    printc
      "Find the differentials of the first kind on curve defined by:";
    mapc(asqrtl,function printsf) >>;
  for each s in asqrtl do <<
    faclist:=for each u in jfactor(s,intvar)
               collect numr u;
    if !*tra then <<
      princ intvar;
      printc " is not a local variable at the roots of:";
      mapc(faclist,function printsf) >>;
    for each uu in faclist do <<
      v:=stt(uu,intvar);
      if 1 neq car v
        then hard!-ones:=uu.hard!-ones
        else <<
          u:=addf(uu,(mksp(intvar,1) .* (negf cdr v)) .+ nil) ./ cdr v;
          % U is now the value at which this SQRT has a zero.
          u:=list(list(intvar,'difference,intvar,prepsq u),
                  list(intvar,'expt,intvar,2));
          for each w in sqrtsign(for each w in union(delete(s,asqrtl),
                                                     delete(uu,faclist))
         conc sqrtsinsq(simpsqrtsq
      multsq(substitutesq(w ./ 1,u),
      1 ./ !*p2f mksp(intvar,2)),
                                      intvar),
                                 intvar)
            do places:=append(u,w).places >> >> >>;
  sqrts!-in!-problem:=nconc(for each u in hard!-ones
                              collect list(intvar.intvar,
                                    (lambda u;u.u) list('sqrt,prepf u)),
                            places);
  basis:=makeinitialbasis sqrts!-in!-problem;
                  % Bodge in any extra SQRTS that we will require later.
%  u:=1 ./ mapply(function multf,
%                for each u in sqrtl collect !*kk2f u);
%  basis:=for each v in basis collect multsq(u,v);
  basis:=integralbasis(mkvec basis,places,mkilist(places,-1),intvar);
  if not !*galois
    then basis:=combine!-sqrts(basis,
                               getsqrtsfromplaces sqrts!-in!-problem);
  if hard!-ones
    then <<
      v:=upbv basis;
      u:=1;
      for each v in hard!-ones do
        u:=multf(u,!*kk2f list('sqrt,prepf v));
      hard!-ones:=1 ./ u;
      for i:=0:v do
        putv(basis,i,multsq(getv(basis,i),hard!-ones)) >>;
  if not !*galois
    then basis:=modify!-sqrts(basis,sqrtl);
  v:=fractional!-degree!-at!-infinity sqrtl;
  if v iequal 1
    then n:=2
    else n:=2*v-1;
    % N  is the degree of the zero we need at INFINITY.
  basis:=normalbasis(basis,intvar,n);
  previousbasis:=nil;
    % it might have been set before, and we have changed its meaning.
  if !*tra or !*trmin then <<
    printc "Differentials are:";
    mapc(basis,function printsq) >>;
  return basis;
  end;

endmodule;


module intbasis;

% Author: James H. Davenport.

fluid '(excoatespoles intvar previousbasis taylorasslist
        taylorvariable);

global '(!*tra !*trmin);

exports completeplaces,completeplaces2,integralbasis;


symbolic procedure deleteplace(a,b);
if null b
  then nil
  else if equalplace(a,car b)
    then cdr b
    else (car b).deleteplace(a,cdr b);


symbolic procedure completeplaces(places,mults);
begin
  scalar current,cp,cm,op,om,ansp,ansm;
  if null places then return nil;       %%% ACH
loop:
  current:=basicplace car places;
  while places do <<
    if current = (basicplace car places)
      then <<
        cp:=(car places).cp;
        cm:=(car mults ).cm >>
      else <<
        op:=(car places).op;
        om:=(car mults ).om >>;
    places:=cdr places;
    mults:=cdr mults >>;
  cp:=completeplaces2(cp,cm,sqrtsinplaces cp);
  ansp:=append(car cp,ansp);
  ansm:=append(cdr cp,ansm);
  places:=op;
  mults:=om;
  cp:=op:=cm:=om:=nil;
  if places
    then go to loop
    else return ansp.ansm
  end;


symbolic procedure completeplaces2(places,mults,sqrts);
  % Adds extra places with multiplicities of 0 as necessary.
begin scalar b,p;
  sqrts:=sqrtsign(sqrts,intvar);
  b:=basicplace car places;
  p:=places;
  while p do <<
    if not(b = (basicplace car p))
      then interr "Multiple places not supported";
    sqrts:=deleteplace(extenplace car p,sqrts);
    p:=cdr p >>;
  mults:=nconc(nlist(0,length sqrts),mults);
  places:=nconc(mappend(sqrts,b),places);
  return places.mults
  end;


symbolic procedure intbasisreduction(zbasis,places,mults);
begin
  scalar i,m,n,v,w,substn,basis;
  substn:=list(intvar.intvar);
    % The X=X substitution.
  n:=upbv zbasis;
  basis:=copyvec(zbasis,n);
  taylorvariable:=intvar;
  v:=sqrtsinplaces places;
  for i:=0:n do
    w:=union(w,sqrtsinsq(getv(basis,i),intvar));
  m:=intersect(v,w);
  v:=purge(m,v);
  w:=purge(m,w);
  for each u in v do <<
    if !*tra or !*trmin then <<
      printc u;
      printc "does not occur in the functions";
      mapvec(basis,function printsq) >>;
    m:=!*q2f simp argof u;
    i:=w;
    while i and not quotf(m,!*q2f simp argof car i)
      do i:=cdr i;
    if null i
      then interr
         "Unable to find equivalent representation of branches";
    i:=car i;
    w:=delete(i,w);
    places:=subst(i,u,places);
    if !*tra or !*trmin then <<
      printc "replaced by";
      printc i >> >>;
  if (length places) neq (iadd1 n) then <<
   if !*tra
      then printc "Too many functions";
    basis := shorten!-basis basis;
    n:=upbv basis >>;
  m:=mkvect n;
  for i:=0:n do
    putv(m,i,cl6roweval(basis.i,places,mults,substn));
reductionloop:
  if !*tra then <<
    printc "Matrix before a reduction step:";
    mapvec(m,function printc) >>;
  v:=firstlinearrelation(m,iadd1 n);
  if null v
    then return replicatebasis(basis,(iadd1 upbv zbasis)/(n+1));
  i:=n;
  while null numr getv(v,i) do
    i:=isub1 i;
  w:=nil ./ 1;
  for j:=0:i do
    w:=!*addsq(w,!*multsq(getv(basis,j),getv(v,j)));
  w:=removecmsq multsq(w,1 ./ !*p2f mksp(intvar,1));
  if null numr w
    then <<
      mapvec(basis,function printsq);
      printc iadd1 i;
      interr "Basis collapses" >>;
  if !*tra then <<
    princ "Element ";
    princ iadd1 i;
    printc " of the basis replaced by ";
    if !*tra then
      printsq w >>;
  putv(basis,i,w);
  putv(m,i,cl6roweval(basis.i,places,mults,substn));
  goto reductionloop
  end;


symbolic procedure integralbasis(basis,places,mults,x);
begin
  scalar z,save,points,p,m,princilap!-part,mm;
  if null places
    then return basis;
  mults:=mapcar(mults,function (lambda u;min(u,0)));
  % this makes sure that we impose constraints only on
  % poles, not on zeroes.
  points:=removeduplicates mapcar(places,function basicplace);
  if points = list(x.x)
    then basis:=intbasisreduction(basis,places,mults)
    else if cdr points
      then go complex
      else <<
        substitutevec(basis,car points);
        if !*tra then <<
          printc "Integral basis reduction at";
          printc car points >>;
        basis:=intbasisreduction(basis,
                                 mapcar(places,function extenplace),
                                 mults);
        substitutevec(basis,antisubs(car points,x)) >>;
join:
  save:=taylorasslist;
  % we will not need te taylorevaluates at gensym.
  z:=gensym();
  places:=mapcons(places,x.list('difference,x,z));
  z:=list(x . z);
%  basis:=intbasisreduction(basis,
%                          places,
%                          nlist(0,length places),
%                          x,z);
  taylorasslist:=save;
  % ***time-hack-2***;
  if not excoatespoles
    then previousbasis:=copyvec(basis,upbv basis);
    % Save only if in COATES/FINDFUNCTION, not if in EXCOATES.
  return basis;
complex:
  while points do <<
    p:=places;
    m:=mults;
    princilap!-part:=mm:=nil;
    while p do <<
    if (car points) = (basicplace car p)
      then <<
        princilap!-part:=(extenplace car p).princilap!-part;
        mm:=(car m).mm >>;
      p:=cdr p;
      m:=cdr m >>;
    substitutevec(basis,car points);
    if !*tra then <<
      printc "Integral basis reduction at";
      printc car points >>;
    basis:=intbasisreduction(basis,princilap!-part,mm);
    substitutevec(basis,antisubs(car points,x));
    points:=cdr points >>;
  go to join
  end;


symbolic procedure cl6roweval(basisloc,places,mults,x!-alpha);
% Evaluates a row of the matrix in coates lemma 6.
begin
  scalar i,v,w,save,basiselement,taysave,mmults,flg;
  i:=isub1 length places;
  v:=mkvect i;
  taysave:=mkvect i;
  i:=0;
  basiselement:=getv(car basisloc,cdr basisloc);
  mmults:=mults;
  while places do <<
    w:=substitutesq(basiselement,car places);
    w:=taylorform substitutesq(w,x!-alpha);
      % The separation of these 2 is essential since the x->x-a
      % must occur after the places are chosen.
    save:=taylorasslist;
    if not flg
      then putv(taysave,i,w);
    w:=taylorevaluate(w,car mmults);
    tayshorten save;
    putv(v,i,w);
    i:=iadd1 i;
    flg:=flg or numr w;
    mmults:=cdr mmults;
    places:=cdr places >>;
  if flg
    then return v;
    % There was a non-zero element in this row.
  save:=0;
loop:
  save:=iadd1 save;
  mmults:=mults;
  i:=0;
  while mmults do <<
    w:=taylorevaluate(getv(taysave,i),save + car mmults);
    flg:=flg or numr w;
    mmults:=cdr mmults;
    putv(v,i,w);
    i:=iadd1 i >>;
  if not flg
    then go to loop;
    % Another zero row.
  putv(car basisloc,cdr basisloc,multsq(basiselement,
                                        1 ./ !*p2f mksp(intvar,save)));
  return v
  end;


symbolic procedure replicatebasis(basis,n);
if n = 1
  then basis
  else if n = 2
    then begin
      scalar b,sqintvar,len;
      len:=upbv basis;
      sqintvar:=!*kk2q intvar;
      b:=mkvect(2*len+1);
      for i:=0:len do <<
        putv(b,i,getv(basis,i));
        putv(b,i+len+1,multsq(sqintvar,getv(basis,i))) >>;
      return b
      end
    else interr "Unexpected replication request";


symbolic procedure shorten!-basis v;
begin
  scalar u,n,sfintvar;
  sfintvar:=!*kk2f intvar;
  n:=upbv v;
  for i:=0:n do begin
    scalar uu;
    uu:=getv(v,i);
    if not quotf(numr uu,sfintvar)
      then u:=uu.u
    end;
  return mkvec u
  end;


endmodule;


module jhddiff;

% Author: James H. Davenport.

fluid '(dw);

% Differentiation routines for algebraic expressions;
symbolic procedure !*diffsq(u,v);
   %U is a standard quotient, V a kernel.
   %Value is the standard quotient derivative of U wrt V.
   %Algorithm: df(x/y,z)= (x'-(x/y)*y')/y;
   !*multsq(!*addsq(!*difff(numr u,v),
                    negsq !*multsq(u,!*difff(denr u,v))),
          1 ./ denr u);

symbolic procedure !*difff(u,v);
   %U is a standard form, V a kernel.
   %Value is the standard quotient derivative of U wrt V;
   if domainp u then nil ./ 1
    else !*addsq(!*addsq(multpq(lpow u,!*difff(lc u,v)),
                        !*multsq(lc u ./ 1,!*diffp(lpow u,v))),
               !*difff(red u,v));

symbolic procedure !*diffp(u,v);
%  Special treatment of SQRT's (JHD is not sure why,
%  but it seems to be necessary);
if atom (car u) then diffp(u,v)
  else if not (caar u) eq 'sqrt then diffp(u,v)
    else begin
           scalar w,dw;
           w:=simp argof car u;
           dw:= !*diffsq(w,v);
           if null numr dw then return dw;
           return !*multsq(!*multsq(dw,invsq w),
                           !*multf(cdr u,mksp(car u,1) .* 1 .+ nil)./ 2)
           end;

endmodule;


module jhdriver;

% Author: James H. Davenport.

fluid '(!*backtrace
        basic!-listofallsqrts
        basic!-listofnewsqrts
        expression
        gaussiani
        intvar
        listofallsqrts
        listofnewsqrts
        previousbasis
        sqrt!-intvar
        sqrtflag
        sqrts!-in!-integrand
        sqrts!-mod!-prime
        taylorasslist
        varlist
        zlist);

global '(!*algint !*coates !*noacn !*tra !*trmin btrlevel tryharder);

switch algint,coates,noacn,tra,trmin;

exports algebraiccase,doalggeom,coates!-multiple;

!*algint := t;   % Assume algebraic integration wanted if this module
                 % is loaded.

symbolic procedure operateon(reslist,x);
begin
  scalar u,v,answer,save;
  scalar sqrts!-mod!-prime;
  u:=zmodule(reslist);
  v:=answer:=nil ./ 1;
  while u and not atom v do <<
    v:=findfunction cdar u;
    if not atom v then <<
      if !*tra or !*trmin then <<
        printc "Extension logarithm is ";
        printsq v >>;
      save:=tryharder;
      tryharder:=x;
      v:= !*multsq(simp!* caar u,
                simplogsq v);
      tryharder:=save;
      answer:=addsq(answer,v);
      u:=cdr u >> >>;
  if atom v
    then return v
    else return answer
  end;


symbolic procedure findfunction divisor;
begin
  scalar v,places,mults,ans,dof1k;
  scalar previousbasis;
  % ***time-hack-2 :::
    % A hack for decreasing the amount of work done in COATES.
  divisor:=for each u in divisor collect
             correct!-mults u;
  if !*coates
    then go to nohack;
  v:=precoates(divisor,intvar,nil);
  if not atom v
    then return v;
nohack:
  for each u in divisor do <<
    places:=(car u).places;
    mults :=(cdr u).mults >>;
  v:=coates(places,mults,intvar);
  if not atom v
    then return v;
  dof1k:=differentials!-1 getsqrtsfromplaces places;
  if null dof1k
    then interr "Must be able to integrate over curves of genus 0";
  if not mazurp(places,dof1k)
    then go to general;
  ans:='provably!-impossible;
  for i:=2:12 do
    if (i neq 11) and
       not atom (ans:=coates!-multiple(places,mults,i))
    then i:=12;   % leave the loop - we have an answer.
  return ans;
general:
  v:=findmaninparm places;
  if null v
     then return algebraic!-divisor(divisor,dof1k);
  if not maninp(divisor,v,dof1k)
    then return 'provably!-impossible;
  v:=1;
loop:
  v:=iadd1 v;
  if not atom (ans:=coates!-multiple(places,mults,v))
    then return ans;
  go to loop
  end;


symbolic procedure correct!-mults u;
begin
  scalar multip;
  multip:=cdr u;
  for each v in car u do
    if (lsubs v eq intvar) and
        eqcar(rsubs v,'expt)
      then multip:=multip * (caddr rsubs v);
    return (car u).multip
  end;


symbolic procedure algebraiccase
    (expression,zlist,varlist);
begin
  scalar rischpart,deriv,w,firstterm;
  scalar sqrtflag;
  sqrtflag:=t;
  sqrtsave(listofallsqrts,listofnewsqrts,list(intvar . intvar));
  rischpart:=errorset('(doalggeom expression),
                      if !*tra or !*trmin then t else btrlevel,
                      !*backtrace);
  newplace list (intvar.intvar);
  if atom rischpart
    then <<
      if !*tra then printc "Inner integration failed";
      deriv:=nil ./ 1;
      % assume no answer.
      rischpart:=deriv >>
    else
      if atom car rischpart
        then <<
          if !*tra or !*trmin then
            printc "The 'logarithmic part' is not elementary";
          return simpint1 list ('int,prepsq expression,intvar) >>
      else <<
        rischpart:=car rischpart;
        deriv:=!*diffsq(rischpart,intvar);
        % deriv := squashsqrt deriv;
        % Should no longer be necessary.
        if !*tra or !*trmin then <<
          printc "Inner working yields";
          printsq rischpart;
          printc "with derivative";
          printsq deriv >> >>;
  deriv:=!*addsq(expression,negsq deriv);
  if null numr deriv
    then return rischpart; % no algebraic part.
  if null involvesq(deriv,intvar)
    then return !*addsq(rischpart,
                !*multsq(deriv,((mksp(intvar,1) .* 1) .+ nil) ./ 1));
                % if the difference is merely a constant.
  varlist:=getvariables deriv;
  zlist:=findzvars(varlist,list intvar,intvar,nil);
  varlist:=purge(zlist,varlist);
  firstterm:=simp!* car zlist; % this may crop up.
  w:=sqrt2top !*multsq(deriv,invsq !*diffsq(firstterm,intvar));
  if null involvesq(w,intvar)
    then return !*addsq(rischpart,!*multsq(w,firstterm));
  if !*noacn then interr "Testing only logarithmic code";
  deriv:=transcendentalcase(deriv,intvar,nil,zlist,varlist);
  return !*addsq(deriv,rischpart)
  end;


symbolic procedure doalggeom(differential);
begin
  scalar reslist,place,placelist,
         savetaylorasslist,sqrts!-in!-integrand,
         taylorasslist;
  placelist:=findpoles(differential,intvar);
  reslist:=nil;
  sqrts!-in!-integrand:=sqrtsinsq (differential,intvar);
  while placelist do <<
    place:=car placelist;
    placelist:=cdr placelist;
    savetaylorasslist:=taylorasslist;
    place:=find!-residue(differential,intvar,place);
    if place
      then reslist:=append(place,reslist)
      else taylorasslist:=savetaylorasslist >>;
  if reslist
    then go to serious;
  if !*tra or !*trmin
    then printc "No residues => no logs";
  return nil ./ 1;
serious:
  placelist:=operateon(reslist,intvar);
  if placelist eq 'failed
    then interr "Divisor operations failed";
  return placelist
  end;


symbolic procedure algebraic!-divisor(divisor,dof1k);
if length dof1k = 1
  then lutz!-nagell(divisor)
  else bound!-torsion(divisor,dof1k);


symbolic procedure coates!-multiple(places,mults,v);
begin
  scalar ans;
  if not atom (ans:=coates(places,
                           for each u in mults collect v*u,
                           intvar))
    then <<
      if !*tra or !*trmin then <<
        princ "Divisor has order ";
        printc v >>;
      return !*kk2q list('nthroot,mk!*sq ans,v) >>
    else return ans
  end;


symbolic procedure mazurp(places,dof1k);
   % Checks to ensure we have an elliptic curve over the rationals.
begin
%  scalar sqrt2,sqrt4,v;
%  sqrt2:=0;
%    % Number of SQRTs of things of degree 1 or 2;
%  sqrt4:=0;
%    % " " " 3 or 4;
%  for each u in getsqrtsfromplaces places do <<
%    v:=!*q2f simp u;
%    if sqrtsinsq(v,intvar)
%      then return nil;
%      % Cannot use nested SQRTs;
%    v:=car stt(v,intvar);
%    if v < 3
%      then if sqrt4>0
%        then return nil
%        else if sqrt2>1
%          then return nil
%          else sqrt2:=iadd1 sqrt2
%      else if v < 5
%        then if sqrt2>0 or sqrt4>0
%          then return nil
%          else sqrt4:=1
%        else return nil >>;
  scalar answer;
  if length dof1k neq 1
    then return nil;
    % Genus = # linearly independent differentials of 1st kind;
    % We know know that it is of genus = 1.
  answer:=t;
  while answer and places do
    if sqrtsintree(basicplace car places,nil,nil)
      then answer:= nil
      else places:=cdr places;
  if null answer then return nil;
  if !*tra then
    <<prin2 "*** We can apply Mazur's bound on the torsion of";
      prin2t "elliptic curves over the rationals">>;
  return t
  end;

endmodule;


module linrel;

% Author: James H. Davenport.

symbolic procedure firstlinearrelation(m,n);
% Returns vector giving first linear relation between
% the rows of n*n matrix m.
begin
  scalar mm,u,uu,v,w,x,xx,i,j,isub1n,ans;
  isub1n:=isub1 n;
  mm:=mkvect(isub1n);
  for i:=0 step 1 until isub1n do
    putv(mm,i,copyvec(getv(m,i),isub1n));
  % mm is a copy of m which we can afford to destroy.
  ans:=mkidenm isub1n;
  i:=0;
outerloop:
  u:=getv(mm,i);
  uu:=getv(ans,i);
  j:=0;
pivotsearch:
  if j iequal n
    then goto zerorow;
  v:=getv(u,j);
  if null numr v then << j:=iadd1 j; goto pivotsearch >>;
  % we now use the j-th element of row i to flatten the j-th
  % element of all later rows.
  if i iequal isub1n then return nil;
    %no further rows to flatten, so no relationships.
  v:=!*invsq negsq v;
  for k:=iadd1 i step 1 until isub1n do <<
    xx:=getv(ans,k);
    x:=getv(mm,k);
    w:=!*multsq(v,getv(x,j));
    for l:=0:isub1n do <<
      putv(x,l,addsq(getv(x,l),!*multsq(w,getv(u,l))));
      putv(xx,l,addsq(getv(xx,l),!*multsq(w,getv(uu,l)))) >> >>;
  i:=iadd1 i;
  if i < n then goto outerloop;
  % no zero rows found at all.
  return nil;
zerorow:
  % the i-t row is all zero, i.e. rows 1...i are dependent.
  return getv(ans,i)
  end;

endmodule;


module maninp;

% Author: James H. Davenport.

fluid '(intvar);

symbolic procedure findmaninparm places;
begin
  scalar sqrts,vars,u;
  sqrts:=sqrtsinplaces places;
loop:
  if null sqrts then return nil;
  vars:=getvariables simp argof car sqrts;
innerloop:
  if null vars
    then <<
      sqrts:=cdr sqrts;
      go to loop >>;
  u:=car vars;
  vars:=cdr vars;
  if u eq intvar
    then go to innerloop;
  if atom u
    then return u;
  if car u eq 'sqrt
    then << u:=simp argof u;
            vars:=varsinsf(numr u,varsinsf(denr u,vars));
            go to innerloop >>;
  interr "Unrecognised differentiation candidate"
  end;

endmodule;


module modify;

% Author: James H. Davenport.

fluid '(intvar);

global '(!*tra);

exports modify!-sqrts,combine!-sqrts;

symbolic procedure modify!-sqrts(basis,sqrtl);
begin
  scalar sqrtl!-in!-sf,n,u,v,f;
  n:=upbv basis;
  sqrtl!-in!-sf:=for each u in sqrtl collect
                    !*q2f simp argof u;
  for i:=0:n do begin
    u:=getv(basis,i);
    v:=sqrtsinsq(u,intvar);
    % We have two tasks to perform,
    % the replacing of SQRT(A)*SQRT(B) by SQRT(A*B)
    % where relevant and the replacing of SQRT(A)
    % by SQRT(A*B) or 1 (depending on whether it occurs in
    % the numerator or the denominator).
    v:=purge(sqrtl,v);
    if null v
      then go to nochange;
    u:=sqrt2top u;
    u:=multsq(modify2(numr u,v,sqrtl!-in!-sf) ./ 1,
              1 ./ modify2(denr u,v,sqrtl!-in!-sf));
    v:=sqrtsinsq(u,intvar);
    v:=purge(sqrtl,v);
    if v then <<
      if !*tra then <<
        printc "Discarding element";
        printsq u >>;
      putv(basis,i,1 ./ 1) >>
      else putv(basis,i,removecmsq u);
    f:=t;
  nochange:
    end;
  basis:=mkuniquevect basis;
  if f and !*tra then <<
    printc "Basis replaced by";
    mapvec(basis,function printsq) >>;
  return basis
  end;


symbolic procedure combine!-sqrts(basis,sqrtl);
begin
  scalar sqrtl!-in!-sf,n,u,v,f;
  n:=upbv basis;
  sqrtl!-in!-sf:=for each u in sqrtl collect
                    !*q2f simp argof u;
  for i:=0:n do begin
    u:=getv(basis,i);
    v:=sqrtsinsq(u,intvar);
    % We have one task to perform,
    % the replacing of SQRT(A)*SQRT(B) by SQRT(A*B)
    % where relevant.
    v:=purge(sqrtl,v);
    if null v
      then go to nochange;
    u:=multsq(modify2(numr u,v,sqrtl!-in!-sf) ./ 1,
              1 ./ modify2(denr u,v,sqrtl!-in!-sf));
    putv(basis,i,u);
    f:=t;
  nochange:
    end;
  if f and !*tra then <<
    printc "Basis replaced by";
    mapvec(basis,function printsq) >>;
  return basis
  end;


symbolic procedure modify2(sf,sqrtsin,realsqrts);
if atom sf
  then sf
  else if atom mvar sf
    then sf
    else if eqcar(mvar sf,'sqrt) and dependsp(mvar sf,intvar)
      then begin
        scalar u,v,w,lcsf,sqrtsin2,w2,lcsf2,temp;
        u:=!*q2f simp argof mvar sf;
        v:=realsqrts;
        while v and null (w:=modify!-quotf(car v,u))
          do v:=cdr v;
        if null v
          then <<
            if !*tra then <<
              printc "Unable to modify (postponed)";
              printsf !*kk2f mvar sf >>;
            return sf >>;
        v:=car v;
        % We must modify SQRT(U) into SQRT(V) if possible.
        lcsf:=lc sf;
        sqrtsin2:=delete(mvar sf,sqrtsin);
        while sqrtsin2 and (w neq 1) do <<
          temp:=!*q2f simp argof car sqrtsin2;
          if (w2:=modify!-quotf(w,temp)) and
             (lcsf2:=modify!-quotf(lcsf,!*kk2f car sqrtsin2))
            then <<
              w:=w2;
              lcsf:=lcsf2 >>;
          sqrtsin2:=cdr sqrtsin2 >>;
        if w = 1
          then return addf(multf(lcsf,formsqrt v),
                           modify2(red sf,sqrtsin,realsqrts));
                           % It is important to use FORMSQRT here since
                           % SIMPSQRT will recreate the factorisation
                           % we are trying to destroy.
          % Satisfactorily explained away.
        return addf(multf(!*p2f lpow sf,
                          modify2(lc sf,sqrtsin,realsqrts)),
                    modify2(red sf,sqrtsin,realsqrts))
        end
      else addf(multf(!*p2f lpow sf,
                      modify2(lc sf,sqrtsin,realsqrts)),
                modify2(red sf,sqrtsin,realsqrts));



%symbolic procedure modifydown(sf,sqrtl);
%if atom sf
%  then sf
%  else if atom mvar sf
%    then sf
%    else if eqcar(mvar sf,'sqrt) and
%            dependsp(mvar sf,intvar) and
%           not member(!*q2f simp argof mvar sf,sqrtl)
%      then addf(modifydown(lc sf,sqrtl),
%                modifydown(red sf,sqrtl))
%      else addf(multf(!*p2f lpow sf,
%                      modifydown(lc sf,sqrtl)),
%                modifydown(red sf,sqrtl));


% symbolic procedure modifyup(sf,sqrtl);
% if atom sf
%   then sf
%   else if atom mvar sf
%     then sf
%     else if eqcar(mvar sf,'sqrt) and
%             dependsp(mvar sf,intvar)
%       then begin
%         scalar u,v;
%         u:=!*q2f simp argof mvar sf;
%         if u member sqrtl
%         then return addf(multf(!*p2f lpow sf,
%                                 modifyup(lc sf,sqrtl)),
%                           modifyup(red sf,sqrtl));
%        v:=sqrtl;
%        while v and not modify!-quotf(car v,u)
%          do v:=cdr v;
%        if null v
%          then interr "No sqrt to upgrade to";
%       return addf(multf(!*kk2f simpsqrt2 car v,
%                          modifyup(lc sf,sqrtl)),
%                    modifyup(red sf,sqrtl))
%        end
%      else addf(multf(!*p2f lpow sf,
%                      modifyup(lc sf,sqrtl)),
%                modifyup(red sf,sqrtl));


symbolic procedure modify!-quotf(u,v);
% Replacement for quotf, in that it gets sqrts right.
if atom v or atom mvar v
  then quotf(u,v)
  else if u=v then 1
  else begin
    scalar sq;
    sq:=sqrt2top(u ./ v);
    if involvesf(denr sq,intvar)
      then return nil;
    if not onep denr sq
      then if not numberp denr sq
        then interr "Gauss' lemma violated in modify"
        else if !*tra
          then <<
            printc "*** Denominator ignored in modify";
            printc denr sq >>;
    return numr sq
    end;

endmodule;


module modlineq;

% Author: James H. Davenport.

fluid '(current!-modulus sqrts!-mod!-prime);

global '(!*tra !*trmin list!-of!-medium!-primes sqrts!-mod!-8);

exports check!-lineq;

list!-of!-medium!-primes:='(101 103 107 109);

sqrts!-mod!-8:=mkvect 7;

putv(sqrts!-mod!-8,0,t);

putv(sqrts!-mod!-8,1,t);

putv(sqrts!-mod!-8,4,t);

symbolic procedure modp!-nth!-root(m,n,p);
begin
  scalar j,p2;
  p2:=p/2;
  for i:=-p2 step 1 until p2 do
    if modular!-expt(i,n) iequal m
      then << j:=i; i:=p2 >>;
  return j
  end;


symbolic procedure modp!-sqrt(n,p);
begin
  scalar p2,s,tt;
  p2:=p/2;
  if n < 0
    then n:=n+p;
  for i:=0:p2 do begin
    tt:=n+p*i;
    if null getv(sqrts!-mod!-8,tt irem 8)
      then return;
      % mod 8 test for perfect squares.
    if (iadd1 tt irem 5) > 2
      then return;
      % squares are -1,0,1 mod 5.
    s:=int!-sqrt tt;
    if fixp s then <<
      p2:=0;
      return >>
    end;
  if (not fixp s) or null s
    then return nil
    else return s
  end;

symbolic procedure subsetp(a,b);
%True if all members of a are also members of b.
    if null a then t
    else if member(car a,b) then subsetp(cdr a,b)
    else nil;

symbolic procedure check!-lineq(m,rightside);
begin
  scalar vlist,n1,n2,u,primelist,mm,v,modp!-subs,atoms;
  n1:=upbv m;
  for i:=0:n1 do <<
    u:=getv(m,i);
    if u
      then for j:=0:(n2:=upbv u) do
        vlist:=varsinsq(getv(u,j),vlist) >>;
  u:=vlist;
  while u do <<
    v:=car u;
    u:=cdr u;
    if atom v
      then atoms:=v.atoms
      else if (car v eq 'sqrt) or (car v eq 'expt)
        then for each w in varsinsf(!*q2f simp argof v,nil) do
             if not (w member vlist)
               then <<
                 u:=w.u;
                 vlist:=w.vlist >>
        else nil
      else interr "Unexpected item" >>;
  if sqrts!-mod!-prime and
     subsetp(vlist,for each u in cdr sqrts!-mod!-prime
                     collect car u)
    then go to end!-of!-loop;
  vlist:=purge(atoms,vlist);
  u:=nil;
  for each v in vlist do
    if car v neq 'sqrt
      then u:=v.u;
  vlist:=nconc(u,sortsqrts(purge(u,vlist),nil));
    % NIL is the variable to measure nesting on:
    % therefore all nesting is being caught.
  primelist:=list!-of!-medium!-primes;
  set!-modulus car primelist;
  atoms:=for each u in atoms collect
           u . modular!-number random car primelist;
  goto try!-prime;
next!-prime:
  primelist:=cdr primelist;
  if null primelist and !*tra
    then printc "Ran out of primes in check!-lineq";
  if null primelist
    then return t;
  set!-modulus car primelist;
try!-prime:
  modp!-subs:=atoms;
  v:=vlist;
loop:
  if null v
    then go to end!-of!-loop;
  u:=modp!-subst(simp argof car v,modp!-subs);
  if caar v eq 'sqrt
    then u:=modp!-sqrt(u,car primelist)
    else if caar v eq 'expt
      then u:=modp!-nth!-root(modular!-expt(u,cadr caddr car v),
        caddr caddr car v,car primelist)
      else interr "Unexpected item";
  if null u
    then go to next!-prime;
  modp!-subs:=(car v . u) . modp!-subs;
  v:=cdr v;
  go to loop;
end!-of!-loop:
  if null primelist
    then <<
      setmod(car sqrts!-mod!-prime);
      modp!-subs:=cdr sqrts!-mod!-prime >>
    else sqrts!-mod!-prime:=(car primelist).modp!-subs;
  mm:=mkvect n1;
  for i:=0:n1 do begin
    u:=getv(m,i);
    if null u
      then return;
    putv(mm,i,v:=mkvect n2);
    for j:=0:n2 do
      putv(v,j,modp!-subst(getv(u,j),modp!-subs))
    end;
  v:=mkvect n1;
  for i:=0:n1 do
    putv(v,i,modp!-subst(getv(rightside,i),modp!-subs));
  u:=mod!-jhdsolve(mm,v);
  if (u eq 'failed) and (!*tra or !*trmin)
    then <<
      princ "Proved insoluble mod ";
      printc car sqrts!-mod!-prime >>;
  return u
  end;


symbolic procedure modp!-subst(sq,slist);
modular!-quotient(modp!-subf(numr sq,slist),
                  modp!-subf(denr sq,slist));


symbolic procedure modp!-subf(sf,slist);
if atom sf
  then if null sf
    then 0
    else modular!-number sf
  else begin
    scalar u;
    u:=assoc(mvar sf,slist);
    if null u
      then interr "Unexpected variable";
    return modular!-plus(modular!-times(modular!-expt(cdr u,ldeg sf),
                                        modp!-subf(lc sf,slist)),
                         modp!-subf(red sf,slist))
    end;


symbolic procedure mod!-jhdsolve(m,rightside);
% Returns answer to m.answer=rightside.
% Matrix m not necessarily square.
begin
  scalar n1,n2,ans,u,row,swapflg,swaps;
  % The SWAPFLG is true if we have changed the order of the
  % columns and need later to invert this via SWAPS.
  n1:=upbv m;
  for i:=0:n1 do
    if (u:=getv(m,i))
      then (n2:=upbv u);
  swaps:=mkvect n2;
  for i:=0:n2 do
    putv(swaps,i,n2-i);
    % We have the SWAPS vector, which should be a vector of indices,
    % arranged like this because VECSORT sorts in decreasing order.
  for i:=0:isub1 n1 do begin
    scalar k,v,pivot;
  tryagain:
    row:=getv(m,i);
    if null row
      then go to interchange;
    % look for a pivot in row.
    k:=-1;
    for j:=0:n2 do
      if not zerop (pivot:=getv(row,j))
        then <<
          k:=j;
          j:=n2 >>;
    if k neq -1
      then goto newrow;
    if not zerop getv(rightside,i)
      then <<
        m:='failed;
        i:=sub1 n1; %Force end of loop.
        go to finished >>;
interchange:
    % now interchange i and last element.
    swap(m,i,n1);
    swap(rightside,i,n1);
    n1:=isub1 n1;
    if i iequal n1
      then goto finished
      else goto tryagain;
  newrow:
    if i neq k
      then <<
        swapflg:=t;
        swap(swaps,i,k);
          % record what we have done.
        for l:=0:n1 do
          swap(getv(m,l),i,k) >>;
        % place pivot on diagonal.
    pivot:=modular!-minus modular!-reciprocal pivot;
    for j:=iadd1 i:n1 do begin
      u:=getv(m,j);
      if null u
        then return;
      v:=modular!-times(getv(u,i),pivot);
      if not zerop v then <<
        putv(rightside,j,
            modular!-plus(getv(rightside,j),
                modular!-times(v,getv(rightside,i))));
        for l:=0:n2 do
          putv(u,l,
             modular!-plus(getv(u,l),
                 modular!-times(v,getv(row,l)))) >>
      end;
  finished:
    end;
  if m eq 'failed
    then go to failed;
    % Equations were inconsistent.
  while null (row:=getv(m,n1)) do
    n1:=isub1 n1;
  u:=nil;
  for i:=0:n2 do
    if not zerop getv(row,i)
      then u:='t;
  if null u
    then if not zerop getv(rightside,n1)
      then go to failed
      else n1:=isub1 n1;
      % Deals with a last equation which is all zero.
  if n1 > n2
    then go to failed;
    % Too many equations to satisfy.
  ans:=mkvect n2;
  for i:=0:n2 do
    putv(ans,i,0);
  % now to do the back-substitution.
  for i:=n1 step -1 until 0 do begin
    row:=getv(m,i);
    if null row
      then return;
    u:=getv(rightside,i);
    for j:=iadd1 i:n2 do
      u:=modular!-plus(u,
         modular!-times(getv(row,j),modular!-minus getv(ans,j)));
    putv(ans,i,modular!-times(u,modular!-reciprocal getv(row,i)))
    end;
  if swapflg
    then vecsort(swaps,list ans);
  return ans;
failed:
  if !*tra
    then printc "Unable to force correct zeroes";
  return 'failed
  end;

endmodule;


module nagell;

% Author: James H. Davenport.

fluid '(intvar);

global '(!*tra !*trmin);

exports lutz!-nagell;

symbolic procedure lutz!-nagell(divisor);
begin
  scalar ans,places,mults,save!*tra;
  for each u in divisor do <<
    places:=(car u).places;
    mults :=(cdr u).mults >>;
  ans:=lutz!-nagell!-2(places,mults);
  save!*tra:=!*tra;
  if !*trmin
    then !*tra:=nil;
  ans:=coates!-multiple(places,mults,ans);
  !*tra:=save!*tra;
  return ans
  end;


symbolic procedure lutz!-nagell!-2(places,mults);
begin
  scalar wst,x,y,equation,point,a;
  wst:=weierstrass!-form getsqrtsfromplaces places;
  x:=car wst;
  y:=cadr wst;
  equation:=caddr wst;
  equation:=!*q2f !*multsq(equation,equation);
  equation:=makemainvar(equation,intvar);
  if ldeg equation = 3
    then equation:=red equation
    else interr "Equation not of correct form";
  if mvar equation eq intvar
    then if ldeg equation = 1
      then <<
        a:=(lc equation) ./ 1;
        equation:=red equation >>
      else interr "Equation should not have a x**2 term"
    else a:=nil ./ 1;
  equation:= a . (equation ./ 1);
  places:=for each u in places collect
            wst!-convert(u,x,y);
  point:=elliptic!-sum(places,mults,equation);
  a:=lutz!-nagell!-bound(point,equation);
  if !*tra or !*trmin then <<
    princ "Point actually is of order ";
    printc a >>;
  return a
  end;


symbolic procedure wst!-convert(place,x,y);
begin
  x:=subzero(xsubstitutesq(x,place),intvar);
  y:=subzero(xsubstitutesq(y,place),intvar);
  return x.y
  end;


symbolic procedure elliptic!-sum(places,mults,equation);
begin
  scalar point;
  point:=elliptic!-multiply(car places,car mults,equation);
  places:=cdr places;
  mults:=cdr mults;
  while places do <<
    point:=elliptic!-add(point,
                         elliptic!-multiply(car places,car mults,
                                            equation),
                         equation);
    places:=cdr places;
    mults:=cdr mults >>;
  return point
  end;


symbolic procedure elliptic!-multiply(point,n,equation);
if n < 0
  then elliptic!-multiply( (car point) . (negsq cdr point),
                           -n,
                           equation)
  else if n = 0
    then interr "N=0 in elliptic!-multiply"
    else if n = 1
      then point
      else begin
        scalar q,r;
        q:=divide(n,2);
        r:=cdr q;
        q:=car q;
        q:=elliptic!-multiply(elliptic!-add(point,point,equation),q,
                                            equation);
        if r = 0
          then return q
          else return elliptic!-add(point,q,equation)
        end;


symbolic procedure elliptic!-add(p1,p2,equation);
begin
  scalar x1,x2,y1,y2,x3,y3,inf,a,b,lhs,rhs;
  a:=car equation;
  b:=cdr equation;
  inf:=!*kk2q 'infinity;
  x1:=car p1;
  y1:=cdr p1;
  x2:=car p2;
  y2:=cdr p2;
  if x1 = x2
    then if y1 = y2
      then <<
        % this is the doubling case.
        x3:=!*multsq(!*addsq(!*addsq(!*multsq(a,a),
                                     !*exptsq(x1,4)),
                             !*addsq(multsq(-8 ./ 1,!*multsq(x1,b)),
                                     !*multsq(!*multsq(x1,x1),
                                              multsq(-2 ./ 1,a)))),
                     !*invsq multsq(4 ./ 1,
                          !*addsq(b,!*multsq(x1,!*addsq(a,
                                                   !*exptsq(x1,2))))));
        y3:=!*addsq(y1,!*multsq(!*multsq(!*addsq(x3,negsq x1),
                                         !*addsq(a,multsq(3 ./ 1,
                                                     !*multsq(x1,x1)))),
                                 !*invsq multsq(2 ./ 1,
                                                y1))) >>
      else x3:=(y3:=inf)
    else if x1 = inf
      then <<
        x3:=x2;
        y3:=y2 >>
      else if x2 = inf
        then <<
          x3:=x1;
          y3:=y1 >>
        else <<
          x3:=!*multsq(!*addsq(!*multsq(a,!*addsq(x1,x2)),
                               !*addsq(multsq(2 ./ 1,b),
                                       !*addsq(!*multsq(!*multsq(x1,x2),
                                                        !*addsq(x1,x2)),
                                               multsq(-2 ./ 1,
                                                    !*multsq(y1,y2))))),
                       !*invsq !*exptsq(!*addsq(x1,negsq x2),2));
          y3:=!*multsq(!*addsq(!*multsq(!*addsq(y2,negsq y1),x3),
                               !*addsq(!*multsq(x2,y1),
                                       !*multsq(x1,negsq y2))),
                       !*invsq !*addsq(x1,negsq x2)) >>;
  if x3 = inf
    then return x3.y3;
  lhs:=!*multsq(y3,y3);
  rhs:=!*addsq(b,!*multsq(x3,!*addsq(a,!*multsq(x3,x3))));
  if numr !*addsq(lhs,negsq rhs) % We can't just compare them
                                  % since they're algebraic numbers.
                                  % JHD Jan 14th. 1987.
    then <<
      prin2t "Point defined by X and Y as follows:";
      printsq x3;
      printsq y3;
      prin2t "on the curve defined by A and B as follows:";
      printsq a;
      printsq b;
      prin2t "gives a consistency check between:";
      printsq lhs;
      printsq rhs;
      interr "Consistency check failed in elliptic!-add" >>;
  return x3.y3
  end;




symbolic procedure infinitep u;
kernp u and (mvar numr u eq 'infinite);


symbolic procedure lutz!-nagell!-bound(point,equation);
begin
  scalar x,y,a,b,lutz!-alist,n,point2,p,l,ans;
    % THE LUTZ!-ALIST is an association list of elements of the form
    % [X-value].([Y-value].[value of N for this point])
    % See thesis, chapter 7, algorithm LUTZ!-NAGELL, step [1].
  x:=car point;
  y:=cdr point;
  if !*tra or !*trmin then <<
    printc "Point to have torsion investigated is";
    printsq x;
    printsq y >>;
  a:=car equation;
  b:=cdr equation;
  if denr y neq 1 then <<
    l:=denr y;
    % we can in fact make l an item whose cube is > denr y.
    y:=!*multsq(y,!*exptf(l,3) ./ 1);
    x:=!*multsq(x,!*exptf(l,2) ./ 1);
    a:=!*multsq(a,!*exptf(l,4) ./ 1);
    b:=!*multsq(b,!*exptf(l,6) ./ 1) >>;
  if denr x neq 1 then <<
    l:=denr x;
    % we can in fact make l an item whose square is > denr x.
    y:=!*multsq(y,!*exptf(l,3) ./ 1);
    x:=!*multsq(x,!*exptf(l,2) ./ 1);
    a:=!*multsq(a,!*exptf(l,4) ./ 1);
    b:=!*multsq(b,!*exptf(l,6) ./ 1) >>;
  % we now have integral co-ordinates for x,y.
  lutz!-alist:=list (x . (y . 0));
  if (x neq car point) and (!*tra or !*trmin) then <<
    printc "Point made integral as ";
    printsq x;
    printsq y;
    printc "on the curve with coefficients";
    printsq a;
    printsq b >>;
  point:=x.y;
  equation:=a.b;
  n:=0;
loop:
  n:=n+1;
  point2:=elliptic!-multiply(x.y,2,equation);
  x:=car point2;
  y:=cdr point2;
  if infinitep x
    then return 2**n;
  if denr x neq 1
    then go to special!-denr;
  if a:=assoc(x,lutz!-alist)
    then if y = cadr a
      then return (ans:=lutz!-reduce(point,equation,2**n-2**(cddr a)))
      else if null numr !*addsq(y,cadr a)
        then return (ans:=lutz!-reduce(point,equation,2**n+2**(cddr a)))
        else interr "Cannot have 3 points here";
  lutz!-alist:=(x.(y.n)).lutz!-alist;
  if ans
    then return ans;
  go to loop;
special!-denr:
  p:=denr x;
  if not jhd!-primep p
    then return 'infinite;
  n:=1;
  n:=1;
loop2:
  point:=elliptic!-multiply(point,p,equation);
  n:=n*p;
  if infinitep car point
    then return n;
  if quotf(p,denr car point)
    then go to loop2;
  return 'infinite
  end;


symbolic procedure lutz!-reduce(point,equation,power);
begin
  scalar n;
  if !*tra or !*trmin then <<
    princ "Point is of order dividing ";
    printc power >>;
  n:=1;
  while evenp power do <<
    power:=power/2;
    n:=n*2;
    point:=elliptic!-add(point,point,equation) >>;
    % we know that all the powers of 2 must appear in the answer.
  if power = 1
    then return n;
  if jhd!-primep power
    then return n*power;
  return n*lutz!-reduce2(point,equation,power,3)
  end;



symbolic procedure lutz!-reduce2(point,equation,power,prime);
if power = 1
  then if infinitep car point
    then 1
    else nil
  else if infinitep car point
    then power
    else begin
      scalar n,prime2,u,ans;
      n:=0;
      while zerop cdr divide(power,prime) do <<
        n:=n+1;
        power:=power/prime >>;
      prime2:=nextprime prime;
      for i:=0:n do <<
        u:=lutz!-reduce2(point,equation,power,prime2);
        if u
          then <<
              ans:=u*prime**i;
              i:=n >>
         else <<
          power:=power*prime;
          point:=elliptic!-multiply(point,prime,equation) >> >>;
      if ans
        then return ans
        else return nil
      end;

endmodule;


module nbasis;

% Author: James H. Davenport.

fluid '(nestedsqrts sqrt!-intvar taylorasslist);

global '(!*tra);

exports normalbasis;
imports substitutesq,taylorform,printsq,newplace,sqrtsinsq,union,
        sqrtsign,interr,vecsort,mapvec,firstlinearrelation,mksp,multsq,
        !*multsq,addsq,removecmsq,antisubs,involvesq;


symbolic procedure normalbasis(zbasis,x,infdegree);
begin
  scalar n,nestedsqrts,sqrts,u,v,w,li,m,lam,i,inf,basis,save;
  save:=taylorasslist;
  inf:=list list(x,'quotient,1,x);
  n:=upbv zbasis;
  basis:=mkvect n;
  lam:=mkvect n;
  m:=mkvect n;
  goto  a;
square:
  sqrts:=nil;
  inf:=append(inf,list list(x,'expt,x,2));
  % we were in danger of getting sqrt(x) where we didnt want it.
a:
  newplace(inf);
  for i:=0:n do <<
    v:=substitutesq(getv(zbasis,i),inf);
    putv(basis,i,v);
    sqrts:=union(sqrts,sqrtsinsq(v,x)) >>;
  if !*tra then <<
    princ "Normal integral basis reduction with the";
    printc " following sqrts lying over infinity:";
    superprint sqrts >>;
  if member(list('sqrt,x),sqrts)
    then goto square;
  sqrts:=sqrtsign(sqrts,x);
  if iadd1 n neq length sqrts
    then interr "Length mismatch in normalbasis";
  for i:=0:n do <<
    v:=cl8roweval(getv(basis,i),sqrts);
    putv(m,i,cdr v);
    putv(lam,i,car v) >>;
reductionloop:
  vecsort(lam,list(basis,m));
  if !*tra then <<
    printc "Matrix before a reduction step at infinity is:";
    mapvec(m,function printc) >>;
  v:=firstlinearrelation(m,iadd1 n);
  if null v
    then goto ret;
  i:=n;
  while null numr getv(v,i) do
    i:=isub1 i;
  li:=getv(lam,i);
  w:=nil ./ 1;
  for j:=0:i do
    w:=addsq(w,!*multsq(getv(basis,j),
                 multsq(getv(v,j),1 ./  !*fmksp(x,-li+getv(lam,j)) )));
           % note the change of sign. my x is coates 1/x at this point!.
  if !*tra then <<
    princ "Element ";
    princ i;
    printc " replaced by the function printed below:" >>;
  w:=removecmsq w;
  putv(basis,i,w);
  w:=cl8roweval(w,sqrts);
  if car w <= li
    then interr "Normal basis reduction did not work";
  putv(lam,i,car w);
  putv(m,i,cdr w);
  goto reductionloop;
ret:
  newplace list (x.x);
  u:= 1 ./ !*p2f mksp(x,1);
  inf:=antisubs(inf,x);
  u:=substitutesq(u,inf);
  m:=nil;
  for i:=0:n do begin
    v:=getv(lam,i)-infdegree;
    if v < 0
      then goto next;
    w:=substitutesq(getv(basis,i),inf);
    for j:=0:v do <<
      if not involvesq(w,sqrt!-intvar)
        then m:=w.m;
      w:=!*multsq(w,u) >>;
  next:
    end;
  tayshorten save;
  return m
  end;


symbolic procedure !*fmksp(x,i);
% sf for x**i.
if i iequal 0
  then 1
  else !*p2f mksp(x,i);


symbolic procedure cl8roweval(basiselement,sqrts);
begin
  scalar lam,row,i,v,minimum,n;
  n:=isub1 length sqrts;
  lam:=mkvect n;
  row:=mkvect n;
  i:=0;
  minimum:=1000000;
  while sqrts do <<
    v:=taylorform substitutesq(basiselement,car sqrts);
    v:=assoc(taylorfirst v,taylorlist v);
    putv(row,i,cdr v);
    v:=car v;
    putv(lam,i,v);
    if v < minimum
      then minimum:=v;
    i:=iadd1 i;
    sqrts:=cdr sqrts >>;
  if !*tra then <<
    princ "Evaluating ";
    printsq basiselement;
    printc lam;
    printc row >>;
  v:=1000000;
  for i:=0:n do <<
    v:=getv(lam,i);
    if v > minimum
      then putv(row,i,nil ./ 1) >>;
  return minimum.row
  end;

endmodule;


module places;
 
% Author: James H. Davenport.
 
fluid '(basic!-listofallsqrts
        basic!-listofnewsqrts
        intvar
        listofallsqrts
        listofnewsqrts
        sqrt!-intvar
        sqrt!-places!-alist
        sqrts!-in!-integrand);
 
exports getsqrtsfromplaces,sqrtsinplaces,get!-correct!-sqrts,basicplace,
        extenplace,equalplace,printplace;
 
 
 
% Function to manipulate places
% a place is stored as a list of substitutions
% substitutions (x.f(x)) define the algrbraic number
% of which this place is an extension,
% while places (f(x).g(x)) define the extension.
%    currently g(x( is list ('minus,f(x))
%       or similar,e.g. (sqrt(sqrt x)).(sqrt(-sqrt x)).

 
 
% Given a list of places, produces a list of all
% the SQRTs in it that depend on INTVAR.
symbolic procedure getsqrtsfromplaces places;
  % The following loop finds all the SQRTs for a basis,
  % taking account of BASICPLACEs.
begin
  scalar basis,v,b,c,vv;
  for each u in places do <<
    v:=antisubs(basicplace u,intvar);
    vv:=sqrtsinsq (substitutesq(!*kk2q intvar,v),intvar);
      % We must go via SUBSTITUTESQ to get parallel
      % substitutions performed correctly.
    if vv
      then vv:=simp argof car vv;
    for each w in extenplace u do <<
      b:=substitutesq(simp lsubs w,v);
      b:=delete(sqrt!-intvar,sqrtsinsq(b,intvar));
      for each u in b do
        for each v in delete(u,b) do
          if dependsp(v,u)
            then b:=delete(u,b);
            % remove all the "inner" items, since they will
            % be accounted for anyway.
      if length b iequal 1
        then b:=car b
 else b:=mvar numr simpsqrtsq mapply(function !*multsq,
                                for each u in b collect simp argof u);
      if vv and not (b member sqrts!-in!-integrand)
        then <<
          c:=numr multsq(simp argof b,vv);
          c:=car sqrtsinsf(simpsqrt2 c,nil,intvar);
   if c member sqrts!-in!-integrand
            then b:=c >>;
      if not (b member basis)
        then basis:=b.basis >> >>;
  % The following loop deals with the annoying case of, say,
  % (X DIFFERENCE X 1) (X EXPT X 2) which should give rise to
  % SQRT(X-1).
  for each u in places do begin
    v:=cdr u;
    if null v or (car rfirstsubs v neq 'expt)
      then return;
    u:=simp!* subst(list('minus,intvar),intvar,rfirstsubs u);
    while v and (car rfirstsubs v eq 'expt) do <<
      u:=simpsqrtsq u;
      v:=cdr v;
      basis:=union(basis,delete(sqrt!-intvar,sqrtsinsq(u,intvar))) >>
    end;
  return remove!-extra!-sqrts basis
  end;
 
 
 
symbolic procedure sqrtsinplaces u;
% Note the difference between this procedure and
% the previous one: this one does not take account
% of the BASICPLACE component (& is pretty useless).
if null u
  then nil
  else sqrtsintree(for each v in car u collect lsubs v,
                   intvar,
                   sqrtsinplaces cdr u);
 
 
 
%symbolic procedure placesindiv places;
% Given a list of places (i.e. a divisor),
% produces a list of all the SQRTs on which the places
% explicitly depend.
%begin scalar v;
%  for each u in places do
%    for each uu in u do
%      if not (lsubs uu member v)
%        then v:=(lsubs uu) . v;
%  return v
%  end;

 
 
symbolic procedure get!-correct!-sqrts u;
% u is a basicplace.
begin
  scalar v;
  v:=assoc(u,sqrt!-places!-alist);
  if v
    then <<
      v:=cdr v;
      listofallsqrts:=cdr v;
      listofnewsqrts:=car v
      >>
    else <<
      listofnewsqrts:=basic!-listofnewsqrts;
      listofallsqrts:=basic!-listofallsqrts
      >>;
  return nil
  end;
 
 
 
%symbolic procedure change!-place(old,new);
%% old and new are basicplaces;
%begin
%  scalar v;
%  v:=assoc(new,sqrt!-places!-alist);
%  if v
%    then sqrtsave(cddr v,cadr v,old)
%    else <<
%      listofnewsqrts:=basic!-listofnewsqrts;
%      listofallsqrts:=basic!-listofallsqrts
%      >>;
%  return nil
%  end;

 
 
symbolic procedure basicplace(u);
% Returns the basic part of a place.
if null u
  then nil
  else if atom caar u
    then (car u).basicplace cdr u
    else nil;
 
 
 
symbolic procedure extenplace(u);
% Returns the extension part of a place.
if u and atom caar u
  then extenplace cdr u
  else u;
 
 
 
symbolic procedure equalplace(a,b);
% Sees if two extension places represent the same place or not.
if null a
  then if null b
    then t
    else nil
  else if null b
    then nil
    else if member(car a,b)
      then equalplace(cdr a,delete(car a,b))
      else nil;
 
 
 
symbolic procedure remove!-extra!-sqrts basis;
begin
  scalar basis2,save;
  save:=basis2:=for each u in basis collect !*q2f simp argof u;
  for each u in basis2 do
    for each v in delete(u,basis2) do
      if quotf(v,u)
        then basis2:=delete(v,basis2);
  if basis2 eq save
    then return basis
    else return for each u in basis2 collect list('sqrt,prepf u)
  end;
 
 
 
symbolic procedure printplace u;
begin
  scalar a,n,v;
  a:=rfirstsubs u;
  princ (v:=lfirstsubs u);
  princ "=";
  if atom a
    then princ "0"
    else if (car a eq 'quotient) and (cadr a=1)
      then princ "infinity"
      else <<
 n:=negsq addsq(!*kk2q v,negsq simp!* a);
% NEGSQ added JHD 22.3.87 - the previous value was wrong.
% If the substitution is (X-v) then this takes -v to 0,
% so the place was at -v.
        if (numberp numr n) and (numberp denr n)
          then <<
            princ numr n;
            if not onep denr n
              then <<
                princ " / ";
                princ denr n >> >>
          else <<
            if degreein(numr n,intvar) > 1
             then printc "Any root of:";
            printsq n;
            if cdr u
              then princ "at the place " >> >>;
  u:=cdr u;
  if null u
    then goto nl!-return;
  n:=1;
  while u and (car rfirstsubs u eq 'expt) do <<
    n:=n * caddr rfirstsubs u;
    u:=cdr u >>;
  if n neq 1 then <<
    terpri!* nil;
    prin2 " ";
    princ v;
    princ "=>";
    princ v;
    princ "**";
    princ n >>;
  while u do <<
    if car rfirstsubs u eq 'minus
      then princ "-"
      else princ "+";
    u:=cdr u >>;
nl!-return:
  terpri();
  return
  end;
 
 
 
symbolic procedure degreein(sf,var);
if atom sf
  then 0
  else if mvar sf eq var
    then ldeg sf
    else max(degreein(lc sf,var),degreein(red sf,var));
 
endmodule;


module precoats;

% Author: James H. Davenport.

fluid '(basic!-listofallsqrts
        basic!-listofnewsqrts
        sqrt!-intvar
        taylorvariable
        thisplace);

global '(!*tra);

exports precoates;
imports mksp,algint!-subf,subzero2,substitutesq,removeduplicates,
        printsq,basicplace,extenplace,interr,get!-correct!-sqrts,
        printplace,simptimes,subzero,negsq,addsq,involvesq,taylorform,
        taylorevaluate,mk!*sq,!*exptsq,!*multsq,!*invsq,sqrt2top,
        jfactor,sqrtsave,antisubs;


symbolic procedure infsubs(w);
if caar w = thisplace
  then (cdar w).(cdr w)
  else (thisplace.(car w)).(cdr w);
% thisplace is (z quotient 1 z) so we are moving to infinity.


symbolic procedure precoates(residues,x,movedtoinfinity);
begin
  scalar answer,placeval,reslist,placelist,placelist2,thisplace;
  reslist:=residues;
  placelist:=nil;
  while reslist do <<
    % car reslist = <substitution list>.<value>;
    placeval:=algint!-subf((mksp(x,1) .* 1) .+ nil,caar reslist);
    if 0 neq cdar reslist
      then if null numr subzero2(denr placeval,x)
        then <<
          if null answer
            then answer:='infinity
            else if answer eq 'finite
              then answer:='mixed;
          if !*tra
            then printc "We have an residue at infinity" >>
        else <<
          if null answer
            then answer:='finite
            else if answer eq 'infinity
              then answer:='mixed;
          placelist:=placeval.placelist;
          if !*tra
            then printc "This is a finite residue" >>;
    reslist:=cdr reslist >>;
  if answer eq 'mixed
    then return answer;
  if answer eq 'infinity
    then <<
      thisplace:=list(x,'quotient,1,x);
      % maps x to 1/x.
      answer:=precoates(for each u in residues collect infsubs u,x,t);
                % derivative of 1/x is -1/x**2.
      if atom answer
        then return answer
        else return substitutesq(answer,list(thisplace)) >>;
  placelist2:=removeduplicates placelist;
  answer := 1 ./ 1;
  % the null divisor.
  if !*tra then <<
    printc "The divisor has elements at:";
    mapcar(placelist2,function printsq) >>;
  while placelist2 do begin
    scalar placelist3,extrasubs,u,bplace;
    % loop over all distinct places.
    reslist:=residues;
    placelist3:=placelist;
    placeval:=nil;
    while reslist do <<
      if car placelist2 = car placelist3
        then <<
          placeval:=(cdar reslist).placeval;
          thisplace:= caar reslist;
          % the substitutions defining car placelist.
          u:=caar reslist;
          bplace:=basicplace u;
          u:=extenplace u;
          extrasubs:=u.extrasubs >>;
      reslist:=cdr reslist;
      placelist3:=cdr placelist3 >>;
    % placeval is a list of all the residues at this place.
    if !*tra then <<
      princ "List of multiplicities at this place:";
      printc placeval;
      princ "with substitutions:";
      superprint extrasubs >>;
    if 0 neq mapply(function plus2,placeval)
      then interr "Divisor not effective";
    get!-correct!-sqrts bplace;
    u:=pbuild(x,extrasubs,placeval);
    sqrtsave(basic!-listofallsqrts,basic!-listofnewsqrts,bplace);
    if atom u
      then <<
        placelist2:=nil;
        % set to terminate loop.
        answer:=u >>
      else <<
        answer:=substitutesq(!*multsq(answer,u),antisubs(thisplace,x));
        placelist2:=cdr placelist2 >>
    end;
    % loaded in pbuild to check for poles at the correct places.
  return answer
  end;



symbolic procedure dlist(u);
% Given a list of lists,converts to a list.
if null u
  then nil
  else if null car u
    then dlist cdr u
    else append(car u,dlist cdr u);


symbolic procedure debranch(extrasubs,reslist);
begin
  scalar substlist;
  % remove spurious substitutions.
  for each u in dlist extrasubs do
    if not ((car u) member substlist)
      then substlist:=(car u).substlist;
  % substlist is a list of all the possible substitutions).
  while substlist do
    begin scalar tsqrt,usqrt;
      scalar with1,with2,without1,without2,wres;
    scalar a1,a2,b1,b2;
    % decide if tsqrt is redundant.
    tsqrt:=car substlist;
    substlist:=cdr substlist;
    wres:=reslist;
    for each place in extrasubs do <<
      usqrt:=assoc(tsqrt,place);
        % usqrt is s.s' or s.(minus s').
      if null usqrt
        then interr "Places not all there";
      if cadr usqrt eq 'sqrt
        then<<
          with2:=(car wres).with2;
          with1:=delete(usqrt,place).with1>>
        else<<
          if not (cadr usqrt eq 'minus)
            then interr "Ramification format error";
          without2:=(car wres).without2;
          without1:=delete(usqrt,place).without1 >>;
      wres:=cdr wres>>;
    % first see if one item appears passim.
    if null with1
      then go to itswithout;
    if null without1
      then go to itswith;
    % Now must see if WITH2 matches WITHOUT2 in order WITH1/WITHOUT1.
    a1:=with1;
    a2:=with2;
  outerloop:
    b1:=without1;
    b2:=without2;
  innerloop:
    if (car a1) = (car b1)
      then << if (car a2) neq (car b2)
           then return;
           else go to outeriterate >>;
    b1:=cdr b1;
    b2:=cdr b2;
    if null b1
      then return
      else go to innerloop;
      % null b1 => lists do not match at all.
  outeriterate:
    a1:=cdr a1;
    a2:=cdr a2;
    if a1
      then go to outerloop;
    if !*tra then <<
      princ "Residues reduce to:";
      printc without2;
      printc "at ";
      mapc(without1,function printplace) >>;
    extrasubs:=without1;
    reslist:=without2;
    return;
  itswithout:
    % everything is in the "without" list.
    with1:=without1;
    with2:=without2;
  itswith:
    % remove usqrt from the with lists.
    extrasubs:=for each u in with1 collect delete(assoc(tsqrt,u),u);
    if !*tra then <<
      printc "The following appears throughout the list ";
      printc tsqrt >>;
    reslist:=with2
    end;
  return extrasubs.reslist
  end;


symbolic procedure pbuild(x,extrasubs,placeval);
begin
  scalar multivals,u,v,answer;
  u:=debranch(extrasubs,placeval);
  extrasubs:=car u;
  placeval:=cdr u;
  % remove spurious entries.
  if (length car extrasubs) > 1
    then return 'difficult;
  % hard cases not allowed for.
  multivals:=mapcar(dlist extrasubs,function car);
  u:=simptimes removeduplicates multivals;
  answer:= 1 ./ 1;
    while extrasubs do <<
      v:=substitutesq(u,car extrasubs);
      v:=addsq(u,negsq subzero(v,x));
      v:=mkord1(v,x);
      if !*tra then <<
        princ "Required component is ";
        printsq v >>;
      answer:=!*multsq(answer,!*exptsq(v,car placeval));
      % place introduced with correct multiplicity.
      extrasubs:=cdr extrasubs;
      placeval:=cdr placeval >>;
  if length jfactor(denr sqrt2top !*invsq  answer,x) > 1
    then return 'many!-poles
    else return answer
  end;


symbolic procedure findord(v,x);
begin
  scalar nord,vd;
  %given v(x) with v(0)=0, makes v'(0) nonzero.
  nord:=0;
  taylorvariable:=x;
  while involvesq(v,sqrt!-intvar) do
    v:=substitutesq(v,list(x.list('expt,x,2)));
  vd:=taylorform v;
loop:
  nord:=nord+1;
  if null numr taylorevaluate(vd,nord)
    then go to loop;
  return nord
  end;


symbolic procedure mkord1(v,x);
begin
  scalar nord;
  nord:=findord(v,x);
  if nord iequal 1
    then return v;
  if !*tra then <<
    princ "Order reduction: ";
    printsq v;
    princ "from order ";
    princ nord;
    printc " to order 1" >>;
  % Note that here we do not need to simplify, since SIMPLOG will
  % remove all these SQRTs or EXPTs later.
  return !*p2q mksp(list('nthroot,mk!*sq v,nord),1)
  end;

endmodule;


module primes;

% Author: James H. Davenport.

exports nextprime,jhd!-primep;

symbolic procedure nextprime p;
% Returns the next prime number bigger than p.
    if p=0 then 1
    else if p=1 then 2
    else begin
        if evenp p then p:=p+1 else p:=p+2;
 test:  if jhd!-primep p then return p;
        p:=p+2;
        go to test end;

symbolic procedure jhd!-primep p;
    if p < 4 then t
    else if evenp p then nil
    else begin
      scalar n;
      n:=3; %trial factor.
 top: if n*n>p then return t
      else if remainder(p,n)=0 then return nil;
      n:=n+2;
      go to top end;

endmodule;


module removecm;  % Routines to remove constant factors from expresions.

% Author: James H. Davenport.

fluid '(intvar);

% New improved REMOVECOMMOMMULTIPLES routines.
% These routines replace a straightforward pair with GCDF instead of
% CMGCDF and its associates.  The saving is large in complicated
% expressions (in the "general point of order 7" calculations, they
% exceeded 90% in some cases, being 1.5 secs as opposed to > 15 secs.).
% They are about 1K larger, but this seems a small price to pay.

exports removecmsq,removeconstantsf;
imports ordop,addf,gcdn,gcdf,gcdk,involvesf,dependsp,makemainvar,quotf;

symbolic procedure removecmsq sq;
(removecmsf numr sq) ./ (removecmsf denr sq);

symbolic procedure removecmsf sf;
if atom sf or not ordop(mvar sf,intvar) or not involvesf(sf,intvar)
  then if sf
    then 1
    else nil
  else if null red sf
    then if dependsp(mvar sf,intvar)
      then (lpow sf .* removecmsf lc sf) .+ nil
      else removecmsf lc sf
    else begin
      scalar u,v;
      % The general principle here is to find a (non-INTVAR-depending)
      % coefficient of a purely INTVAR-depending monomial, and then
      % perform a g.c.d. to discover that factor of this which is a CM.
      u:=sf;
      while (v:=involvesf(u,intvar)) do u:=lc makemainvar(u,v);
      if u iequal 1
        then return sf;
      return quotf(sf,cmgcdf(sf,u))
      end;

symbolic procedure cmgcdf(sf,u);
if numberp u
  then if atom sf
    then if null sf
      then u
      else gcdn(sf,u)
    else if u = 1
      then 1
      else cmgcdf(red sf,cmgcdf(lc sf,u))
  else if atom sf
    then gcdf(sf,u)
    else if mvar u eq mvar sf
      then if ordop(intvar,mvar u)
        then gcdf(sf,u)
        else cmgcdf2(sf,u)
      else if ordop(mvar sf,mvar u)
        then cmgcdf(red sf,cmgcdf(lc sf,u))
        else cmgcdf(u,sf);

symbolic procedure remove!-maxdeg(sf,var);
if atom sf
  then 0
  else if mvar sf eq var
    then ldeg sf
    else if ordop(var,mvar sf)
      then 0
      else max(remove!-maxdeg(lc sf,var),remove!-maxdeg(red sf,var));

symbolic procedure cmgcdf2(sf,u);
% SF and U have the same MVAR, but INTVAR comes somewhere
% down in SF.  Therefore we can do better than a straight
% GCDK, or even a straight MAKEMAINVAR.
begin
  scalar n;
  n:=remove!-maxdeg(sf,intvar);
  if n = 0
    then return gcdf(sf,u);
    % Doesn't actually depend on INTVAR.
loop:
  if u = 1
    then return 1;
  u:=gcdf(u,collectterms(sf,intvar,n));
  n:=isub1 n;
  if n < 0
    then return u
    else go loop
  end;

symbolic procedure collectterms(sf,var,n);
if atom sf
  then if n = 0
    then sf
    else nil
  else if mvar sf eq var
    then if ldeg sf = n
      then lc sf
      else if ldeg sf > n
        then collectterms(red sf,var,n)
        else nil
    else if ordop(var,mvar sf)
      then if n = 0
        then sf
        else nil
      else begin
        scalar v,w;
        v:=collectterms(lc sf,var,n);
        w:=collectterms(red sf,var,n);
        if null v
          then return w
          else return addf(w,(lpow sf .* v) .+ nil)
        end;

symbolic procedure removeconstantsf sf;
% Very simple version for now.
begin
  scalar u;
  if null sf
    then return nil
    else if atom sf
      then return 1;
  while (null red sf) and (remove!-constantp mvar sf) do
    sf:=lc sf;
  u:=remove!-const!-content sf;
  if u = 1
    then return sf
    else return quotf!*(sf,u)
  end;

symbolic procedure remove!-constantp pf;
if numberp pf
  then t
  else if atom pf
    then nil
    else if car pf eq 'sqrt
      then remove!-constantp argof pf
      else if (car pf eq 'expt) or (car pf eq 'quotient)
        then (remove!-constantp argof pf)
             and (remove!-constantp caddr pf)
        else nil;

symbolic procedure remove!-const!-content sf;
if numberp sf
  then sf
  else if null red sf
    then if remove!-constantp mvar sf
      then (lpow sf .* remove!-const!-content lc sf) .+ nil
      else remove!-const!-content lc sf
    else begin
      scalar u;
      u:=remove!-const!-content lc sf;
      if u = 1
        then return u;
      return gcdf(u,remove!-const!-content red sf)
      end;

endmodule;


module sqfrnorm;

% Author: James H. Davenport.

fluid '(!*pvar listofallsqrts);

global '(modevalcount);

modevalcount:=1;

exports sqfr!-norm2,res!-sqrt;

%symbolic procedure resultant(u,v);
%begin
%  scalar maxdeg,zeroes,ldegu,ldegv,m;
%  % we can have gone makemainvar on u and v;
%  ldegu:=ldeg u;
%  ldegv:=ldeg v;
%  maxdeg:=isub1 max2(ldegu,ldegv);
%  zeroes:=nlist(nil,maxdeg);
%  u:=remake(u,mvar u,ldegu);
%  v:=remake(v,mvar v,ldegv);
%  m:=nil;
%  ldegu:=isub1 ldegu;
%  ldegv:=isub1 ldegv;
%  for i:=0 step 1 until ldegv do
%    m:=append(ncdr(zeroes,maxdeg-ldegv+i),
%              append(u,ncdr(zeroes,maxdeg-i))).m;
%  for i:=0 step 1 until ldegu do
%    m:=append(ncdr(zeroes,maxdeg-ldegu+i),
%              append(v,ncdr(zeroes,maxdeg-i))).m;
%  return detqf m
%  end;


%symbolic procedure remake(u,v,w);
%% remakes u into a list of sf's representing its coefficients;
%if w iequal 0 then list u
%  else if (pairp u) and (mvar u eq v) and (ldeg u iequal w)
%    then (lc u).remake(red u,v,isub1 w)
%    else (nil ).remake(    u,v,isub1 w);

%fluid '(n); %needed for the mapcar;

%symbolic procedure detqf u;
%   %u is a square matrix standard form.
%%  %value is the determinant of u.
%%  %algorithm is expansion by minors of first row/column;
%   begin integer n;
%   scalar x,y,z;
%        if length u neq length car u then rederr "Non square matrix"
%         else if null cdr u then return caar u;
%        if length u < 3
%          then go to noopt;
%        % try to remove a row with only one non-zero in it;
%        z:=1;
%        x:=u;
%      loop:
%        n:=posnnonnull car x;
%        if n eq t
%          then return nil;
%        % special test for all null;
%        if n then <<
%          y:=nth(car x,n);
%          % next line is equivalent to:
%%           onne of n,z is even;
%          if evenp (n+z-1)
%            then y:=negf y;
%          u:=remove(u,z);
%          return !*multf(y,detqf remove2 u) >>;
%       x:=cdr x;
%       z:=z+1;
%       if x
%         then go to loop;
%     noopt:
%        x := u;
%        n := 1;                 %number of current row/column;
%        z := nil;
%        if nonnull car u < nonnullcar u
%         then go to row!-expand;
%        u:=mapcar(u,function cdr);
%    a:  if null x then return z;
%        y := caar x;
%        if null y then go to b
%         else if evenp n then y := negf y;
%        z := addf(!*multf(y,detqf remove(u,n)),z);
%    b:  x := cdr x;
%        n := iadd1 n;
%        go to a;
%      row!-expand:
%        u:=cdr u;
%        x:=car x;
%      aa:
%        if null x then return z;
%        y:=car x;
%        if null y
%          then go to bb
%          else if evenp n then y:=negf y;
%        z:=addf(!*multf(y,detqf remove2 u),z);
%      bb:
%        x:=cdr x;
%        n:=iadd1 n;
%        go to aa
%   end;
%
%
%symbolic procedure remove2 u;
%mapcar(u,function (lambda x;
%                    remove(x,n)));
%
%unfluid '(n);
%
%symbolic procedure nonnull u;
%if null u
%  then 0
%  else if null car u
%    then nonnull cdr u
%    else iadd1 (nonnull cdr u);
%
%
%symbolic procedure nonnullcar u;
%if null u
%  then 0
%  else if null caar u
%    then nonnullcar cdr u
%    else iadd1 (nonnullcar cdr u);
%
%
%
%symbolic procedure posnnonnull u;
%% returns t if u has no non-null elements
%% nil if more than one
%% else position of the first;
%begin
%  scalar n,x;
%  n:=1;
%loop:
%  if null u
%    then return
%      if x
%        then x
%        else t;
%  if car u
%    then if x
%      then return nil
%      else x:=n;
%  n:=iadd1 n;
%  u:=cdr u;
%  go to loop
%  end;


symbolic procedure res!-sqrt(u,a);
% Evaluates resultant of u ( as a poly in its mvar) and x**-a.
begin
  scalar x,n,v,k,l;
  x:=mvar u;
  n:=ldeg u;
  n:=quotient(n,2);
  v:=mkvect n;
  putv(v,0,1);
  for i:=1:n do
    putv(v,i,!*multf(a,getv(v,i-1)));
  % now substitute for x**2 in u leaving k*x+l.
  k:=l:=nil;
  while u do
    if mvar u neq x
      then <<
        l:=addf(l,u);
        u:=nil >>
      else <<
        if evenp ldeg u
          then l:=addf(l,!*multf(lc u,getv(v,(ldeg u)/2)))
          else k:=addf(k,!*multf(lc u,getv(v,(ldeg u -1)/2)));
        u:=red u >>;
  % now have k*x+l,x**2-a, giving l*l-a*k*k.
  return addf(!*multf(l,l),!*multf(negf a,multf(k,k)))
  end;


symbolic procedure sqfr!-norm2 (f,mvarf,a);
begin
  scalar u,w,aa,ff,resfn;
  resfn:='resultant;
  if eqcar(a,'sqrt)
    then <<
      resfn:='res!-sqrt;
      aa:=!*q2f simp argof a >>
    else rederr "Norms over transcendental extensions";
  f:=pvarsub(f,a,'! gerbil);
  w:=nil;
  if involvesf(f,'! gerbil) then goto l1;
increase:
  w:=addf(w,!*p2f mksp(a,1));
  f:=!*q2f algint!-subf(f,list(mvarf . list('plus,mvarf,
                                            list('minus,'! gerbil))));
l1:
  u:=apply(resfn,list(makemainvar(f,'! gerbil),aa));
  ff:=nsqfrp(u,mvarf);
  if ff
    then go to increase;
  f:=!*q2f algint!-subf(f,list('! gerbil.a));
  % cannot use pvarsub since want to squash higher powers.
  return list(u,w,f)
  end;

symbolic procedure nsqfrp(u,v);
begin
  scalar w;
  w:=modeval(u,v);
  if w eq 'failed
    then go to normal;
  if atom w
    then go to normal;
  if ldegvar(w,v) neq ldegvar(u,v)
    then go to normal;
%  printc "Modular image is:";
%  printsf w;
  w:=gcdf(w,partialdiff(w,v));
%  printc "Answer is:";
%  printsf w;
  if w iequal 1
    then return nil;
normal;
  w:=gcdf(u,partialdiff(u,v));
  if involvesf(w,v)
    then return w
    else return nil
  end;

symbolic procedure ldegvar(u,v);
if atom u
  then 0
  else if mvar u eq v
    then ldeg u
    else if ordop(v,mvar u)
      then 0
      else max2(ldegvar(lc u,v),ldegvar(red u,v));


symbolic procedure modeval(u,v);
if atom u
  then u
  else if v eq mvar u
    then begin
      scalar w,x;
      w:=modeval(lc u,v);
      if w eq 'failed
        then return w;
      x:=modeval(red u,v);
      if x eq 'failed
        then return x;
      if null w
        then return x
        else return (lpow u .* w) .+ x
      end
    else begin
      scalar w,x;
      x:=mvar u;
      if not atom x
        then if dependsp(x,v)
          then return 'failed;
      x:=modevalvar x;
      if x eq 'failed
        then return x;
      w:=modeval(lc u,v);
      if w eq 'failed
        then return w;
      if x
        then w:=multf(w,exptf(x,ldeg u));
      x:=modeval(red u,v);
      if x eq 'failed
        then return x;
      return addf(w,x)
      end;


symbolic procedure modevalvar v;
begin
  scalar w,x;
  if not atom v
    then go to alg;
  w:=get(v,'modvalue);
  if w
    then return w;
  put(v,'modvalue,modevalcount);
  modevalcount:=modevalcount+1;
  return modevalcount-1;
alg:
  if car v neq 'sqrt
    then rederr "Unexpected algebraic";
  if numberp argof v
    then return (mksp(v,1) .* 1) .+ nil;
  w:=modeval(!*q2f simp argof v,!*pvar);
  w:=assoc(w,listofallsqrts);
  % the variable does not matter, since we know that it does not depend.
  if w
    then return cdr w
    else return 'failed
  end;

% unglobal '(modevalcount);

endmodule;


module substns;

% Author: James H. Davenport.

exports xsubstitutep,xsubstitutesq,substitutevec,substitutesq,subzero,
        subzero2,pvarsub;


symbolic procedure xsubstitutep(pf,slist);
simp xsubstitutep2(pf,slist);


symbolic procedure xsubstitutep2(pf,slist);
if null slist
  then pf
  else xsubstitutep2(subst(rfirstsubs slist,
                           lfirstsubs slist,
                           pf),
                     cdr slist);


symbolic procedure xsubstitutesq(sq,slist);
substitutesq(substitutesq(sq,basicplace slist),extenplace slist);


symbolic procedure substitutevec(v,slist);
for i:=0:upbv v do
  putv(v,i,substitutesq(getv(v,i),slist));


symbolic procedure substitutesq(sq,slist);
begin
  scalar list2,nm;
  list2:=nil;
  while slist do <<
    if cdar slist iequal 0
      then <<
        if list2
          then sq:=substitutesq(sq,reversewoc list2);
        list2:=nil;
        sq:=subzero(sq,caar slist) >>
      else if not (caar slist = cdar slist)
        then if assoc(caar slist,list2)
          then list2:=for each u in list2 collect
                  (car u).subst(cdar slist,caar slist,cdr u)
          else list2:=(car slist).list2;
        % don't bother with the null substitution.
    slist:=cdr slist >>;
  list2:=reversewoc list2;
  if null list2
    then return sq;
  nm:=algint!-subf(numr sq,list2);
  if numr nm
    then nm:=!*multsq(nm,invsq algint!-subf(denr sq,list2));
  return nm
  end;

% standard interface.
symbolic procedure subzero(exprn,var);
begin
  scalar top;
  top:=subzero2(numr exprn,var);
  if null numr top
    then return nil ./ 1;
  return !*multsq(top,!*invsq subzero2(denr exprn,var))
  end;


symbolic procedure subzero2(sf,var);
if not involvesf(sf,var)
  then sf ./ 1
  else if var eq mvar sf
    then subzero2(red sf,var)
    else if ordop(var,mvar sf)
      then sf ./ 1
      else begin
        scalar u,v;
        if dependsp(mvar sf,var)
          then <<
            u:=simp subst(0,var,mvar sf);
            if numr u
              then u:=!*exptsq(u,ldeg sf) >>
          else u:=((lpow sf .* 1) .+ nil) ./ 1;
        if null numr u
          then return subzero2(red sf,var);
        v:=subzero2(lc sf,var);
        if null numr v
          then return subzero2(red sf,var);
        return !*addsq(subzero2(red sf,var),
                       !*multsq(u,v))
        end;



symbolic procedure pvarsub(f,u,v);
% Changes u to v in polynomial f. No proper substitutions at all.
if atom f
  then f
  else if mvar f equal u
    then addf(multf(lc f,!*p2f mksp(v,ldeg f)),
              pvarsub(red f,u,v))
    else if ordop(u,mvar f)
      then f
      else addf(multf(pvarsub(lc f,u,v),!*p2f lpow f),
                pvarsub(red f,u,v));

endmodule;


module taylor;

% Author: James H. Davenport.

fluid '(const taylorasslist taylorvariable);

exports taylorform,taylorformp,taylorevaluate,return0,taylorplus,
         initialtaylorplus,taylorminus,initialtaylorminus,
         tayloroptminus,tayloroptplus,taylorctimes,initialtaylortimes,
         tayloroptctimes,taylorsqrtx,initialtaylorsqrtx,
         taylorquotient,initialtaylorquotient,taylorformersqrt,
         taylorbtimes,taylorformertimes,taylorformerexpt;

 symbolic procedure taylorform sq;
 if involvesf(denr sq,taylorvariable)
   then taylorformp list('quotient,tayprepf numr sq,tayprepf denr sq)
   else if 1 iequal denr sq
     then taylorformp tayprepf numr sq
     else taylorformp list('constanttimes,
                           tayprepf numr sq,
                           mk!*sq(1 ./ (denr sq)));
 % get division by a constant right.


 symbolic procedure taylorformp pf;
 if null pf
   then nil
   else if not dependsp(pf,taylorvariable)
     then taylorconst simp pf
     else begin
       scalar fn,initial,args;
       if atom pf
         then if pf eq taylorvariable
           then return taylorformp list ('expt,pf,1)
           else interr "False atom in taylorformp";
       % get 'x right as reduce shorthand for x**1.
       if taylorp pf
         then return pf;
       % cope with pre-expressed cases.
       % ***store-hack-1***
       % remove the (car pf eq 'sqrt) if more store is available.
       if (car pf eq 'sqrt) and
          (fn:=assoc(pf,taylorasslist))
         then go to lookupok;
       % look it up first.
       fn:=get(car pf,'taylorformer);
       if null fn
         then go to ordinary;
       fn:=apply(fn,list cdr pf);
       % ***store-hack-1***
       % remove the test if more store is available.
       if car pf eq 'sqrt
         then taylorasslist:=(pf.fn).taylorasslist;
       return fn;
       % cope with the special cases.
     ordinary:
       args:=mapcar(cdr pf,function taylorformp);
       fn:=get(car pf,'tayloropt);
       if null fn
         then go to nooptimisation;
       fn:=apply(fn,list args);
       if fn
         then go to ananswer;
       % an optimisation has been made.
     nooptimisation:
       fn:=get(car pf,'taylorfunction);
       if null fn
         then interr "No Taylor function provided";
       fn:=fn.args;
       % fn is now the "how to compute" code.
       initial:=get(car pf,'initialtaylorfunction);
       if null initial
         then interr "No initial Taylor function";
       initial:=apply(initial,
                      list for each u in cdr fn collect firstterm u);
       % the first term in the expansion.
       fn:=list(fn,(car initial).(car initial),initial);
     ananswer:
       % ***store-hack-1***
       % uncomment this if more store is available;
       % taylorasslist:=(pf.fn).taylorasslist;
       return fn;
     lookupok:
       % These PRINT statements can be enabled in order to test the
       % efficacy of the association list
 %      printc "Taylor lookup succeeded";
 %      superprint car fn;
 %      printc length taylorasslist;
       return cdr fn
       end;


 symbolic procedure taylorevaluate(texpr,n);
 if n<taylorfirst texpr
   then nil ./ 1
   else if n>taylorlast texpr
     then tayloreval2(texpr,n)
     else begin
       scalar u;
       u:=assoc(n,taylorlist texpr);
       if u
         then return cdr u
         else return tayloreval2(texpr,n)
       end;


 symbolic procedure tayloreval2(texpr,n);
 begin
   scalar u;
   % actually evaluates from scratch.
   u:=apply(taylorfunction texpr, list(n,texpr,cdr taylordefn texpr));
   if 'return0 eq taylorfunction texpr
     then return u;
   % no need to update with trivial zeroes.
   rplacd(cdr texpr,(n.u).taylorlist texpr);
   % update the association list.
   if n>taylorlast texpr
     then rplacd(taylornumbers texpr,n);
   % update the first/last pointer.
   return u
   end;


 symbolic procedure taylorconst sq;
 list('return0 . nil,0 . 0,0 . sq);


 symbolic procedure return0 (a,b,c);
 nil ./ 1;

 flag('(return0),'taylor);


 symbolic procedure firstterm texpr;
 begin
   scalar n,i;
   i:=taylorfirst texpr;
 trynext:
   n:=taylorevaluate(texpr,i);
   if numr n
     then return i.n;
   if i > 50
     then interr "Potentially zero Taylor series";
   i:=iadd1 i;
   rplaca(taylornumbers texpr,i);
   go to trynext
   end;


 symbolic procedure tayloroneterm u;
 % See if a Taylor expression has only one term.
  'return0 eq taylorfunction u and taylorfirst u=taylorlast u;


 % ***store-hack-1***;
 % uncomment this procedure if more store is available;
 % there is a smacro for this at the start of the file
 % for use if no store can be spared;
 %symbolic procedure tayshorten(save);
 %begin
 %  scalar z;
 %  % shortens the association list back to save,
 %    removing all the non-sqrts from it;
 %  while taylorasslist neq save do <<
 %    if caar taylorasslist eq 'sqrt
 %      then z:=(car taylorasslist).z;
 %    taylorasslist:=cdr taylorasslist >>;
 %  taylorasslist:=nconc(z,taylorasslist);
 %  return nil
 %  end;


 symbolic procedure tayprepf sf;
 if atom sf
   then sf
   else if atom mvar sf
     then taylorpoly makemainvar(sf,taylorvariable)
     else if null red sf
       then tayprept lt sf
       else list('plus,tayprept lt sf,tayprepf red sf);


 symbolic procedure tayprept term;
 if tdeg term = 1
   then if tc term = 1
     then tvar term
     else list('times,tvar term,tayprepf tc term)
   else if tc term = 1
     then list ('expt,tvar term,tdeg term)
     else list('times,list('expt,tvar term,tdeg term),
                    tayprepf tc term);


 symbolic procedure taylorpoly sf;
 % SF is a poly with MVAR = TAYLORVARIABLE.
 begin
   scalar tmax,tmin,u;
   tmax:=tmin:=ldeg sf;
   while sf do
     if atom sf or (mvar sf neq taylorvariable)
       then <<
         tmin:=0;
         u:=(0 . !*f2q sf).u;
         sf:=nil >>
       else <<
         u:=((tmin:=ldeg sf) . !*f2q lc sf) . u;
         sf:=red sf >>;
   return (list 'return0) . ((tmin.tmax).u)
   end;


 symbolic procedure taylorplus(n,texpr,args);
 mapply(function addsq,
        for each u in args collect taylorevaluate(u,n));


 symbolic procedure initialtaylorplus slist;
 begin
   scalar n,numlst;
   n:=mapply(function min2,mapcar(slist,function car));
   % the least of the degrees.
   numlst:=nil;
   while slist do <<
     if caar slist iequal n
       then numlst:=(cdar slist).numlst;
     slist:=cdr slist >>;
   return n.mapply(function addsq,numlst)
   end;


 put ('plus,'taylorfunction,'taylorplus);
 put ('plus,'initialtaylorfunction,'initialtaylorplus);


 symbolic procedure taylorminus(n,texpr,args);
 negsq taylorevaluate(car args,n);


 symbolic procedure initialtaylorminus slist;
 (caar slist).(negsq cdar slist);


 put('minus,'taylorfunction,'taylorminus);
 put('minus,'initialtaylorfunction,'initialtaylorminus);


 flag('(taylorplus taylorminus),'taylor);


 symbolic procedure tayloroptminus(u);
 if 'return0 eq taylorfunction car u
   then taylormake(taylordefn car u,
                   taylornumbers car u,
                   taylorneglist taylorlist car u)
   else if 'taylorctimes eq taylorfunction car u
     then begin
       scalar const;
       u:=car u;
       const:=caddr taylordefn u;
       % the item to be negated.
       const:=taylormake(taylordefn const,
                         taylornumbers const,
                         taylorneglist taylorlist const);
       return taylormake(list(taylorfunction u,
                              argof taylordefn u,
                              const),
                         taylornumbers u,
                         taylorneglist taylorlist u)
       end
     else nil;
 put('minus,'tayloropt,'tayloroptminus);


 symbolic procedure taylorneglist u;
 mapcar(u,function (lambda v;
                    (car v).(negsq cdr v)));



 symbolic procedure tayloroptplus args;
 begin
   scalar ret,hard,u;
   u:=args;
   while u do <<
     if 'return0 eq taylorfunction car u
       then ret:=(car u).ret
       else hard:=(car u).hard;
     u:=cdr u >>;
   if null ret or
       null cdr ret
     then return nil;
   ret:=mapply(function joinret,ret);
   if null hard
     then return ret;
   rplaca(args,ret);
   rplacd(args,hard);
    return nil
   end;
 put('plus,'tayloropt,'tayloroptplus);


 symbolic procedure joinret(u,v);
 begin
   scalar nums,a,b,al;
   nums:=(min2(taylorfirst u,taylorfirst v).
          max2(taylorlast u,taylorlast v));
   al:=nil;
   u:=taylorlist u;
   v:=taylorlist v;
   for i:=(car nums) step 1 until (cdr nums) do <<
     a:=assoc(i,u);
     b:=assoc(i,v);
     if a
       then if b
         then al:=(i.addsq(cdr a,cdr b)).al
         else al:=a.al
       else if b
         then al:=b.al  >>;
   return taylormake(list 'return0,nums,al)
   end;




 % the operator constanttimes
 % has two arguments (actually a list)
 % 1) a form dependent on the taylorvariable
 % 2) a form which is not.


 % the operator binarytimes has two arguments (actually a list)
   % but behaves like times otherwise.


 symbolic procedure taylorctimes(n,texpr,args);
 !*multsq(taylorevaluate(car args,n-(taylorfirst cadr args)),
        taylorevaluate(cadr args,taylorfirst cadr args));


 symbolic procedure initialtaylortimes slist;
 % Multiply the variable by the constant.
 ((caar slist)+(caadr slist)). !*multsq(cdar slist,cdadr slist);


 symbolic procedure tayloroptctimes u;
 if 'taylorctimes eq taylorfunction car u
   then begin
     scalar reala,const,iconst,degg;
     % we have nested multiplication.
     reala:=argof taylordefn car u;
     % the thing to be multiplied by the two constants.
     const:=car taylorlist cadr u;
     %the actual outer constant: deg.sq.
     iconst:=caddr taylordefn car u;
     %the inner constant.
     degg:=(taylorfirst iconst)+(car const);
     iconst:=list(taylordefn iconst,
                   degg.degg,
                   degg.!*multsq(cdar taylorlist iconst,cdr const));
     return list('taylorctimes,reala,iconst).
                 ((((taylorfirst car u) + (car const)).
                         ((taylorlast car u) + (car const))).
                  mapcar(taylorlist car u,function multconst))
     end
   else if 'return0 eq taylorfunction car u
     then begin
       scalar const;
       const:=car taylorlist cadr u;
       % the actual constant:deg.sq.
       u:=car u;
       return (taylordefn u).
                   ((((taylorfirst u)+car const).
                         ((taylorlast u)+car const)).
                   mapcar(taylorlist u,function multconst))
       end
     else nil;


 symbolic procedure multconst v;
 % Multiplies v by const in deg.sq form.
 ((car v)+(car const)) . !*multsq(cdr v,cdr const);


 put('constanttimes,'tayloropt,'tayloroptctimes);
 put('constanttimes,'simpfn,'simptimes);
 put('constanttimes,'taylorfunction,'taylorctimes);
 put('constanttimes,'initialtaylorfunction,'initialtaylortimes);


 symbolic procedure taylorbtimes(n,texpr,args);
 begin
   scalar answer,i,n1,n2;
   answer:= nil ./ 1;
   n1:=car firstterm car args;
   % the first term in one argument.
   n2:=car firstterm cadr args;
   % the first term in the other.
   for i:=n1 step 1 until (n-n2) do
     answer:=addsq(answer,!*multsq(taylorevaluate(cadr args,n-i),
                                       taylorevaluate(car args,i)));
   return answer
   end;




 put('binarytimes,'taylorfunction,'taylorbtimes);
 put('binarytimes,'initialtaylorfunction,'initialtaylortimes);
 put('binarytimes,'simpfn,'simptimes);


symbolic procedure taylorformertimes arglist;
begin
  scalar const,var,degg,wsqrt,negcount,u;
  negcount:=0;
  degg:=0;% the deggrees of any solitary x we may meet.
  const:=nil;
  var:=nil;
  wsqrt:=nil;
  while arglist do <<
    if dependsp(car arglist,taylorvariable)
      then if and(eqcar(car arglist,'expt),
                        cadar arglist eq taylorvariable,
                        numberp caddar arglist)
        then degg:=degg+caddar arglist
% removed JHD 21.8.86 - while it is anoptimisation,
% it runs the risk of proving that -1 = +1 by ignoring the
% number of "i" needed - despite the attempts we went to.
%        else if eqcar(car arglist,'sqrt)
%          then <<
%            u:=argof car arglist;
%            wsqrt:=u.wsqrt;
%            if minusq cdr firstterm taylorformp u
%              then negcount:=1+negcount >>
          else if car arglist eq taylorvariable
            then degg:=degg + 1
            else var:=(car arglist).var
      else const:=(car arglist).const;
    arglist:=cdr arglist >>;
  if wsqrt
    then if cdr wsqrt
      then var:=list('sqrt,prepsq simptimes wsqrt).var
      else var:=('sqrt.wsqrt).var;
  if var
    then var:=mapply(function (lambda u,v;
                               list('binarytimes,u,v)),var);
  % insert binary multiplications.
  negcount:=negcount/2;
  if onep cdr divide(negcount,2)
    then const:= (-1).const;
  % we had an odd number of (-1) from i*i.
  if const or (degg neq 0)
    then <<
      if const
        then const:=simptimes const
        else const:=1 ./ 1;
      const:=taylormake(list 'return0,degg.degg,list(degg.const));
      if null var
        then var:=const
        else var:=list('constanttimes,var,const) >>;
  return taylorformp var
  end;

put('times,'taylorformer,'taylorformertimes);




flag('(taylorbtimes taylorctimes taylorquotient),'taylor);
symbolic procedure taylorformerexpt arglist;
begin
  scalar base,expon;
  base:=car arglist;
  expon:=simpcar cdr arglist;
  if (denr expon neq 1) or
     (not numberp numr expon)
    then interr "Hard exponent";
  expon:=numr expon;
  if base neq taylorvariable
    then interr "Hard base";
  return list('return0 . nil,expon.expon,expon.(1 ./ 1))
  end;
put ('expt,'taylorformer,'taylorformerexpt);


symbolic procedure initialtaylorquotient slist;
(caar slist - caadr slist).!*multsq(cdar slist,!*invsq cdadr slist);


symbolic procedure taylorquotient(n,texpr,args);
begin
  % problem is texpr=b/c or c*texpr=b.
  scalar sofar,b,c,cfirst;
  b:=car args;
  c:=cadr args;
  cfirst:=taylorfirst c;
  sofar:=taylorevaluate(b,n+cfirst);
  for i:=taylorfirst texpr step 1 until n-1 do
    sofar:=addsq(sofar,!*multsq(taylorevaluate(texpr,i),
                              negsq taylorevaluate(c,n+cfirst-i)));
  return !*multsq(sofar,!*invsq taylorevaluate(c,cfirst))
  end;


put('quotient,'taylorfunction,'taylorquotient);
put('quotient,'initialtaylorfunction,'initialtaylorquotient);


symbolic procedure minusq sq;
if null sq
  then nil
  else if minusf numr sq
    then not minusf denr sq
     else minusf denr sq;



% This is wrapped round TAYLORFORMERSQRT2 in order to
% remove the innards of the SQRT from the asslist.
% note the precautions for nested SQRTs.

symbolic procedure taylorformersqrt arglist;
% ***store-hack-1***;
% Uncomment these lines if more store is available.
%begin
%  scalar z;
%  z:=taylorasslist;
%  if sqrtsintree(car arglist,taylorvariable)
%    then return taylorformersqrt2 arglist;
%  arglist:=taylorformersqrt2 arglist;
%  taylorasslist:=z;
%  return arglist
%  end;
%
%
%symbolic procedure taylorformersqrt2 arglist;
begin
  scalar f,realargs,ff,realsqrt;
  realargs:=taylorformp carx(arglist,'taylorformersqrt2);
  f:=firstterm realargs;
  if not evenp car f
    then interr "Extra sqrt substitution needed";
  if and(0 iequal car f,
         1 iequal numr cdr f,
         1 iequal denr cdr f)
    then return taylorformp list('sqrtx,realargs);
  % if it starts with 1 already then it is easy.
  ff:=- car f;
  ff:=list(list 'return0,ff.ff,ff.(!*invsq cdr f));
  % ff is the leading term in the expansion of realargs.
  realsqrt:=list('sqrtx,list('constanttimes,realargs,ff));
  ff:=(car f)/2;
  return taylorformp list('constanttimes,
                          realsqrt,
                          list(list 'return0,
                               ff.ff,
                               ff.(simpsqrtsq cdr f)))
  end;


put('sqrt,'taylorformer,'taylorformersqrt);


symbolic procedure initialtaylorsqrtx slist;
0 . (1 ./ 1);
% sqrt(1+ ...) = 1+....


symbolic procedure taylorsqrtx(n,texpr,args);
begin
  scalar sofar,i;
  sofar:=taylorevaluate(car args,n);
  % (1+.....+a(n)*x**n)**2
  % = ....+x**n*(2*a(n)+sum(0<i<n,a(i)*a(n-i))).
  % So a(n)=(coeff(x**n)-sum) /2.
  for i:=1 step 1 until (n-1) do
    sofar:=addsq(sofar,negsq !*multsq(taylorevaluate(texpr,i),
                                    taylorevaluate(texpr,n-i)));
  return multsq(sofar,1 ./ 2)
  end;


flag('(taylorsqrtx),'taylor);
put('sqrtx,'taylorfunction,'taylorsqrtx);
put('sqrtx,'initialtaylorfunction,'initialtaylorsqrtx);

endmodule;


module torsionb;

% Author: James H. Davenport.

fluid '(intvar nestedsqrts);

global '(!*tra !*trmin);

exports bound!-torsion;

symbolic procedure bound!-torsion(divisor,dof1k);
% Version 1 (see Trinity Thesis for difference).
begin
  scalar field,prime1,prime2,prime3,minimum,places;
  scalar non!-p1,non!-p2,non!-p3,curve,curve2,nestedsqrts;
  places:=for each u in divisor
            collect car u;
  curve:=getsqrtsfromplaces places;
  if nestedsqrts
    then rederr "Not yet implemented"
    else curve2:=curve;
  for each u in places do begin
    u:=rfirstsubs u;
    if eqcar(u,'quotient) and cadr u = 1
      then return;
    u:=substitutesq(simp u,list(intvar . 0));
    field:=union(field,sqrtsinsq(u,nil));
    u:=list(intvar . prepsq u);
    for each v in curve2 do
      field:=union(field,sqrtsinsq(substitutesq(v,u),nil));
    end;
  prime1:=2;
  while null (non!-p1:=good!-reduction(curve,dof1k,field,prime1)) do
    prime1:=nextprime prime1;
  prime2:=nextprime prime1;
  while null (non!-p2:=good!-reduction(curve,dof1k,field,prime2)) do
    prime2:=nextprime prime2;
  prime3:=nextprime prime2;
  while null (non!-p3:=good!-reduction(curve,dof1k,field,prime3)) do
    prime3:=nextprime prime3;
  minimum:=fix sqrt float(non!-p1*non!-p2*non!-p3);
  minimum:=min(minimum,non!-p1*max!-power(prime1,min(non!-p2,non!-p3)));
  minimum:=min(minimum,non!-p2*max!-power(prime2,min(non!-p1,non!-p3)));
  minimum:=min(minimum,non!-p3*max!-power(prime3,min(non!-p2,non!-p1)));
  if !*tra or !*trmin then <<
    princ "Torsion is bounded by ";
    printc minimum >>;
  return minimum
  end;


symbolic procedure max!-power(p,n);
% Greatest power of p not greater than n.
begin scalar ans;
  ans:=1;
  while ans<=n do
    ans:=ans*p;
  ans:=ans/p;
  end;


symbolic procedure good!-reduction(curve,dof1k,field,prime);
begin
  scalar u;
  u:=algebraic!-factorise(prime,field);
  interr "Good reduction not finished";
  end;

endmodule;


module wstrass;

% Author: James H. Davenport.

fluid '(!*backtrace
        intvar
        listofallsqrts
        listofnewsqrts
        magiclist
        previousbasis
        sqrt!-intvar
        sqrtflag
        sqrts!-in!-integrand
        taylorasslist
        taylorvariable
        thisplace
        zlist);

global '(!*tra !*trmin coates!-fdi);

exports simpwstrass,weierstrass!-form,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;
  tt:=readclock();
  sqrtflag:=t;
  taylorvariable:=intvar:=car u;
  sqrt!-intvar:=mvar !*q2f simpsqrti intvar;
  u:=for each v in cdr u
            collect simp!* v;
  sqrts!-in!-integrand:=sqrtsinsql(u,intvar);
  u:=errorset('(weierstrass!-form sqrts!-in!-integrand),
              t,!*backtrace);
  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,vec,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 >>;
  vec:=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(vec,i,!*q2f !*multsq(u,getv(vec,i)));
  if !*tra then <<
    printc "List of seven functions in weierstrass!-form:";
    mapvec(vec,function printsf) >>;
  vec:=wstrass!-lineq vec;
% printsq(addsq(getv(vec,0),addsq(!*multsq(getv(vec,1),x1),
%               addsq(!*multsq(getv(vec,2),x2),
%                     addsq(!*multsq(getv(vec,3),!*multsq(x1,x1)),
%                           addsq(!*multsq(getv(vec,4),!*multsq(x2,x2)),
%                              addsq(!*multsq(getv(vec,5),exptsq(x1,3)),
%                                       !*multsq(getv(vec,6),
%                                               !*multsq(x1,x2)))))))));
  x2:=!*addsq(!*multsq(!*multsq(2 ./ 1,getv(vec,4)),x2),
              addsq(!*multsq(x1,getv(vec,6)),
                    getv(vec,2)));
  putv(vec,4,!*multsq(-4 ./ 1,getv(vec,4)));
  a:=!*multsq(getv(vec,4),getv(vec,5));
  b:=!*addsq(!*multsq(getv(vec,6),getv(vec,6)),
             !*multsq(getv(vec,3),getv(vec,4)));
  c:=!*addsq(!*multsq(2 ./ 1,!*multsq(getv(vec,2),getv(vec,6))),
             !*multsq(getv(vec,1),getv(vec,4)));
  d:=!*addsq(!*multsq(getv(vec,2),getv(vec,2)),
             !*multsq(getv(vec,0),getv(vec,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 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 zerop cdr (uu:=divide(u,b)) do <<
      nc:=nc+1;
      u:=car uu >>;
    if nc < 2
      then go to next;
    uu:=algint!-numeric!-content numr d;
    nd:=0;
    while zerop cdr (uu:=divide(uu,b)) 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 mapcar(places,function basicplace);
%  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();
  places2:=mapcar(places2,function fdi!-upgrade);
  save:=taylorasslist;
  u:=wstrassmodule(places2,
                  mapcar(mults2,function (lambda u;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 vec;
begin
  scalar zlist,powlist,m,rightside,v;
  scalar zero,one;
  zero:=nil ./ 1;
  one:=1 ./ 1;
  for i:=0:6 do
    zlist:=varsinsf(getv(vec,i),zlist);
  zlist:=intvar . findzvars(zlist,nil,intvar,nil);
  for i:=0:6 do
    putv(vec,i,f2df getv(vec,i));
  for i:=0:6 do
    for each u in getv(vec,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(vec,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;


module zmodule;

% Author: James H. Davenport.

fluid '(!*galois
        basic!-listofallsqrts
        basic!-listofnewsqrts
        commonden
        gaussiani
        listofallsqrts
        listofnewsqrts
        sqrt!-places!-alist
        taylorasslist);

global '(!*tra !*trfield !*trmin);

exports zmodule;
imports !*multf,sqrtsinsql,sortsqrts,simp,!*q2f,actualsimpsqrt,printsf;
imports prepf,substitutesq,printsq,mapply,!*multsq,mkilist;
imports mkvecf2q,mkvec,mkidenm,invsq,multsq,negsq,addsq,gcdn;
imports !*invsq,prepsq;

symbolic procedure zmodule(w);
begin
  scalar reslist,denlist,u,commonden,basis,p1,p2,hcf;
  % w is a list of elements (place.residue)=sq.
  for each v in w do <<
    u:=cdr v;
    reslist:=u.reslist;
    denlist:=(denr u).denlist >>;
  basis:=sqrtsinsql(reslist,nil);
  if null u or null cdr u or !*galois
    then go to nochange;
  reslist:=check!-sqrts!-dependence(reslist,basis);
  denlist:=for each u in reslist
             collect denr u;
nochange:
 commonden:=mapply(function(lambda u,v;
                      multf(u,quotf(v,gcdf(u,v)))),denlist)./1;
  u:=nil;
  for each v in reslist do
    u:=(numr !*multsq(v,commonden)).u;
  reslist:=u;
    % We have effectively reserves RESLIST twice,
    % so it is in the corect order.
  u:=bexprn(reslist);
  basis:=car u;
  reslist:=cdr u;
  denlist:=nil;
  while basis do <<
    p1:=reslist;
    p2:=w;
    u:=nil;
    hcf:=0;
    while p1 do <<
      if 0 neq caar p1
        then  <<
          u:=((caar p2).(caar p1)).u;
          hcf:=gcdn(hcf,caar p1) >>;
      p1:=cdr p1;
      p2:=cdr p2 >>;
    if hcf neq 1
     then u:=for each uu in u collect
               (car uu). ( (cdr uu) / hcf);
    denlist:=(prepsq !*multsq(car basis,
                              multsq(!*f2q hcf,!*invsq commonden))
                  .u).denlist;
    basis:=cdr basis;
    reslist:=mapcar(reslist,function cdr) >>;
  return denlist
  end;


symbolic procedure bexprn(wlist);
begin
  scalar basis,replist,w,w2,w3,p1,p2;
  % wlist is a list of sf.
  w:=reverse wlist;
  replist:=nil;
  while w do <<
    w2:=sf2df car w;
    % now ensure that all elements of w2 are in the basis.
    w3:=w2;
    while w3 do <<
      % caar is the sf,cdar a its coefficient.
      if not member(caar w3,basis)
        then <<
          basis:=(caar w3).basis;
          replist:=mapcons(replist,0) >>;
          % adds car w3 to basis.
      w3:=cdr w3 >>;
    replist:=mkilist(basis,0).replist;
    % builds a new zero representation.
    w3:=w2;
    while w3 do <<
      p1:=basis;
      p2:=car replist;
      %the list for this element.
      while p1 do <<
        if caar w3 = car p1
          then rplaca(p2,cdar w3);
        p1:=cdr p1;
        p2:=cdr p2 >>;
      w3:=cdr w3 >>;
    w:=cdr w >>;
  return mkbasis(basis,replist)
  end;


symbolic procedure mkbasis(basis,reslist);
begin
  scalar row,nbasis,nreslist,u,v;
  basis:=for each u in basis collect !*f2q u;
  % basis is a list of sq's
  % reslist is a list of representations in the form
  % ( (coeff1 coeff2 ...)    ...).
  nreslist:=mkilist(reslist,nil);
  % initialise our list-of-lists.
  trynewloop:
  row:=mapcar(reslist,function car);
  reslist:=mapcar(reslist,function cdr);
  if obvindep(row,nreslist)
    then u:=nil
    else u:=lindep(row,nreslist);
  if u
    then <<
      % u contains the numbers with which to add this new item into the
      % basis.
      v:=nil;
      while nbasis do <<
        v:=addsq(car nbasis,!*multsq(car basis,car u)).v;
        nbasis:=cdr nbasis;
        u:=cdr u >>;
      nbasis:=reversewoc v >>
    else <<
      nreslist:=pair(row,nreslist);
      nbasis:=(car basis).nbasis
      >>;
  basis:=cdr basis;
  if basis
   then go to trynewloop;
  return nbasis.nreslist
  end;


symbolic procedure obvindep(row,matrx);
  % True if row is obviously linearly independent of the
  % Rows of the matrix.
begin scalar u;
  if null car matrx
    then return t;
  % no matrix => no dependence.
nexttry:
  if null row
    then return nil;
  if 0 iequal car row
    then go to nouse;
  u:=car matrx;
testloop:
  if 0 neq car u
    then go to nouse;
  u:=cdr u;
  if u
    then go to testloop;
  return t;
nouse:
  row:=cdr row;
  matrx:=cdr matrx;
  go to nexttry
  end;


symbolic procedure sf2df sf;
if null sf
   then nil
   else if numberp sf
    then (1 . sf).nil
    else begin
      scalar a,b,c;
      a:=sf2df lc sf;
      b:=(lpow sf .* 1) .+ nil;
      while a do <<
        c:=(!*multf(caar a,b).(cdar a)).c;
        a :=cdr a >>;
      return nconc(c,sf2df red sf)
      end;





symbolic procedure check!-sqrts!-dependence(sql,sqrtl);
% Resimplifies the list of SQs SQL,
% allowing for all dependencies among the
% sqrts in SQRTl.
begin
  scalar !*galois,sublist,sqrtsavelist,changeflag;
  sqrtsavelist:=listofallsqrts.listofnewsqrts;
  listofnewsqrts:=list mvar gaussiani;
  listofallsqrts:=list((argof mvar gaussiani) . gaussiani);
  !*galois:=t;
  for each u in sortsqrts(sqrtl,nil) do begin
    scalar v,uu;
    uu:=!*q2f simp argof u;
    v:=actualsimpsqrt uu;
    listofallsqrts:=(uu.v).listofallsqrts;
    if domainp v or mvar v neq u
      then <<
        if !*tra or !*trfield then <<
           printc u;
           printc "re-expressed as";
           printsf v >>;
        v:=prepf v;
        sublist:=(u.v) . sublist;
        changeflag:=t >>
    end;
  if changeflag then <<
    sql:=for each u in sql collect
           substitutesq(u,sublist);
    taylorasslist:=nil;
    sqrt!-places!-alist:=nil;
    basic!-listofallsqrts:=listofallsqrts;
    basic!-listofnewsqrts:=listofnewsqrts;
    if !*tra or !*trmin then <<
      printc "New set of residues are";
      mapc(sql,function printsq) >> >>
    else <<
      listofallsqrts:=car sqrtsavelist;
      listofnewsqrts:=cdr sqrtsavelist >>;
  return sql
  end;



symbolic procedure lindep(row,matrx);
  begin
    scalar m,mm,n,i,j,k,u,v,inverse,rowsinuse,failure;
    % Inverse is the answer from the "gaussian elimination"
    % we are doing.
    % Rowsinuse has nil for rows with no "awkward" non-zero entries.
    mm:=length car matrx;
    m:=isub1 mm;
    n:=isub1 length matrx;
    % n=length row.
    row:=mkvecf2q row;
    matrx:=mkvec mapcar(matrx,function mkvecf2q);
    inverse:=mkidenm mm;
    rowsinuse:=mkvect m;
    failure:=t;
    % initialisation complete.
    for i:=0 step 1 until n do begin
      % try to kill off i'th elements in each row.
      u:=nil;
      for j:=0 step 1 until m do <<
        % try to find a  pivot element.
        if  (null u) and
            (null getv(rowsinuse,j)) and
            (numr getv(getv(matrx,i),j))
          then u:=j >>;
      if null u
        then go to nullu;
      putv(rowsinuse,u,t);
      % it is no use trying this again ---
      % u is our pivot element.
      if u iequal m
        then go to nonetokill;
      for j:=iadd1 u step 1 until m do
        if numr getv(getv(matrx,i),j)
          then <<
            v:=negsq multsq(getv(getv(matrx,i),j),
                            invsq getv(getv(matrx,i),u));
            for k:=0 step 1 until mm do
              putv(getv(inverse,k),j,
                addsq(getv(getv(inverse,k),j),
                  multsq(v,getv(getv(inverse,k),u))));
            for k:=0 step 1 until n do
              putv(getv(matrx,k),j,
                addsq(getv(getv(matrx,k),j),
                  multsq(v,getv(getv(matrx,k),u)))) >>;
      %we have now pivoted throughout matrx.
    nonetokill:
      % now do the same in row if necessary.
      if null numr getv(row,i)
        then go to norowop;
      v:=negsq multsq(getv(row,i),
                 invsq getv(getv(matrx,i),u));
      for k:=0 step 1 until mm do
        putv(getv(inverse,k),mm,
          addsq(getv(getv(inverse,k),mm),
            multsq(v,getv(getv(inverse,k),u))));
      for k:=0 step 1 until n do
        putv(row,k,addsq(getv(row,k),
                     multsq(v,getv(getv(matrx,k),u))));
      u:=nil;
      for k:=0 step 1 until n do
        if numr getv(row,k)
          then u:=t;
      % if u is null then row is all 0.
      if null u
        then <<
          n:=-1;
          failure:=nil >>;
    norowop:
      if !*tra then <<
        princ "At end of cycle";
        printc row;
        printc matrx;
        printc inverse >>;
      return;
    nullu:
      % there is no pivot for this u.
      if numr getv(row,i)
        then n:=-1;
        % this element cannot be killed.
      end;
    if failure
      then return nil;
    v:=nil;
    for i:=0 step 1 until m do
      v:=(negsq getv(getv(inverse,m-i),mm)).v;
    return v
    end;

endmodule;


end;