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