reduce-historical/r38/packages/poly/gint.red

module gint;  % Support for gaussian integers (complex numbers).

% Author: Eberhard Schruefer.

global '(domainlist!*);

fluid '(!*complex !*gcd);

switch complex;

domainlist!* := union('(!:gi!:),domainlist!*);

symbolic procedure setcmpxmode(u,bool);
   % Sets polynomial domain mode in complex case.
   begin scalar x,y;
      x := get(u,'tag);
      if u eq 'complex
	then if null dmode!*
	       then return if null bool then nil
                else <<put('i,'idvalfn,'mkdgi);
                       setdmode1('complex,bool)>>
	 else if null bool
		then return if null !*complex then nil
                else if get(dmode!*,'dname) eq 'complex
                      then <<remprop('i,'idvalfn);
                             setdmode1('complex,nil)>>
                     else <<remprop('i,'idvalfn);
                           setdmode1(get(get(dmode!*,'realtype),'dname),
                                       t)>>
         else if dmode!* eq '!:gi!: or get(dmode!*,'realtype)
            then return nil
         else if not (y := get(dmode!*,'cmpxtype))
               then dmoderr(dmode!*,x)
         else <<put('i,'idvalfn,get(y,'ivalue));
                     return setdmode1(get(y,'dname),bool)>>
      else if null bool then
         if u eq (y := get(get(dmode!*,'realtype),'dname)) then
            <<put('i,'idvalfn,'mkdgi); return setdmode1('complex,t)>>
	    else if null y then return nil
            else offmoderr(u,y)
      else <<u := get(u,'tag);
             y := get(u,'cmpxtype);
             if null y then dmoderr(u,'!:gi!:);
             put('i,'idvalfn,get(y,'ivalue));
	     return setdmode1(get(y,'dname),bool)>>
   end;

% Used by gcdk.

symbolic procedure intgcdd(u,v);
   if null u then v
   else if atom u then
      if atom v then gcdn(u,v) else gcdn(cadr v,gcdn(cddr v,u))
   else if atom v then intgcdd(v,u)
      else intgcdd(cadr u,intgcdd(cddr u,v));

put('complex,'tag,'!:gi!:);
put('!:gi!:,'dname,'complex);
put('!:gi!:,'i2d,'!*i2gi);
put('!:gi!:,'minusp,'giminusp!:);
put('!:gi!:,'zerop,'gizerop!:);
put('!:gi!:,'onep,'gionep!:);
put('!:gi!:,'plus,'giplus!:);
put('!:gi!:,'difference,'gidifference!:);
put('!:gi!:,'times,'gitimes!:);
put('!:gi!:,'quotient,'giquotient!:);
put('!:gi!:,'divide,'gidivide!:);
put('!:gi!:,'gcd,'gigcd!:);
put('!:gi!:,'factorfn,'gifactor!:);
% put('!:gi!:,'rationalizefn,'girationalize!:);
put('!:gi!:,'prepfn,'giprep!:);
put('!:gi!:,'intequivfn,'gintequiv!:);
put('!:gi!:,'specprn,'giprn!:);
put('!:gi!:,'prifn,'giprn!:);
put('!:gi!:,'cmpxfn,'mkgi);

put('!:gi!:,'unitsfn,'!:gi!:unitconv);

symbolic procedure !:gi!:unitconv(u,v);
   unitconv(u,v,get('!:gi!:,'units));

put('!:gi!:,'units,'(((!:gi!: 0 . 1) . (!:gi!: 0 . -1))
	     ((!:gi!: 0 . -1) . (!:gi!: 0 . 1))));

symbolic procedure unitconv(u,v,w);
   begin scalar z;
   a: if null w then return u;
      z := quotf1(v,caar w);
      if null z or not fixp z then <<w := cdr w; go to a>>;
      z := multf(denr u,cdar w);
      w := multf(numr u,cdar w);
      if minusf z then <<w := negf w; z := negf z>>;
      return w ./ z
   end;

symbolic procedure !*i2gi u; '!:gi!: . (u . 0);

symbolic procedure giminusp!: u;
   %*** this is rather a test for u being in a canonical form! ***;
   if cadr u = 0 then minusp cddr u else minusp cadr u;

symbolic procedure gizerop!: u;
   cadr u = 0 and cddr u = 0;

symbolic procedure gionep!: u;
   cadr u=1 and cddr u=0;

symbolic procedure gintequiv!: u;
   if cddr u=0 then cadr u else nil;

symbolic procedure mkdgi u;
   ('!:gi!: . (0 . 1)) ./ 1;

symbolic procedure mkgi(re,im);
   '!:gi!: . (re . im);

symbolic procedure giplus!:(u,v);
   mkgi(cadr u+cadr v,cddr u+cddr v);

symbolic procedure gidifference!:(u,v);
   mkgi(cadr u-cadr v,cddr u-cddr v);

symbolic procedure gitimes!:(u,v);
   (lambda r1,i1,r2,i2;
       mkgi(r1*r2-i1*i2,r1*i2+r2*i1))
    (cadr u,cddr u,cadr v,cddr v);

symbolic procedure giquotient!:(u,v);
   begin integer r1,i1,r2,i2,d; scalar rr,ii;
     r1 := cadr u; i1 := cddr u;
     r2 := cadr v; i2 := cddr v;
     d := r2*r2+i2*i2;
     rr := divide(r1*r2+i1*i2,d);
     ii := divide(i1*r2-i2*r1,d);
     return if cdr ii=0 and cdr rr=0 then mkgi(car rr,car ii)
             else '!:gi!: . (0 . 0)
   end;

symbolic procedure gidivide!:(u,v);
   begin integer r1,i1,r2,i2,d,rr,ir,rq,iq;
     r1 := cadr u; i1 := cddr u;
     r2 := cadr v; i2 := cddr v;
     d := r2*r2+i2*i2;
     rq := r1*r2+i1*i2;
     iq := i1*r2-i2*r1;
     rq := car divide(2*rq+if rq<0 then -d else d,2*d);
     iq := car divide(2*iq+if iq<0 then -d else d,2*d);
     rr := r1-(rq*r2-iq*i2);
     ir := i1-(iq*r2+rq*i2);
     return mkgi(rq,iq) . mkgi(rr,ir)
   end;

symbolic procedure giremainder(u,v);
   begin integer r1,i1,r2,i2,d,rr,ir,rq,iq;
     r1 := cadr u; i1 := cddr u;
     r2 := cadr v; i2 := cddr v;
     d := r2*r2+i2*i2;
     rq := r1*r2+i1*i2;
     iq := i1*r2-i2*r1;
     rq := car divide(2*rq+if rq<0 then -d else d,2*d);
     iq := car divide(2*iq+if iq<0 then -d else d,2*d);
     rr := r1-(rq*r2-iq*i2);
     ir := i1-(iq*r2+rq*i2);
     return '!:gi!: . (rr . ir)
   end;

symbolic procedure gigcd!:(u,v);
   % Straightforward Euclidean algorithm.
   if gizerop!: v then fqa u else gigcd!:(v,giremainder(u,v));

symbolic procedure fqa u;
   %calculates the unique first-quadrant associate of u;
  if cddr u=0 then abs cadr u
   else if cadr u=0 then '!:gi!: . (abs cddr u . 0)
   else if (cadr u*cddr u)>0 then
		  '!:gi!: . (abs cadr u . abs cddr u)
   else '!:gi!: . (abs cddr u . abs cadr u);

symbolic procedure gifactor!: u;
   % Trager's modified version of Van der Waerdens algorithm.
   begin scalar x,y,norm,aftrs,ftrs,mvu,dmode!*,!*exp,w,z,l,bool;
         integer s;
     !*exp := t;
     if realp u
        then u := cdr factorf u
      else u := list(u . 1);
     w := 1;
     for each f in u do begin
     u := car f;
     dmode!* := '!:gi!:;
     mvu := power!-sort powers u;
     bool := contains!-oddpower mvu;
     if realp u and bool
        then <<u := normalize!-lcgi u;
               w := multd(car u,w);
               aftrs := (cdr u . 1) . aftrs;
               return>>;
     mvu := caar mvu;
     y := u;
     go to b;
     a: l := list(mvu . prepf addf(!*k2f mvu,multd(s,!*k2f 'i)));
        u := numr subf1(y,l);
     b: if realp u
           then if bool
                   then <<y := normalize!-lcgi y;
                          w := multd(car y,w);
                          aftrs := (cdr y . 1) . aftrs;
                          return>>
                 else <<s := s-1; go to a>>;
	norm := multf(u,conjgd u);
        if not sqfrp norm then <<s := s-1; go to a>>;
     dmode!* := nil;
     ftrs := factorf norm;
     dmode!* := '!:gi!:;
     if null cddr ftrs then <<y := normalize!-lcgi y;
                              w := multd(car y,w);
                              aftrs := (cdr y . 1) . aftrs;
                              return>>;
     w := car ftrs;
     l := if s=0 then nil
           else list(mvu . prepf addf(!*k2f mvu,
                                      negf multd(s,!*k2f 'i)));
     for each j in cdr ftrs do
         <<x := gcdf!*(car j,u);
           u := quotf!*(u,x);
           z := if l then numr subf1(x,l) else x;
           z := normalize!-lcgi z;
           w := multd(car z,w);
           aftrs := (cdr z . 1) . aftrs>>;
     w := multf(u,w) end;
     return w . aftrs
   end;

symbolic procedure normalize!-lcgi u;
   % Normalize lnc by using units as canonsq would do it.
   begin scalar l,x,y;
     l := lnc u;
     if numberp l then return if minusp l then (-1) . negf u
                               else 1 . u;
     x := get('!:gi!:,'units);
     a: if null x then return 1 . u;
        y := quotf1(l,caar x);
        if null y or null fixp y then <<x := cdr x; go to a>>;
        u := multd(cdar x,u);
        return if minusf u then negf caar x . negf u
                else caar x . u
   end;

symbolic procedure contains!-oddpower u;
   if null u then nil
    else if evenp cdar u then contains!-oddpower cdr u
    else t;

symbolic procedure power!-sort u;
   begin scalar x,y;
     while u do
       <<x := maxdeg(cdr u,car u);
	 u := delallasc(car x,u);
         y := x . y>>;
     return y
   end;

symbolic procedure sqfrp u;
   % Square free test for poly u over the integers.
   % It works best with ezgcd on.
   begin scalar !*ezgcd, dmode!*;
     % Make sure ezgcd loaded.
     if null getd 'ezgcdf1 then load_package ezgcd;
     !*ezgcd := t;
     return domainp gcdf!*(u,diff(u,mvar u))
   end;

symbolic procedure realp u;
   if domainp u
     then atom u
        or not get(car u,'cmpxfn)
        or cddr u = cddr apply1(get(car u,'i2d),1)
    else realp lc u and realp red u;

symbolic procedure fd2f u;
   if atom u then u
    else if car u eq '!:gi!:
        then addf(!*n2f cadr u,multf(!*k2f 'i,!*n2f cddr u))
    else addf(multf(!*p2f lpow u,fd2f lc u),fd2f red u);


% symbolic procedure giprep!: u;  %giprep1 cdr u;
%   prepsq!* addsq(!*n2f cadr u ./ 1,
%                  multsq(!*n2f cddr u ./ 1, !*k2q 'i));

% symbolic procedure giprep1 u;  %not used now;
%    if cdr u=0 then car u
%     else if car u=0 then retimes list(cdr u,'i)
%     else begin scalar gn;
%            gn := gcdn(car u,cdr u);
%            return retimes list(gn,
%                       replus list(car u/gn,retimes list(cdr u/gn,'i)))
%          end;

symbolic procedure giprep!: u;
   <<if im=0 then rl else if rl=0 then giprim im
     else if im<0 then list('difference,rl,giprim(minus im))
     else list('plus,rl,giprim im)>>
   where rl=cadr u,im=cddr u;

symbolic procedure giprim im; if im=1 then 'i else list('times,im,'i);

symbolic procedure giprn!: v;
   (lambda u;
    if atom u or (car u eq 'times) then maprin u
     else <<prin2!* "("; maprin u; prin2!* ")" >>) giprep!: v;

% symbolic procedure girationalize!: u;
%    %Rationalizes standard quotient u over the gaussian integers.
%    begin scalar x,y,z;
%       y := denr u;
%       z := conjgd y;
%       if y=z then return u;
%       x := multf(numr u,z);
%       y := multf(y,z);
%       return x ./ y
%    end;

symbolic procedure girationalize!: u;
   % Rationalizes standard quotient u over the gaussian integers.
   begin scalar y,z,!*gcd;
      !*gcd := t;
      if y=(z := conjgd(y := denr u)) then return u;
      % Remove from z any real polynomial factors of y and z.
      z := quotf(z,quotf(
	 gcdf(addf(y,z),multf(addf(z,negf y),'!:gi!: . (0 . 1))),2));
      % The following subs2 can undo the above if !*match is non-NIL.
%     return subs2 gigcdsq(multf(numr u,z),multf(y,z))
      return gigcdsq(multf(numr u,z),multf(y,z))
   end;

symbolic procedure gigcdsq(x,y); % remove integer common factor.
   <<if x then
        (<<if d neq 1 and (d := giintgcd(x,d)) neq 1 then
              <<x := quotf(x,d); y := quotf(y,d)>> >>
         where d=giintgcd(y,0)); x ./ y >>;

symbolic procedure giintgcd(u,d);
   if d=1 then 1 else if null u then d else if atom u then gcdn(u,d)
   else if eqcar(u,'!:gi!:) then gcdn(cadr u,gcdn(cddr u,d))
   else giintgcd(lc u,giintgcd(red u,d));

symbolic procedure conjgd u;
   begin scalar x;
      return if atom u then u
	      else if domainp u and (x := get(car u,'cmpxfn))
	       then
		  apply2(x,cadr u,
		   if numberp cddr u then -cddr u
		       % Allow for tagged parts of complex object.
		       else if domainp cddr u and not numberp caddr u
			  then !:minus cddr u
		       else cdr !:minus (get(car u,'realtype).cddr u))
	      else if domainp u then u   % Should be a real number now.
	      else addf(multpf(lpow u,conjgd lc u),conjgd red u)
   end;

initdmode 'complex;

endmodule;

end;