Artifact c8f90d0d50b5a9ddd7425eac58da020fdaa35f6331cc2eb2e71a89fcd639380a:
- Executable file
r38/packages/poly/specfac.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: 10905) [annotate] [blame] [check-ins using] [more...]
module specfac; % Splitting of low degree polynomials. % Author: Anthony C. Hearn. % Copyright (c) 1991 RAND. All rights reserved. fluid '(!*keepsqrts !*sub2 !*surds knowndiscrimsign kord!* zlist); % switch surds; exports cubicf,quadraticf,quarticf; symbolic procedure coeffs pol; % Extract coefficients of polynomial wrt its main variable and leading % degree. Result is a list of coefficients. begin integer degree,deg1; scalar cofs,mv; mv := mvar pol; degree := ldeg pol; while not domainp pol and mvar pol eq mv do <<deg1 := ldeg pol; for i:= 1:(degree-deg1-1) do cofs := nil . cofs; cofs := lc pol . cofs; pol := red pol; degree := deg1>>; for i:=1:degree-1 do cofs := nil . cofs; return reversip(pol . cofs) end; symbolic procedure shift!-pol pol; % Shifts main variable, mv, of square free nth degree polynomial pol so % that coefficient of mv**(n-1) is zero. % Does not assume pol is univariate. begin scalar lc1,ld,mv,pol1,redp,shift,x; mv := mvar pol; ld := ldeg pol; redp := red pol; if domainp redp or not(mvar redp eq mv) or ldeg redp<(ld-1) then return list(pol,1,nil ./ 1); lc1 := lc pol; x := lc redp; shift := quotsq(!*f2q x,!*f2q multd(ld,lc1)); pol1 := subf1(pol,list(mv . mk!*sq addsq(!*k2q mv,negsq shift))); return list(numr pol1,denr pol1,shift) end; symbolic procedure quadraticf!*(pol,var); if domainp pol then errach "invalid quadratic to factr" else if mvar pol = var then quadraticf pol else begin scalar kord,w; kord := kord!*; kord!* := list var; w := coeffs !*q2f resimp(pol ./ 1); kord!* := kord; w := quadraticf1(car w,cadr w,caddr w); if w eq 'failed then return list(1,pol); var := !*k2f var; return list(if car w neq 1 then mkrn(1,car w) else 1, addf(multf(var,cadr w),caddr w), addf(multf(var,cadddr w),car cddddr w)) end; symbolic procedure quadraticf pol; % Finds factors of square free quadratic polynomial pol (if they % exist). Does not assume pol is univariate. (if x eq 'failed then list(1,pol) else if not domainp car x then list(1,pol) % Answer would be rational. else list(if car x neq 1 then mkrn(1,car x) else 1, y .* cadr x .+ caddr x,y .* cadddr x .+ car cddddr x) where y = (mvar pol .** 1)) where x = quadraticf1(car w,cadr w,caddr w) where w = coeffs pol; symbolic procedure quadraticf1(a,b,c); begin scalar a1,denom,discrim,w; if null b and minusf c and not minusf a then <<a := rootxf(a,2); c := rootxf(negf c,2); return if a eq 'failed or c eq 'failed then 'failed else list(1,a,c,a,negf c)>>; discrim := powsubsf addf(exptf(b,2),multd(-4,multf(a,c))); % A null discriminator can arise from a polynomial such as % 16x^2+(32i-8)*x-8i-15; if null discrim then nil else <<if knowndiscrimsign then <<if knowndiscrimsign eq 'negative then return 'failed>> % else if not clogflag and minusf discrim % then return 'failed; else if minusf discrim then return 'failed; discrim:=rootxf(discrim,2); if discrim='failed then return discrim>>; denom := multd(4,a); a := a1 := multd(2,a); w := addf(b,discrim); c := addf(b,negf discrim); b := w; if (w := gcdf(a,b)) neq 1 then <<a1 := quotf(a,w); b := quotf(b,w); denom := quotf(denom,w)>>; if (w := gcdf(a,denom)) neq 1 and (w := gcdf(c,denom)) then <<a := quotf(a,w); c := quotf(c,w); denom := quotf(denom,w)>>; return list(denom,a1,b,a,c) end; symbolic procedure rootxf(u,n); % Return either polynomial nth root of u or "failed". begin scalar x,y,z,w; if domainp u then return if minusf u then 'failed else if atom u and (y := irootn(u,n))**n=u then y else if not atom u and (x := get(car u,'rootfn)) then apply2(x,u,n) else if !*surds and not(u member zlist) then nrootn!*(u,n) else 'failed; x := comfac u; u := quotf(u,comfac!-to!-poly x); z := 1; if car x then if cdr(y := divide(cdar x,n)) = 0 then z := multpf(caar x .** car y,z) else if !*surds then <<z := multf(mkrootf(caar x,n,cdr y),z); if car y neq 0 then z := multpf(caar x .** car y,z)>> else return 'failed; x := cdr x; if domainp x then if minusf x then return 'failed else if fixp x and (y := irootn(x,n))**n=x then z := multd(y,z) else if !*surds and fixp x then z := multf(nrootn!*(x,n),z) else if not atom x and (w := get(car x,'rootfn)) then apply2(w,x,n) else return 'failed else if (y := rootxf(x,n)) eq 'failed then return y else z := multf(y,z); if u=1 then return z; x := sqfrf u; c: if null x then return z else if cdr(y := divide(cdar x,n)) = 0 then <<z := multf(exptf(caar x,car y),z); x := cdr x>> else if !*surds then <<z := multf(mkrootf(prepf caar x,n,cdr y), multf(exptf(caar x,car y),z)); x := cdr x>> else return 'failed; go to c end; symbolic procedure mkrootf(u,m,n); if m neq 2 or null !*keepsqrts then !*p2f mksp(list('expt,u,list('quotient,1,m)),n) else if n neq 1 then errach 'mkrootf else !*q2f simpsqrt list u; symbolic procedure nrootn!*(u,n); % Returns a standard form representation of the nth root of u. begin scalar x; if null u then return nil; u := nrootn(u,n); x := cdr u; % surd part. u := car u; % rational part. if x=1 then return x; x := mkrootf(prepf x,n,1); return powsubsf multf(u,x) end; symbolic procedure cubicf pol; % Split the cubic pol if a change of origin puts it in the form % (x-a)**3-b=0. begin scalar a,a0,a1,b,neg,p; p := shift!-pol pol; a := coeffs car p; if cadr a then return list(1,pol) % Cadr a non nil probably means there are some surds in the % coefficients that don't reduce to 0. else if caddr a then return list(1,pol); % Factorization not possible by this method. a0 := cadddr a; a := car a; if minusf a0 then <<neg := t; a0 := negf a0>>; if (a := rootxf(a,3)) eq 'failed or (a0 := rootxf(a0,3)) eq 'failed then return list(1,pol); if neg then a0 := negf a0; a := !*f2q a; a0 := !*f2q a0; p := addsq(!*k2q mvar pol,caddr p); % Now numr (a*(mv+shift)+a0) is a factor of pol. a1 := numr addsq(multsq(a,p),a0); % quotf(pol,a) is quadratic factor. However, the surd division may % not work properly, so we calculate factor directly. b := multsq(a0,a0); b := addsq(b,multsq(negsq multsq(a,a0),p)); b := numr addsq(b,multsq(multsq(a,a),exptsq(p,2))); return aconc!*(quadraticf b,a1) end; symbolic procedure powsubsf u; % We believe that the result of this operation must be a polynomial. % If subs2q returns a rational, it must be because there are % unsimplified surds. Hopefully rationalizesq can fix those. begin scalar !*sub2; u := subs2q !*f2q u; if denr u neq 1 then <<u := rationalizesq u; if denr u neq 1 then errach list('powsubsf,u)>>; return numr u end; symbolic procedure quarticf pol; % Splits quartics that can be written in the form % (x-a)**4+b*(x-a)**2+c. % Note that any call of rootxf can lead to a result "failed." begin scalar !*sub2,a,a2,a0,b,dsc,p,p1,p2,q,shift,var; var := mvar pol; p := shift!-pol pol; a := coeffs car p; shift := caddr p; if cadr a % pol not correctly shifted, possibly due to sqrt. % e.g., 729para^4*be^4 - 81para^3*sqrt(27*be^2*para^2 - 8cte1^3)* % sqrt(3)*be^3 - 216para^2*be^2*cte1^3 + 12para*sqrt(27be^2*para^2 % - 8*cte1^3)*sqrt(3) *be*cte1^3 + 8*cte1^6. or cadddr a then return list(1,pol); % Factorization not possible by this method. a2 := cddr a; a0 := caddr a2; a2 := car a2; a := car a; q := quadraticf1(a,a2,a0); if not(q eq 'failed) then <<a2 := car q; q := cdr q; a := exptsq(addsq(!*k2q mvar pol,shift),2); b := numr subs2q quotsq(addsq(multsq(!*f2q car q,a), !*f2q cadr q), !*f2q cadr p); a := numr subs2q quotsq(addsq(multsq(!*f2q caddr q,a), !*f2q cadddr q), !*f2q cadr p); a := quadraticf!*(a,var); b := quadraticf!*(b,var); return multf(a2,multf(car a,car b)) . nconc!*(cdr a,cdr b)>> else if null !*surds or denr shift neq 1 then return list(1,pol); % Factorization not possible by this method. shift := numr shift; if knowndiscrimsign eq 'negative then go to complex; dsc := powsubsf addf(exptf(a2,2),multd(-4,multf(a,a0))); p2 := minusf a0; if not p2 and minusf dsc then go to complex; p1 := not a2 or minusf a2; if not p1 then if p2 then p1 := t else p2 := t; p1 := if p1 then 'positive else 'negative; p2 := if p2 then 'negative else 'positive; a := rootxf(a,2); if a eq 'failed then return list(1,pol); dsc := rootxf(dsc,2); if dsc eq 'failed then return list(1,pol); p := invsq !*f2q addf(a,a); q := multsq(!*f2q addf(a2,negf dsc),p); p := multsq(!*f2q addf(a2,dsc),p); b := multf(a,exptf(addf(!*k2f mvar pol,shift),2)); a := powsubsf addf(b,q); b := powsubsf addf(b,p); knowndiscrimsign := p1; a := quadraticf!*(a,var); knowndiscrimsign := p2; b := quadraticf!*(b,var); knowndiscrimsign := nil; return multf(car a,car b) . nconc!*(cdr a,cdr b); % Complex case. complex: a := rootxf(a,2); if a eq 'failed then return list(1,pol); a0 := rootxf(a0,2); if a0 eq 'failed then return list(1,pol); a2 := powsubsf addf(multf(2,multf(a,a0)),negf a2); a2 := rootxf(a2,2); if a2 eq 'failed then return list(1,pol); % Now a*(x+shift)**2 (+/-) b*(x+shift) + c is a factor. p := addf(!*k2f mvar pol,shift); q := addf(multf(a,exptf(p,2)),a0); p := multf(a2,p); a := powsubsf addf(q,p); b := powsubsf addf(q,negf p); knowndiscrimsign := 'negative; a := quadraticf!*(a,var); b := quadraticf!*(b,var); knowndiscrimsign := nil; return multf(car a,car b) . nconc!*(cdr a,cdr b) end; endmodule; end;