Artifact c11d6aeccb7b50378242703ea1c89e738da0c6add4ea03a757dafee2da3ce402:

• Executable file r37/packages/arnum/arinv.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5399) [annotate] [blame] [check-ins using] [more...]

• Executable file r38/packages/arnum/arinv.red — part of check-in [f2fda60abd] at 2011-09-02 18:13:33 on branch master — Some historical releases purely for archival purposes

  git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 5399) [annotate] [blame] [check-ins using]

module arinv;   % Routines for inversion of algebraic numbers.

% Author: Eberhard Schruefer.

fluid '(dmode!*);

global '(arbase!* curdefpol!*);

symbolic procedure arquot(u,v);
   % U is an ar domain element, result is a matrix form which
   % needs to be inverted for calculating the inverse of ar.
   begin scalar mv,r,sgn,x,y,z,w,dmode!*;
         integer n;
     mv := mvar curdefpol!*;
     x := u;
     w := for each k in cdr arbase!* collect
             (x := reducepowers multf(!*k2f mv,x));
     x := nil;
     w := negf v . (u . w);
     for j := (ldeg curdefpol!* - 1) step -1 until 0 do
       <<y := nil;
         z := 1;
         n := -2;
         w := for each k in w collect
                if (degr(k,mv) = j) and k then
                   <<y := list(n := n+1)
                          . * (r := if domainp k then k else lc k)
                          .+ y;
                     if eqcar(r,'!:rn!:) then z := ilcm(z,cddr r);
                     if null domainp k then red k>>
                 else <<n := n + 1; k>>;
         if z neq 1 then
            y := for each j on y collect
                     lpow j .* if eqcar(lc j,'!:rn!:) then
                                  cadr lc j*z/cddr lc j
                                else lc j*z;
         if null x then x := y
          else x := b!:extmult(y,x)>>;
     sgn := evenp length lpow x;
     z := nil;
     for each j in lpow x do
         z := addf(if j = 0 then arnum!-mkglsol(0,x,sgn := not sgn,-1)
                    else multpf(mv to j,
				arnum!-mkglsol(j,x,sgn:=not sgn,-1)),z);
     return z
   end;

symbolic procedure arquot1 u;
   % U is an ar domain element, result is a matrix form which
   % needs to be inverted for calculating the inverse of ar.
   begin scalar mv,r,sgn,x,y,z,w,dmode!*;
         integer n;
     mv := mvar curdefpol!*;
     x := u;
     w := for each k in cdr arbase!* collect
             (x := reducepowers multf(!*k2f mv,x));
     x := nil;
     w := -1 . (u . w);
     for j := (ldeg curdefpol!* - 1) step -1 until 0 do
       <<y := nil;
         z := 1;
         n := -2;
         w := for each k in w collect
                if (degr(k,mv) = j) and k then
                   <<y := list(n := n+1)
                          . * (r := if domainp k then k else lc k)
                          .+ y;
                     if eqcar(r,'!:rn!:) then z := ilcm(z,cddr r);
                     if null domainp k then red k>>
                 else <<n := n + 1; k>>;
         if z neq 1 then
            y := for each j on y collect
                     lpow j .* if eqcar(lc j,'!:rn!:) then
                                  cadr lc j*z/cddr lc j
                                else lc j*z;
         if null x then x := y
          else x := b!:extmult(y,x)>>;
     sgn := evenp length lpow x;
     z := nil;
     for each j in lpow x do
         z := addf(if j = 0 then arnum!-mkglsol(0,x,sgn := not sgn,-1)
                    else multpf(mv to j,
                             arnum!-mkglsol(j,x,sgn := not sgn,-1)),z);
     return z
   end;

symbolic procedure arinv u;
   % Sort of half-extended gcd. No B-technique applied yet.
   % Performance is pretty bad.
   begin scalar mv,sgn,x,z,v,dmode!*;
         integer k,m,n;
     m := ldeg curdefpol!*;
     n := ldeg u;
     v := curdefpol!*;
     mv := mvar curdefpol!*;
     x := list(m-1) .* lc u .+ (list(-1) .* lc v .+ nil);
     for j := 2:(n+m) do
       begin scalar y;
         if j=(n+m) then y := list(-n-1) .* -1 .+ nil
          else nil;
         if (n + m - j - degr(v,mv) + 1) = 0 then v := red v;
         if (n + m - j - degr(u,mv) + 1) = 0 then u := red u;
             z := u;
         a:  if z and ((k := m - j + n - degr(z,mv))<m) then
               y := list k .* (if domainp z then z else lc z) .+ y
              else go to b;
             if null domainp z then z := red z
              else z := nil;
             go to a;
         b:  z := v;
         c:  if z and ((k := - j + m - degr(z,mv))<0) then
               y := list k .* (if domainp z then z else lc z) .+ y
              else go to d;
             if null domainp z then z := red z
              else z := nil;
             go to c;
         d:  x := b!:extmult(y,x)
       end;
     sgn := evenp length lpow x;
     z := nil;
     for each j in lpow x do
       if j > -1 then
         z := addf(if j = 0 then arnum!-mkglsol(0,x,sgn := not sgn,-n-1)
                    else multpf(mv to j,
			     arnum!-mkglsol(j,x,sgn:=not sgn,-n-1)),z);
     return z
   end;

symbolic procedure arnum!-mkglsol(u,v,sgn,n);
   begin scalar s,x,y,dmode!*;
     dmode!* := '!:rn!:;
     y := lpow v;
     for each j on red v do
       if s := arnum!-glsolterm(u,y,j,n)
          then x := s;
     return int!-equiv!-chk mkrn(if sgn then -x else x,lc v)
   end;

symbolic procedure arnum!-glsolterm(u,v,w,n);
   begin scalar x,y,sgn;
     x := lpow w;
     a: if null x then return
           if car y = n then lc w;
        if car x = u then return nil
         else if car x member v then <<x := cdr x;
                                     if y then sgn := not sgn>>
         else if y then return nil
               else <<y := list car x; x := cdr x>>;
        go to a
   end;

endmodule;

end;