Artifact 37296c12393883c5c7c205756fdc029a07eeac3ba6c8c462e8b9fa554753333a:
- Executable file
r37/packages/rataprx/contfrac.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: 14617) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/rataprx/contfrac.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: 14617) [annotate] [blame] [check-ins using]
module contfrac; % Continued fractions. % Author: Lisa Temme % Date: August 1995. % Code to check for rational polynomials. % polynomials and rational functions % by Winfried Neun symbolic procedure PolynomQQQ (x); (if fixp xx then 1 else if not onep denr (xx := cadr xx) then NIL else begin scalar kerns,kern,aa,var,fform,mvv,degg; fform := sfp mvar numr xx; var := reval cadr x; if fform then << xx := numr xx; while (xx neq 1) do << mvv := mvar xx; degg := ldeg xx; xx := lc xx; if domainp mvv then <<if not freeof(mvv,var) then << xx := 1 ; kerns := list list('sin,var) >> >> else kerns := append ( append (kernels mvv,kernels degg),kerns) >> >> else kerns := kernels !*q2f xx; aa: if null kerns then return 1; kern := first kerns; kerns := cdr kerns; if not(eq (kern, var)) and depends(kern,var) then return NIL else go aa; end) where xx = aeval(car x); put('PolynomQQ,'psopfn,'polynomQQQ); symbolic procedure ttttype_ratpoly(u); ( if fixp xx then 1 else if not eqcar (xx , '!*sq) then nil else and(polynomQQQ(list(mk!*sq (numr cadr xx ./ 1), reval cadr u)), polynomQQQ(list(mk!*sq (denr cadr xx ./ 1), reval cadr u))) ) where xx = aeval(car u); flag ('(type_ratpoly),'boolean); put('type_ratpoly,'psopfn,'ttttype_ratpoly); symbolic procedure type_ratpoly(f,z); ttttype_ratpoly list(f,z); %% To combine number, rational and non-rational approaches %% (including truncated versions) include the following %% boolean returns and the cfracrules rulelist. flag ('(vari),'boolean); symbolic procedure vari(x); idp x; procedure polynomialp(u,x); if den u = 1 and (freeof (u,x) or deg(u,x) >= 1 ) then t else nil; flag ('(polynomialp),'boolean); algebraic; operator cfrac; operator contfrac; procedure a_constant (x); lisp constant_exprp (x); cfracrules := { cfrac (~x) => (begin scalar cf, pt2, q, res; cf := continued_fraction x; pt2 := part(cf,2); res := for q := 2:(length pt2) collect append({1},{part(pt2,q)}); return contfrac(part(cf,1), append({part(pt2,1)},res)); end) when a_constant(x), cfrac (~x,~s) => (begin scalar kk, cf, cf1, pt2, cf2, cf3, bs, m, p, q, res; cf := continued_fraction(x); pt2 := part(cf,2); if s>=length(part(cf,2)) then << cf1 := for q:=2:(length pt2) collect append({1},{part(pt2,q)}); res := contfrac(part(cf,1), append({part(pt2,1)},cf1)); >> else << cf2 := for kk:=1:s+1 collect part(pt2,kk); bs := part(cf2,s+1); for m:= s step -1 until 1 do bs := part(cf2,m)+1/bs; cf3 := for p:=2:(length cf2) collect append({1},{part(cf2,p)}); res:=contfrac(bs,append({part(cf2,1)},cf3)) >>; %% res := continued_fraction(x,s); return res; end) when a_constant(x) and numberp s, cfrac (~x,~s) => (begin scalar cf, pt2, q, r, res; cf := cfracall(x,s); pt2 := part(cf,2); if type_ratpoly(x,s) then <<res := for r:=2:(length pt2) collect append({1},{part(pt2,r)}) >> else <<res := for q:=2:(length pt2) collect list(num(part(pt2,q)), den(part(pt2,q))) >>; return contfrac(part(cf,1), append({part(pt2,1)},res)); end) when not numberp x and vari s, cfrac(~a,~b,~c) => (begin scalar cf, pt2, q, res; cf := cfrac_ratpoly(a,b,c); pt2 := part(cf,2); res := for q:=2:(length pt2) collect append({1},{part(pt2, q)}); return contfrac(part(cf,1), append({part(pt2,1)},res)); end) when numberp c and vari b and type_ratpoly(a,b), cfrac(~a,~b,~c) => (begin scalar cf, pt2, q, res; cf := cfrac_nonratpoly(a,b,c); pt2 := part(cf, 2); res := for q:=2:length(pt2) collect list(num(part(pt2,q)), den(part(pt2,q))); return contfrac(part(cf,1), append({part(pt2,1)},res)); end) when numberp c and vari b and NOT(type_ratpoly(a,b))%, }; let cfracrules; % LOAD Taylor Package for non-rationals load taylor; %INPUT my code for rational polynomials procedure cfracall(rat_poly,var); begin scalar top_poly, bot_poly, euclidslist, ld_return; if type_ratpoly(rat_poly,var) then << top_poly := num rat_poly; bot_poly := den rat_poly; euclidslist := {}; while part(longdiv(top_poly, bot_poly, var),2) neq 0 do << ld_return := longdiv(top_poly, bot_poly, var); top_poly := bot_poly; bot_poly := part(ld_return,2); euclidslist := part(ld_return,1).euclidslist; >>; euclidslist := part(longdiv(top_poly, bot_poly, var),1) . euclidslist; return list(inv_cfracall(reverse(euclidslist)), reverse(euclidslist)); >> else << return cfrac_nonratpoly(rat_poly,var,5) >>; end; %************ %INPUT my code for rational polynomials (truncated) procedure cfrac_ratpoly(rat_poly,var,number); begin scalar top_poly, bot_poly, euclidslist, ld_return, k; if type_ratpoly(rat_poly,var) then << top_poly := num rat_poly; bot_poly := den rat_poly; euclidslist := {}; k:=number; %-1; while part(longdiv(top_poly, bot_poly, var),2) neq 0 and k neq 0 do << ld_return := longdiv(top_poly, bot_poly, var); top_poly := bot_poly; bot_poly := part(ld_return,2); euclidslist := part(ld_return,1).euclidslist; k := k-1; >>; euclidslist := part(longdiv(top_poly, bot_poly, var),1) . euclidslist; return list(inv_cfracall(reverse(euclidslist)), reverse(euclidslist)); >> else << return cfrac_nonratpoly(rat_poly,var,number) >>; end; procedure longdiv(poly1, poly2,x); begin scalar numer, denom, div, div_list, elmt, flag, rem, answer; %longdiv called by cfracall so poly2 will never be zero. %on rounded; numer := poly1; denom := poly2; div_list := {}; div := 0; flag := 0; answer := 0; if longdivdeg(numer,x) < longdivdeg(denom,x) then rem := numer else << while (longdivdeg(numer,x) >= longdivdeg(denom,x)) AND flag neq 1 do << if longdivlterm(numer,x) = 0 then << div := numer/denom; rem :=0; flag :=1; >> else << div := longdivlterm(numer,x)/longdivlterm(denom,x); numer := numer - denom*div; rem := numer; >>; div_list := div.div_list; >>; answer := for each elmt in div_list sum elmt; >>; return list(answer,rem) end; procedure longdivdeg(i_p,i_p_var); begin scalar a; a:= if numberp(den(i_p)) then deg(i_p*den(i_p),i_p_var); return a end; procedure longdivlterm(i_p,i_p_var); begin scalar b; b := if numberp(den(i_p)) then lterm(den(i_p)*i_p,i_p_var)/den(i_p); return b end; %**************** %Check for a polynomial %% flag ('(type_poly),'boolean); %% put('type_poly,'psopfn,'PolynomQQQ); %INPUT my code for non-rationals procedure cfrac_nonratpoly(nonrat,x,n); begin scalar hh,g, a_0, a_1, coeff_list, flag1, flag2, k, j, h, oneplus, xover; g := taylor(nonrat,x,0,2*n); h := 1; k:=n; if taylorp(taylortostandard g) then rederr "not yet implemented" else << %%CHANGE TO: if not type_poly then ERROR %Include error here so that COEFF can be used in while condition if not type_ratpoly(taylortostandard g,x) or (type_ratpoly(taylortostandard g,x) and not(freeof(den(taylortostandard g),x))) then rederr "not yet implemented"; while (length(coeff(taylortostandard g, x)) >1 and k>=0) do %0) do << %%CHANGE TO: if not type_poly then ERROR %Include error here so that each time a new "g" is generated % it will be checked to see if it is a polynomial if not type_ratpoly(taylortostandard g,x) or (type_ratpoly(taylortostandard g,x) and not(freeof(den(taylortostandard g),x))) then rederr "not yet implemented"; a_0 := first coeff(taylortostandard g, x); a_1 := second coeff(taylortostandard g, x); if flag1 =0 then << coeff_list := {a_0}; flag1 := 1; >> else << if a_1 neq 0 then << g := taylorcombine(a_1*taylor(x^h,x,0,2*n)/(g - a_0)); coeff_list := (a_1*x^h).coeff_list; >> else << j := 2; while j <= length coeff(taylortostandard g, x) and flag2=0 do << if coeffn(taylortostandard g, x, j) neq 0 then << a_n := coeffn(taylortostandard g, x, j); flag2 := 1; >> else j := j+1; >>; coeff_list := (a_n*x^j).coeff_list; g := taylorcombine(a_n*taylor(x^j,x,0,2*n)/(g - a_0)); flag2 := 0; h := j >>; >>; k := k-1; >>; %% %"1+" form %% oneplus := list(inv_cfrac_nonratpoly1(reverse(coeff_list)), %% reverse(coeff_list)); %"x/" form xover:= list(inv_cfrac_nonratpoly2(adaptcfrac(reverse(coeff_list))), adaptcfrac(reverse(coeff_list))); return xover %% list(oneplus,xover) >>; end; %*************** %INPUT my code for different representation of cfrac_nonratpoly procedure adaptcfrac(l_list); begin scalar h, l, k, n, m, new_list; new_list := {}; if length l_list < 3 then return l_list else << h := first l_list; l := second l_list; k := 2; while length l_list >= k do << n := num l; d := den l; new_list := (n/d).new_list; k := k+1; if length l_list >= k then << l := part(l_list, k); l := d*l >>; >>; >>; return h.reverse(new_list) end; procedure inv_cfrac_nonratpoly1(c_list); begin scalar ans, j, expan; j := length c_list; if j < 3 then << ans := for each m in c_list sum m; return ans; >> else << for k:=j step -1 until 2 do << if k=j then expan := part(c_list,k) else expan := part(c_list,k) / (1 + expan); >>; expan := part(c_list,1) + expan; return expan >>; end; procedure inv_cfrac_nonratpoly2(c_list); begin scalar ans, j, expan; j := length c_list; if j < 3 then << ans := for each m in c_list sum m; return ans; >> else << for k:=j step -1 until 2 do << if k=j then expan := part(c_list,k) else expan := num(part(c_list,k)) / (den(part(c_list,k)) + expan); >>; expan := part(c_list,1) + expan; return expan >>; end; procedure inv_cfracall(c_list); begin scalar ans, j; j := length c_list; if j=0 then return {} else << if j=1 then ans := part(c_list,1) else << ans := part(c_list,j); for k:=j-1 step -1 until 1 do << ans := part(c_list,k) + 1/ans >>; >>; >>; return ans end; symbolic procedure print!-contfract(x); % printing continued fractions begin scalar xx,xxx; if null !*nat or atom x or length x < 3 or not eqcar(caddr x,'list) then return 'failed; xx := reverse cddr caddr x; if length xx > 12 then return 'failed; if xx then <<xxx := list('quotient ,cadr first xx,caddr first xx); for each tt in rest xx do xxx := list('quotient ,cadr tt,list('plus,caddr tt,xxx)); if cadr caddr x = 0 then maprin list('list,cadr x,xxx) else maprin list('list,cadr x,list ('plus,cadr caddr x ,xxx)); >> else maprin list('list,cadr x,cadr caddr x); return t; end; put('contfrac,'prifn,'print!-contfract); endmodule; end;