1on rounded;
2
3
4
5% This is part of ongoing work... the main Reduce sources in arith/rounded.red
6% define safe!_fp!-plus etc but both CSL and PSL defined their own
7% "better" versions. The CSL one in cslbase/arith08.c and the PSL one
8% is packages/support/psl.red.
9
10% As at least a temporary measure the CSL code has been adjusted to try to
11% match the portable code. The PSL version should be trivial to do the
12% same to by commenting out the versions in psl.red. However a better
13% long term solution will be to end up with compatible nice versions.
14
15% The code here checks against a copy of the reference code. At present I
16% believe that CSL matches that. PSL is not expected to - but I *really*
17% do not understand the nature of some of the discrepancies and they look
18% a bit like bugs to me. If nothing else on one Linux I see results that
19% that suggest "1.0 = 1.0" sometimes says "no"...
20
21lisp;
22
23
24nil
25
26
27on echo,comp;
28
29
30nil
31
32load compiler;
33
34
35nil
36
37off pwrds;
38
39
40nil
41
42
43errs := 0;
44
45
460
47
48
49maxx := 1.797e+308;
50
51
521.797e+308
53
54
55
56ll := v := nil;
57
58
59nil
60
61for i := 0:20 do <<
62   w := expt(1.12345, i);
63   v := w . (1.0/w) . v;
64   v := (!!plumin*w) . (!!plumin/w) . v;
65   if w < maxx/!!plumax then
66     v := (!!plumax*w) . (!!plumax/w) . v;
67   v := (!!timmin*w) . (!!timmax/w) . v
68   >>;
69
70
71nil
72
73
74v;
75
76
77(2.16402e-153 4.62102e+152 4.56509e-295 4.3381e-297 10.2583 0.0974821
781.92623e-153 5.19148e+152 4.06346e-295 4.87364e-297 9.13106 0.109516
791.71457e-153 5.83237e+152 3.61694e-295 5.47529e-297 8.1277 0.123036 1.52616e-153
806.55238e+152 1.6257e+308 3.10607e+306 3.2195e-295 6.15121e-297 7.23459 0.138225
811.35846e-153 7.36127e+152 1.44706e+308 3.48952e+306 2.86572e-295 6.91058e-297
826.43962 0.155289 1.20919e-153 8.27002e+152 1.28805e+308 3.9203e+306 2.55082e-295
837.76369e-297 5.732 0.174459 1.07632e-153 9.29095e+152 1.14651e+308 4.40426e+306
842.27053e-295 8.72212e-297 5.10214 0.195996 9.58045e-154 1.04379e+153
851.02053e+308 4.94797e+306 2.02103e-295 9.79886e-297 4.54149 0.220192
868.52771e-154 1.17265e+153 9.08386e+307 5.55879e+306 1.79895e-295 1.10085e-296
874.04245 0.247375 7.59064e-154 1.31741e+153 8.08568e+307 6.24503e+306
881.60127e-295 1.23675e-296 3.59825 0.277913 6.75655e-154 1.48005e+153
897.19719e+307 7.01598e+306 1.42532e-295 1.38943e-296 3.20286 0.312221 6.0141e-154
901.66276e+153 6.40633e+307 7.8821e+306 1.2687e-295 1.56096e-296 2.85091 0.350765
915.35325e-154 1.86803e+153 5.70237e+307 8.85514e+306 1.12929e-295 1.75366e-296
922.53764 0.394067 4.76501e-154 2.09863e+153 5.07577e+307 9.94831e+306 1.0052e-295
931.97014e-296 2.25879 0.442715 4.2414e-154 2.35771e+153 4.51802e+307 1.11764e+307
948.9474e-296 2.21336e-296 2.01059 0.497368 3.77534e-154 2.64877e+153 4.02156e+307
951.25562e+307 7.96422e-296 2.4866e-296 1.78965 0.558768 3.36049e-154 2.97576e+153
963.57965e+307 1.41062e+307 7.08907e-296 2.79357e-296 1.593 0.627748 2.99122e-154
973.34312e+153 3.1863e+307 1.58476e+307 6.31009e-296 3.13844e-296 1.41795 0.705243
982.66253e-154 3.75582e+153 2.83618e+307 1.7804e+307 5.61671e-296 3.52588e-296
991.26214 0.792305 2.36996e-154 4.21948e+153 2.52452e+307 2.00019e+307
1004.99952e-296 3.96114e-296 1.12345 0.890115 2.10954e-154 4.74038e+153
1012.24712e+307 2.24712e+307 4.45015e-296 4.45015e-296 1.0 1.0)
102
103
104for each x in v do
105  ll := x . (-x) . ll;
106
107
108nil
109
110
111% ll is now a list of critical values
112
113length ll;
114
115
116324
117
118
119fluid '(errs);
120
121
122nil
123
124
125symbolic procedure badcase();
126  errs := errs + 1;
127
128
129badcase
130
131
132
133symbolic procedure portable!-fp!-plus(x,y);
134   if zerop x then y
135    else if zerop y then x
136    else if x>0.0 and y>0.0
137     then if x<!!plumax and y<!!plumax then plus2(x,y) else nil
138    else if x<0.0 and y<0.0
139     then if -x<!!plumax and -y<!!plumax then plus2(x,y) else nil
140    else if abs x<!!plumin and abs y<!!plumin then nil
141    else (if u=0.0 then u else if abs u<!!fleps1*abs x then 0.0 else u)
142         where u = plus2(x,y);
143
144
145portable!-fp!-plus
146
147
148symbolic procedure portable!-fp!-times(x,y);
149 if zerop x or zerop y then 0.0
150 else if x=1.0 then y else if y=1.0 then x else
151   begin scalar u,v; u := abs x; v := abs y;
152      if u>=1.0 and u<=!!timmax then
153         if v<=!!timmax then go to ret else return nil;
154      if u>!!timmax then if v<=1.0 then go to ret else return nil;
155      if u<1.0 and u>=!!timmin then
156         if v>=!!timmin then go to ret else return nil;
157      if u<!!timmin and v<1.0 then return nil;
158 ret: return times2(x,y) end;
159
160
161portable!-fp!-times
162
163
164symbolic procedure portable!-fp!-quot(x,y);
165 if zerop y then rdqoterr()
166 else if zerop x then 0.0 else if y=1.0 then x else
167   begin scalar u,v; u := abs x; v := abs y;
168      if u>=1.0 and u<=!!timmax then
169         if v>=!!timmin then go to ret else return nil;
170      if u>!!timmax then if v>=1.0 then go to ret else return nil;
171      if u<1.0 and u>=!!timmin then
172         if v<=!!timmax then go to ret else return nil;
173      if u<!!timmin and v>1.0 then return nil;
174 ret: return quotient(x,y) end;
175
176
177portable!-fp!-quot
178
179
180symbolic procedure tab_to n;
181  while posn() < n do prin2 " ";
182
183
184tab_to
185
186
187for each x in ll do
188  for each y in ll do <<
189     if errs < 20 then <<
190        a1 := safe!-fp!-plus(x, y);
191        a2 := portable!-fp!-plus(x, y);
192        if not eqn(a1, a2) then <<
193            terpri();
194            prin2t "safe-fp-plus incorrect";
195            prin2 "x: "; prin2 x; tab_to 40; prin2t hexfloat1 x;
196            prin2 "y: "; prin2 y; tab_to 40; prin2t hexfloat1 y;
197            prin2 "new: "; prin2 a1; tab_to 40; prin2t hexfloat1 a1;
198            prin2 "ref: "; prin2 a2; tab_to 40; prin2t hexfloat1 y;
199            terpri();
200            badcase() >>;
201        a1 := safe!-fp!-times(x, y);
202        a2 := portable!-fp!-times(x, y);
203        if not eqn(a1, a2) then <<
204            terpri();
205            prin2t "safe-fp-times incorrect";
206            prin2 "x: "; prin2 x; tab_to 40; prin2t hexfloat1 x;
207            prin2 "y: "; prin2 y; tab_to 40; prin2t hexfloat1 y;
208            prin2 "new: "; prin2 a1; tab_to 40; prin2t hexfloat1 a1;
209            prin2 "ref: "; prin2 a2; tab_to 40; prin2t hexfloat1 y;
210            terpri();
211            badcase() >>;
212        a1 := safe!-fp!-quot(x, y);
213        a2 := portable!-fp!-quot(x, y);
214        if not eqn(a1, a2) then <<
215            terpri();
216            prin2t "safe-fp-quot incorrect";
217            prin2 "x: "; prin2 x; tab_to 40; prin2t hexfloat1 x;
218            prin2 "y: "; prin2 y; tab_to 40; prin2t hexfloat1 y;
219            prin2 "new: "; prin2 a1; tab_to 40; prin2t hexfloat1 a1;
220            prin2 "ref: "; prin2 a2; tab_to 40; prin2t hexfloat1 y;
221            terpri();
222            badcase() >> >> >>;
223
224
225safe-fp-quot incorrect
226x: 1.0                                  10_0000_0000_0000_B1
227y: 4.45015e-296                         1d_1a94_a200_0000_B-981
228new: 2.24712e+295                       11_9799_812d_ea11_B982
229ref: nil                                1d_1a94_a200_0000_B-981
230
231
232safe-fp-quot incorrect
233x: 1.0                                  10_0000_0000_0000_B1
234y: -4.45015e-296                        -1d_1a94_a200_0000_B-981
235new: -2.24712e+295                      -11_9799_812d_ea11_B982
236ref: nil                                -1d_1a94_a200_0000_B-981
237
238
239safe-fp-quot incorrect
240x: 1.0                                  10_0000_0000_0000_B1
241y: 4.45015e-296                         1d_1a94_a200_0000_B-981
242new: 2.24712e+295                       11_9799_812d_ea11_B982
243ref: nil                                1d_1a94_a200_0000_B-981
244
245
246safe-fp-quot incorrect
247x: 1.0                                  10_0000_0000_0000_B1
248y: -4.45015e-296                        -1d_1a94_a200_0000_B-981
249new: -2.24712e+295                      -11_9799_812d_ea11_B982
250ref: nil                                -1d_1a94_a200_0000_B-981
251
252
253safe-fp-plus incorrect
254x: 1.0                                  10_0000_0000_0000_B1
255y: 2.24712e+307                         10_0000_0000_0000_B1022
256new: 2.24712e+307                       10_0000_0000_0000_B1022
257ref: nil                                10_0000_0000_0000_B1022
258
259
260safe-fp-quot incorrect
261x: 1.0                                  10_0000_0000_0000_B1
262y: 2.24712e+307                         10_0000_0000_0000_B1022
263new: nil                                nil
264ref: 4.45015e-308                       10_0000_0000_0000_B1022
265
266
267safe-fp-quot incorrect
268x: 1.0                                  10_0000_0000_0000_B1
269y: -2.24712e+307                        -10_0000_0000_0000_B1022
270new: nil                                nil
271ref: -4.45015e-308                      -10_0000_0000_0000_B1022
272
273
274safe-fp-plus incorrect
275x: 1.0                                  10_0000_0000_0000_B1
276y: 2.24712e+307                         10_0000_0000_0000_B1022
277new: 2.24712e+307                       10_0000_0000_0000_B1022
278ref: nil                                10_0000_0000_0000_B1022
279
280
281safe-fp-quot incorrect
282x: 1.0                                  10_0000_0000_0000_B1
283y: 2.24712e+307                         10_0000_0000_0000_B1022
284new: nil                                nil
285ref: 4.45015e-308                       10_0000_0000_0000_B1022
286
287
288safe-fp-quot incorrect
289x: 1.0                                  10_0000_0000_0000_B1
290y: -2.24712e+307                        -10_0000_0000_0000_B1022
291new: nil                                nil
292ref: -4.45015e-308                      -10_0000_0000_0000_B1022
293
294
295safe-fp-quot incorrect
296x: 1.0                                  10_0000_0000_0000_B1
297y: 3.96114e-296                         19_e7e0_24a4_ee94_B-981
298new: 2.52452e+295                       13_c391_aa12_e0e2_B982
299ref: nil                                19_e7e0_24a4_ee94_B-981
300
301
302safe-fp-quot incorrect
303x: 1.0                                  10_0000_0000_0000_B1
304y: -3.96114e-296                        -19_e7e0_24a4_ee94_B-981
305new: -2.52452e+295                      -13_c391_aa12_e0e2_B982
306ref: nil                                -19_e7e0_24a4_ee94_B-981
307
308
309safe-fp-quot incorrect
310x: 1.0                                  10_0000_0000_0000_B1
311y: 4.99952e-296                         10_592d_6928_0000_B-980
312new: 2.00019e+295                       1f_5172_1293_117a_B981
313ref: nil                                10_592d_6928_0000_B-980
314
315
316safe-fp-quot incorrect
317x: 1.0                                  10_0000_0000_0000_B1
318y: -4.99952e-296                        -10_592d_6928_0000_B-980
319new: -2.00019e+295                      -1f_5172_1293_117a_B981
320ref: nil                                -10_592d_6928_0000_B-980
321
322
323safe-fp-plus incorrect
324x: 1.0                                  10_0000_0000_0000_B1
325y: 2.52452e+307                         11_f9a6_b50b_0f28_B1022
326new: 2.52452e+307                       11_f9a6_b50b_0f28_B1022
327ref: nil                                11_f9a6_b50b_0f28_B1022
328
329
330safe-fp-quot incorrect
331x: 1.0                                  10_0000_0000_0000_B1
332y: 2.52452e+307                         11_f9a6_b50b_0f28_B1022
333new: nil                                nil
334ref: 3.96114e-308                       11_f9a6_b50b_0f28_B1022
335
336
337safe-fp-quot incorrect
338x: 1.0                                  10_0000_0000_0000_B1
339y: -2.52452e+307                        -11_f9a6_b50b_0f28_B1022
340new: nil                                nil
341ref: -3.96114e-308                      -11_f9a6_b50b_0f28_B1022
342
343
344safe-fp-quot incorrect
345x: 1.0                                  10_0000_0000_0000_B1
346y: 3.52588e-296                         17_0f22_3a63_c0dd_B-981
347new: 2.83618e+295                       16_342c_3c44_2242_B982
348ref: nil                                17_0f22_3a63_c0dd_B-981
349
350
351safe-fp-quot incorrect
352x: 1.0                                  10_0000_0000_0000_B1
353y: -3.52588e-296                        -17_0f22_3a63_c0dd_B-981
354new: -2.83618e+295                      -16_342c_3c44_2242_B982
355ref: nil                                -17_0f22_3a63_c0dd_B-981
356
357
358safe-fp-quot incorrect
359x: 1.0                                  10_0000_0000_0000_B1
360y: 5.61671e-296                         12_5dd6_68a2_4001_B-980
361new: 1.7804e+295                        1b_e073_6466_9b36_B981
362ref: nil                                12_5dd6_68a2_4001_B-980
363
364
365nil
366
367
368
369end;
370
371nil
372Tested on x86_64-pc-windows CSL
373Time (counter 1): 16 ms
374
375End of Lisp run after 0.01+0.06 seconds
376real 0.21
377user 0.03
378sys 0.04
379