module dmodeop; % Generic operators for domain arithmetic.
% Author: Anthony C. Hearn.
% Copyright (c) 1991 The RAND Corporation. All rights reserved.
% internal dividef;
fluid '(!*noequiv); % !*convert
% switch convert;
% !*convert := t;
symbolic procedure !:difference(u,v);
if null u then !:minus v else if null v then u
else if u=v then nil
else if atom u and atom v then u-v
else dcombine(u,v,'difference);
symbolic procedure !:divide(u,v);
% Returns a dotted pair of quotient and remainder of non-invertable
% domain element U divided by non-invertable domain element V.
% Note that a zero is returned as NIL.
if null u then nil . nil
else if null v then rerror(poly,202,"zero divisor")
else if atom u and atom v then dividef(u,v)
else dcombine(u,v,'divide);
symbolic procedure dividef(m,n);
((if car x=0 then nil else car x) . if cdr x=0 then nil else cdr x)
where x=divide(m,n);
symbolic procedure !:expt(u,n);
% Raises domain element U to integer power N. Value is a domain
% element.
if null u then if n=0 then rerror(poly,11,"0/0 formed") else nil
else if n=0 then 1
else if n=1 then u
else if u=1 then 1
else if n<0
then !:recip !:expt(if not fieldp u then mkratnum u else u,-n)
else if atom u then u**n
% Moved into the exponentiation method of !:mod!:
% else if car u eq '!:mod!:
% then (lambda x; if x=0 then nil else if x=1 then 1 else car u . x)
% general!-modular!-expt(cdr u,n)
else begin scalar v,w,x;
if x := get(car u,'expt)
then return apply2(x,u,n);
% There was a special exponentiation method.
v := apply1(get(car u,'i2d),1); % unit element.
x := get(car u,'times);
a: w := n;
if w-2*(n := n/2) neq 0 then v := apply2(x,u,v);
if n=0 then return v;
u := apply2(x,u,u);
go to a
end;
symbolic procedure !:gcd(u,v);
if null u then v else if null v then u
else if atom u and atom v then gcdn(u,v)
else if fieldp u or fieldp v then 1
else dcombine(u,v,'gcd);
% symbolic procedure !:i2d u;
symbolic procedure !:minus u;
% U is a domain element. Value is -U.
if null u then nil
else if atom u then -u
else (if x then apply1(x,u) else dcombine(u,-1,'times))
where x=get(car u,'minus);
symbolic procedure !:minusp u;
if atom u then minusp u else apply1(get(car u,'minusp),u);
symbolic procedure !:onep u;
if atom u then onep u else apply1(get(car u,'onep),u);
symbolic procedure !:plus(u,v);
if null u then v else if null v then u
else if atom u and atom v
then (if w=0 then nil else w) where w=u+v
else dcombine(u,v,'plus);
% symbolic procedure !:prep u;
% symbolic procedure !:print u;
symbolic procedure !:quotient(u,v);
if null u or u=0 then nil
else if null v or v=0 then rerror(poly,12,"Zero divisor")
else if atom u and atom v
% We might also check that remainder is zero in integer case.
then if null dmode!* then u/v else
(if atom recipv then u*recipv else dcombine(u,recipv,'times))
where recipv=!:recip v
else dcombine(u,v,'quotient);
symbolic procedure !:recip u;
% U is an invertable domain element. Value is 1/U.
begin
if numberp u
then if abs u=1 then return u
else if null dmode!* or dmode!* memq '(!:rd!: !:cr!:)
then return !:rn2rd mkrn(1,u)
else u := apply1(get(dmode!*,'i2d),u);
return (if not atom x and car x='!:rn!: then !:rn2rd x else x)
where x=dcombine(1,u,'quotient)
end;
symbolic procedure !:rn2rd x;
% Convert rn to rd in dmodes rd and cr if roundall is on.
if !*roundall and !*rounded then !*rn2rd x else x;
symbolic procedure !:times(u,v);
% We assume neither u nor v can be 0.
if null u or null v then nil
else if atom u and atom v then u*v
else dcombine(u,v,'times);
symbolic procedure !:zerop u;
if null u or u=0 then t
else if atom u then nil
else apply1(get(car u,'zerop),u);
symbolic procedure fieldp u;
% U is a domain element. Value is T if U is invertable, NIL
% otherwise.
not atom u and flagp(car u,'field);
symbolic procedure gettransferfn(u,v);
% This may be unnecessary. If dmodechk has been called, then all
% transfer functions should be defined.
(if x then x else dmoderr(u,v)) where x=get(u,v);
symbolic procedure dcombine(u,v,fn);
% U and V are domain elements, but not both atoms (integers).
% FN is a binary function on domain elements;
% Value is the domain element representing FN(U,V)
% or pair of domain elements representing divide(u,v).
<<u := if atom u
then apply2(get(car v,fn),apply1(get(car v,'i2d),u),v)
else if atom v
then apply2(get(car u,fn),u,apply1(get(car u,'i2d),v))
else if car u eq car v then apply2(get(car u,fn),u,v) else
% convert anything to :ps: but not the reverse;
% convert real to complex, never the reverse;
% also convert rn or crn to rd or cr but not the reverse:
% hence crn or gi plus rd requires that *both* convert to cr.
(<<if (not(x and atom x)
or (get(du,'cmpxfn) and not get(dv,'cmpxfn)
or du memq dl and not(dv memq dl)) and dv neq '!:ps!:)
% extra test added above by Alan Barnes to ensure
% result is :ps: if either operand is a :ps:
and not flagp(dv,'noconvert) then
% convert v -> u but may first have to convert u.
<<if du memq dml and dv eq '!:rd!:
or du eq '!:rd!: and dv memq dml then
<<u := apply1(get(du,'!:cr!:),u); du := '!:cr!:>>
else if du eq '!:rn!: and dv eq '!:gi!:
or du eq '!:gi!: and dv eq '!:rn!: then
<<u := apply1(get(du,'!:crn!:),u); du := '!:crn!:>>;
v := apply1(get(dv,du),v); x := get(du,fn)>>
else <<u := apply1(x,u); x := get(dv,fn)>>;
apply2(x,u,v)>>
where x=get(du,dv),dml='(!:crn!: !:gi!:),dl='(!:rd!: !:cr!:))
where du=car u,dv=car v;
if !*rounded and !*roundall and not atom u then
% atom test added by Alan Barnes in case a power series
% operation has already produced an integer.
int!-equiv!-chk
if (v := car u) eq '!:rn!: and cddr u neq 1 then !*rn2rd u
else if v eq '!:crn!: and (cdadr u neq 1 or cdddr u neq 1)
then !*crn2cr u
else u
else if fn eq 'divide then % Modified by Francis Wright.
int!-equiv!-chk car u . int!-equiv!-chk cdr u
else int!-equiv!-chk u>>;
symbolic procedure int!-equiv!-chk u;
% U is a domain element. If U can be converted to 0, result is NIL,
% if U can be converted to 1, result is 1,
% if U is a rational or a complex rational and can be converted to
% an integer, result is that integer,
% if *convert is on and U can be converted to an integer, result
% is that integer. Otherwise, U is returned.
% In most cases, U will be structured.
if !*noequiv then u
else begin scalar x;
if atom u then return if u=0 then nil else u
else if apply1(get(car u,'zerop),u) then return nil
else if apply1(get(car u,'onep),u) then return 1
% else if null !*convert then return u
else if (x := get(car u,'intequivfn)) and (x := apply1(x,u))
then return if x=0 then nil else x
else return u
end;
% symbolic procedure minuschk u;
% if eqcar(u,'minus)
% and (numberp cadr u
% or not atom cadr u and idp caadr u and get(caadr u,'dname))
% then !:minus cadr u
% else u;
endmodule;
end;