1module arith;  % Header module for real arith package.
2
3% Redistribution and use in source and binary forms, with or without
4% modification, are permitted provided that the following conditions are met:
5%
6%    * Redistributions of source code must retain the relevant copyright
7%      notice, this list of conditions and the following disclaimer.
8%    * Redistributions in binary form must reproduce the above copyright
9%      notice, this list of conditions and the following disclaimer in the
10%      documentation and/or other materials provided with the distribution.
11%
12% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
13% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
14% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
15% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR
16% CONTRIBUTORS
17% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
18% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
19% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
20% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
21% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
22% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
23% POSSIBILITY OF SUCH DAMAGE.
24%
25
26
27% Last updated Dec 14, 1992
28
29% Assumptions being made on the underlying arithmetic:
30%
31% (1) The integer arithmetic is binary.
32%
33% (2) It is possible to convert any lisp float into an integer
34%     by applying fix, and this yields the result with the full
35%     precision of the float.
36%
37
38
39create!-package('(arith smlbflot bfauxil paraset math rounded comprd
40                  rdelem crelem bfelem), nil);
41
42flag('(arith),'core_package);
43
44exports ashift, bfloat, bfminusp, bfnzp, bfp!:, bfzp, crim, crrl,
45        divbf, ep!:, gfzerop, i2bf!:, lshift, make!:cr, make!:ibf,
46        make!:rd, msd!:, mt!:, oddintp, preci!:, rdp, retag, rndpwr,
47        sgn, tagrl, tagim, timbf;
48
49imports eqcar, round!:mt;
50
51global '(!!minnorm !!minnegnorm);
52
53fluid '(!*noconvert !:bprec!: dmode!*);
54
55
56switch noconvert;
57
58symbolic inline procedure mt!: u;
59   % This function selects the mantissa of U, a binary bigfloat
60   % representation of a number.
61   cadr u;
62
63symbolic inline procedure ep!: u;
64   % This function selects the exponent of U, a binary bigfloat
65   % representation of a number.
66   cddr u;
67
68symbolic inline procedure make!:ibf (mt, ep);
69   '!:rd!: . (mt . ep);
70
71symbolic inline procedure i2bf!: u; make!:ibf (u, 0);
72
73symbolic inline procedure make!:rd u;
74   '!:rd!: . u;
75
76symbolic inline procedure rdp x;
77   % This function returns true if X is a rounded number
78   % representation, else NIL.  X is any Lisp entity.
79   eqcar(x,'!:rd!:);
80
81symbolic inline procedure float!-bfp x; atom cdr x;
82
83symbolic inline procedure rd2fl x; cdr x;
84
85symbolic inline procedure fl2rd x; make!:rd x;
86
87symbolic inline procedure bfp!:(x);
88   % This function returns true if X is a binary bigfloat
89   % representation, else NIL.  X is any Lisp entity.
90   rdp x and not float!-bfp x;
91
92symbolic inline procedure retag u;
93   if atom u then u else '!:rd!: . u;
94
95symbolic inline procedure rndpwr j; normbf round!:mt(j,!:bprec!:);
96
97symbolic procedure msd!: m;
98  % returns the position n of the most significant (binary) digit
99  % of a positive binary integer m, i.e. floor(log2 m) + 1
100  begin integer i,j,k;
101    j := m;
102    while (j := ((k := j) / 65536)) neq 0 do i := i + 16;
103    j := k;
104    while (j := ((k := j) / 256)) neq 0 do i := i + 8;
105    j := k;
106    while (j := ((k := j) / 16)) neq 0 do i := i + 4;
107    j := k;
108    while (j := ((k := j) / 2)) neq 0 do i := i + 1;
109    return (i + 1);
110  end;
111
112% For both PSL and CSL this is a built-in, so making it inline here is not
113% useful and at one time it cause me pain.
114
115symbolic procedure ashift(m,d);
116   % This procedure resembles loosely an arithmetic shift.
117   %  It returns m*2**d
118   if d=0 then m
119    else if d<0 then m/2**(-d)
120    else m*2**d;
121
122symbolic inline procedure lshift(m,d);
123   % Variant of ashift that is called ONLY when m>=0.
124   %  This should be redefined for Lisp systems that provide
125   %  an efficient logical shift.
126   ashift(m,d);
127
128% And a shifts for use on signed inputs. For positive d this shifts left and
129% will be equivalent to multiplying by a power of 1. For negative d it
130% shifts right. It behaves live division by a power of two but rounding
131% towards -infinity. A different way to express that is that it takes the
132% binary representation of m and shifts it right just discarding any bits
133% that fall off the right hand end. For negative m it behaves as if there
134% was an infinite stream of "1" bits in high order positions, so those
135% come in to fill any vacated spaces. It is provided because the lshift
136% as per the specification above should not be used whem m<0.
137
138symbolic inline procedure shift(m,d);
139   if minusp m and minusp d then lnot lshift(lnot m,d)
140   else lshift(m,d);
141
142symbolic inline procedure oddintp n; not evenp n;
143
144symbolic inline procedure preci!: nmbr;
145   % This function counts the precision of a number "n". NMBR is a
146   % binary bigfloat representation of "n".
147   msd!: abs mt!: nmbr;
148
149
150symbolic inline procedure divbf(u,v); normbf divide!:(u,v,!:bprec!:);
151
152symbolic inline procedure timbf(u,v); rndpwr times!:(u,v);
153
154symbolic inline procedure bfminusp u;
155  if atom u then minusp u else minusp!: u;
156
157% The following function is not intended for very general use - it should
158% only be needed within code that evaluates elementary functions of
159% complex arguments. It differs from the regulat bfminusp function in that
160% a native floating point value of -0.0 is reported as being negative.
161% checking for this may be expensive!
162
163symbolic procedure special_bfminusp u;
164  if atom u then
165    if floatp u and u = 0.0 then
166% If I am checking of a floating point zero is "negative" in thie special
167% manner I want to detect an IEEE "-0.0" value. At present the Lisp systems
168% I use do not provide direct ways of achieving this, however using atan2
169% does the job. I hope that this case arises infrequently so that the
170% apparent costs are unimportant.
171      atan2(1.0, u) < 0.0      % atan2(1.0, -0.0) => -pi/2, while
172                               % atan2(1.0, +0.0) => pi/2
173    else minusp u
174  else minusp!: u;
175
176symbolic inline procedure bfzp u;
177  if atom u then zerop u else mt!: u=0;
178
179symbolic inline procedure bfnzp u; not bfzp u;
180
181symbolic inline procedure bfloat x;
182  if floatp x then fl2bf x
183   else normbf(if not atom x then x
184                else if fixp x then i2bf!: x
185                else read!:num x);
186
187symbolic inline procedure rdfl2rdbf x; fl2bf rd2fl x;
188
189symbolic inline procedure rd!:forcebf x;
190   % forces rounded number x to binary bigfloat representation
191   if float!-bfp x then rdfl2rdbf x else x;
192
193symbolic inline procedure crrl x; cadr x;
194
195symbolic inline procedure crim x; cddr x;
196
197symbolic inline procedure make!:cr (re,im);
198   '!:cr!: . (re . im);
199
200symbolic inline procedure crp x;
201   % This function returns true if X is a complex rounded number
202   % representation, else NIL.  X is any Lisp entity.
203   eqcar(x,'!:cr!:);
204
205symbolic inline procedure tagrl x; make!:rd crrl x;
206
207symbolic inline procedure tagim x; make!:rd crim x;
208
209symbolic inline procedure gfrl u; car u;
210
211symbolic inline procedure gfim u; cdr u;
212
213symbolic inline procedure mkgf (rl,im); rl . im;
214
215symbolic inline procedure gfzerop u;
216  if not atom gfrl u then mt!: gfrl u = 0 and mt!: gfim u = 0
217   else u = '(0.0 . 0.0);
218
219% symbolic inline procedure sgn x;
220%   if x>0 then 1 else if x<0 then -1 else 0;
221
222
223global '(bfz!* bfhalf!* bfone!* bftwo!* bfthree!* bffive!* bften!*
224         !:bf60!* !:180!* !:bf1!.5!*);
225
226global '(!:bf!-0!.25        %0.25
227         !:bf!-0!.0625      %0.0625
228         !:bf0!.419921875   %0.419921875
229        );
230
231%Miscellaneous constants
232
233bfz!* := make!:ibf(0,0);
234
235bfhalf!* := make!:ibf(1,-1);
236
237bfone!* := make!:ibf(1,0);
238
239!:bf1!.5!* := make!:ibf (3, -1);
240
241bftwo!* := make!:ibf (2, 0);
242
243bfthree!* := make!:ibf (3, 0);
244
245bffive!* := make!:ibf (5, 0);
246
247bften!* := make!:ibf (5, 1);
248
249!:bf60!* := make!:ibf (15, 2);
250
251!:180!* := make!:ibf (45, 2);
252
253!:bf!-0!.25 := make!:ibf(1,-2);
254
255!:bf!-0!.0625 := make!:ibf (1, -4);
256
257!:bf0!.419921875 := make!:ibf(215, -9);
258
259% These need to be added to other modules.
260
261symbolic procedure dn!:simp u;
262  if car u = 0 then nil ./ 1
263   else if u = '(10 . -1) and null !*noconvert then 1 ./ 1
264   else if dmode!* memq '(!:rd!: !:cr!:)
265    then rd!:simp cdr decimal2internal (car u, cdr u)
266   else if cdr u >= 0 then !*f2q (car u * 10**cdr u)
267   else simp {'quotient, car u, 10**(-cdr u)};
268
269put ('!:dn!:, 'simpfn, 'dn!:simp);
270
271symbolic procedure dn!:prin u;
272  bfprin0x (cadr u, cddr u);
273
274put ('!:dn!:, 'prifn, 'dn!:prin);
275
276endmodule;
277
278end;
279