Artifact 53ca062a1481b6093e377fea042db9eec13a89dea839aaf778cea23e9a48f9c3:
- Executable file
r37/packages/ratint/ratint.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: 20611) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/ratint/ratint.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: 20611) [annotate] [blame] [check-ins using]
% ------------------------------------------------------------------------- % 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;