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