Artifact bb1341c3d62cf7e4317e3a736c73eaec548df90a12657af3d7a8557e1eae2a88:
- Executable file
r37/packages/int/tidysqrt.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: 5767) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/int/tidysqrt.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: 5767) [annotate] [blame] [check-ins using]
module tidysqrt; % General tidying up of square roots. % Authors: Mary Ann Moore and Arthur C. Norman. % Modifications by J.H. Davenport. exports sqrt2top,tidysqrt; %symbolic procedure tidysqrtdf a; % if null a then nil % else begin scalar tt,r; % tt:=tidysqrt lc a; % r:=tidysqrtdf red a; % if null numr tt then return r; % return ((lpow a) .* tt) .+ r % end; % symbolic procedure tidysqrt q; begin scalar n1,dd; n1:=tidysqrtf numr q; if null n1 then return nil ./ 1; %answer is zero. dd:=tidysqrtf denr q; return multsq(n1,invsq dd) end; symbolic procedure tidysqrtf p; %Input - standard form. %Output - standard quotient. %Simplifies sqrt(a)**n with n>1. if domainp p then p ./ 1 else begin scalar v,w; v:=lpow p; if car v='i then v:=mksp('(sqrt -1),cdr v); %I->sqrt(-1); if eqcar(car v,'sqrt) and not onep cdr v then begin scalar x; %here we have a reduction to apply. x:=divide(cdr v,2); %halve exponent. w:=exptsq(simp cadar v,car x); %rational part of answer. if cdr x neq 0 then w := multsq(w,((mksp(car v,1) .* 1) .+ nil) ./ 1); %the next line allows for the horrors of nested sqrts. w:=tidysqrt w end else w:=((v .* 1) .+ nil) ./ 1; v:=multsq(w,tidysqrtf lc p); return addsq(v,tidysqrtf red p) end; symbolic procedure multoutdenr q; % Move sqrts in a sq to the numerator. begin scalar n,d,root,conj; n:=numr q; d:=denr q; while (root:=findsquareroot d) do << conj:=conjugatewrt(d,root); n:=!*multf(n,conj); d:=!*multf(d,conj) >>; while (root:=findnthroot d) do << conj:=conjugateexpt(d,root,kord!*); n:=!*multf(n,conj); d:=!*multf(d,conj) >>; return (n . d); end; symbolic procedure conjugateexpt(d,root,kord!*); begin scalar ord,ans,repl,xi; ord:=caddr caddr root; % the denominator of the exponent; ans:=1; kord!*:= (xi:=gensym()) . kord!*; % XI is an ORD'th root of unity; for i:=1:ord-1 do << ans:=!*multf(ans,numr subf(d, list(root . list('times,root,list('explt,xi,i))))); while (mvar ans eq xi) and ldeg ans > ord do ans:=addf(red ans,(xi) to (ldeg ans - ord) .* lc ans .+ nil); if (mvar ans eq xi) and ldeg ans = ord then ans:=addf(red ans,lc ans) >>; if (mvar ans eq xi) and ldeg ans = ord-1 then << repl:=-1; for i:=1:ord-2 do repl:=(xi) to i .* -1 .+ repl; ans:=addf(red ans,!*multf(lc ans,repl)) >>; if not domainp ans and mvar ans eq xi then interr "Conjugation failure"; return ans; end; symbolic procedure sqrt2top q; begin scalar n,d; n:=multoutdenr q; d:=denr n; n:=numr n; if d eq denr q then return q;%no change. if d iequal 1 then return (n ./ 1); q:=gcdcoeffsofsqrts n; if q iequal 1 then if minusf d then return (negf n ./ negf d) else return (n ./ d); q:=gcdf(q,d); n:=quotf(n,q); d:=quotf(d,q); if minusf d then return (negf n ./ negf d) else return (n ./ d) end; %symbolic procedure denrsqrt2top q; %begin % scalar n,d; % n:=multoutdenr q; % d:=denr n; % n:=numr n; % if d eq denr q % then return d; % no changes; % if d iequal 1 % then return 1; % q:=gcdcoeffsofsqrts n; % if q iequal 1 % then return d; % q:=gcdf(q,d); % if q iequal 1 % then return d % else return quotf(d,q) % end; symbolic procedure findsquareroot p; % Locate a sqrt symbol in poly p. if domainp p then nil else begin scalar w; w:=mvar p; %check main var first. if atom w then return nil; %we have passed all sqrts. if eqcar(w,'sqrt) then return w; w:=findsquareroot lc p; if null w then w:=findsquareroot red p; return w end; symbolic procedure findnthroot p; nil; % Until corrected. % symbolic procedure x!-findnthroot p; % % Locate an n-th root symbol in poly p. % if domainp p then nil % else begin scalar w; % w:=mvar p; %check main var first. % if atom w % then return nil; %we have passed all sqrts. % if eqcar(w,'expt) and eqcar(caddr w,'quotient) then return w; % w:=findnthroot lc p; % if null w then w:=findnthroot red p; % return w % end; symbolic procedure conjugatewrt(p,var); % Var -> -var in form p. if domainp p then p else if mvar p=var then begin scalar x,c,r; x:=tdeg lt p; %degree c:=lc p; %coefficient r:=red p; %reductum x:=remainder(x,2); %now just 0 or 1. if x=1 then c:=negf c; %-coefficient. return (lpow p .* c) .+ conjugatewrt(r,var) end else if ordop(var,mvar p) then p else (lpow p .* conjugatewrt(lc p,var)) .+ conjugatewrt(red p,var); symbolic procedure gcdcoeffsofsqrts u; if atom u then if numberp u and minusp u then -u else u else if eqcar(mvar u,'sqrt) then begin scalar v; v:=gcdcoeffsofsqrts lc u; if v iequal 1 then return v else return gcdf(v,gcdcoeffsofsqrts red u) end else begin scalar root; root:=findsquareroot u; if null root then return u; u:=makemainvar(u,root); root:=gcdcoeffsofsqrts lc u; if root iequal 1 then return 1 else return gcdf(root,gcdcoeffsofsqrts red u) end; endmodule; end;