1 /*  arith08.c                         Copyright (C) 1990-2008 Codemist Ltd */
2 
3 /*
4  * Arithmetic functions.
5  */
6 
7 /**************************************************************************
8  * Copyright (C) 2008, Codemist Ltd.                     A C Norman       *
9  *                                                                        *
10  * Redistribution and use in source and binary forms, with or without     *
11  * modification, are permitted provided that the following conditions are *
12  * met:                                                                   *
13  *                                                                        *
14  *     * Redistributions of source code must retain the relevant          *
15  *       copyright notice, this list of conditions and the following      *
16  *       disclaimer.                                                      *
17  *     * Redistributions in binary form must reproduce the above          *
18  *       copyright notice, this list of conditions and the following      *
19  *       disclaimer in the documentation and/or other materials provided  *
20  *       with the distribution.                                           *
21  *                                                                        *
22  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS    *
23  * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT      *
24  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS      *
25  * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE         *
26  * COPYRIGHT OWNERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,   *
27  * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,   *
28  * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS  *
29  * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
30  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR  *
31  * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF     *
32  * THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH   *
33  * DAMAGE.                                                                *
34  *************************************************************************/
35 
36 
37 
38 /* Signature: 4e04c032 22-Aug-2010 */
39 
40 #include "headers.h"
41 
42 
43 #ifdef COMMON
44 
Lboole(Lisp_Object nil,int nargs,...)45 static Lisp_Object MS_CDECL Lboole(Lisp_Object nil, int nargs, ...)
46 {
47     Lisp_Object r, op, a, b;
48     va_list aa;
49     argcheck(nargs, 3, "boole");
50     va_start(aa, nargs);
51     op = va_arg(aa, Lisp_Object);
52     a = va_arg(aa, Lisp_Object);
53     b = va_arg(aa, Lisp_Object);
54     va_end(aa);
55     switch (op)
56     {
57 case fixnum_of_int(boole_clr):
58         return onevalue(fixnum_of_int(0));
59 case fixnum_of_int(boole_and):
60         r = logand2(a, b);
61         break;
62 case fixnum_of_int(boole_andc2):
63         push(a);
64         b = lognot(b);
65         pop(a);
66         errexit();
67         r = logand2(a, b);
68         break;
69 case fixnum_of_int(boole_1):
70         return onevalue(a);
71 case fixnum_of_int(boole_andc1):
72         push(b);
73         a = lognot(a);
74         pop(b);
75         errexit();
76         r = logand2(a, b);
77         break;
78 case fixnum_of_int(boole_2):
79         return onevalue(b);
80 case fixnum_of_int(boole_xor):
81         r = logxor2(a, b);
82         break;
83 case fixnum_of_int(boole_ior):
84         r = logior2(a, b);
85         break;
86 case fixnum_of_int(boole_nor):
87         a = logior2(a, b);
88         errexit();
89         r = lognot(a);
90         break;
91 case fixnum_of_int(boole_eqv):
92         r = logeqv2(a, b);
93         break;
94 case fixnum_of_int(boole_c2):
95         r = lognot(b);
96         break;
97 case fixnum_of_int(boole_orc2):
98         b = lognot(b);
99         errexit();
100         r = logior2(a, b);
101         break;
102 case fixnum_of_int(boole_c1):
103         r = lognot(a);
104         break;
105 case fixnum_of_int(boole_orc1):
106         push(b);
107         a = lognot(a);
108         pop(b);
109         errexit();
110         r = logior2(a, b);
111         break;
112 case fixnum_of_int(boole_nand):
113         a = logand2(a, b);
114         errexit();
115         r = lognot(a);
116         break;
117 case fixnum_of_int(boole_set):
118         return onevalue(fixnum_of_int(-1));
119 default:
120         return aerror1("bad arg for boole",  op);
121     }
122     errexit();
123     return onevalue(r);
124 }
125 
Lbyte(Lisp_Object nil,Lisp_Object a,Lisp_Object b)126 static Lisp_Object Lbyte(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
127 {
128     a = cons(a, b);
129     errexit();
130     return onevalue(a);
131 }
132 
Lbyte_position(Lisp_Object nil,Lisp_Object a)133 static Lisp_Object Lbyte_position(Lisp_Object nil, Lisp_Object a)
134 {
135     if (!consp(a)) return aerror1("byte-position", a);
136     else return onevalue(qcdr(a));
137 }
138 
Lbyte_size(Lisp_Object nil,Lisp_Object a)139 static Lisp_Object Lbyte_size(Lisp_Object nil, Lisp_Object a)
140 {
141     if (!consp(a)) return aerror1("byte-size", a);
142     else return onevalue(qcar(a));
143 }
144 
Lcomplex_2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)145 static Lisp_Object Lcomplex_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
146 {
147 /* /* Need to coerce a and b to the same type here... */
148     a = make_complex(a, b);
149     errexit();
150     return onevalue(a);
151 }
152 
Lcomplex_1(Lisp_Object nil,Lisp_Object a)153 static Lisp_Object Lcomplex_1(Lisp_Object nil, Lisp_Object a)
154 {
155 /* /* Need to make zero of same type as a */
156     a = make_complex(a, fixnum_of_int(0));
157     errexit();
158     return onevalue(a);
159 }
160 
Lconjugate(Lisp_Object nil,Lisp_Object a)161 static Lisp_Object Lconjugate(Lisp_Object nil, Lisp_Object a)
162 {
163     if (!is_number(a)) return aerror1("conjugate", a);
164     if (is_numbers(a) && is_complex(a))
165     {   Lisp_Object r = real_part(a),
166                     i = imag_part(a);
167         push(r);
168         i = negate(i);
169         pop(r);
170         errexit();
171         a = make_complex(r, i);
172         errexit();
173         return onevalue(a);
174     }
175     else return onevalue(a);
176 }
177 
Ldenominator(Lisp_Object nil,Lisp_Object a)178 static Lisp_Object Ldenominator(Lisp_Object nil, Lisp_Object a)
179 {
180     CSL_IGNORE(nil);
181     if (!is_number(a)) return aerror1("denominator", a);
182     if (is_numbers(a) && is_ratio(a))
183         return onevalue(denominator(a));
184     else return onevalue(fixnum_of_int(1));
185 }
186 
Ldeposit_field(Lisp_Object nil,int nargs,...)187 static Lisp_Object MS_CDECL Ldeposit_field(Lisp_Object nil, int nargs, ...)
188 {
189     va_list aa;
190     Lisp_Object a, b, c;
191     CSL_IGNORE(nil);
192     CSL_IGNORE(nargs);
193     va_start(aa, nargs);
194     a = va_arg(aa, Lisp_Object);
195     b = va_arg(aa, Lisp_Object);
196     c = va_arg(aa, Lisp_Object);
197     va_end(aa);
198     return aerror("deposit-field");
199 }
200 
Ldpb(Lisp_Object nil,int nargs,...)201 static Lisp_Object MS_CDECL Ldpb(Lisp_Object nil, int nargs, ...)
202 {
203     va_list aa;
204     Lisp_Object a, b, c;
205     CSL_IGNORE(nil);
206     CSL_IGNORE(nargs);
207     va_start(aa, nargs);
208     a = va_arg(aa, Lisp_Object);
209     b = va_arg(aa, Lisp_Object);
210     c = va_arg(aa, Lisp_Object);
211     va_end(aa);
212     return aerror("dpb");
213 }
214 
Lffloor(Lisp_Object nil,Lisp_Object a,Lisp_Object b)215 static Lisp_Object Lffloor(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
216 {
217     CSL_IGNORE(nil);
218     CSL_IGNORE(a);
219     CSL_IGNORE(b);
220     return aerror("ffloor");
221 }
222 
Lgcd_n(Lisp_Object nil,int nargs,...)223 Lisp_Object MS_CDECL Lgcd_n(Lisp_Object nil, int nargs, ...)
224 {
225     va_list a;
226     int i;
227     Lisp_Object r;
228     if (nargs == 0) return fixnum_of_int(0);
229     va_start(a, nargs);
230     push_args(a, nargs);
231 /*
232  * The actual args have been passed a C args - I can not afford to
233  * risk garbage collection until they have all been moved somewhere safe,
234  * and here that safe place is the Lisp stack.  I have to delay checking for
235  * overflow on same until all args have been pushed.
236  */
237     stackcheck0(nargs);
238     pop(r);
239     for (i = 1; i<nargs; i++)
240     {   Lisp_Object w;
241         if (r == fixnum_of_int(1))
242         {   popv(nargs-i);
243             break;
244         }
245         pop(w);
246         r = gcd(r, w);
247         errexitn(nargs-i-1);
248     }
249     return onevalue(r);
250 }
251 
Lgcd(Lisp_Object nil,Lisp_Object a,Lisp_Object b)252 Lisp_Object Lgcd(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
253 {
254     a = gcd(a, b);
255     errexit();
256     return onevalue(a);
257 }
258 
Lgcd_1(Lisp_Object nil,Lisp_Object a)259 Lisp_Object Lgcd_1(Lisp_Object nil, Lisp_Object a)
260 {
261     CSL_IGNORE(nil);
262     return onevalue(a);
263 }
264 
Limagpart(Lisp_Object nil,Lisp_Object a)265 static Lisp_Object Limagpart(Lisp_Object nil, Lisp_Object a)
266 {
267     CSL_IGNORE(nil);
268     if (!is_number(a)) return aerror1("imagpart", a);
269     if (is_numbers(a) && is_complex(a))
270         return onevalue(imag_part(a));
271 /* /* the 0.0 returned here ought to be the same type as a has */
272     else return onevalue(fixnum_of_int(0));
273 }
274 
Lldb(Lisp_Object nil,Lisp_Object a,Lisp_Object b)275 static Lisp_Object Lldb(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
276 {
277     CSL_IGNORE(nil);
278     CSL_IGNORE(a);
279     CSL_IGNORE(b);
280     return aerror("ldb");
281 }
282 
Llcm_n(Lisp_Object nil,int nargs,...)283 Lisp_Object MS_CDECL Llcm_n(Lisp_Object nil, int nargs, ...)
284 {
285     va_list a;
286     int i;
287     Lisp_Object r;
288     if (nargs == 0) return onevalue(fixnum_of_int(1));
289     va_start(a, nargs);
290     push_args(a, nargs);
291     stackcheck0(nargs);
292     pop(r);
293     for (i = 1; i<nargs; i++)
294     {   Lisp_Object w;
295         pop(w);
296         r = lcm(r, w);
297         errexitn(nargs-i-1);
298     }
299     return onevalue(r);
300 }
301 
Llcm(Lisp_Object nil,Lisp_Object a,Lisp_Object b)302 Lisp_Object Llcm(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
303 {
304     a = lcm(a, b);
305     errexit();
306     return onevalue(a);
307 }
308 
Llcm_1(Lisp_Object nil,Lisp_Object a)309 Lisp_Object Llcm_1(Lisp_Object nil, Lisp_Object a)
310 {
311     CSL_IGNORE(nil);
312     return onevalue(a);
313 }
314 
Lldb_test(Lisp_Object nil,Lisp_Object a,Lisp_Object b)315 static Lisp_Object Lldb_test(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
316 {
317     CSL_IGNORE(nil);
318     CSL_IGNORE(a);
319     CSL_IGNORE(b);
320     return aerror("ldb-test");
321 }
322 
Lnumerator(Lisp_Object nil,Lisp_Object a)323 static Lisp_Object Lnumerator(Lisp_Object nil, Lisp_Object a)
324 {
325     CSL_IGNORE(nil);
326     if (!is_number(a)) return aerror1("numerator", a);
327     if (is_numbers(a) && is_ratio(a))
328         return onevalue(numerator(a));
329     else return onevalue(a);
330 }
331 
Lrealpart(Lisp_Object nil,Lisp_Object a)332 static Lisp_Object Lrealpart(Lisp_Object nil, Lisp_Object a)
333 {
334     CSL_IGNORE(nil);
335     if (!is_number(a)) return aerror1("realpart", a);
336     if (is_numbers(a) && is_complex(a))
337         return onevalue(real_part(a));
338     else return onevalue(a);
339 }
340 
341 #else /* COMMON */
342 
Lgcd(Lisp_Object nil,Lisp_Object a,Lisp_Object b)343 Lisp_Object Lgcd(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
344 {
345     a = gcd(a, b);
346     errexit();
347     return onevalue(a);
348 }
349 
Llcm(Lisp_Object nil,Lisp_Object a,Lisp_Object b)350 Lisp_Object Llcm(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
351 {
352     Lisp_Object g;
353     push2(b, a);
354     g = gcd(a, b);
355     errexitn(2);
356     pop(a);
357     a = quot2(a, g);
358     pop(b);
359     errexit();
360     a = times2(a, b);
361     errexit();
362     return onevalue(a);
363 }
364 
365 #endif /* COMMON */
366 
Ldecode_float(Lisp_Object nil,Lisp_Object a)367 static Lisp_Object Ldecode_float(Lisp_Object nil, Lisp_Object a)
368 {
369     double d = float_of_number(a), neg = 1.0;
370     int x;
371     Lisp_Object sign;
372     if (!is_float(a)) return aerror("decode-float");
373     if (d < 0.0) d = -d, neg = -1.0;
374     if (d == 0.0) x = 0;
375     else
376     {   d = frexp(d, &x);
377         if (d == 1.0) d = 0.5, x++;
378     }
379 #ifdef COMMON
380     if (is_sfloat(a)) sign = make_sfloat(neg);
381     else
382 #endif
383         sign = make_boxfloat(neg, type_of_header(flthdr(a)));
384     errexit();
385     push(sign);
386 #ifdef COMMON
387     if (is_sfloat(a)) a = make_sfloat(d);
388     else
389 #endif
390         a = make_boxfloat(d, type_of_header(flthdr(a)));
391     pop(sign);
392     errexit();
393 #ifdef COMMON
394     mv_2 = fixnum_of_int(x);
395     mv_3 = sign;
396     return nvalues(a, 3);
397 #else
398     return list3(sign, fixnum_of_int(x), a);
399 #endif
400 }
401 
402 /*
403  * The next two functions depend on IEEE-format floating point numbers.
404  * They are (thus?) potentially a portability trap, but may suffice for
405  * MOST systems. I need to test the handling of double precision values
406  * on computers that store the two words of a double in each of the
407  * possible orders. If the exponent field in floats was not stored in the
408  * position that would be 0x7f800000 in an integer my treatment of short
409  * floats would fail, so I have already assumed that that is so. I think.
410  */
411 
Lfloat_denormalized_p(Lisp_Object nil,Lisp_Object a)412 static Lisp_Object Lfloat_denormalized_p(Lisp_Object nil, Lisp_Object a)
413 {
414     int x = 0;
415     switch ((int)a & TAG_BITS)
416     {
417 #ifdef COMMON
418 case TAG_SFLOAT:
419         if ((a & 0x7fffffff) == TAG_SFLOAT) return onevalue(nil);  /* 0.0 */
420         x = (int32_t)a & 0x7f800000;
421         return onevalue(x == 0 ? lisp_true : nil);
422 #endif
423 case TAG_BOXFLOAT:
424         switch (type_of_header(flthdr(a)))
425         {
426 #ifdef COMMON
427     case TYPE_SINGLE_FLOAT:
428             if (single_float_val(a) == 0.0) return onevalue(nil);
429             x = ((Single_Float *)((char *)a-TAG_BOXFLOAT))->f.i & 0x7f800000;
430             return onevalue(x == 0 ? lisp_true : nil);
431     case TYPE_LONG_FLOAT:
432             if (long_float_val(a) == 0.0) return onevalue(nil);
433             x = ((int32_t *)long_float_addr(a))[
434                     current_fp_rep&FP_WORD_ORDER ? 0 : 1] & 0x7ff00000;
435             return onevalue(x == 0 ? lisp_true : nil);
436 #endif
437     case TYPE_DOUBLE_FLOAT:
438             if (double_float_val(a) == 0.0) return onevalue(nil);
439             x = ((int32_t *)double_float_addr(a))[
440                     current_fp_rep&FP_WORD_ORDER ? 0 : 1] & 0x7ff00000;
441             return onevalue(x == 0 ? lisp_true : nil);
442         }
443 default:
444         break;
445     }
446     return onevalue(nil);
447 }
448 
Lfloat_infinity_p(Lisp_Object nil,Lisp_Object a)449 static Lisp_Object Lfloat_infinity_p(Lisp_Object nil, Lisp_Object a)
450 {
451 /*
452  * The issue of NaN values is one I do not want to worry about at present.
453  * I deem anything with the largest possibl exponent value to be infinite.
454  */
455     int x = 0;
456     switch ((int)a & TAG_BITS)
457     {
458 #ifdef COMMON
459 case TAG_SFLOAT:
460         x = (int32_t)a & 0x7f800000;
461         return onevalue(x == 0x7f800000 ? lisp_true : nil);
462 #endif
463 case TAG_BOXFLOAT:
464         switch (type_of_header(flthdr(a)))
465         {
466 #ifdef COMMON
467     case TYPE_SINGLE_FLOAT:
468             if (single_float_val(a) == 0.0) return onevalue(nil);
469             x = ((Single_Float *)((char *)a-TAG_BOXFLOAT))->f.i & 0x7f800000;
470             return onevalue(x == 0x7f800000 ? lisp_true : nil);
471     case TYPE_LONG_FLOAT:
472             if (long_float_val(a) == 0.0) return onevalue(nil);
473             x = ((int32_t *)long_float_addr(a))[
474                     current_fp_rep&FP_WORD_ORDER ? 0 : 1] & 0x7ff00000;
475             return onevalue(x == 0x7ff00000 ? lisp_true : nil);
476 #endif
477     case TYPE_DOUBLE_FLOAT:
478             if (double_float_val(a) == 0.0) return onevalue(nil);
479             x = ((int32_t *)double_float_addr(a))[
480                     current_fp_rep&FP_WORD_ORDER ? 0 : 1] & 0x7ff00000;
481             return onevalue(x == 0x7ff00000 ? lisp_true : nil);
482         }
483 default:
484         break;
485     }
486     return onevalue(nil);
487 }
488 
489 /*
490  * The functions such as float-digits, float-precision, float-radix are
491  * coded here assuming that IEEE-style arithmetic is being used. If this is
492  * not so then all the code in this file needs review.  Furthermore I will
493  * be lazy about NaNs and denormalised numbers for now and come back to
494  * worry about them later on if I am really forced to.
495  */
496 
Lfloat_digits(Lisp_Object nil,Lisp_Object a)497 static Lisp_Object Lfloat_digits(Lisp_Object nil, Lisp_Object a)
498 {
499     int tag = (int)a & TAG_BITS;
500     CSL_IGNORE(nil);
501     switch (tag)
502     {
503 #ifdef COMMON
504 case TAG_SFLOAT:
505         return onevalue(fixnum_of_int(20));
506 #endif
507 case TAG_BOXFLOAT:
508         switch (type_of_header(flthdr(a)))
509         {
510 #ifdef COMMON
511     case TYPE_SINGLE_FLOAT:
512             return onevalue(fixnum_of_int(24));
513 #endif
514     default:
515             return onevalue(fixnum_of_int(53));
516         }
517 default:
518         return aerror("float-digits");
519     }
520 }
521 
Lfloat_precision(Lisp_Object nil,Lisp_Object a)522 static Lisp_Object Lfloat_precision(Lisp_Object nil, Lisp_Object a)
523 {
524     int tag = (int)a & TAG_BITS;
525     double d = float_of_number(a);
526     CSL_IGNORE(nil);
527     if (d == 0.0) return onevalue(fixnum_of_int(0));
528 /* /* I do not cope with de-normalised numbers here */
529     switch (tag)
530     {
531 #ifdef COMMON
532 case TAG_SFLOAT:
533         return onevalue(fixnum_of_int(20));
534 #endif
535 case TAG_BOXFLOAT:
536         switch (type_of_header(flthdr(a)))
537         {
538 #ifdef COMMON
539     case TYPE_SINGLE_FLOAT:
540             return onevalue(fixnum_of_int(24));
541 #endif
542     default:
543             return onevalue(fixnum_of_int(53));
544         }
545 default:
546         return aerror("float-precision");
547     }
548 }
549 
Lfloat_radix(Lisp_Object nil,Lisp_Object a)550 static Lisp_Object Lfloat_radix(Lisp_Object nil, Lisp_Object a)
551 {
552     CSL_IGNORE(nil);
553     CSL_IGNORE(a);
554     return onevalue(fixnum_of_int(FLT_RADIX));
555 }
556 
Lfloat_sign2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)557 static Lisp_Object Lfloat_sign2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
558 {
559     double d = float_of_number(b);
560     CSL_IGNORE(nil);
561     if (float_of_number(a) < 0.0) d = -d;
562 #ifdef COMMON
563     if (is_sfloat(b)) return onevalue(make_sfloat(d));
564     else
565 #endif
566     if (!is_bfloat(b)) return aerror1("bad arg for float-sign",  b);
567     else return onevalue(make_boxfloat(d, type_of_header(flthdr(b))));
568 }
569 
Lfloat_sign1(Lisp_Object nil,Lisp_Object a)570 static Lisp_Object Lfloat_sign1(Lisp_Object nil, Lisp_Object a)
571 {
572     double d = float_of_number(a);
573     CSL_IGNORE(nil);
574     if (d < 0.0) d = -1.0; else d = 1.0;
575 #ifdef COMMON
576     if (is_sfloat(a)) return onevalue(make_sfloat(d));
577     else
578 #endif
579     if (!is_bfloat(a)) return aerror1("bad arg for float-sign",  a);
580     else return onevalue(make_boxfloat(d, type_of_header(flthdr(a))));
581 }
582 
Lfround(Lisp_Object nil,Lisp_Object a,Lisp_Object b)583 static Lisp_Object Lfround(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
584 {
585     CSL_IGNORE(nil);
586     CSL_IGNORE(a);
587     CSL_IGNORE(b);
588     return aerror("fround");
589 }
590 
Lftruncate(Lisp_Object nil,Lisp_Object a,Lisp_Object b)591 static Lisp_Object Lftruncate(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
592 {
593     CSL_IGNORE(nil);
594     CSL_IGNORE(a);
595     CSL_IGNORE(b);
596     return aerror("ftruncate");
597 }
598 
Linteger_decode_float(Lisp_Object nil,Lisp_Object a)599 static Lisp_Object Linteger_decode_float(Lisp_Object nil, Lisp_Object a)
600 {
601     double d = float_of_number(a);
602 #ifdef COMMON
603     int tag = (int)a & TAG_BITS;
604 #endif
605     int x, neg = 0;
606     int32_t a1, a2;
607     CSL_IGNORE(nil);
608     if (!is_float(a)) return aerror("integer-decode-float");
609     if (d == 0.0)
610 #ifdef COMMON
611     {   mv_2 = fixnum_of_int(0);
612         mv_3 = fixnum_of_int(d<0 ? -1 : 1);
613         nvalues(fixnum_of_int(0), 3);
614     }
615 #else
616         return list3(fixnum_of_int(0), fixnum_of_int(0),
617                      fixnum_of_int(d<0 ? -1 : 1));
618 #endif
619     if (d < 0.0) d = -d, neg = 1;
620     d = frexp(d, &x);
621     if (d == 1.0) d = 0.5, x++;
622 #ifdef COMMON
623     if (tag == TAG_SFLOAT)
624     {   d *= TWO_20;
625         x -= 20;
626         a1 = (int32_t)d;
627         a = fixnum_of_int(a1);
628     }
629     else if (tag == TAG_BOXFLOAT &&
630              type_of_header(flthdr(a)) == TYPE_SINGLE_FLOAT)
631     {   d *= TWO_24;
632         x -= 24;
633         a1 = (int32_t)d;
634         a = fixnum_of_int(a1);
635     }
636     else
637 #endif
638     {   d *= TWO_22;
639         a1 = (int32_t)d;
640         d -= (double)a1;
641         a2 = (int32_t)(d*TWO_31);  /* This conversion should be exact */
642         x -= 53;
643         a = make_two_word_bignum(a1, a2);
644         errexit();
645     }
646 #ifdef COMMON
647     {   mv_2 = fixnum_of_int(x);
648         mv_3 = neg ? fixnum_of_int(-1) : fixnum_of_int(1);
649         return nvalues(a, 3);
650     }
651 #else
652         return list3(a, fixnum_of_int(x),
653                      neg ? fixnum_of_int(-1) : fixnum_of_int(1));
654 #endif
655 }
656 
Linteger_length(Lisp_Object nil,Lisp_Object a)657 static Lisp_Object Linteger_length(Lisp_Object nil, Lisp_Object a)
658 {
659     a = Labsval(nil, a);
660     errexit();
661     return Lmsd(nil, a);
662 }
663 
Llogbitp(Lisp_Object nil,Lisp_Object a,Lisp_Object b)664 static Lisp_Object Llogbitp(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
665 {
666     CSL_IGNORE(nil);
667     CSL_IGNORE(a);
668     CSL_IGNORE(b);
669     return aerror("logbitp");
670 }
671 
Llogcount(Lisp_Object nil,Lisp_Object a)672 static Lisp_Object Llogcount(Lisp_Object nil, Lisp_Object a)
673 {
674     CSL_IGNORE(nil);
675     CSL_IGNORE(a);
676     return aerror("logcount");
677 }
678 
Llogtest(Lisp_Object nil,Lisp_Object a,Lisp_Object b)679 static Lisp_Object Llogtest(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
680 {
681     CSL_IGNORE(nil);
682     CSL_IGNORE(a);
683     CSL_IGNORE(b);
684     return aerror("logtest");
685 }
686 
Lmask_field(Lisp_Object nil,Lisp_Object a,Lisp_Object b)687 static Lisp_Object Lmask_field(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
688 {
689     CSL_IGNORE(nil);
690     CSL_IGNORE(a);
691     CSL_IGNORE(b);
692     return aerror("mask-field");
693 }
694 
Lscale_float(Lisp_Object nil,Lisp_Object a,Lisp_Object b)695 static Lisp_Object Lscale_float(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
696 {
697     double d = float_of_number(a);
698     CSL_IGNORE(nil);
699     if (!is_fixnum(b)) return aerror("scale-float");
700     d = ldexp(d, int_of_fixnum(b));
701 #ifdef COMMON
702     if (is_sfloat(a)) return onevalue(make_sfloat(d));
703     else
704 #endif
705     if (!is_bfloat(a)) return aerror1("bad arg for scale-float",  a);
706     else return onevalue(make_boxfloat(d, type_of_header(flthdr(a))));
707 }
708 
709 #define FIX_TRUNCATE 0
710 #define FIX_ROUND    1
711 #define FIX_FLOOR    2
712 #define FIX_CEILING  3
713 
lisp_fix_sub(Lisp_Object a,int roundmode)714 static Lisp_Object lisp_fix_sub(Lisp_Object a, int roundmode)
715 /*
716  * This converts from a double to a Lisp integer, which will
717  * quite often have to be a bignum.  No overflow is permitted - the
718  * result can always be accurate.
719  */
720 {
721     int32_t a0, a1, a2, a3;
722     int x, x1;
723     CSLbool negative;
724     double d = float_of_number(a);
725     if (M2_31_1 < d && d < TWO_31)
726     {   int32_t a = (int32_t)d;
727 /*
728  * If my floating point value is in the range -(2^31+1) to +2^31 (exclusive)
729  * then when C truncates it I will get an integer in the inclusive range
730  * from -2^31 to 2^31-1, i.e. the whole range of C 32-bit integers.
731  * If I then apply a rounding mode other than truncation I may just have
732  * overflow, which I have to detect specially.
733  */
734         int32_t w;
735         double f;
736         switch (roundmode)
737         {
738     case FIX_TRUNCATE:
739             break;          /* C casts truncate, so this is OK */
740     case FIX_ROUND:
741             f = d - (double)a;
742             if (f > 0.5)
743             {   if (a == 0x7fffffff) return make_two_word_bignum(1, 0);
744                 else a++;
745             }
746             else if (f < -0.5)
747             {   if (a == ~0x7fffffff)
748                     return make_two_word_bignum(-2, 0x7fffffff);
749                 else a--;
750             }
751    /* The rounding rule on the next lines show what I do in 1/2ulp cases */
752             else if (f == 0.5) a = (a+1) & ~1;
753             else if (f == -0.5)
754             {   if (a == ~0x7fffffff)
755                     return make_two_word_bignum(-2, 0x7fffffff);
756                 else a = a & ~1;
757             }
758             break;
759     case FIX_FLOOR:
760             f = d - (double)a;
761             if (f < 0.0)
762             {   if (a == ~0x7fffffff)
763                     return make_two_word_bignum(-2, 0x7fffffff);
764                 else a--;
765             }
766             break;
767     case FIX_CEILING:
768             f = d - (double)a;
769             if (f > 0.0)
770             {   if (a == 0x7fffffff) return make_two_word_bignum(1, 0);
771                 else a++;
772             }
773             break;
774         }
775         w = a & fix_mask;
776         if (w == 0 || w == fix_mask) return fixnum_of_int(a);
777         else if (!signed_overflow(a)) return make_one_word_bignum(a);
778         else if (a > 0) return make_two_word_bignum(0, a);
779         else return make_two_word_bignum(-1, clear_top_bit(a));
780     }
781     if (d < 0.0) d = -d, negative = YES;
782     else negative = NO;
783     d = frexp(d, &x); /* 0.5 <= abs(d) < 1.0, x = the (binary) exponent */
784     if (d == 1.0) d = 0.5, x++;
785     d *= TWO_31;
786     a1 = (int32_t)d;      /* 2^31 > d >= 2^30 */
787     d -= (double)a1;
788     a2 = (uint32_t)(d*TWO_31);  /* This conversion should be exact */
789     if (negative)
790     {   if (a2 == 0) a1 = -a1;
791         else a2 = clear_top_bit(-a2), a1 = ~a1;
792     }
793     x -= 62;
794     if (x < 0)              /* Need to shift right x places */
795     {   x = -x;             /* The shift amount here can be 31 at most... */
796         a3 = a2 << (32 - x);
797         a2 = clear_top_bit((a2 >> x) | (a1 << (31 - x)));
798 #ifdef SIGNED_SHIFTS_ARE_LOGICAL
799         if (a1 < 0) a1 = (a1 >> x) | (((int32_t)-1) << (31 - x));
800         else a1 = a1 >> x;
801 #else
802         a1 = a1 >> x;
803 #endif
804         switch (roundmode)
805         {
806     case FIX_TRUNCATE:
807             if (a1 < 0 && a3 != 0)  /* proper rounding on -ve values */
808             {   a2++;
809                 if (a2 < 0) /* carry */
810                 {   a2 = 0;
811                     a1++;
812                 }
813             }
814             break;
815     case FIX_ROUND:
816             if ((a3 & 0x80000000U) != 0 &&
817                (a3 != ~0x7fffffff || (a2 & 1) != 0))
818             {   a2++;
819                 if (a2 < 0) /* carry */
820                 {   a2 = 0;
821                     a1++;
822                 }
823             }
824             break;
825     case FIX_FLOOR: /* Comes out in wash of 2's complement */
826             break;
827     case FIX_CEILING:
828             if (a3 != 0)
829             {   a2++;
830                 if (a2 < 0) /* carry */
831                 {   a2 = 0;
832                     a1++;
833                 }
834             }
835             break;
836         }
837         return make_two_word_bignum(a1, a2);
838     }
839 /* Now the double-length value (a1,a2) is correct and exact, and it
840  * needs shifting left by x bits.  This will give a 3 or more word bignum.
841  * Note that no rounding etc is needed at all here since there are no
842  * fractional bits in the representation.
843  */
844     if (a1 < 0)
845     {   a0 = -1;
846         a1 = clear_top_bit(a1);
847     }
848     else a0 = 0;
849     x1 = x / 31;
850     x = x % 31;
851     a0 = (a0 << x) | (a1 >> (31-x));
852     a1 = clear_top_bit(a1 << x) | (a2 >> (31-x));
853     a2 = clear_top_bit(a2 << x);
854     return make_n_word_bignum(a0, a1, a2, x1);
855 }
856 
857 #ifdef COMMON
858 
lisp_fix_ratio(Lisp_Object a,int roundmode)859 static Lisp_Object lisp_fix_ratio(Lisp_Object a, int roundmode)
860 /*
861  * This converts from a ratio to a Lisp integer.  It has to apply
862  * the specified rounding regime.
863  */
864 {
865     Lisp_Object p, q, r, nil = C_nil;;
866     p = numerator(a);
867     q = denominator(a);
868     push2(q, p);
869     r = quot2(p, q);
870     errexitn(2);
871     p = stack[0];
872     stack[0] = r;
873     p = Cremainder(p, stack[-1]);
874     pop2(r, q);
875     errexit();
876     switch (roundmode)
877     {
878 case FIX_TRUNCATE:
879         break;
880 case FIX_ROUND:
881         /* /* This case unfinished at present */
882         break;
883 case FIX_FLOOR:
884         if (minusp(p))
885         {   push(r);
886             p = plus2(p, q);
887             pop(r);
888             errexit();
889             r = sub1(r);
890             errexit();
891         }
892         break;
893 case FIX_CEILING:
894         if (plusp(p))
895         {   push(r);
896             p = difference2(p, q);
897             pop(r);
898             errexit();
899             r = add1(r);
900             errexit();
901         }
902         break;
903     }
904     mv_2 = p;
905     return nvalues(r, 2);
906 }
907 #endif
908 
909 
lisp_fix(Lisp_Object a,int roundmode)910 static Lisp_Object lisp_fix(Lisp_Object a, int roundmode)
911 {
912     Lisp_Object r, nil;
913     push(a);
914     r = lisp_fix_sub(a, roundmode);
915     errexitn(1);
916     a = stack[0];
917     stack[0] = r;
918     a = difference2(a, r);
919     pop(r);
920     errexit();
921     mv_2 = a;
922     return nvalues(r, 2);
923 }
924 
lisp_ifix(Lisp_Object a,Lisp_Object b,int roundmode)925 static Lisp_Object lisp_ifix(Lisp_Object a, Lisp_Object b, int roundmode)
926 {
927     Lisp_Object q, r, nil;
928     if (is_float(a) || is_float(b))
929     {   double p = float_of_number(a), q = float_of_number(b), d = p/q;
930         a = make_boxfloat(d, TYPE_DOUBLE_FLOAT);
931         errexit();
932         r = lisp_fix(a, roundmode);
933         errexit();
934         push(r);
935         a = make_boxfloat(float_of_number(mv_2)*q, TYPE_DOUBLE_FLOAT);
936         pop(r);
937         errexit();
938         mv_2 = a;
939         return nvalues(r, 2);
940     }
941     push2(a, b);
942     q = quot2(a, b);
943     errexitn(2);
944     a = stack[-1];
945     b = stack[0];
946     push(q);
947     r = Cremainder(a, b);
948     errexitn(3);
949     switch (roundmode)
950     {
951 case FIX_TRUNCATE:
952         break;
953 case FIX_ROUND:
954         popv(3); return aerror("two-arg ROUND");
955 case FIX_FLOOR:
956         if (!minusp(r)) break;
957         r = plus2(r, stack[-1]);
958         errexitn(3);
959         q = stack[0];
960         push(r);
961         q = sub1(q);
962         pop(r);
963         errexitn(3);
964         stack[0] = q;
965         break;
966 case FIX_CEILING:
967         if (!plusp(r)) break;
968         r = difference2(r, stack[-1]);
969         errexitn(3);
970         q = stack[0];
971         push(r);
972         q = add1(q);
973         pop(r);
974         errexitn(3);
975         stack[0] = q;
976         break;
977     }
978     pop3(q, b, a);
979     mv_2 = r;
980     return nvalues(q, 2);
981 }
982 
983 /*
984  * So far I have not implemented support for rational numbers in the 2-arg
985  * versions of these functions. /*
986  */
987 
Lceiling_2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)988 static Lisp_Object Lceiling_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
989 {
990     CSL_IGNORE(nil);
991     if (!is_number(a) || !is_number(b)) return aerror1("ceiling", a);
992     return lisp_ifix(a, b, FIX_CEILING);
993 }
994 
Lfloor_2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)995 static Lisp_Object Lfloor_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
996 {
997     CSL_IGNORE(nil);
998     if (!is_number(a) || !is_number(b)) return aerror1("floor", a);
999     return lisp_ifix(a, b, FIX_FLOOR);
1000 }
1001 
Lround_2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)1002 static Lisp_Object Lround_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
1003 {
1004     CSL_IGNORE(nil);
1005     if (!is_number(a) || !is_number(b)) return aerror1("round", a);
1006     return lisp_ifix(a, b, FIX_ROUND);
1007 }
1008 
Ltruncate_2(Lisp_Object nil,Lisp_Object a,Lisp_Object b)1009 Lisp_Object Ltruncate_2(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
1010 {
1011     CSL_IGNORE(nil);
1012     if (!is_number(a) || !is_number(b)) return aerror1("truncate", a);
1013     return lisp_ifix(a, b, FIX_TRUNCATE);
1014 }
1015 
Lceiling(Lisp_Object nil,Lisp_Object a)1016 static Lisp_Object Lceiling(Lisp_Object nil, Lisp_Object a)
1017 {
1018     CSL_IGNORE(nil);
1019     if (!is_number(a)) return aerror1("ceiling", a);
1020 #ifdef COMMON
1021     if (is_numbers(a) && is_ratio(a)) return lisp_fix_ratio(a, FIX_CEILING);
1022 #endif
1023     if (!is_float(a)) { mv_2 = fixnum_of_int(0); return nvalues(a, 2); }
1024     return lisp_fix(a, FIX_CEILING);
1025 }
1026 
Lfloor(Lisp_Object nil,Lisp_Object a)1027 static Lisp_Object Lfloor(Lisp_Object nil, Lisp_Object a)
1028 {
1029     CSL_IGNORE(nil);
1030     if (!is_number(a)) return aerror1("floor", a);
1031 #ifdef COMMON
1032     if (is_numbers(a) && is_ratio(a)) return lisp_fix_ratio(a, FIX_FLOOR);
1033 #endif
1034     if (!is_float(a)) { mv_2 = fixnum_of_int(0); return nvalues(a, 2); }
1035     return lisp_fix(a, FIX_FLOOR);
1036 }
1037 
Lround(Lisp_Object nil,Lisp_Object a)1038 static Lisp_Object Lround(Lisp_Object nil, Lisp_Object a)
1039 {
1040     CSL_IGNORE(nil);
1041     if (!is_number(a)) return aerror1("round", a);
1042 #ifdef COMMON
1043     if (is_numbers(a) && is_ratio(a)) return lisp_fix_ratio(a, FIX_ROUND);
1044 #endif
1045     if (!is_float(a)) { mv_2 = fixnum_of_int(0); return nvalues(a, 2); }
1046     return lisp_fix(a, FIX_ROUND);
1047 }
1048 
Ltruncate(Lisp_Object nil,Lisp_Object a)1049 Lisp_Object Ltruncate(Lisp_Object nil, Lisp_Object a)
1050 {
1051     CSL_IGNORE(nil);
1052     if (!is_number(a)) return aerror1("fix", a);
1053 #ifdef COMMON
1054     if (is_numbers(a) && is_ratio(a)) return lisp_fix_ratio(a, FIX_TRUNCATE);
1055 #endif
1056     if (!is_float(a)) { mv_2 = fixnum_of_int(0); return nvalues(a, 2); }
1057     return lisp_fix(a, FIX_TRUNCATE);
1058 }
1059 
1060 /*
1061  * The following macro is expected to select out the 32-bit word within a
1062  * floating point number that has the exponent field packed within it.
1063  * On IEEE-format machines I expect to find the exponent in the
1064  * 0x7ff00000 bits within this word. It sets i to the top half of the
1065  * floating value d, and worries somewhat about GCCs strict aliasing
1066  * optimisation!! Specifically it avoids the use of pointers in the
1067  * punning that has to go on.
1068  */
1069 
1070 #define top_half(ii, d)                          \
1071   { union { double f; int32_t i[2]; } w;         \
1072     w.f = d;                                     \
1073     ii = w.i[(~current_fp_rep) & FP_WORD_ORDER]; \
1074   }
1075 
1076 /*
1077  *  symbolic procedure safe!-fp!-plus(x,y);
1078  *     if zerop x then y else if zerop y then x else
1079  *     begin scalar u;
1080  *        if x>0.0 and y>0.0 then
1081  *           if x<!!plumax and y<!!plumax then go to ret else return nil;
1082  *        if x<0.0 and y<0.0 then
1083  *           if -x<!!plumax and -y<!!plumax then go to ret else return nil;
1084  *        if abs x<!!plumin and abs y<!!plumin then return nil;
1085  *   ret: return
1086  *          if (u := plus2(x,y))=0.0
1087  *             or x>0.0 and y>0.0 or x<0.0 and y<0.0 then u
1088  *          else if abs u<(abs x)*!!fleps1 then 0.0 else u end;
1089  */
1090 
Lsafe_fp_plus(Lisp_Object env,Lisp_Object a,Lisp_Object b)1091 static Lisp_Object Lsafe_fp_plus(Lisp_Object env, Lisp_Object a, Lisp_Object b)
1092 /*
1093  * safe!-fp!-plus must always be given floating point arguments.  In
1094  * most reasonable cases it just returns the floating point sum.  If
1095  * there was some chance that the sum might either overflow or underflow
1096  * then the value NIL is returned instead.  The exact place where this
1097  * function starts to report overflow is not precisely defined - all that
1098  * is guaranteed is that if a result is returned then it will be of full
1099  * precision.
1100  */
1101 {
1102     Lisp_Object nil = C_nil;
1103     double aa, bb, cc;
1104     int32_t ah, bh;
1105     if (!is_float(a) || !is_float(b)) return aerror2("safe-fp-plus", a, b);
1106 /*
1107  * I accept here that adding 0.0 to anything is always possible without
1108  * problem.
1109  */
1110     if ((aa = double_float_val(a)) == 0.0) return onevalue(b);
1111     if ((bb = double_float_val(b)) == 0.0) return onevalue(a);
1112     top_half(ah, aa);
1113     top_half(bh, bb);
1114 /*
1115  * Here I am going to assume IEEE-format numbers. I hope that I have
1116  * picked out the word containing the exponent and that it is positioned in
1117  * the word where I expect.
1118  */
1119 /*
1120  * I deem overflow POSSIBLE if both numbers have the same sign and if at
1121  * least one of them has an exponent 0x7fe (i.e. the highest exponent that
1122  * a non-infinite number can have).
1123  */
1124     if (ah >= 0)
1125     {   if (bh >= 0)
1126         {   if (ah >= 0x7fe00000 || bh >= 0x7fe00000) return onevalue(nil);
1127             cc = aa + bb;
1128         }
1129         else
1130 /*
1131  * I deem underflow POSSIBLE if the numbers have opposite signs and BOTH
1132  * of them have exponents less than 0x035 (which is 53 in decimal, and
1133  * IEEE-format numbers have 53-bit mantissas.  The critical case would
1134  * be the subtraction
1135  *   (53) 1 00000000000000000000 00000000000000000000000000000001
1136  * - (53) 1 00000000000000000000 00000000000000000000000000000000
1137  * where you see the LSB (which is all that is left after the subtraction)
1138  * has to be shifted left 52 bits to get to the position of the implied
1139  * leading 1 bit in the mantissa.  As that happens 52 gets subtracted
1140  * from the exponent, leaving it as 1, the smallest exponent for a normalised
1141  * non-zero value.
1142  */
1143         {   double eps, ra;
1144             bh &= 0x7fffffff;
1145             if (ah < 0x03500000 && bh < 0x03500000) return onevalue(nil);
1146 /*
1147  * The next few lines check to see if addition has led to cancellation of
1148  * leading digits to an extent greater than !!fleps1.  The environent cell
1149  * of safe!-fp!-plus must be set to !!fleps, and then the value cell of
1150  * this symbol is accessed to obtain the limit.
1151  */
1152             eps = 0.0;
1153             if (is_symbol(env))
1154             {   Lisp_Object ve = qvalue(env);
1155                 if (is_float(ve)) eps = double_float_val(ve);
1156             }
1157             cc = aa + bb;
1158             ra = cc/aa;
1159             if (ra < 0.0) ra = -ra;
1160             if (ra < eps) cc = 0.0; /* Force to zero if too small */
1161         }
1162     }
1163     else
1164     {   if (bh >= 0)
1165         {   double eps, ra;
1166             ah &= 0x7fffffff;
1167             if (ah < 0x03500000 && bh < 0x03500000) return onevalue(nil);
1168             eps = 0.0;
1169             if (is_symbol(env))
1170             {   Lisp_Object ve = qvalue(env);
1171                 if (is_float(ve)) eps = double_float_val(ve);
1172             }
1173             cc = aa + bb;
1174             ra = cc/aa;
1175             if (ra < 0.0) ra = -ra;
1176             if (ra < eps) cc = 0.0;
1177         }
1178         else
1179         {   ah &= 0x7fffffff;
1180             bh &= 0x7fffffff;
1181             if (ah >= 0x7fe00000 || bh >= 0x7fe00000) return onevalue(nil);
1182             cc = aa + bb;
1183         }
1184     }
1185     a = make_boxfloat(cc, TYPE_DOUBLE_FLOAT);
1186     errexit();
1187     return onevalue(a);
1188 }
1189 
1190 /*
1191  *  symbolic procedure safe!-fp!-times(x,y);
1192  *   if zerop x or zerop y then 0.0
1193  *   else if x=1.0 then y else if y=1.0 then x else
1194  *     begin scalar u,v; u := abs x; v := abs y;
1195  *        if u>=1.0 and u<=!!timmax then
1196  *           if v<=!!timmax then go to ret else return nil;
1197  *        if u>!!timmax then if v<=1.0 then go to ret else return nil;
1198  *        if u<1.0 and u>=!!timmin then
1199  *           if v>=!!timmin then go to ret else return nil;
1200  *        if u<!!timmin and v<1.0 then return nil;
1201  *   ret: return times2(x,y) end;
1202  */
1203 
Lsafe_fp_times(Lisp_Object nil,Lisp_Object a,Lisp_Object b)1204 static Lisp_Object Lsafe_fp_times(Lisp_Object nil,
1205                                   Lisp_Object a, Lisp_Object b)
1206 /*
1207  * safe!-fp!-times must always be given floating point arguments.  In
1208  * most reasonable cases it just returns the floating point product.  If
1209  * there was some chance that the product might either overflow or underflow
1210  * then the value NIL is returned instead.  REDUCE requires that if this
1211  * function returns a non-overflow result than two such returned values
1212  * can be added witjout fear of overflow, as in a*b+c*d.  Hence I ought to
1213  * report trouble slightly earlier than I might otherwise. - but REDUCE is
1214  * being changed to remove this oddity, and it seems it could only cause
1215  * trouble in (1.0,max)*(max, 1.0) complex multiplication, so I am not
1216  * going to worry...
1217  */
1218 {
1219     double aa, bb, cc;
1220     int32_t ah, bh;
1221     if (!is_float(a) || !is_float(b)) return aerror2("safe-fp-times", a, b);
1222 /*
1223  * Multiplication by zero is handled as a special case here - doing so
1224  * means that I do not have to worry about the special case of zero
1225  * exponents later on, and it also avoids allocating fresh space to hold
1226  * a new floating point zero value.
1227  */
1228     if ((aa = double_float_val(a)) == 0.0) return onevalue(a);
1229     if ((bb = double_float_val(b)) == 0.0) return onevalue(b);
1230     top_half(ah, aa);
1231     top_half(bh, bb);
1232 /*
1233  * Here I am going to assume IEEE-format numbers. I hope that I have
1234  * picked out the word containing the exponent and that it is positioned in
1235  * the word where I expect.
1236  */
1237     ah = (ah >> 20) & 0x7ff;
1238     bh = (bh >> 20) & 0x7ff;
1239     ah = ah + bh;
1240 /*
1241  * I can estimate the value of the product by adding the exponents of the
1242  * two operands.
1243  */
1244     if (ah < 0x400 || ah >= 0xbfd) return onevalue(nil);
1245     cc = aa * bb;
1246     a = make_boxfloat(cc, TYPE_DOUBLE_FLOAT);
1247     errexit();
1248     return onevalue(a);
1249 }
1250 
1251 /*
1252  *  symbolic procedure safe!-fp!-quot(x,y);
1253  *   if zerop y then rdqoterr()
1254  *   else if zerop x then 0.0 else if y=1.0 then x else
1255  *     begin scalar u,v; u := abs x; v := abs y;
1256  *        if u>=1.0 and u<=!!timmax then
1257  *           if v>=!!timmin then go to ret else return nil;
1258  *        if u>!!timmax then if v>=1.0 then go to ret else return nil;
1259  *        if u<1.0 and u>=!!timmin then
1260  *           if v<=!!timmax then go to ret else return nil;
1261  *        if u<!!timmin and v>1.0 then return nil;
1262  *   ret: return quotient(x,y) end;
1263  */
1264 
Lsafe_fp_quot(Lisp_Object nil,Lisp_Object a,Lisp_Object b)1265 static Lisp_Object Lsafe_fp_quot(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
1266 {
1267     double aa, bb, cc;
1268     int32_t ah, bh;
1269     if (!is_float(a) || !is_float(b)) return aerror2("safe-fp-quot", a, b);
1270     if ((aa = double_float_val(a)) == 0.0) return onevalue(a);
1271 /*
1272  * I pass division by zero back to the general case here.
1273  */
1274     if ((bb = double_float_val(b)) == 0.0) return onevalue(nil);
1275     top_half(ah, aa);
1276     top_half(bh, bb);
1277     ah = (ah >> 20) & 0x7ff;
1278     bh = (bh >> 20) & 0x7ff;
1279     ah = ah - bh;
1280     if (ah <= -0x3fe || ah >= 0x3fe) return onevalue(nil);
1281     cc = aa / bb;
1282     a = make_boxfloat(cc, TYPE_DOUBLE_FLOAT);
1283     errexit();
1284     return onevalue(a);
1285 }
1286 
1287 /*
1288  *  symbolic procedure safe!-fp!-pl(x,y);
1289  *    % floating plus protect from under- and over-flows but no zero
1290  *    % result check.
1291  *     if zerop x then y else if zerop y then x
1292  *     else if x>0 and y>0 then
1293  *        if x<!!plumax and y<!!plumax then plus2(x,y) else nil
1294  *     else if x<0 and y<0 then
1295  *        if (-x<!!plumax and -y<!!plumax) then plus2(x,y) else nil
1296  *     else if abs x<!!plumin or abs y<!!plumin then nil else plus2(x,y);
1297  */
1298 
Lsafe_fp_pl(Lisp_Object nil,Lisp_Object a,Lisp_Object b)1299 static Lisp_Object Lsafe_fp_pl(Lisp_Object nil, Lisp_Object a, Lisp_Object b)
1300 {
1301     double aa, bb, cc;
1302     int32_t ah, bh;
1303     if (!is_float(a) || !is_float(b)) return aerror2("safe-fp-pl", a, b);
1304     if ((aa = double_float_val(a)) == 0.0) return onevalue(b);
1305     if ((bb = double_float_val(b)) == 0.0) return onevalue(a);
1306     top_half(ah, aa);
1307     top_half(bh, bb);
1308     if (ah >= 0)
1309     {   if (bh >= 0)
1310         {   if (ah >= 0x7fe00000 || bh >= 0x7fe00000) return onevalue(nil);
1311         }
1312         else
1313         {   bh &= 0x7fffffff;
1314             if (ah < 0x03500000 && bh < 0x03500000) return onevalue(nil);
1315         }
1316     }
1317     else
1318     {   if (bh >= 0)
1319         {   ah &= 0x7fffffff;
1320             if (ah < 0x03500000 && bh < 0x03500000) return onevalue(nil);
1321         }
1322         else
1323         {   ah &= 0x7fffffff;
1324             bh &= 0x7fffffff;
1325             if (ah >= 0x7fe00000 || bh >= 0x7fe00000) return onevalue(nil);
1326         }
1327     }
1328     cc = aa + bb;
1329     a = make_boxfloat(cc, TYPE_DOUBLE_FLOAT);
1330     errexit();
1331     return onevalue(a);
1332 }
1333 
1334 /*
1335  *  symbolic procedure safe!-fp!-pl0(x,y);
1336  *    % protects floating plus against under-flow only.
1337  *     if zerop x then y else if zerop y then x
1338  *     else if abs x<!!plumin and abs y<!!plumin then nil else plus2(x,y);
1339  *
1340  * implement as safe!-fp!-pl without MUCH loss of speed.
1341  */
1342 
Llose_precision(Lisp_Object nil,Lisp_Object a,Lisp_Object n)1343 static Lisp_Object Llose_precision(Lisp_Object nil,
1344                                    Lisp_Object a, Lisp_Object n)
1345 {
1346     double aa;
1347     char buffer[64];
1348     int32_t nn;
1349     if (!is_float(a) || !is_fixnum(n)) return aerror2("lose_precision", a, n);
1350     nn = int_of_fixnum(n);
1351     if (nn <= 0 || nn >= 20) return onevalue(a);
1352     sprintf(buffer, "%.*g", (int)nn, double_float_val(a));
1353     if (sscanf(buffer, "%lg", &aa) != 1) return aerror("lose-precision");
1354     a = make_boxfloat(aa, TYPE_DOUBLE_FLOAT);
1355     errexit();
1356     return onevalue(a);
1357 }
1358 
1359 setup_type const arith08_setup[] =
1360 {
1361 /*
1362  * The next few are JUST for REDUCE, but they are expected to speed up
1363  * rounded-mode arithmetic rather a lot.
1364  */
1365     {"lose-precision",          too_few_2, Llose_precision, wrong_no_2},
1366     {"safe-fp-plus",            too_few_2, Lsafe_fp_plus, wrong_no_2},
1367     {"safe-fp-times",           too_few_2, Lsafe_fp_times, wrong_no_2},
1368     {"safe-fp-quot",            too_few_2, Lsafe_fp_quot, wrong_no_2},
1369     {"safe-fp-pl",              too_few_2, Lsafe_fp_pl, wrong_no_2},
1370     {"safe-fp-pl0",             too_few_2, Lsafe_fp_pl, wrong_no_2},
1371 /* End of REDUCE specialities */
1372     {"ceiling",                 Lceiling, Lceiling_2, wrong_no_1},
1373     {"floor",                   Lfloor, Lfloor_2, wrong_no_1},
1374     {"round",                   Lround, Lround_2, wrong_no_1},
1375     {"truncate",                Ltruncate, Ltruncate_2, wrong_no_1},
1376     {"decode-float",            Ldecode_float, too_many_1, wrong_no_1},
1377     {"float-denormalized-p",    Lfloat_denormalized_p, too_many_1, wrong_no_1},
1378     {"float-infinity-p",        Lfloat_infinity_p, too_many_1, wrong_no_1},
1379     {"integer-decode-float",    Linteger_decode_float, too_many_1, wrong_no_1},
1380     {"integer-length",          Linteger_length, too_many_1, wrong_no_1},
1381     {"float-digits",            Lfloat_digits, too_many_1, wrong_no_1},
1382     {"float-precision",         Lfloat_precision, too_many_1, wrong_no_1},
1383     {"float-radix",             Lfloat_radix, too_many_1, wrong_no_1},
1384     {"float-sign",              Lfloat_sign1, Lfloat_sign2, wrong_no_2},
1385     {"fround",                  too_few_2, Lfround, wrong_no_2},
1386     {"ftruncate",               too_few_2, Lftruncate, wrong_no_2},
1387     {"logbitp",                 too_few_2, Llogbitp, wrong_no_2},
1388     {"logcount",                Llogcount, too_many_1, wrong_no_1},
1389     {"logtest",                 too_few_2, Llogtest, wrong_no_2},
1390     {"mask-field",              too_few_2, Lmask_field, wrong_no_2},
1391     {"scale-float",             too_few_2, Lscale_float, wrong_no_2},
1392 #ifdef COMMON
1393     {"boole",                   wrong_no_na, wrong_no_nb, Lboole},
1394     {"byte",                    too_few_2, Lbyte, wrong_no_2},
1395     {"byte-position",           Lbyte_position, too_many_1, wrong_no_1},
1396     {"byte-size",               Lbyte_size, too_many_1, wrong_no_1},
1397     {"complex",                 Lcomplex_1, Lcomplex_2, wrong_no_2},
1398     {"conjugate",               Lconjugate, too_many_1, wrong_no_1},
1399     {"decode-float",            Ldecode_float, too_many_1, wrong_no_1},
1400     {"float-denormalized-p",    Lfloat_denormalized_p, too_many_1, wrong_no_1},
1401     {"float-infinity-p",        Lfloat_infinity_p, too_many_1, wrong_no_1},
1402     {"denominator",             Ldenominator, too_many_1, wrong_no_1},
1403     {"deposit-field",           wrong_no_na, wrong_no_nb, Ldeposit_field},
1404     {"dpb",                     wrong_no_na, wrong_no_nb, Ldpb},
1405     {"ffloor",                  too_few_2, Lffloor, wrong_no_2},
1406     {"gcd",                     Lgcd_1, Lgcd, Lgcd_n},
1407     {"imagpart",                Limagpart, too_many_1, wrong_no_1},
1408     {"ldb",                     too_few_2, Lldb, wrong_no_2},
1409     {"ldb-test",                too_few_2, Lldb_test, wrong_no_2},
1410     {"lcm",                     Llcm_1, Llcm, Llcm_n},
1411     {"numerator",               Lnumerator, too_many_1, wrong_no_1},
1412     {"realpart",                Lrealpart, too_many_1, wrong_no_1},
1413 #else
1414     {"fix",                     Ltruncate, too_many_1, wrong_no_1},
1415     {"gcdn",                    too_few_2, Lgcd, wrong_no_2},
1416     {"lcmn",                    too_few_2, Llcm, wrong_no_2},
1417 #endif
1418     {NULL,                      0,0,0}
1419 };
1420 
1421 /* end of arith08.c */
1422 
1423