reduce-historical/r37/packages/ratint/ratint.red

% -------------------------------------------------------------------------
% Neil Langmead
% ZIB Berlin, December 1996 / January 1997
% 
% Package to integrate rational functions. Uses the Hermite Horowitz Rothstein
% Trager algorithms to determine firstly, the reduction of the rational fn
% into its polynomial and logarithmic parts, then the integration  of these
% parts seperately.

create!-package('(ratint convert),nil);

global '(traceratint!*); % for the tracing facility

switch traceratint;

off traceratint; % set to off by default

algebraic;   

% We first need a few utility functions, given below
% 
% -------------------------------------------------------------------------
% routine to return the square free factorisation of a polynomial
% outputs factors in a list, coupled with their exponents.
% However, such a function is already defined in poly/facform.

% expr procedure square_free(poly,x);
% begin  scalar  !*EXP, !*FACTOR, !*DIV, !*RATARG, !*GCD, !*RATIONAL,
% l,k,c,w,z,y,output, answer, result;
% on exp; on div; on ratarg; off factor; on rational;
% k:=1; output:={}; answer:= df(poly,x);
% poly:=poly;
% c:=gcd(poly,answer);
% w:=poly/c;
% while (c neq 1) do
%        <<
%          y:=gcd(w,c);
%          z:=w/y;
%          if (z neq 1) then
%           output:=append({{z,k}},output);
%           k:=k+1;
%           w:=y; c:=c/y;
%         >>;
% output:=append({{w,k}},output);
% return output;
% end;

% expr procedure eval_square_free(list,var);
% begin scalar output;
% off exp;
% output:= for l:=1:arglength(list)
% product(part(part(list,l),1)^(part(part(list,l),2)));
% return output;
% end;

expr procedure make_mon(li,var);

begin scalar current, li2;
li2:={};
for k:=1:arglength(li) do
  <<
    current:=part(li,k);
    current:=monic(current,var);
    li2:=append(li2,{current});
   >>;
return(li2);
end;

expr procedure monic(exp,var);
begin scalar lecof, temp;
lecof:=lcof(exp,var);
exp:=exp/lecof;
return exp;
end;
%x^8-2*x^6+2*x^2-1;
%z:=square_free(x^8-2*x^6+2*x^2-1,x);
%res:= for l:=1:2 product(part(part(z,l),1)^(part(part(z,l),2)));

% -------------------------------------------------------------------------
% An implementation of the extended Eucledian algorithm. Given polynomials
% a and b in the variable x in a Euclidean domain D, this computes elements
% s and t in D such that g:=gcd(a,b)=sa+tb

%off factor;
% input is polynomials a and b, in variable var
on rational;
%expr procedure u(exp,var);
%lcof(exp,var);

expr procedure norm(exp,var);
exp/lcof(exp,var);

expr procedure gcd_ex(a,b,var);
begin scalar c, c1, c2,d,d1,d2,q,r,m,r1,r2,g,s,b;
on rational;
c:=norm(a,var); d:=norm(b,var);
c1:=1; d1:=0; c2:=0; d2:=1;

while(d neq 0) do
  <<
    on rational;
    m:=pseudorem(c,d,var); 
    q:=part(pseudorem(c,d,var),3)/part(pseudorem(c,d,var),2);
    %q:=part(pf(c/d,var),1); q should be quo(c,d)
    %off rational;
    r:=c-q*d;
    %r:=part(m,1);
    r1:=c1-q*d1; r2:=c2-q*d2;
    c:=d; c1:=d1; c2:=d2;
    d:=r; d1:=r1; d2:=r2;
  >>;
s:=c1/(u(a,var)*u(c,var));
b:=c2/(u(b,var)*u(c,var));
return({s,b});
end;

expr procedure u(exp,var);
if(numberp(exp)) then exp else
lcof(exp,var);

%l:=3*x^3+x^2+x+5;% p:=5*x^2-3*x+1
%pseudorem(l,p,x);
%quote:=part(pseudorem(o,p,x),3)/part(pseudorem(o,p,x),2);
%gcd_ex(x^2-1,x+1,x);
%a:=x^2+(7/6)*x+(1/3); b:=2*x+(7/6);
%gcd_ex(a,b,x);
%f:=48*x^3-84*x^2+42*x-36; g:=-4*x^3-10*x^2+44*x-30;
%gcd_ex(f,g,x);

% -------------------------------------------------------------------------
% routine to remove elements from a list which are zero

expr procedure rem_zero(li);
begin scalar k,j,l,li1,li2;
for k:=1:arglength(li) do
   <<
     j:=part(li,k);
     if(j neq 0) then nil else
       <<
         li1:=for l:=1:k-1 collect part(li,l);
         li2:=for l:=k+1:arglength(li) collect part(li,l);
         li:=append(li1,li2);
       >>;
   >>;
return li;
end;

% -------------------------------------------------------------------------- 
% This routine takes as input a rational function p/q^(expt), returning p, q
% and expt seperately in a list

expr procedure hr_monic_den(li,x);
begin scalar !*EXP, !*FACTOR, q, lc;
on EXP;
li:= (for each r in li collect << lc:=lcof(den(r),x);
				  {num(r)/lc, den(r)/lc} >> );
on FACTOR;
li:= (for each r in li collect 
     <<  q:=part(r,2);
         if(arglength(q) > -1 and part(q,0)=expt) then
         {part(r,1),part(q,1),part(q,2)} else
         {part(r,1),part(q,1),1} >> ) ; return li;
end;

%q:={3/(x+3)^2,4/(x+1)^5};
%hr_monic_den(q,x);
%in "monic";

% -------------------------------------------------------------------------
% The implementation of the Rostein Trager algorithm
% Takes as input a/b in x, and returns a two element list, with the polynomial
% and logarithmic parts of the integral. For aesthetic reasons in REDUCE, the
% values aren't added. This should be done manually by the user, but is a 
% trivial task.
% in "mkmonic.red";
operator c, rtof,v, alpha;
load_package arnum;

load_package assist;

expr procedure not_numberp(x);
if (not numberp(x)) then t else nil;

expr procedure rt(a,b,x);
begin scalar vv, j,k,i,current,sol,res,cc,b_prime,extra_term,
      current1,vvv,integral, eqn,d,v_list, sol1,sol2, temp, temp2;
b_prime:=df(b,x);
v_list:={};
on rational;
res:=resultant(a-z*b_prime,b,x);
on rational; on ifactor;
res:= old_factorize(res);
res:=extractlist(res, not_numberp);
res:=make_mon(res,z);
res:=mkset(res); % removes duplicates by turning list into a set
%write "res is ", res;
   
integral:=0;

for k:=1:arglength(res) do
<<
  current:=part(res,k);% write "current is ", current;
  d:=deg(current,z); %write "d is ", d;
  if(d=1) then
    <<
      sol:=solve(current=0,z);
      sol:=part(sol,1);
      cc:=part(sol,2);% write "cc is ", cc; 
      vv:=gcd(a-cc*b_prime,b);% write "vv is " , vv;
      vv:=vv/(lcof(vv,x));
      extra_term:=append({cc},{log(vv)});% write extra_term;
      extra_term:=part(extra_term,1)*part(extra_term,2); %write extra_term;
      integral:=extra_term+integral; %write "integral is ", integral;
    >>
   else
    <<
      current:=sub(z=alpha,current);% write "current is ", current;
      current1:=sub(alpha=alp,current);% write "current1 is ", current1;
      defpoly(current);
      %write "alpha is ", alpha;
      a:=sub(x=z,a); b:=sub(x=z,b); b_prime:=sub(x=z,b_prime);
      vv:=gcd(a-alpha*b_prime,b);% write "vv is ", vv; % OK up to here
      
      off arnum;
      on fullroots;
      vv:=sub(a1=alpha*8,vv); 
      vv:=sub(z=x,vv);
      vvv:=solve(vv=0,x);
      vvv:=sub(a1=1/alp,vvv);% write "vvv is ", vvv;
      eqn:=part(part(vvv,1),1)-part(part(vvv,1),2); 
      %write "eqn is ", eqn;
      if(d=2) then 
         <<
             sol:=solve(current1=0,alp);
             sol1:=part(sol,1); sol2:=part(sol,2);
             %write "sol1, 2 are ", sol1, sol2;
             c(1):=part(sol1,2); c(2):=part(sol2,2);
             %write "c(1), c(2) are ", c(1), c(2);
             for j:=1:2 do
             <<    
               v(j):=sub(alp=c(j),eqn);
               integral:=integral+c(j)*log(v(j));
               %write "integral is ", integral;
             >>;
           >>
       else
        <<
           k:=1;
           %write "d is ", d;
           while (k<=d) do
           %for k:=1:3 do
           <<
              c(k):=rtof(current1);% write "c(k) is ", c(k);
              v(k):=sub({alp=c(k)},eqn);% write "v(k) is ", v(k);
              integral:=integral+c(k)*log(v(k));
              %write "integral is ", integral;
              k:=k+1;
           >>;
        >>;
    >>;
lisp null remprop ('alpha,'currep);
lisp null remprop ('alpha,'idvalfn);
>>;
return(integral);
end;

% -------------------------------------------------------------------------
% This piece of code was written by Matt Rebeck. Input are the functions
% p, q and variable x. It returns the pseudo remainder of p and q, and the
% quotient.

symbolic procedure prem(r,v,var);
  begin
    scalar d,dr,dv,l,n,tt,rule_list,m,q,input1,input2,rr,vv;
    on rational; off factor;
    rr := r;
    vv := v;
    dr := deg(r,var);
    dv := deg(v,var);
    if dv <= dr then 
    <<
      l := reval coeffn(v,var,dv);
      v := reval{'plus,v,{'minus,{'times,l,{'expt,var,dv}}}};
    >>
    else l := 1;
    d := dr-dv+1;
    n := 0;
    while dv<=dr and r neq 0 do
    <<
      tt := reval{'times,{'expt,var,(dr-dv)},v,coeffn(r,var,dr)};
      if dr = 0 then r := 0 
      else 
      <<
        rule_list := {'expt,var,dr}=>0;
        let rule_list;
        r := reval r;
        clearrules rule_list;
      >>;
      r := reval{'plus,{'times,l,r},{'minus,tt}};
      dr := deg(r,var);
      n := n+1;
    >>;
    r := reval{'times,{'expt,l,(d-n)},r};
    m := reval{'expt,l,d};
    input1 := reval{'plus,{'times,{'expt,l,d},rr},{'minus,r}};
    input2 := vv;
    q := reval{'quotient,input1,input2};
    return {r,m,q};
  end;

procedure pseudorem(x,y,var); lisp ('list . prem(x,y,var));
%e.g.
 %pseudorem(3x^5+4,x^2+1,x);
%exp1:=441*x^7+780*x^6-2861*x^5+4085*x^4+7695*x^3+3713*x^2-43253*x+24500;
%exp2:=9*x^6+6*x^5-65*x^4+20*x^3+135*x^2-154*x+49;
%pseudorem(exp1,exp2,x);
%a:=x^8+x^6-3*x^4-3*x^3+8*x^2+2*x-5;
%b:=3*x^6+5*x^4-4*x^2-9*x+21;
%pseudorem(a,b,x);
%r:=-15*x^4+3*x^2-9;
%rr:=
%operator a,c, neil;
%in "rem";

% -------------------------------------------------------------------------
% this routine is the implementation of Horowitz' method of reducing the 
% rational function into a polynomial and logarithmic part.
 
operator a,c, neil ; 
expr procedure howy(p,q,x);

begin
scalar pseudo, quo,rem,pp, poly_part,d,mm,b,nn,j,k,aa,cc,
pseudo3,i,quo3,r,pseudo2,eqn,l,neil1,sol, var1, temp,var2,var3,p,test,output;
pseudo:=pseudorem(p,q,x);    
quo:=part(pseudo,3)/part(pseudo,2);
rem:=part(pseudo,1)/part(pseudo,2);

poly_part:=quo;
pp:=rem;% write "pp is ", pp;   
d:=gcd(q,df(q,x));
pseudo2:=pseudorem(q,d,x);
b:=part(pseudo2,3)/part(pseudo2,2);

mm:=deg(b,x); 
nn:=deg(d,x);

aa:=for k:=0:(mm-1) sum (a(k))*(x^(k)); 
cc:=for j:=0:(nn-1) sum (c(j))*(x^(j));

var1:=for i:=0:(mm-1) collect a(i); 
var2:=for k:=0:(nn-1) collect c(k);
var3:=append(var1,var2);
%write var3;
on rational;
pseudo3:=pseudorem(b*df(d,x),d,x);
quo3:=part(pseudo3,3)/part(pseudo3,2);

temp:=b*df(d,x)/d; 

temp:=pseudorem(num(temp),den(temp),x);
temp:=part(temp,3)/part(temp,2);
%write "temp is ", temp;
r:=b*df(cc,x)-cc*temp+d*aa;% write "r is: ", r;

for k:=0:(mm+nn-1) do
     <<
       %on factor;
       neil(k):=coeffn(pp,x,k)-coeffn(r,x,k);
       %write "neil(k)= ", neil(k);
     >>;

neil1:=for k:=0:(mm+nn-1) collect neil(k)=0;
%write "neil1= ", neil1;

sol:=solve(neil1,var3);% write "sol= ", sol;
sol:=first(sol);% write "sol= ", sol;
aa:=sub(sol,aa);
%write "aa= ", aa;
%aa:=for k:=1:mm sum(part(sol,k));
cc:=sub(sol,cc);% write "cc is ", cc;
ans1:=cc/d;
ans2:=int(poly_part,x);
ans3:=(aa/b);
output:={ans1,ans2,ans3};
return output;
end;

% -------------------------------------------------------------------------
%in "eea"; in "rem"; in "phi";

expr procedure newton(a,p,u1,w1,B);
begin scalar alpha,gamma,eea_result,s,tt,u,w,ef,modulus,c,sigma,
sigma_tilde,tau, tau_tilde,re,r,quo;

alpha:=lcof(a,x);
a:=alpha*a;
gamma:=alpha;
%a:=gamma*a;

%u1:=n(u1,x);
u1:=phi(u1,x,p); write "u1 is ", u1;

off modular;
w1:=phi(alpha*w1,x,p); off modular;% write "w1 is ", w1;
%w1:=phi(alpha*w1,x,p); off modular;% write "w1 is ",w1;

eea_result:=gcd_ex(u1,w1,x);

on modular; setmod p;
s:=part(eea_result,1);   tt:=part(eea_result,2);

   
on modular; setmod p;
u:=replace_lc(u1,x,gamma);   w:=replace_lc(w1,x,alpha);

off modular; 
ef:=a-u*w; off modular;  modulus:=p;
%write "ef is ", ef; write "u,w are  ", u,w;
% iterate until either the factorisation in Z[x] is obtained, or else
% the bound on modulus is reached

on modular; setmod p;
while(ef neq 0 and modulus<2*B*gamma) do
    <<

       c:=ef/modulus;% write "c is ", c; off modular;
      sigma_tilde:=phi(s*c,x,p); off modular;
      tau_tilde:=phi(tt*c,x,p); off modular;
      %
       re:=pseudorem(sigma_tilde,w1,x);
      
      r:=part(re,1)/part(re,2); quo:=part(re,3)/part(re,2);
     sigma:=re; tau:=phi(tau_tilde+quo*u1,x,p); off modular;
      % update the factors and compute the error
      u:=u+tau*modulus;  w:=w+sigma*modulus;
     % write "u is ",u; write "w is ",w; write "ef is ", ef;
      ef:=a-u*w; modulus:=modulus*p; 
   >>; 

% check termination status
	if(ef=0) then << u:=u; w:=w/gamma;   >> 
                          else rederr "nsfe";
return {u,w};
end; 

%trst newton;
 %newton(12*x^3+10*x^2-36*x+35,5,x,x^2+3,10000);
% in "phi.red";
% in "eea.red"; in "rem.red"; %in "replace_lc";
clear p;
%trst newton;trst newton;
%newton(12*x^3+10*x^2-36*x+35,5,2*x,x^2+2,10000)

% -------------------------------------------------------------------------
expr procedure replace_lc(exp,var,val);
begin scalar lead_term, new_lead_term,red;
lead_term:=lterm(exp,var);
red:=reduct(exp,var);
new_lead_term:=lead_term/lcof(exp,var);
new_lead_term:=new_lead_term*val;
new_exp:=new_lead_term+red;

return new_exp;
end;

% -------------------------------------------------------------------------
% in "rem"; in "eea";
% routine to solve the polynomial diophantine equation
% s(x)a(x)+t(x)b(x)=c(x) for the unknown polynomials s and t

expr procedure polydi(a,b,c,x);
begin scalar q,r, sigma,tau, s,tt,sol,sigma_tilde,tau_tilde,g;
on rational;
g:=gcd(a,b);

s:=part(gcd_ex(a,b,x),1);
tt:=part(gcd_ex(a,b,x),2);
sol:=(s*c/g)*a+(tt*c/g)*b;
% here, sol=c(x), our right hand side

sigma_tilde:=s*c/g;
tau_tilde:=tt*c/g;

result:=pseudorem(sigma_tilde,b/g,x);
q:=part(result,3)/part(result,2);
r:=part(result,1)/part(result,2);

sigma:=r;
tau:=tau_tilde+(q*(a/g));

return {sigma,tau};
end;
%in "rem"; in "eea";
%trst polydi;
%polydi(x+(7/3),1,294,x);
%polydi(x^2+(7/6)*x+(1/3),2*x+(7/6),-(4425/2)*x-(5525/4),x);

% -------------------------------------------------------------------------
expr procedure phi(exp,var,p);
 begin scalar prime;
 prime:=p;
 if(primep p) then << on modular;
                      setmod p;
                      exp:=exp mod p;% off modular;
                   >> else rederr "p sahould be prime";
 return exp;
 end;


expr procedure nn(exp,var,p);
 begin scalar lcoef;
 lcoef:=lcof(exp,var);
 if(primep p) then <<   on modular; setmod p;
                        exp:=exp/lcoef;
             
                  >> else rederr "p should be prime";

return exp;
 end;
off modular;
%in "ratint" %in "examples"
%in "make_monic"
%operator c, rtof,v, alpha;
load_package arnum;

%load_package assist

expr procedure not_numberp(x);
if (not numberp(x)) then t else nil;

operator log_sum;
expr procedure rt(a,b,x);
begin scalar vv, j,k,i,current,sol,res,cc,b_prime,extra_term,
      current1,vvv,integral, eqn,d,v_list, sol1,sol2, temp, temp2;
b_prime:=df(b,x);
v_list:={};
on rational;
res:=resultant(a-z*b_prime,b,x);
on rational; on ifactor;
res:= old_factorize(res);
res:=extractlist(res, not_numberp);
res:=make_mon(res,z);
res:=mkset(res); % removes duplicates by turning list into a set
%write "res is ", res;
   
integral:=0;

for k:=1:arglength(res) do
<<
  current:=part(res,k); %write "current is ", current;
  d:=deg(current,z); %write "d is ", d;
  if(d=1) then
    <<
      sol:=solve(current=0,z);
      sol:=part(sol,1);
      cc:=part(sol,2);% write "cc is ", cc; 
      vv:=gcd(a-cc*b_prime,b);% write "vv is " , vv;
      vv:=vv/(lcof(vv,x));
      extra_term:=append({cc},{log(vv)});% write extra_term;
      extra_term:=part(extra_term,1)*part(extra_term,2);%write extra_term;
      integral:=extra_term+integral; %write "integral is ", integral
  if(lisp !*traceratint) then write "integral in Rothstein T is ", integral;
    >>
   else
    <<
      current:=sub(z=alpha,current);% write "current is ", current;
      current1:=sub(alpha=alp,current);% write "current1 is ", current1;
      off mcd; current:=current;
      defpoly(current);
      on mcd;
      %write "alpha is ", alpha; %write part(alpha,1);
      a:=sub(x=z,a); b:=sub(x=z,b); b_prime:=sub(x=z,b_prime);
      vv:=gcd(a-alpha*b_prime,b); % write "vv is ", vv; % OK up to here
      off arnum;
      on fullroots;
      %vv:=sub(a1=alpha*(1/part(alpha,1)),vv); 
      vv:=sub(z=x,vv); % vv:=sub(a1=part(alpha,1)*alpha,vv);
      %write "vv is ", vv;
      %write "deg is ", deg(vv,x);
      on rational; on ratarg;
      % write "deg is ", deg(vv,x);
      if(deg(vv,x)>2) then 
<<
% we want to give the answer not in terms of a complete pf decomposition,
% but without splitting the field

integral:=integral+log_sum(alpha_a,current1,0,alpha_a*log(vv),x);
integral:=sub(alpha_a=alpha,integral);
%integral:=sub(a1=part(alpha,1)*alpha,integral)
integral:=sub(alp=alpha,integral);
if(lisp !*traceratint) then write "integral in Rothstein T is ", integral;
>> 
else % degree less than or eq to 2, so no problem solving vv=0
<< 
      % write "current is ", current;
      current:=sub(alpha=beta,current); 
      vv:=sub(alpha=beta,vv);
      integral:=integral+log_sum(beta,current,0,beta*log(vv));
      % vvv:=solve(vv=0,x);
      %vvv:=sub(a1=1/alp,vvv); write "vvv is ", vvv;
      %eqn:=part(part(vvv,1),1)-part(part(vvv,1),2); 
      %write "eqn is ", eqn;
      if(d=2) then 
         <<
             sol:=solve(current1=0,alp);
             sol1:=part(sol,1); sol2:=part(sol,2);
             %write "sol1, 2 are ", sol1, sol2;
             c(1):=part(sol1,2); c(2):=part(sol2,2);
             %write "c(1), c(2) are ", c(1), c(2)
             for j:=1:2 do
             <<    
               v(j):=sub(alp=c(j),eqn);
               %integral:=integral+c(j)*log(v(j));
               % write "integral is ", integral;
             >>;
           >>
       else
        <<
           k:=1;
           %write "d is ", d;
           while (k<=d) do
           %for k:=1:3 do
           <<
              c(k):=rtof(current1); % write "c(k) is ", c(k);
              v(k):=sub({alp=c(k)},eqn); % write "v(k) is ", v(k);
              %integral:=integral+c(k)*log(v(k));
              %write "integral is ", integral;
              k:=k+1;
           >>;
        >>;
    >>;
>>;
lisp null remprop ('alpha,'currep);
lisp null remprop ('alpha,'idvalfn);
>>;
return(integral);
end;

expr procedure dependp(exp,x);

  if(freeof(exp,x) and not numberp(exp)) then nil else t ;
% -------------------------------------------------------------------------
% procedure to integrate any rational function, using the implementations of
% the above algorithms.

expr procedure ratint(p,q,x);
begin scalar s_list,first_term, second_term, r_part, answer;
% check input carefully

if(not dependp(p,x) and not dependp(q,x)) then
   return (p/q)*x
else <<
   if(not dependp(p,x) and dependp(q,x)) then return
      p*ratint(1,q,x)
 else <<
   if( dependp(p,x) and not dependp(q,x)) then return
       (1/q)*int(p,x);
      >>;
     >>;
if(numberp p and numberp q) then return (p/q)*x;
%if(not polynomp p or not polynomp q) then rederr "input must be polynomials"
if(lisp !*traceratint) then write "performing Howoritz reduction on ", p/q;
s_list:=howy(p,q,x); 
if(lisp !*traceratint) then write "Howoritz gives: ", s_list;
first_term:=part(s_list,1);
second_term:=part(s_list,2);
r_part:=part(s_list,3);% write "r_part is ", r_part;
if(lisp !*traceratint) then write "computing Rothstein Trager on ", r_part;
r_part:=rt(num(r_part),den(r_part),x);
answer:={first_term+second_term,r_part};
return (answer);
end;

% examples

%exp1:=441*x^7+780*x^6-2861*x^5+4085*x^4+7695*x^3+3713*x^2-43253*x+24500;

%exp2:=9*x^6+6*x^5-65*x^4+20*x^3+135*x^2-154*x+49;

%k:=36*x^6+126*x^5+183*x^4+(13807/6)*x^3-407*x^2-(3242/5)*x+(3044/15);
%l:=(x^2+(7/6)*x+(1/3))^2*(x-(2/5))^3;

%ratint(k,l,x);

%aa:=7*x^13+10*x^8+4*x^7-7*x^6-4*x^3-4*x^2+3*x+3;
%bb:=x^14-2*x^8-2*x^7-2*x^4-4*x^3-x^2+2*x+1;
%trst ratint;
%ratint(aa,bb,x);
%-----------------------------------------------------------------------------

end;