reduce-historical/r38/packages/algint/afactor.red

module afactor;

% Author: James H. Davenport.

fluid '(!*galois !*noextend !*sqfree !*trfield afactorvar listofnewsqrts
	monicpart varlist zlist sqrtlist sqrtflag indexlist);

switch trfield;   % traces the algebraic number field manipluation

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

% internal!-fluid '(monicpart);

% This routine, and its subsidiaries, do factorization over algebraic
% extensions, by the Trager modification of van der Waerden's algorithm
% See SYMSAC '76.

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 "factorizes as ";
    mapc(v,function printsf) >>;
  return v
  end;


symbolic procedure algfactor2(f,a);
if null a
  then for each u in cdr factorf f collect %car factorf is a constant
    if cdr u = 1
      then car u
      else interr "repeated factors found while processing algebraics"
  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 w and 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 null ff then return list f         % does not factor.
   else 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 factorization 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 for each j in a collect demonise j;
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 for each j in z collect demonise j;
  % 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 gcdinonevar(u,v,x);
% Gcd of u and v, regarded as polynomials in x.
if null u
  then v
  else if null v
    then u
    else begin
      scalar u1,v1,z,w;
      u1:=stt(u,x);
      v1:=stt(v,x);
    loop:
      if (car u1) > (car v1)
        then go to ok;
      z:=u1;u1:=v1;v1:=z;
      z:=u;u:=v;v:=z;
    ok:
      if car v1 iequal 0
        then interr "Coprimeness in gcd";
      z:=gcdf(cdr u1,cdr v1);
      w:=quotf(cdr u1,z);
      if (car u1) neq (car v1)
        then w:=!*multf(w,!*p2f mksp(x,(car u1)-(car v1)));
      u:=addf(!*multf(v,w),
              !*multf(u,negf quotf(cdr v1,z)));
      if null u
        then return v;
      u1:=stt(u,x);
      go to loop
      end;

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;

symbolic procedure jfactor(sf,var);
   % Algebraic integrator's sole interface to factorization.
   % except for a direct call to the true factoriser from
   % inside afactor
begin
  scalar varlist,zlist,indexlist,sqrtlist,sqrtflag;
  scalar prim,res,answer,u,v,x,y; % scalar var2
  prim:=jsqfree(sf,var);
  indexlist:=createindices zlist;
  while not null prim do <<
    x:=car prim;
    while not null x do <<
        y:=facbypp(car x,varlist);
        while not null y do <<
            res:=append(int!-fac car y,res);
            y:=cdr y >>;
        x:=cdr x >>;
    prim:=cdr prim >>;
  while res do <<
    if caar res eq 'log then <<
        u:=cdar res;
        u:=!*multsq(numr u ./ 1,1 ./ cdr stt(numr u,var));
        v:=denr u;
        while not atom v do v:=lc v;
        if  (numberp v) and (0> v)
          then u:=(negf numr u) ./ (negf denr u);
        if u neq '(1 . 1) then answer := u . answer>>
      else if caar res eq 'atan then nil
      % We can ignore this, since we also get LOG (U**2+1) as another
      % term.
      else interr "Unexpected term in jfactor";
    res:=cdr res >>;
  return answer
  end;

% unfluid '(monicpart);

endmodule;

end;