1 /* Part of SWI-Prolog
2
3 Author: Jan Wielemaker
4 E-mail: J.Wielemaker@vu.nl
5 WWW: http://www.swi-prolog.org
6 Copyright (c) 1985-2020, University of Amsterdam
7 VU University Amsterdam
8 CWI, Amsterdam
9 All rights reserved.
10
11 Redistribution and use in source and binary forms, with or without
12 modification, are permitted provided that the following conditions
13 are met:
14
15 1. Redistributions of source code must retain the above copyright
16 notice, this list of conditions and the following disclaimer.
17
18 2. Redistributions in binary form must reproduce the above copyright
19 notice, this list of conditions and the following disclaimer in
20 the documentation and/or other materials provided with the
21 distribution.
22
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
31 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
32 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
33 ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34 POSSIBILITY OF SUCH DAMAGE.
35 */
36
37 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38 The arithmetic module defines a small set of logical integer predicates
39 as well as the evaluation of arbitrary arithmetic expressions.
40 Arithmetic can be interpreted or compiled (see -O flag). Interpreted
41 arithmetic is supported by the built-in predicates is/2, >/2, etc.
42 These functions call valueExpression() to evaluate a Prolog term holding
43 an arithmetic expression.
44
45 For compiled arithmetic, the compiler generates WAM codes that execute a
46 stack machine. This module maintains an array of arithmetic functions.
47 These functions are addressed by the WAM instructions using their index
48 in this array.
49 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50
51 #ifdef __MINGW32__
52 #include <winsock2.h>
53 #include <windows.h>
54 #include <wincrypt.h>
55 #endif
56
57 /*#define O_DEBUG 1*/
58 #include "pl-incl.h"
59 #include "pl-arith.h"
60 #include <math.h>
61 #include <limits.h>
62 #ifdef HAVE_FLOAT_H
63 #include <float.h>
64 #ifdef _MSC_VER
65 #ifndef isnan
66 #define isnan(x) _isnan(x)
67 #endif
68 #define copysign(x,y) _copysign(x,y)
69 #endif
70 #endif
71 #ifdef HAVE_IEEEFP_H
72 #include <ieeefp.h>
73 #endif
74 #include <fenv.h>
75
76 #ifndef DBL_MAX
77 #define DBL_MAX 1.7976931348623157e+308
78 #endif
79 #ifndef DBL_MIN
80 #define DBL_MIN 2.2250738585072014e-308
81 #endif
82 #ifndef DBL_EPSILON
83 #define DBL_EPSILON 0.00000000000000022204
84 #endif
85
86
87 #ifdef fpclassify
88 #define HAVE_FPCLASSIFY 1
89 #endif
90
91 #undef LD
92 #define LD LOCAL_LD
93
94 #ifndef M_PI
95 #define M_PI (3.14159265358979323846)
96 #endif
97 #ifndef M_E
98 #define M_E (2.7182818284590452354)
99 #endif
100
101 static double const_nan;
102 static double const_inf;
103 static double const_neg_inf;
104
105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 On some machines, notably FreeBSD upto version 3.x, floating point
107 operations raise signals rather then leaving an error condition and this
108 behaviour can be changed to be IEEE754 using fpsetmask() and friends.
109 Here we test whether this interface is present and set it up
110 accordingly.
111
112 With many thanks to NIDE Naoyuki for the clear explanation of the
113 problem.
114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115
116 #if defined(HAVE_FLOATINGPOINT_H) && defined(HAVE_FPSETMASK) && defined(HAVE_FPRESETSTICKY)
117 #define O_INHIBIT_FP_SIGNALS
118 #include <floatingpoint.h>
119 #ifndef FP_X_DZ
120 #define FP_X_DZ 0
121 #endif
122 #ifndef FP_X_INV
123 #define FP_X_INV 0
124 #endif
125 #ifndef FP_X_OFL
126 #define FP_X_OFL 0
127 #endif
128 #endif
129
130 static int ar_minus(Number n1, Number n2, Number r);
131 static int mul64(int64_t x, int64_t y, int64_t *r);
132 static int notLessThanZero(const char *f, int a, Number n);
133 static int mustBePositive(const char *f, int a, Number n);
134 static int set_roundtoward(Word p, Number old ARG_LD);
135
136
137 /********************************
138 * LOGICAL INTEGER FUNCTIONS *
139 *********************************/
140
141 static inline void
clearInteger(Number n)142 clearInteger(Number n)
143 {
144 #ifdef O_GMP
145 if ( n->type == V_MPZ && n->value.mpz->_mp_alloc )
146 mpz_clear(n->value.mpz);
147 #endif
148 }
149
150
151 typedef struct between_state
152 { number low;
153 number high;
154 int hinf;
155 } between_state;
156
157
158 static
159 PRED_IMPL("between", 3, between, PL_FA_NONDETERMINISTIC)
160 { PRED_LD
161 between_state *state;
162 term_t low = A1;
163 term_t high = A2;
164 term_t n = A3;
165 int rc = TRUE;
166
167 switch( CTX_CNTRL )
168 { case FRG_FIRST_CALL:
169 { number l, h, i;
170 int hinf = FALSE;
171
172 if ( !PL_get_number(low, &l) || !intNumber(&l) )
173 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, low);
174 if ( !PL_get_number(high, &h) || !intNumber(&h) )
175 { if ( PL_is_inf(high) )
176 { h.type = V_INTEGER; /* make clearInteger() safe */
177 hinf = TRUE;
178 } else
179 { return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, high);
180 }
181 }
182
183 /* between(+,+,+) */
184 if ( PL_get_number(n, &i) && intNumber(&i) )
185 { int rc;
186
187 if ( hinf )
188 { rc = cmpNumbers(&i, &l) >= 0;
189 } else
190 { rc = cmpNumbers(&i, &l) >= 0 && cmpNumbers(&i, &h) <= 0;
191 }
192
193 clearInteger(&l);
194 clearInteger(&i);
195 if ( !hinf )
196 clearInteger(&h);
197
198 return rc;
199 }
200
201 /* between(+,+,-) */
202 if ( !PL_is_variable(n) )
203 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, n);
204 if ( hinf == FALSE && cmpNumbers(&h, &l) < 0 )
205 { clearInteger(&l);
206 clearInteger(&h);
207 fail;
208 }
209
210 if ( !PL_unify(n, low) )
211 fail;
212 if ( hinf == FALSE && cmpNumbers(&l, &h) == 0 )
213 { clearInteger(&l);
214 clearInteger(&h);
215 succeed;
216 }
217
218 state = allocForeignState(sizeof(*state));
219 cpNumber(&state->low, &l);
220 cpNumber(&state->high, &h);
221 state->hinf = hinf;
222 clearInteger(&l);
223 clearInteger(&h);
224 ForeignRedoPtr(state);
225 /*NOTREACHED*/
226 }
227 case FRG_REDO:
228 { state = CTX_PTR;
229
230 if ( !ar_add_ui(&state->low, 1) ||
231 !PL_unify_number(n, &state->low) )
232 { rc = FALSE;
233 goto cleanup;
234 }
235 if ( !state->hinf &&
236 cmpNumbers(&state->low, &state->high) == 0 )
237 goto cleanup;
238 ForeignRedoPtr(state);
239 /*NOTREACHED*/
240 }
241 case FRG_CUTTED:
242 { state = CTX_PTR;
243 cleanup:
244 clearInteger(&state->low);
245 clearInteger(&state->high);
246 freeForeignState(state, sizeof(*state));
247 /*FALLTHROUGH*/
248 }
249 default:;
250 return rc;
251 }
252 }
253
254 static
255 PRED_IMPL("succ", 2, succ, 0)
256 { GET_LD
257 Word p1, p2;
258 number i1, i2, one;
259 int rc;
260
261 p1 = valTermRef(A1); deRef(p1);
262
263 one.type = V_INTEGER;
264 one.value.i = 1;
265
266 if ( isInteger(*p1) )
267 { get_integer(*p1, &i1);
268 if ( ar_sign_i(&i1) < 0 )
269 return PL_error(NULL, 0, NULL, ERR_DOMAIN,
270 ATOM_not_less_than_zero, A1);
271 rc = ( pl_ar_add(&i1, &one, &i2) &&
272 PL_unify_number(A2, &i2)
273 );
274 } else if ( !canBind(*p1) )
275 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, A1);
276
277 p2 = valTermRef(A2); deRef(p2);
278
279 if ( isInteger(*p2) )
280 { get_integer(*p2, &i2);
281 switch( ar_sign_i(&i2) )
282 { case 1:
283 rc = ( ar_minus(&i2, &one, &i1) &&
284 PL_unify_number(A1, &i1)
285 );
286 break;
287 case 0:
288 fail;
289 case -1:
290 default:
291 return PL_error(NULL, 0, NULL, ERR_DOMAIN,
292 ATOM_not_less_than_zero, A2);
293 }
294 } else if ( !canBind(*p2) )
295 { return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, A2);
296 } else
297 return PL_error(NULL, 0, NULL, ERR_INSTANTIATION);
298
299 clearInteger(&i1);
300 clearInteger(&i2);
301 clearInteger(&one);
302
303 return rc;
304 }
305
306
307 static int
var_or_integer(term_t t,number * n,int which,int * mask ARG_LD)308 var_or_integer(term_t t, number *n, int which, int *mask ARG_LD)
309 { Word p = valTermRef(t);
310
311 deRef(p);
312 if ( isInteger(*p) )
313 { get_integer(*p, n);
314 *mask |= which;
315 succeed;
316 }
317 if ( canBind(*p) )
318 succeed;
319
320 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer, t);
321 }
322
323
324 static
325 PRED_IMPL("plus", 3, plus, 0)
326 { GET_LD
327 number m, n, o;
328 int mask = 0;
329 int rc;
330
331 if ( !var_or_integer(A1, &m, 0x1, &mask PASS_LD) ||
332 !var_or_integer(A2, &n, 0x2, &mask PASS_LD) ||
333 !var_or_integer(A3, &o, 0x4, &mask PASS_LD) )
334 fail;
335
336 switch(mask)
337 { case 0x7: /* +, +, + */
338 case 0x3: /* +, +, - */
339 pl_ar_add(&m, &n, &o);
340 rc = PL_unify_number(A3, &o);
341 break;
342 case 0x5: /* +, -, + */
343 ar_minus(&o, &m, &n);
344 rc = PL_unify_number(A2, &n);
345 break;
346 case 0x6: /* -, +, + */
347 ar_minus(&o, &n, &m);
348 rc = PL_unify_number(A1, &m);
349 break;
350 default:
351 return PL_error(NULL, 0, NULL, ERR_INSTANTIATION);
352 }
353
354 clearInteger(&m);
355 clearInteger(&n);
356 clearInteger(&o);
357
358 return rc;
359 }
360
361 /********************************
362 * LOGICAL NUMBER FUNCTION *
363 *********************************/
364
365 static
366 PRED_IMPL("bounded_number", 3, bounded_number, 0)
367 { PRED_LD
368 number n, lo, hi;
369 int rc;
370
371 if ( PL_get_number(A3, &n) )
372 { switch(n.type)
373 {
374 #ifdef O_GMP
375 case V_MPZ:
376 #endif
377 case V_INTEGER:
378 { cpNumber(&lo, &n);
379 cpNumber(&hi, &n);
380 ar_add_ui(&lo, -1);
381 ar_add_ui(&hi, 1);
382 break;
383 }
384 #if O_GMP
385 case V_MPQ:
386 promoteToFloatNumber(&n);
387 /*FALLTHROUGH*/
388 #endif
389 case V_FLOAT:
390 { if ( isfinite(n.value.f) )
391 { lo.type = V_FLOAT;
392 lo.value.f = nexttoward(n.value.f,-INFINITY);
393 hi.type = V_FLOAT;
394 hi.value.f = nexttoward(n.value.f, INFINITY);
395 } else
396 { clearNumber(&n);
397 return FALSE;
398 }
399 break;
400 }
401 }
402
403 rc = ( ((PL_get_number(A1, &lo)) ? (cmpNumbers(&lo, &n) == -1)
404 : PL_unify_number(A1, &lo)) &&
405 ((PL_get_number(A2, &hi)) ? (cmpNumbers(&n, &hi) == -1)
406 : PL_unify_number(A2, &hi))
407 );
408
409 } else
410 { rc = PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_number, A1);
411 }
412 clearNumber(&n);
413 clearNumber(&lo);
414 clearNumber(&hi);
415
416 return rc;
417 }
418
419 /*******************************
420 * BIGNUM FUNCTIONS *
421 *******************************/
422
423 #ifdef O_GMP
424
425 static int
get_mpz(term_t t,Number n ARG_LD)426 get_mpz(term_t t, Number n ARG_LD)
427 { Word p = valTermRef(t);
428
429 deRef(p);
430 if ( isInteger(*p) )
431 { get_integer(*p, n);
432 promoteToMPZNumber(n);
433
434 return TRUE;
435 }
436
437 return PL_type_error("integer", t);
438 }
439
440
441 /**
442 * divmod(+Dividend, +Divisor, -Quotient, -Remainder)
443 *
444 * Defined as
445 *
446 * - Quotient is div(Dividend, Divisor)
447 * - Remainder is mod(Dividend, Divisor)
448 */
449
450 static
451 PRED_IMPL("divmod", 4, divmod, 0)
452 { PRED_LD
453 number N = {V_INTEGER}, D = {V_INTEGER};
454 int rc = FALSE;
455
456 if ( get_mpz(A1, &N PASS_LD) &&
457 get_mpz(A2, &D PASS_LD) )
458 { if ( mpz_sgn(D.value.mpz) != 0 )
459 { number Q = {V_MPZ}, R = {V_MPZ};
460
461 mpz_init(R.value.mpz);
462 mpz_init(Q.value.mpz);
463 mpz_fdiv_qr(Q.value.mpz, R.value.mpz, N.value.mpz, D.value.mpz);
464 rc = ( PL_unify_number(A3, &Q) &&
465 PL_unify_number(A4, &R)
466 );
467 clearNumber(&R);
468 clearNumber(&Q);
469 } else
470 { rc = PL_error("divmod", 2, NULL, ERR_DIV_BY_ZERO);
471 }
472 }
473
474 clearNumber(&N);
475 clearNumber(&D);
476
477 return rc;
478 }
479
480 /**
481 * nth_integer_root_and_remainder(+N, +I, -Root, -Remainder)
482 */
483
484 static
485 PRED_IMPL("nth_integer_root_and_remainder", 4,
486 nth_integer_root_and_remainder, 0)
487 { PRED_LD
488 number N = {V_INTEGER};
489 long I;
490 int rc = FALSE;
491
492 if ( PL_get_long_ex(A1, &I) &&
493 get_mpz(A2, &N PASS_LD) )
494 { if ( I >= 1 )
495 { number root = {V_MPZ};
496 number rem = {V_MPZ};
497
498 if ( mpz_sgn(N.value.mpz) < 0 &&
499 I % 2 == 0 )
500 { rc = PL_error(NULL, 0, NULL, ERR_AR_UNDEF);
501 goto out;
502 }
503
504 mpz_init(root.value.mpz);
505 mpz_init(rem.value.mpz);
506 mpz_rootrem(root.value.mpz, rem.value.mpz,
507 N.value.mpz, (unsigned long)I);
508 rc = ( PL_unify_number(A3, &root) &&
509 PL_unify_number(A4, &rem)
510 );
511 clearNumber(&root);
512 clearNumber(&rem);
513 } else
514 { rc = PL_domain_error("not_less_than_one", A1);
515 }
516 }
517
518 out:
519 clearNumber(&N);
520
521 return rc;
522 }
523
524
525 static
526 PRED_IMPL("rational", 3, rational, 0)
527 { PRED_LD
528 Word p = valTermRef(A1);
529
530 deRef(p);
531 if ( isRational(*p) )
532 { if ( isMPQNum(*p) )
533 { number n, num, den;
534 int rc;
535
536 get_rational(*p, &n);
537 assert(n.type == V_MPQ);
538
539 num.type = V_MPZ;
540 den.type = V_MPZ;
541 mpz_init(num.value.mpz);
542 mpz_init(den.value.mpz);
543 mpz_set(num.value.mpz, mpq_numref(n.value.mpq));
544 mpz_set(den.value.mpz, mpq_denref(n.value.mpq));
545
546 rc = ( PL_unify_number(A2, &num) &&
547 PL_unify_number(A3, &den) );
548
549 clearNumber(&num);
550 clearNumber(&den);
551
552 return rc;
553 } else
554 { return ( PL_unify(A1, A2) &&
555 PL_unify_integer(A3, 1) );
556 }
557 }
558
559 return FALSE;
560 }
561
562
563 #endif /*O_GMP*/
564
565 static
566 PRED_IMPL("float_parts", 4, float_parts, 0)
567 { PRED_LD
568 double d;
569
570 if ( PL_get_float_ex(A1, &d) )
571 { double m;
572 int e;
573
574 m = frexp(d, &e);
575 return ( PL_unify_float(A2, m) &&
576 PL_unify_integer(A3, 2) &&
577 PL_unify_integer(A4, e) );
578 }
579
580 return FALSE;
581 }
582
583
584 /********************************
585 * COMPARISON *
586 *********************************/
587
588 /* implements <, =<, >, >=, =:= and =\=
589 */
590
591 int
ar_compare(Number n1,Number n2,int what)592 ar_compare(Number n1, Number n2, int what)
593 { int diff = cmpNumbers(n1, n2); /* nan compares CMP_NOTEQ */
594
595 switch(what)
596 { case LT: return diff == CMP_LESS;
597 case GT: return diff == CMP_GREATER;
598 case LE: return (diff == CMP_LESS) || (diff == CMP_EQUAL);
599 case GE: return (diff == CMP_GREATER) || (diff == CMP_EQUAL);
600 case NE: return diff != CMP_EQUAL;
601 case EQ: return diff == CMP_EQUAL;
602 default:
603 assert(0);
604 return FALSE;
605 }
606 }
607
608
609 static word
compareNumbers(term_t n1,term_t n2,int what ARG_LD)610 compareNumbers(term_t n1, term_t n2, int what ARG_LD)
611 { AR_CTX
612 number left, right;
613 int rc;
614
615 AR_BEGIN();
616
617 if ( valueExpression(n1, &left PASS_LD) &&
618 valueExpression(n2, &right PASS_LD) )
619 { rc = ar_compare(&left, &right, what);
620
621 clearNumber(&left);
622 clearNumber(&right);
623 AR_END();
624 } else
625 { AR_CLEANUP();
626 rc = FALSE;
627 }
628
629 return rc;
630 }
631
632 static
633 PRED_IMPL("<", 2, lt, PL_FA_ISO)
634 { PRED_LD
635 return compareNumbers(A1, A2, LT PASS_LD);
636 }
637
638 static
639 PRED_IMPL(">", 2, gt, PL_FA_ISO)
640 { PRED_LD
641 return compareNumbers(A1, A2, GT PASS_LD);
642 }
643
644 static
645 PRED_IMPL("=<", 2, leq, PL_FA_ISO)
646 { PRED_LD
647 return compareNumbers(A1, A2, LE PASS_LD);
648 }
649
650 static
651 PRED_IMPL(">=", 2, geq, PL_FA_ISO)
652 { PRED_LD
653 return compareNumbers(A1, A2, GE PASS_LD);
654 }
655
656 static
657 PRED_IMPL("=\\=", 2, neq, PL_FA_ISO)
658 { PRED_LD
659 return compareNumbers(A1, A2, NE PASS_LD);
660 }
661
662 static
663 PRED_IMPL("=:=", 2, eq, PL_FA_ISO)
664 { PRED_LD
665 return compareNumbers(A1, A2, EQ PASS_LD);
666 }
667
668
669 /*******************************
670 * ARITHMETIC STACK *
671 *******************************/
672
673 Number
growArithStack(ARG1_LD)674 growArithStack(ARG1_LD)
675 { Number n;
676
677 if ( LD->arith.stack.top == LD->arith.stack.max )
678 { size_t size;
679
680 if ( LD->arith.stack.base )
681 { size = (size_t)(LD->arith.stack.max - LD->arith.stack.base);
682 LD->arith.stack.base = PL_realloc(LD->arith.stack.base,
683 size*sizeof(number)*2);
684 LD->arith.stack.top = LD->arith.stack.base+size;
685 size *= 2;
686 } else
687 { size = 16;
688 LD->arith.stack.base = PL_malloc(size*sizeof(number));
689 LD->arith.stack.top = LD->arith.stack.base;
690 }
691
692 LD->arith.stack.max = LD->arith.stack.base+size;
693 }
694
695 n = LD->arith.stack.top;
696 LD->arith.stack.top++;
697
698 return n;
699 }
700
701
702 void
freeArithLocalData(PL_local_data_t * ld)703 freeArithLocalData(PL_local_data_t *ld)
704 { if ( ld->arith.stack.base )
705 PL_free(ld->arith.stack.base);
706 #ifdef O_GMP
707 if ( ld->arith.random.initialised )
708 { DEBUG(0, { GET_LD
709 assert(ld == LD);
710 });
711
712 ld->gmp.persistent++;
713 gmp_randclear(ld->arith.random.state);
714 ld->gmp.persistent--;
715 ld->arith.random.initialised = FALSE;
716 }
717 #endif
718 }
719
720
721 /********************************
722 * FUNCTIONS *
723 *********************************/
724
725 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
726 isCurrentArithFunction(functor_t f)
727 Find existing arithmetic function definition for f.
728 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
729
730 static inline ArithF
isCurrentArithFunction(functor_t f)731 isCurrentArithFunction(functor_t f)
732 { size_t index = indexFunctor(f);
733
734 if ( index < GD->arith.functions_allocated )
735 { return GD->arith.functions[index];
736 }
737
738 return NULL;
739 }
740
741
742 int
check_float(Number n)743 check_float(Number n)
744 { PL_error_code code = ERR_NO_ERROR;
745 #ifdef HAVE_FPCLASSIFY
746 switch(fpclassify(n->value.f))
747 { case FP_NAN:
748 code = ERR_AR_UNDEF;
749 break;
750 case FP_SUBNORMAL:
751 code = ERR_AR_UNDERFLOW;
752 break;
753 case FP_INFINITE:
754 code = ERR_AR_OVERFLOW;
755 break;
756 }
757 #else
758 #ifdef HAVE_FPCLASS
759 switch(fpclass(n->value.f))
760 { case FP_SNAN:
761 case FP_QNAN:
762 code = ERR_AR_UNDEF;
763 break;
764 case FP_NINF:
765 case FP_PINF:
766 code = ERR_AR_OVERFLOW;
767 break;
768 case FP_NDENORM: /* pos/neg denormalized non-zero */
769 case FP_PDENORM:
770 code = ERR_AR_UNDERFLOW;
771 break;
772 case FP_NNORM: /* pos/neg normalized non-zero */
773 case FP_PNORM:
774 case FP_NZERO: /* pos/neg zero */
775 case FP_PZERO:
776 break;
777 }
778 #else
779 #ifdef HAVE__FPCLASS
780 switch(_fpclass(n->value.f))
781 { case _FPCLASS_SNAN:
782 case _FPCLASS_QNAN:
783 code = ERR_AR_UNDEF;
784 break;
785 case _FPCLASS_NINF:
786 case _FPCLASS_PINF:
787 code = ERR_AR_OVERFLOW;
788 break;
789 }
790 #else
791 #ifdef HAVE_ISNAN
792 if ( isnan(n->value.f) )
793 code = ERR_AR_UNDEF;
794 #endif
795 #ifdef HAVE_ISINF
796 if ( isinf(n->value.f) )
797 code = ERR_AR_OVERFLOW;
798 #endif
799 #endif /*HAVE__FPCLASS*/
800 #endif /*HAVE_FPCLASS*/
801 #endif /*HAVE_FPCLASSIFY*/
802
803 if ( code != ERR_NO_ERROR )
804 { GET_LD
805
806 switch(code)
807 { case ERR_AR_OVERFLOW:
808 if ( LD->arith.f.flags & FLT_OVERFLOW )
809 return TRUE;
810 break;
811 case ERR_AR_UNDERFLOW:
812 if ( LD->arith.f.flags & FLT_UNDERFLOW )
813 return TRUE;
814 break;
815 case ERR_AR_UNDEF:
816 n->value.f = const_nan;
817 if ( LD->arith.f.flags & FLT_UNDEFINED )
818 return TRUE;
819 break;
820 default:
821 assert(0);
822 }
823
824 return PL_error(NULL, 0, NULL, code);
825 }
826
827 return TRUE;
828 }
829
830 static int
check_zero_div(int sign_n,Number r,char * func,int arity)831 check_zero_div(int sign_n, Number r, char *func, int arity)
832 { GET_LD
833
834 if ( LD->arith.f.flags & FLT_ZERO_DIV )
835 { r->type = V_FLOAT;
836 r->value.f = copysign(const_inf,sign_n);
837 return TRUE;
838 } else
839 { return PL_error(func, arity, NULL, ERR_DIV_BY_ZERO);
840 }
841 }
842
843
844 #ifdef O_GMP
845 static int
check_mpq(Number r)846 check_mpq(Number r)
847 { GET_LD
848 size_t sz;
849
850 if ( (sz=LD->arith.rat.max_rational_size) != (size_t)-1 )
851 { int szn = mpq_numref(r->value.mpq)->_mp_size;
852 int szd = mpq_denref(r->value.mpq)->_mp_size;
853
854 if ( szn < 0 ) szn = -szn;
855 if ( szd < 0 ) szd = -szd;
856
857 if ( ( szn + szd ) * sizeof(mp_limb_t) > sz )
858 { atom_t action = LD->arith.rat.max_rational_size_action;
859
860 if ( action == ATOM_float )
861 promoteToFloatNumber(r);
862 else if ( action == ATOM_error )
863 return PL_error(NULL, 0, "requires more than max_rational_size bytes",
864 ERR_AR_TRIPWIRE, ATOM_max_rational_size, r);
865 else
866 assert(0);
867 }
868 }
869
870 return TRUE;
871 }
872 #endif
873
874 /*******************************
875 * EVALULATE *
876 *******************************/
877
878 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
879 valueExpression() evaluates an `evaluable term'.
880
881 This new implementation avoids using the C-stack to be able to process
882 more deeply nested terms and to be able to recover in the unlikely case
883 that terms are still too deeply nested.
884
885 If it finds a term, it starts processing at the last argument, working back
886 to the start. If it finds the functor itself, it evaluates the pushed
887 arguments. Using this technique we push as few as possible arguments on
888 terms that are nested on the left (as in (1+2)+3, while we only push a
889 single pointer for each recursion level in the evaluable term.
890 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
891
892 static int
pushForMark(segstack * stack,Word p,int wr)893 pushForMark(segstack *stack, Word p, int wr)
894 { word w = ((word)p)|wr;
895
896 return pushSegStack(stack, w, word);
897 }
898
899 static void
popForMark(segstack * stack,Word * pp,int * wr)900 popForMark(segstack *stack, Word *pp, int *wr)
901 { word w = 0;
902
903 popSegStack(stack, &w, word);
904 *wr = w & (word)0x1;
905 *pp = (Word)(w & ~(word)0x1);
906 }
907
908
909 int
valueExpression(term_t expr,number * result ARG_LD)910 valueExpression(term_t expr, number *result ARG_LD)
911 { segstack term_stack;
912 segstack arg_stack;
913 Word term_buf[16];
914 number arg_buf[16];
915 number *n = result;
916 number n_tmp;
917 int walk_ref = FALSE;
918 Word p = valTermRef(expr);
919 Word start;
920 int known_acyclic = FALSE;
921 int pushed = 0;
922 functor_t functor;
923 int old_round_mode = fegetround();
924
925 deRef(p);
926 start = p;
927 LD->in_arithmetic++;
928
929 for(;;)
930 { switch(tag(*p))
931 { case TAG_INTEGER:
932 get_rational(*p, n);
933 break;
934 case TAG_FLOAT:
935 n->value.f = valFloat(*p);
936 n->type = V_FLOAT;
937 break;
938 case TAG_VAR:
939 PL_error(NULL, 0, NULL, ERR_INSTANTIATION);
940 goto error;
941 case TAG_REFERENCE:
942 { if ( !pushForMark(&term_stack, p, walk_ref) )
943 { PL_no_memory();
944 goto error;
945 }
946 walk_ref = TRUE;
947 deRef(p);
948 continue;
949 }
950 case TAG_ATOM:
951 { ArithF f;
952
953 functor = lookupFunctorDef(*p, 0);
954 arity0:
955 if ( (f = isCurrentArithFunction(functor)) )
956 { if ( (*f)(n) != TRUE )
957 goto error;
958 } else
959 { if ( isTextAtom(*p) )
960 { PL_error(NULL, 0, NULL, ERR_NOT_EVALUABLE, functor);
961 } else
962 { PL_error(NULL, 0, NULL, ERR_TYPE,
963 ATOM_evaluable, pushWordAsTermRef(p));
964 popTermRef();
965 }
966 goto error;
967 }
968 break;
969 }
970 case TAG_STRING:
971 if ( getCharExpression(p, n PASS_LD) != TRUE )
972 goto error;
973 break;
974 case TAG_COMPOUND:
975 { Functor term = valueTerm(*p);
976 size_t arity = arityFunctor(term->definition);
977
978 if ( term->definition == FUNCTOR_dot2 )
979 { if ( getCharExpression(p, n PASS_LD) != TRUE )
980 goto error;
981 break;
982 }
983
984 if ( arity == 0 )
985 { functor = term->definition;
986 goto arity0;
987 }
988
989 if ( p == start )
990 { initSegStack(&term_stack, sizeof(Word), sizeof(term_buf), term_buf);
991 initSegStack(&arg_stack, sizeof(number), sizeof(arg_buf), arg_buf);
992 }
993
994 if ( !pushForMark(&term_stack, p, walk_ref) )
995 { PL_no_memory();
996 goto error;
997 }
998 if ( ++pushed > 100 && !known_acyclic )
999 { int rc;
1000
1001 if ( (rc=is_acyclic(start PASS_LD)) == TRUE )
1002 { known_acyclic = TRUE;
1003 } else
1004 { if ( rc == MEMORY_OVERFLOW )
1005 PL_error(NULL, 0, NULL, ERR_NOMEM);
1006 else
1007 PL_error(NULL, 0, "cyclic term", ERR_TYPE, ATOM_expression, expr);
1008 goto error;
1009 }
1010 }
1011 if ( term->definition == FUNCTOR_roundtoward2 )
1012 { number crnd;
1013
1014 if ( !set_roundtoward(&term->arguments[1], &crnd PASS_LD) )
1015 goto error;
1016 if ( !pushSegStack(&arg_stack, crnd, number) )
1017 { PL_no_memory();
1018 goto error;
1019 }
1020 p = &term->arguments[0];
1021 } else
1022 { p = &term->arguments[arity-1];
1023 }
1024 walk_ref = FALSE;
1025 n = &n_tmp;
1026 continue;
1027 }
1028 default:
1029 PL_error(NULL, 0, NULL, ERR_PTR_TYPE, ATOM_number, p);
1030 goto error;
1031 }
1032
1033 if ( p == start )
1034 { LD->in_arithmetic--;
1035 assert(n == result);
1036
1037 return TRUE;
1038 }
1039
1040 if ( walk_ref )
1041 popForMark(&term_stack, &p, &walk_ref);
1042
1043 if ( !pushSegStack(&arg_stack, n_tmp, number) )
1044 { PL_no_memory();
1045 goto error;
1046 }
1047
1048 while ( tagex(*--p) == (TAG_ATOM|STG_GLOBAL) )
1049 { functor_t functor = *p;
1050 ArithF f;
1051
1052 DEBUG(1, Sdprintf("Eval %s/%d\n",
1053 stringAtom(nameFunctor(functor)),
1054 arityFunctor(functor)));
1055
1056 if ( (f = isCurrentArithFunction(functor)) )
1057 { size_t arity = arityFunctor(functor);
1058
1059 switch(arity)
1060 { case 1:
1061 { int rc;
1062 number *a0 = topOfSegStack(&arg_stack);
1063
1064 rc = (*f)(a0, n);
1065 clearNumber(a0);
1066 if ( rc == TRUE )
1067 { *a0 = *n;
1068 } else
1069 { popTopOfSegStack(&arg_stack);
1070 goto error;
1071 }
1072
1073 break;
1074 }
1075 case 2:
1076 { int rc;
1077 void *a[2];
1078
1079 topsOfSegStack(&arg_stack, 2, a);
1080 rc = (*f)((number*)a[0], (number*)a[1], n);
1081 clearNumber((number*)a[0]);
1082 clearNumber((number*)a[1]);
1083 popTopOfSegStack(&arg_stack);
1084
1085 if ( rc == TRUE )
1086 { number *n1 = a[1];
1087 *n1 = *n;
1088 } else
1089 { popTopOfSegStack(&arg_stack);
1090 goto error;
1091 }
1092
1093 break;
1094 }
1095 case 3:
1096 { int rc;
1097 void *a[3];
1098
1099 topsOfSegStack(&arg_stack, 3, a);
1100 rc = (*f)((number*)a[0], (number*)a[1], (number*)a[2], n);
1101 clearNumber((number*)a[0]);
1102 clearNumber((number*)a[1]);
1103 clearNumber((number*)a[2]);
1104 popTopOfSegStack(&arg_stack);
1105 popTopOfSegStack(&arg_stack);
1106
1107 if ( rc == TRUE )
1108 { number *n2 = a[2];
1109 *n2 = *n;
1110 } else
1111 { popTopOfSegStack(&arg_stack);
1112 goto error;
1113 }
1114
1115 break;
1116 }
1117 default:
1118 assert(0);
1119 }
1120
1121 popForMark(&term_stack, &p, &walk_ref);
1122 if ( p == start )
1123 { LD->in_arithmetic--;
1124 *result = *n;
1125
1126 return TRUE;
1127 }
1128 if ( walk_ref )
1129 popForMark(&term_stack, &p, &walk_ref);
1130 } else
1131 { PL_error(NULL, 0, NULL, ERR_NOT_EVALUABLE, functor);
1132 goto error;
1133 }
1134 }
1135 }
1136
1137 error:
1138 if ( p != start )
1139 { number n;
1140
1141 clearSegStack(&term_stack);
1142 while( popSegStack(&arg_stack, &n, number) )
1143 clearNumber(&n);
1144 fesetround(old_round_mode);
1145 }
1146 LD->in_arithmetic--;
1147
1148 return FALSE;
1149 }
1150
1151
1152 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1153 int arithChar(Word p)
1154 Handle arithmetic argument "x", normally appearing as [X], where X
1155 is an integer or one-character atom.
1156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1157
1158 int
arithChar(Word p ARG_LD)1159 arithChar(Word p ARG_LD)
1160 { deRef(p);
1161
1162 if ( isInteger(*p) )
1163 { intptr_t chr = valInt(*p);
1164
1165 if ( chr >= 0 && chr <= PLMAXWCHAR )
1166 return (int)chr;
1167 } else if ( isAtom(*p) )
1168 { PL_chars_t txt;
1169
1170 if ( get_atom_text(*p, &txt) && txt.length == 1 )
1171 { if ( txt.encoding == ENC_WCHAR )
1172 return txt.text.w[0];
1173 else
1174 return txt.text.t[0]&0xff;
1175 }
1176 }
1177
1178 PL_error(NULL, 0, NULL, ERR_TYPE,
1179 ATOM_character, pushWordAsTermRef(p));
1180 popTermRef();
1181
1182 return EOF;
1183 }
1184
1185
1186 int
getCharExpression(Word p,Number r ARG_LD)1187 getCharExpression(Word p, Number r ARG_LD)
1188 { word w = *p;
1189
1190 switch(tag(w))
1191 { case TAG_STRING:
1192 { size_t len;
1193
1194 if ( isBString(w) )
1195 { char *s = getCharsString(w, &len);
1196
1197 if ( len == 1 )
1198 { r->value.i = s[0]&0xff;
1199 r->type = V_INTEGER;
1200 return TRUE;
1201 }
1202 } else
1203 { pl_wchar_t *ws = getCharsWString(w, &len);
1204
1205 if ( len == 1 )
1206 { r->value.i = ws[0];
1207 r->type = V_INTEGER;
1208 return TRUE;
1209 }
1210 }
1211
1212 len_not_one:
1213 PL_error(NULL, 0, "\"x\" must hold one character", ERR_TYPE,
1214 ATOM_nil, pushWordAsTermRef(p));
1215 popTermRef();
1216 return FALSE;
1217 }
1218 case TAG_COMPOUND:
1219 { Word a = argTermP(w, 0);
1220 int chr;
1221
1222 if ( (chr = arithChar(a PASS_LD)) == EOF )
1223 fail;
1224
1225 a = argTermP(w, 1);
1226 if ( !isNil(*a) )
1227 goto len_not_one;
1228
1229 r->value.i = chr;
1230 r->type = V_INTEGER;
1231
1232 return TRUE;
1233 }
1234 default:
1235 assert(0);
1236 return FALSE;
1237 }
1238 }
1239
1240
1241
1242
1243 /*******************************
1244 * CONVERSION *
1245 *******************************/
1246
1247 static int
double_in_int64_range(double x)1248 double_in_int64_range(double x)
1249 { int k;
1250 double y = frexp(x, &k);
1251
1252 if ( k < 8*(int)sizeof(int64_t) ||
1253 (y == -0.5 && k == 8*(int)sizeof(int64_t)) )
1254 return TRUE;
1255
1256 return FALSE;
1257 }
1258
1259
1260 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1261 toIntegerNumber(Number n, int flags)
1262
1263 Convert a number to an integer. Default, only rationals that happen to
1264 be integer are converted. If TOINT_CONVERT_FLOAT is present, floating
1265 point numbers are converted if they represent integers. If also
1266 TOINT_TRUNCATE is provided non-integer floats are truncated to integers.
1267
1268 Note that if a double is out of range for int64_t, it never has a
1269 fractional part.
1270 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1271
1272 int
toIntegerNumber(Number n,int flags)1273 toIntegerNumber(Number n, int flags)
1274 { switch(n->type)
1275 { case V_INTEGER:
1276 succeed;
1277 #ifdef O_GMP
1278 case V_MPZ:
1279 succeed;
1280 case V_MPQ: /* never from stacks iff integer */
1281 if ( mpz_cmp_ui(mpq_denref(n->value.mpq), 1L) == 0 )
1282 { mpz_clear(mpq_denref(n->value.mpq));
1283 n->value.mpz[0] = mpq_numref(n->value.mpq)[0];
1284 n->type = V_MPZ;
1285 succeed;
1286 }
1287 fail;
1288 #endif
1289 case V_FLOAT:
1290 if ( !check_float(n) )
1291 return FALSE;
1292 if ( (flags & TOINT_CONVERT_FLOAT) )
1293 { if ( double_in_int64_range(n->value.f) )
1294 { int64_t l = (int64_t)n->value.f;
1295
1296 if ( (flags & TOINT_TRUNCATE) ||
1297 (double)l == n->value.f )
1298 { n->value.i = l;
1299 n->type = V_INTEGER;
1300
1301 return TRUE;
1302 }
1303 return FALSE;
1304 #ifdef O_GMP
1305 } else
1306 { mpz_init_set_d(n->value.mpz, n->value.f);
1307 n->type = V_MPZ;
1308
1309 return TRUE;
1310 #endif
1311 }
1312 }
1313 return FALSE;
1314 }
1315
1316 assert(0);
1317 fail;
1318 }
1319
1320
1321 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1322 promoteIntNumber() promotes a number of type V_INTEGER to a number with
1323 larger capacity.
1324 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1325
1326 static int
promoteIntNumber(Number n)1327 promoteIntNumber(Number n)
1328 {
1329 #ifdef O_GMP
1330 promoteToMPZNumber(n);
1331 #else
1332 GET_LD
1333
1334 if ( truePrologFlag(PLFLAG_ISO) )
1335 return PL_error(NULL, 0, NULL, ERR_EVALUATION, ATOM_int_overflow);
1336
1337 return promoteToFloatNumber(n);
1338 #endif
1339
1340 succeed;
1341 }
1342
1343
1344
1345 /********************************
1346 * ARITHMETIC FUNCTIONS *
1347 *********************************/
1348
1349 static int ar_u_minus(Number n1, Number r);
1350
1351 int
ar_add_ui(Number n,intptr_t add)1352 ar_add_ui(Number n, intptr_t add)
1353 { switch(n->type)
1354 { case V_INTEGER:
1355 { int64_t r = n->value.i + add;
1356
1357 if ( (r < 0 && add > 0 && n->value.i > 0) ||
1358 (r > 0 && add < 0 && n->value.i < 0) )
1359 { if ( !promoteIntNumber(n) )
1360 fail;
1361 } else
1362 { n->value.i = r;
1363 succeed;
1364 }
1365 }
1366 /*FALLTHROUGH*/
1367 #ifdef O_GMP
1368 case V_MPZ:
1369 { if ( add > 0 )
1370 mpz_add_ui(n->value.mpz, n->value.mpz, (unsigned long)add);
1371 else
1372 mpz_sub_ui(n->value.mpz, n->value.mpz, (unsigned long)-add);
1373
1374 succeed;
1375 }
1376 case V_MPQ:
1377 { if ( add > 0 )
1378 mpz_addmul_ui(mpq_numref(n->value.mpq), mpq_denref(n->value.mpq),
1379 (unsigned long)add);
1380 else
1381 mpz_submul_ui(mpq_numref(n->value.mpq), mpq_denref(n->value.mpq),
1382 (unsigned long)-add);
1383
1384 return check_mpq(n);
1385 }
1386 #endif
1387 case V_FLOAT:
1388 { n->value.f += (double)add;
1389
1390 return check_float(n);
1391 }
1392 default:
1393 ;
1394 }
1395
1396 assert(0);
1397 fail;
1398 }
1399
1400 #define SAME_SIGN(i1, i2) (((i1) ^ (i2)) >= 0)
1401
1402 int
pl_ar_add(Number n1,Number n2,Number r)1403 pl_ar_add(Number n1, Number n2, Number r)
1404 { if ( !same_type_numbers(n1, n2) )
1405 return FALSE;
1406
1407 switch(n1->type)
1408 { case V_INTEGER:
1409 { if ( SAME_SIGN(n1->value.i, n2->value.i) )
1410 { if ( n2->value.i < 0 ) /* both negative */
1411 { if ( n1->value.i < PLMININT - n2->value.i )
1412 goto overflow;
1413 } else /* both positive */
1414 { if ( PLMAXINT - n1->value.i < n2->value.i )
1415 goto overflow;
1416 }
1417 }
1418 r->value.i = n1->value.i + n2->value.i;
1419 r->type = V_INTEGER;
1420 succeed;
1421 overflow:
1422 if ( !promoteIntNumber(n1) ||
1423 !promoteIntNumber(n2) )
1424 fail;
1425 }
1426 /*FALLTHROUGH*/
1427 #ifdef O_GMP
1428 case V_MPZ:
1429 { r->type = V_MPZ;
1430 mpz_init(r->value.mpz);
1431 mpz_add(r->value.mpz, n1->value.mpz, n2->value.mpz);
1432 succeed;
1433 }
1434 case V_MPQ:
1435 { r->type = V_MPQ;
1436 mpq_init(r->value.mpq);
1437 mpq_add(r->value.mpq, n1->value.mpq, n2->value.mpq);
1438 return check_mpq(r);
1439 }
1440 #endif
1441 case V_FLOAT:
1442 { r->value.f = n1->value.f + n2->value.f;
1443 r->type = V_FLOAT;
1444
1445 return check_float(r);
1446 }
1447 }
1448
1449 assert(0);
1450 fail;
1451 }
1452
1453
1454 static int
ar_minus(Number n1,Number n2,Number r)1455 ar_minus(Number n1, Number n2, Number r)
1456 { if ( !same_type_numbers(n1, n2) )
1457 return FALSE;
1458
1459 switch(n1->type)
1460 { case V_INTEGER:
1461 { r->value.i = (uint64_t)n1->value.i - (uint64_t)n2->value.i;
1462
1463 if ( (n1->value.i >= 0 && n2->value.i < 0 && r->value.i <= 0) ||
1464 (n1->value.i < 0 && n2->value.i > 0 && r->value.i >= 0) )
1465 { /* overflow */
1466 if ( !promoteIntNumber(n1) ||
1467 !promoteIntNumber(n2) )
1468 fail;
1469 } else
1470 { r->type = V_INTEGER;
1471 succeed;
1472 }
1473 }
1474 #ifdef O_GMP
1475 case V_MPZ:
1476 { r->type = V_MPZ;
1477 mpz_init(r->value.mpz);
1478 mpz_sub(r->value.mpz, n1->value.mpz, n2->value.mpz);
1479 succeed;
1480 }
1481 case V_MPQ:
1482 { r->type = V_MPQ;
1483 mpq_init(r->value.mpq);
1484 mpq_sub(r->value.mpq, n1->value.mpq, n2->value.mpq);
1485 return check_mpq(r);
1486 succeed;
1487 }
1488 #endif
1489 case V_FLOAT:
1490 { r->value.f = n1->value.f - n2->value.f;
1491 r->type = V_FLOAT;
1492
1493 return check_float(r);
1494 }
1495 }
1496
1497 assert(0);
1498 fail;
1499 }
1500
1501
1502 static int
ar_even(Number i)1503 ar_even(Number i)
1504 { switch(i->type)
1505 { case V_INTEGER:
1506 return i->value.i % 2 == 0;
1507 #ifdef O_GMP
1508 case V_MPZ:
1509 return mpz_fdiv_ui(i->value.mpz, 2) == 0;
1510 #endif
1511 default:
1512 assert(0);
1513 fail;
1514 }
1515 }
1516
1517
1518 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1519 mod(X, Y) = X - (floor(X/Y) * Y)
1520 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1521
1522 static inline int64_t
mod(int64_t x,int64_t y)1523 mod(int64_t x, int64_t y)
1524 { int64_t r = x % y;
1525
1526 if ( r != 0 && (r<0) != (y<0) )
1527 r += y;
1528
1529 return r;
1530 }
1531
1532
1533 static int
ar_mod(Number n1,Number n2,Number r)1534 ar_mod(Number n1, Number n2, Number r)
1535 { if ( !toIntegerNumber(n1, 0) )
1536 return PL_error("mod", 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
1537 if ( !toIntegerNumber(n2, 0) )
1538 return PL_error("mod", 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
1539
1540 if ( !same_type_numbers(n1, n2) )
1541 return FALSE;
1542
1543 switch(n1->type)
1544 { case V_INTEGER:
1545 if ( n2->value.i == 0 )
1546 return PL_error("mod", 2, NULL, ERR_DIV_BY_ZERO);
1547
1548 if ( n2->value.i != -1 || n1->value.i != INT64_MIN )
1549 r->value.i = mod(n1->value.i, n2->value.i);
1550 else
1551 r->value.i = 0;
1552 r->type = V_INTEGER;
1553 break;
1554 #ifdef O_GMP
1555 case V_MPZ:
1556 if ( mpz_sgn(n2->value.mpz) == 0 )
1557 return PL_error("mod", 2, NULL, ERR_DIV_BY_ZERO);
1558
1559 r->type = V_MPZ;
1560 mpz_init(r->value.mpz);
1561 mpz_fdiv_r(r->value.mpz, n1->value.mpz, n2->value.mpz);
1562 break;
1563 #endif
1564 default:
1565 assert(0);
1566 }
1567
1568 succeed;
1569 }
1570
1571
1572 static int
int_too_big(void)1573 int_too_big(void)
1574 { GET_LD
1575 return (int)outOfStack((Stack)&LD->stacks.global, STACK_OVERFLOW_RAISE);
1576 }
1577
1578
1579 static int
shift_to_far(Number shift,Number r,int dir)1580 shift_to_far(Number shift, Number r, int dir)
1581 { if ( ar_sign_i(shift) * dir < 0 ) /* << */
1582 { return int_too_big();
1583 } else
1584 { r->value.i = 0;
1585 r->type = V_INTEGER;
1586
1587 return TRUE;
1588 }
1589 }
1590
1591
1592 static int
ar_shift(Number n1,Number n2,Number r,int dir)1593 ar_shift(Number n1, Number n2, Number r, int dir)
1594 { long shift;
1595 const char *plop = (dir < 0 ? "<<" : ">>");
1596
1597 if ( !toIntegerNumber(n1, 0) )
1598 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
1599 if ( !toIntegerNumber(n2, 0) )
1600 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
1601
1602 if ( ar_sign_i(n1) == 0 ) /* shift of 0 is always 0 */
1603 { r->value.i = 0;
1604 r->type = V_INTEGER;
1605 }
1606
1607 switch(n2->type) /* amount to shift */
1608 { case V_INTEGER:
1609 if ( n2->value.i < LONG_MIN ||
1610 n2->value.i > LONG_MAX )
1611 return shift_to_far(n2, r, dir);
1612 else
1613 shift = (long)n2->value.i;
1614 break;
1615 #ifdef O_GMP
1616 case V_MPZ:
1617 if ( mpz_cmp_si(n2->value.mpz, LONG_MIN) < 0 ||
1618 mpz_cmp_si(n2->value.mpz, LONG_MAX) > 0 )
1619 return shift_to_far(n2, r, dir);
1620 else
1621 shift = mpz_get_si(n2->value.mpz);
1622 break;
1623 #endif
1624 default:
1625 assert(0);
1626 fail;
1627 }
1628
1629 if ( shift < 0 )
1630 { shift = -shift;
1631 dir = -dir;
1632 }
1633
1634 switch(n1->type)
1635 { case V_INTEGER:
1636 if ( dir < 0 ) /* shift left (<<) */
1637 {
1638 #ifdef O_GMP /* msb() is 0..63 */
1639 int bits = shift;
1640
1641 if ( n1->value.i >= 0 )
1642 bits += MSB64(n1->value.i);
1643 else if ( n1->value.i == PLMININT )
1644 bits += sizeof(int64_t)*8;
1645 else
1646 bits += MSB64(-n1->value.i);
1647
1648 if ( bits >= (int)(sizeof(int64_t)*8-1) )
1649 { promoteToMPZNumber(n1);
1650 goto mpz;
1651 } else
1652 #endif
1653 { r->value.i = n1->value.i << shift;
1654 }
1655 } else
1656 { if ( shift >= (long)sizeof(int64_t)*8 )
1657 r->value.i = (n1->value.i >= 0 ? 0 : -1);
1658 else
1659 r->value.i = n1->value.i >> shift;
1660 }
1661 r->type = V_INTEGER;
1662 succeed;
1663 #ifdef O_GMP
1664 case V_MPZ:
1665 mpz:
1666 r->type = V_MPZ;
1667 mpz_init(r->value.mpz);
1668 if ( dir < 0 )
1669 {
1670 #ifdef O_GMP_PRECHECK_ALLOCATIONS
1671 GET_LD
1672 uint64_t msb = mpz_sizeinbase(n1->value.mpz, 2)+shift;
1673
1674 if ( (msb/sizeof(char)) > (uint64_t)globalStackLimit() )
1675 { mpz_clear(r->value.mpz);
1676 return int_too_big();
1677 }
1678 #endif /*O_GMP_PRECHECK_ALLOCATIONS*/
1679 mpz_mul_2exp(r->value.mpz, n1->value.mpz, shift);
1680 } else
1681 { mpz_fdiv_q_2exp(r->value.mpz, n1->value.mpz, shift);
1682 }
1683 succeed;
1684 #endif
1685 default:
1686 assert(0);
1687 fail;
1688 }
1689 }
1690
1691
1692 static int
ar_shift_left(Number n1,Number n2,Number r)1693 ar_shift_left(Number n1, Number n2, Number r)
1694 { return ar_shift(n1, n2, r, -1);
1695 }
1696
1697
1698 static int
ar_shift_right(Number n1,Number n2,Number r)1699 ar_shift_right(Number n1, Number n2, Number r)
1700 { return ar_shift(n1, n2, r, 1);
1701 }
1702
1703
1704 static int
same_positive_ints(const char * fname,Number n1,Number n2)1705 same_positive_ints(const char *fname, Number n1, Number n2)
1706 { if ( !toIntegerNumber(n1, 0) )
1707 return PL_error(fname, 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
1708 if ( !toIntegerNumber(n2, 0) )
1709 return PL_error(fname, 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
1710
1711 if ( !same_type_numbers(n1, n2) )
1712 return FALSE;
1713
1714 switch(n1->type)
1715 { case V_INTEGER:
1716 { int64_t a = n1->value.i;
1717 int64_t b = n2->value.i;
1718
1719 if ( a < 0 )
1720 { a = -(uint64_t)a;
1721 if ( a < 0 )
1722 { promote:
1723 #ifdef O_GMP
1724 promoteToMPZNumber(n1);
1725 promoteToMPZNumber(n2);
1726 goto case_gmp;
1727 #else
1728 return PL_error(fname, 2, NULL, ERR_EVALUATION, ATOM_int_overflow);
1729 #endif
1730 }
1731 }
1732 if ( b < 0 )
1733 { b = -(uint64_t)b;
1734 if ( b < 0 )
1735 goto promote;
1736 }
1737
1738 n1->value.i = a;
1739 n2->value.i = b;
1740 break;
1741 }
1742 #ifdef O_GMP
1743 case V_MPZ:
1744 case_gmp:
1745 /* we don't really need to make absolute here as the GMP functions
1746 * ignore the sign anyway
1747 */
1748 break;
1749 #endif
1750 default:
1751 assert(0);
1752 }
1753
1754 return TRUE;
1755 }
1756
1757
1758 static int64_t
i64_gcd(int64_t a,int64_t b)1759 i64_gcd(int64_t a, int64_t b)
1760 { int64_t t;
1761
1762 if ( a == 0 )
1763 return b;
1764 if ( b == 0 )
1765 return a;
1766
1767 while(b != 0)
1768 { t = b;
1769 b = a % b;
1770 a = t;
1771 }
1772
1773 return a;
1774 }
1775
1776
1777 static int
ar_gcd(Number n1,Number n2,Number r)1778 ar_gcd(Number n1, Number n2, Number r)
1779 { if ( !same_positive_ints("gcd", n1, n2) )
1780 return FALSE;
1781
1782 switch(n1->type)
1783 { case V_INTEGER:
1784 { r->type = V_INTEGER;
1785 r->value.i = i64_gcd(n1->value.i, n2->value.i);
1786 break;
1787 }
1788 #ifdef O_GMP
1789 case V_MPZ:
1790 r->type = V_MPZ;
1791 mpz_init(r->value.mpz);
1792 mpz_gcd(r->value.mpz, n1->value.mpz, n2->value.mpz);
1793 break;
1794 #endif
1795 default:
1796 assert(0);
1797 }
1798
1799 return TRUE;
1800 }
1801
1802 static int
ar_lcm(Number n1,Number n2,Number r)1803 ar_lcm(Number n1, Number n2, Number r)
1804 { if ( !same_positive_ints("lcm", n1, n2) )
1805 return FALSE;
1806
1807 switch(n1->type)
1808 { case V_INTEGER:
1809 { int64_t prod;
1810
1811 if ( mul64(n1->value.i, n2->value.i, &prod) )
1812 { r->type = V_INTEGER;
1813 if ( prod != 0 )
1814 r->value.i = prod/i64_gcd(n1->value.i, n2->value.i);
1815 else
1816 r->value.i = 0;
1817 return TRUE;
1818 }
1819 }
1820 #ifndef O_GMP
1821 return PL_error("lcm", 2, NULL, ERR_EVALUATION, ATOM_int_overflow);
1822 #else
1823 promoteToMPZNumber(n1);
1824 promoteToMPZNumber(n2);
1825 case V_MPZ:
1826 r->type = V_MPZ;
1827 mpz_init(r->value.mpz);
1828 mpz_lcm(r->value.mpz, n1->value.mpz, n2->value.mpz);
1829 break;
1830 #endif
1831 default:
1832 assert(0);
1833 }
1834
1835 return TRUE;
1836 }
1837
1838
1839 /* Unary functions requiring double argument */
1840
1841 #define UNAIRY_FLOAT_FUNCTION(name, op) \
1842 static int \
1843 name(Number n1, Number r) \
1844 { if ( !promoteToFloatNumber(n1) ) return FALSE; \
1845 r->value.f = op(n1->value.f); \
1846 r->type = V_FLOAT; \
1847 return check_float(r); \
1848 }
1849
1850 /* Binary functions requiring integer argument */
1851
1852 #ifdef O_GMP
1853 #define BINAIRY_INT_FUNCTION(name, plop, op, mpop) \
1854 static int \
1855 name(Number n1, Number n2, Number r) \
1856 { if ( !toIntegerNumber(n1, 0) ) \
1857 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n1); \
1858 if ( !toIntegerNumber(n2, 0) ) \
1859 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n2); \
1860 if ( !same_type_numbers(n1, n2) ) \
1861 return FALSE; \
1862 switch(n1->type) \
1863 { case V_INTEGER: \
1864 r->value.i = n1->value.i op n2->value.i; \
1865 r->type = V_INTEGER; \
1866 succeed; \
1867 case V_MPZ: \
1868 r->type = V_MPZ; \
1869 mpz_init(r->value.mpz); \
1870 mpop(r->value.mpz, n1->value.mpz, n2->value.mpz); \
1871 succeed; \
1872 default: \
1873 assert(0); \
1874 fail; \
1875 } \
1876 }
1877
1878 #else /*O_GMP*/
1879
1880 #define BINAIRY_INT_FUNCTION(name, plop, op, mpop) \
1881 static int \
1882 name(Number n1, Number n2, Number r) \
1883 { if ( !toIntegerNumber(n1, 0) ) \
1884 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n1); \
1885 if ( !toIntegerNumber(n2, 0) ) \
1886 return PL_error(plop, 2, NULL, ERR_AR_TYPE, ATOM_integer, n2); \
1887 if ( !same_type_numbers(n1, n2) ) \
1888 return FALSE; \
1889 switch(n1->type) \
1890 { case V_INTEGER: \
1891 r->value.i = n1->value.i op n2->value.i; \
1892 r->type = V_INTEGER; \
1893 succeed; \
1894 default: \
1895 assert(0); \
1896 fail; \
1897 } \
1898 }
1899 #endif /*O_GMP*/
1900
1901 #define BINAIRY_FLOAT_FUNCTION(name, func) \
1902 static int \
1903 name(Number n1, Number n2, Number r) \
1904 { if ( !promoteToFloatNumber(n1) || \
1905 !promoteToFloatNumber(n2) ) return FALSE; \
1906 r->value.f = func(n1->value.f, n2->value.f); \
1907 r->type = V_FLOAT; \
1908 return check_float(r); \
1909 }
1910
UNAIRY_FLOAT_FUNCTION(ar_sin,sin)1911 UNAIRY_FLOAT_FUNCTION(ar_sin, sin)
1912 UNAIRY_FLOAT_FUNCTION(ar_cos, cos)
1913 UNAIRY_FLOAT_FUNCTION(ar_tan, tan)
1914 UNAIRY_FLOAT_FUNCTION(ar_sinh, sinh)
1915 UNAIRY_FLOAT_FUNCTION(ar_cosh, cosh)
1916 UNAIRY_FLOAT_FUNCTION(ar_tanh, tanh)
1917 UNAIRY_FLOAT_FUNCTION(ar_asinh, asinh)
1918 UNAIRY_FLOAT_FUNCTION(ar_acosh, acosh)
1919 UNAIRY_FLOAT_FUNCTION(ar_atanh, atanh)
1920 UNAIRY_FLOAT_FUNCTION(ar_atan, atan)
1921 UNAIRY_FLOAT_FUNCTION(ar_exp, exp)
1922 UNAIRY_FLOAT_FUNCTION(ar_erf, erf)
1923 UNAIRY_FLOAT_FUNCTION(ar_erfc, erfc)
1924
1925 BINAIRY_FLOAT_FUNCTION(ar_atan2, atan2)
1926
1927 BINAIRY_INT_FUNCTION(ar_disjunct, "\\/", |, mpz_ior)
1928 BINAIRY_INT_FUNCTION(ar_conjunct, "/\\", &, mpz_and)
1929 BINAIRY_INT_FUNCTION(ar_xor, "xor", ^, mpz_xor)
1930
1931 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1932 ar_pow() is exponentiation. We do this in integers if possible. However,
1933 GMP crashes the entire process by calling abort() if it discovers that
1934 the resulting value will not fit in the address space. Therefore we
1935 estimage the size and verify that it will in on the global stack limit.
1936
1937 I doubt that the computation is accurate, but it is highly unlikely we
1938 won't run out of memory if we create an integer that requires almost the
1939 complete stack size. It is also not a problem if we underestimate a bit
1940 as long as the result fits in the address space. In that case, the
1941 normal overflow handling will nicely generate a resource error.
1942 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1943
1944 #ifdef O_GMP
1945 static void
1946 mpz_set_num(mpz_t mpz, Number n)
1947 { switch ( n->type )
1948 { case V_MPZ:
1949 mpz_set(mpz, n->value.mpz);
1950 break;
1951 case V_INTEGER:
1952 mpz_init_set_si64(mpz, n->value.i);
1953 break;
1954 default:
1955 assert(0);
1956 }
1957 }
1958
1959 static int
get_int_exponent(Number n,unsigned long * expp)1960 get_int_exponent(Number n, unsigned long *expp)
1961 { long exp;
1962 int64_t i;
1963
1964 switch(n->type)
1965 { case V_INTEGER:
1966 i = n->value.i;
1967 break;
1968 case V_MPZ:
1969 if ( !mpz_to_int64(n->value.mpz, &i) )
1970 return int_too_big();
1971 break;
1972 default:
1973 assert(0);
1974 return FALSE;
1975 }
1976
1977 exp = (long)i;
1978 #if SIZEOF_LONG < 8
1979 if ( (int64_t)exp != i )
1980 return int_too_big();
1981 #endif
1982
1983 if ( exp >= 0 )
1984 *expp = (unsigned long)exp;
1985 else if ( -exp != exp )
1986 *expp = (unsigned long)-exp;
1987 else
1988 return int_too_big();
1989
1990 return TRUE;
1991 }
1992
1993 /* cond_minus_pow() handles rounding mode issues calculating pow with
1994 negative base float have to reverse to_positive and to_negative.
1995 */
1996
1997 static double
cond_minus_pow(double base,double exp)1998 cond_minus_pow(double base, double exp)
1999 { double res;
2000
2001 if ( base < 0 )
2002 { switch( fegetround() )
2003 { case FE_UPWARD:
2004 fesetround(FE_DOWNWARD);
2005 res = -pow(-base,exp);
2006 fesetround(FE_UPWARD);
2007 break;
2008 case FE_DOWNWARD:
2009 fesetround(FE_UPWARD);
2010 res = -pow(-base,exp);
2011 fesetround(FE_DOWNWARD);
2012 break;
2013 default:
2014 res = -pow(-base,exp);
2015 }
2016 } else
2017 { res = pow(base,exp);
2018 }
2019
2020 return res;
2021 }
2022
2023
2024 static int
ar_smallint(Number n,int * i)2025 ar_smallint(Number n, int *i)
2026 { switch(n->type)
2027 { case V_INTEGER:
2028 if ( n->value.i >= -1 && n->value.i <= 1 )
2029 { *i = n->value.i;
2030 return TRUE;
2031 }
2032 return FALSE;
2033 case V_MPZ:
2034 if ( mpz_cmp_si(n->value.mpz, -1L) >= 0 &&
2035 mpz_cmp_si(n->value.mpz, 1L) <= 0 )
2036 { *i = mpz_get_si(n->value.mpz);
2037 return TRUE;
2038 }
2039 return FALSE;
2040 default:
2041 assert(0);
2042 return FALSE;
2043 }
2044 }
2045
2046 #endif /*O_GMP*/
2047
2048 static inline int
sign_f(double f)2049 sign_f(double f)
2050 { return
2051 f < 0 ? -1 :
2052 f > 0 ? 1 :
2053 0 ; /* sign_f(NaN) = 0 */
2054 }
2055
2056 static int
ar_pow(Number n1,Number n2,Number r)2057 ar_pow(Number n1, Number n2, Number r)
2058 { int zero_div_sign;
2059 int exp_sign;
2060 #ifdef O_GMP
2061 unsigned long exp;
2062 int exp_nan;
2063 int n1_val;
2064
2065 if ( n2->type == V_FLOAT )
2066 { exp_nan = isnan(n2->value.f);
2067 exp_sign = sign_f(n2->value.f);
2068 } else
2069 { exp_nan = FALSE;
2070 exp_sign = ar_sign_i(n2);
2071 }
2072 r->type = V_INTEGER; /* for all special cases */
2073
2074 if ( exp_sign == 0 && !exp_nan ) /* test for X**0 */
2075 { r->value.i = 1;
2076 return TRUE;
2077 }
2078
2079 if ( intNumber(n1) && ar_smallint(n1, &n1_val) )
2080 { if ( n1_val == 1 ) /* 1**X => 1 */
2081 { r->value.i = 1;
2082 return TRUE;
2083 }
2084 if ( n1_val == 0 && !exp_nan ) /* n1==0, non-zero(n2) */
2085 { if ( exp_sign > 0)
2086 { r->value.i = 0; /* positive exp => 0 */
2087 return TRUE;
2088 } else /* negative exp => zero_div */
2089 { return check_zero_div(ar_sign_i(n1), r, "**", 2);
2090 }
2091 }
2092 if ( n1_val == -1 && intNumber(n2) ) /* check n1 == -1 */
2093 { r->value.i = ar_even(n2) ? 1 : -1;
2094 return TRUE;
2095 }
2096 }
2097
2098 if ( intNumber(n1) && intNumber(n2) )
2099 { if ( !get_int_exponent(n2, &exp) )
2100 return FALSE;
2101
2102 if ( exp_sign < 0 )
2103 { GET_LD
2104
2105 if ( truePrologFlag(PLFLAG_RATIONAL) )
2106 { promoteToMPQNumber(n1);
2107 goto int_pow_neg_int;
2108 }
2109 goto doreal;
2110 }
2111
2112 { GET_LD /* estimate the size, see above */
2113 size_t op1_bits;
2114 int64_t r_bits;
2115
2116 switch(n1->type)
2117 { case V_INTEGER:
2118 op1_bits = MSB64(n1->value.i);
2119 break;
2120 case V_MPZ:
2121 op1_bits = mpz_sizeinbase(n1->value.mpz, 2);
2122 break;
2123 default:
2124 assert(0);
2125 fail;
2126 }
2127
2128 if ( !( mul64(op1_bits, exp, &r_bits) &&
2129 r_bits/8 < (int64_t)globalStackLimit()
2130 ) )
2131 return int_too_big();
2132 }
2133
2134 r->type = V_MPZ;
2135 mpz_init(r->value.mpz);
2136
2137 switch(n1->type)
2138 { case V_INTEGER:
2139 if ( n1->value.i >= 0L && n1->value.i <= LONG_MAX )
2140 { mpz_ui_pow_ui(r->value.mpz, (unsigned long)n1->value.i, exp);
2141 succeed;
2142 } else
2143 { promoteToMPZNumber(n1);
2144 /*FALLTHROUGH*/
2145 }
2146 case V_MPZ:
2147 mpz_pow_ui(r->value.mpz, n1->value.mpz, exp);
2148 succeed;
2149 default:
2150 assert(0);
2151 fail;
2152 }
2153 } /* end if ( intNumber(n1) && intNumber(n2) ) */
2154
2155 if ( n1->type == V_MPQ && intNumber(n2) )
2156 { number nr, nd, nrp, ndp, nexp;
2157
2158 if ( !get_int_exponent(n2, &exp) )
2159 return FALSE;
2160
2161 if ( exp_sign == 0 )
2162 { r->type = V_INTEGER;
2163 r->value.i = 1;
2164 return TRUE;
2165 }
2166
2167 int_pow_neg_int:
2168 nexp.type = V_INTEGER;
2169 nexp.value.i = exp;
2170
2171 nr.type = V_MPZ;
2172 nr.value.mpz[0] = mpq_numref(n1->value.mpq)[0];
2173 nr.value.mpz->_mp_alloc = 0; /* read-only */
2174 nd.type = V_MPZ;
2175 nd.value.mpz[0] = mpq_denref(n1->value.mpq)[0];
2176 nd.value.mpz->_mp_alloc = 0;
2177
2178 if ( ar_pow(&nr, &nexp, &nrp) &&
2179 ar_pow(&nd, &nexp, &ndp) )
2180 { r->type = V_MPQ;
2181 mpq_init(r->value.mpq);
2182 if ( exp_sign > 0 )
2183 { mpz_set_num(mpq_numref(r->value.mpq), &nrp);
2184 mpz_set_num(mpq_denref(r->value.mpq), &ndp);
2185 } else
2186 { mpz_set_num(mpq_numref(r->value.mpq), &ndp);
2187 mpz_set_num(mpq_denref(r->value.mpq), &nrp);
2188 }
2189 mpq_canonicalize(r->value.mpq);
2190
2191 clearNumber(&nrp);
2192 clearNumber(&ndp);
2193
2194 return check_mpq(r);
2195 }
2196
2197 clearNumber(&nrp);
2198 clearNumber(&ndp);
2199 } /* end MPQ^int */
2200
2201 if ( n2->type == V_MPQ ) /* X ^ rat */
2202 { long r_den;
2203
2204 if ( exp_sign == -1 )
2205 mpz_neg(mpq_numref(n2->value.mpq), mpq_numref(n2->value.mpq));
2206
2207 r_den = mpz_get_ui(mpq_denref(n2->value.mpq));
2208
2209 switch (n1->type)
2210 { case V_INTEGER:
2211 { mpz_init_set_si(r->value.mpz,n1->value.i);
2212 goto int_to_rat;
2213 }
2214 case V_MPZ:
2215 { mpz_init_set(r->value.mpz,n1->value.mpz);
2216
2217 int_to_rat:
2218 r->type = V_MPZ; /* int ^ rat */
2219 /* neg ^ int/even is undefined */
2220 if ( mpz_sgn(r->value.mpz) == -1 && !(r_den & 1))
2221 { mpz_clear(r->value.mpz);
2222 r->type = V_FLOAT;
2223 r->value.f = const_nan;
2224 return check_float(r);
2225 }
2226
2227 if ( mpz_root(r->value.mpz,r->value.mpz,r_den))
2228 { unsigned long r_num = mpz_get_ui(mpq_numref(n2->value.mpq));
2229
2230 if ( r_num > LONG_MAX ) /* numerator exceeds mpz_pow_ui range */
2231 { mpz_clear(r->value.mpz);
2232 if ( promoteToFloatNumber(n1) )
2233 goto doreal_mpq;
2234 else return FALSE;
2235 } else
2236 { mpz_pow_ui(r->value.mpz,r->value.mpz,r_num);
2237
2238 if (exp_sign == -1) /* create mpq=1/r->value */
2239 { mpz_t tempz;
2240
2241 mpz_init_set(tempz,r->value.mpz);
2242 mpz_clear(r->value.mpz);
2243 r->type = V_MPQ;
2244 mpq_init(r->value.mpq);
2245 mpq_set_z(r->value.mpq,tempz);
2246 mpq_inv(r->value.mpq,r->value.mpq);
2247 mpz_clear(tempz);
2248 return check_mpq(r);
2249 } else
2250 { return TRUE;
2251 }
2252 }
2253 } else /* root inexact */
2254 { mpz_clear(r->value.mpz);
2255 if ( promoteToFloatNumber(n1) )
2256 goto doreal_mpq;
2257 else return FALSE;
2258 }
2259 break;
2260 }
2261 case V_MPQ:
2262 { int rat_result;
2263 unsigned long r_num;
2264
2265 r->type = V_MPQ;
2266 mpq_init(r->value.mpq);
2267 mpq_set(r->value.mpq, n1->value.mpq);
2268
2269 /* neg ^ int / even */
2270 if ( (mpq_sgn(r->value.mpq) == -1 ) && !(r_den & 1))
2271 { mpq_clear(r->value.mpq);
2272 r->type = V_FLOAT;
2273 r->value.f = const_nan;
2274 return check_float(r);
2275 }
2276
2277 rat_result = ( mpz_root(mpq_numref(r->value.mpq),
2278 mpq_numref(r->value.mpq),r_den) &&
2279 mpz_root(mpq_denref(r->value.mpq),
2280 mpq_denref(r->value.mpq),r_den)
2281 );
2282
2283 r_num = mpz_get_ui(mpq_numref(n2->value.mpq));
2284 if ( rat_result && (r_num < LONG_MAX) ) /* base = base^P */
2285 { mpz_pow_ui(mpq_numref(r->value.mpq),mpq_numref(r->value.mpq),r_num);
2286 mpz_pow_ui(mpq_denref(r->value.mpq),mpq_denref(r->value.mpq),r_num);
2287
2288 if ( exp_sign == -1 )
2289 mpq_inv(r->value.mpq,r->value.mpq);
2290
2291 return check_mpq(r);
2292 } else /* exponent out of range for mpz_pow_ui */
2293 { mpq_clear(r->value.mpq);
2294
2295 if ( promoteToFloatNumber(n1) )
2296 goto doreal_mpq;
2297 else
2298 return FALSE;
2299 }
2300 assert(0);
2301 }
2302 case V_FLOAT:
2303 { if ( n1->value.f == 0.0 ) goto doreal; /* general case of 0.0**X */
2304 if ( n1->value.f < 0 && !( r_den & 1 ))
2305 { r->value.f = const_nan; /* negative base, even denominator */
2306 } else
2307 {
2308 doreal_mpq:
2309 mpq_init(r->value.mpq);
2310 mpq_set_ui(r->value.mpq,1,mpz_get_ui(mpq_denref(n2->value.mpq)));
2311 double dexp = mpq_get_d(r->value.mpq); /* float(1/n2.den) */
2312 mpq_clear(r->value.mpq);
2313 r->value.f = pow(cond_minus_pow(n1->value.f, dexp),
2314 mpz_get_ui(mpq_numref(n2->value.mpq)));
2315 if ( exp_sign == -1 )
2316 r->value.f = 1.0/r->value.f;
2317 }
2318
2319 r->type = V_FLOAT;
2320 return check_float(r);
2321 }
2322 } /* end switch (n1->type) */
2323 assert(0);
2324 } /* end if ( n2->type == V_MPQ ) */
2325
2326 doreal:
2327 #endif /*O_GMP*/
2328 zero_div_sign = ( (n2->type == V_INTEGER) && (!ar_even(n2)) &&
2329 signbit(n1->value.f) ) ? -1 : 1;
2330
2331 if ( !promoteToFloatNumber(n1) ||
2332 !promoteToFloatNumber(n2) )
2333 return FALSE;
2334
2335 #ifndef O_GMP
2336 exp_sign = sign_f(n2->value.f);
2337 #endif
2338
2339 if ( n1->value.f == 0.0 && exp_sign == -1 )
2340 return check_zero_div(zero_div_sign, r, "**", 2);
2341 if ( n1->value.f == 1.0 )
2342 r->value.f = 1.0; /* broken on e.g., mipsel */
2343 else
2344 r->value.f = pow(n1->value.f, n2->value.f);
2345 r->type = V_FLOAT;
2346
2347 return check_float(r);
2348 }
2349
2350 static int
ar_powm(Number base,Number exp,Number mod,Number r)2351 ar_powm(Number base, Number exp, Number mod, Number r)
2352 {
2353 if ( !intNumber(base) )
2354 PL_error("powm", 3, NULL, ERR_AR_TYPE, ATOM_integer, base);
2355 if ( !intNumber(exp) )
2356 PL_error("powm", 3, NULL, ERR_AR_TYPE, ATOM_integer, exp);
2357 if ( !intNumber(exp) )
2358 PL_error("powm", 3, NULL, ERR_AR_TYPE, ATOM_integer, exp);
2359
2360 #ifdef O_GMP
2361 promoteToMPZNumber(base);
2362 promoteToMPZNumber(exp);
2363 promoteToMPZNumber(mod);
2364
2365 if ( ar_sign_i(base) < 0 ) return notLessThanZero("powm", 3, base);
2366 if ( ar_sign_i(exp) < 0 ) return notLessThanZero("powm", 3, exp);
2367 if ( ar_sign_i(mod) <= 0 ) return mustBePositive("powm", 3, mod);
2368
2369 r->type = V_MPZ;
2370 mpz_init(r->value.mpz);
2371 mpz_powm(r->value.mpz, base->value.mpz, exp->value.mpz, mod->value.mpz);
2372
2373 succeed;
2374 #else
2375 return PL_error("powm", 3, "requires unbounded arithmetic (GMP) support",
2376 ERR_NOT_IMPLEMENTED, "powm/3");
2377 #endif
2378 }
2379
2380 #if 0
2381 /* These tests originate from the days that float errors used
2382 * to be signalling on many systems. Nowadays this is no longer
2383 * the case. We leave the code in for just-in-case.
2384 */
2385 #define AR_UNDEFINED_IF(func, arity, test, r) \
2386 if ( test ) \
2387 { GET_LD \
2388 if ( LD->arith.f.flags & FLT_UNDEFINED ) \
2389 { r->type = V_FLOAT; \
2390 r->value.f = const_nan; \
2391 return TRUE; \
2392 } else \
2393 { return PL_error(func, arity, NULL, ERR_AR_UNDEF); \
2394 } \
2395 }
2396 #define AR_DIV_ZERO_IF(func, arity, n, d, r) \
2397 if ( d == 0.0 ) \
2398 { GET_LD \
2399 if ( LD->arith.f.flags & FLT_ZERO_DIV ) \
2400 { r->type = V_FLOAT; \
2401 r->value.f = signbit(n) == signbit(d) \
2402 ? const_inf \
2403 : const_neg_inf; \
2404 return TRUE; \
2405 } else \
2406 { return PL_error(func, arity, NULL, ERR_DIV_BY_ZERO);\
2407 } \
2408 }
2409 #else
2410 #define AR_UNDEFINED_IF(func, arity, test, r) (void)0
2411 #define AR_DIV_ZERO_IF(func, arity, n, d, r) (void)0
2412 #endif
2413
2414 static int
ar_sqrt(Number n1,Number r)2415 ar_sqrt(Number n1, Number r)
2416 { if ( !promoteToFloatNumber(n1) )
2417 return FALSE;
2418 AR_UNDEFINED_IF("sqrt", 1, n1->value.f < 0, r);
2419 r->value.f = sqrt(n1->value.f);
2420 r->type = V_FLOAT;
2421
2422 return check_float(r);
2423 }
2424
2425
2426 static int
ar_asin(Number n1,Number r)2427 ar_asin(Number n1, Number r)
2428 { if ( !promoteToFloatNumber(n1) )
2429 return FALSE;
2430 AR_UNDEFINED_IF("asin", 1, n1->value.f < -1.0 || n1->value.f > 1.0, r);
2431 r->value.f = asin(n1->value.f);
2432 r->type = V_FLOAT;
2433
2434 return check_float(r);
2435 }
2436
2437
2438 static int
ar_acos(Number n1,Number r)2439 ar_acos(Number n1, Number r)
2440 { if ( !promoteToFloatNumber(n1) )
2441 return FALSE;
2442 AR_UNDEFINED_IF("ascos", 1, n1->value.f < -1.0 || n1->value.f > 1.0, r);
2443 r->value.f = acos(n1->value.f);
2444 r->type = V_FLOAT;
2445
2446 return check_float(r);
2447 }
2448
2449
2450 static int
ar_log(Number n1,Number r)2451 ar_log(Number n1, Number r)
2452 { if ( !promoteToFloatNumber(n1) )
2453 return FALSE;
2454 AR_UNDEFINED_IF("log", 1, n1->value.f <= 0.0 , r);
2455 r->value.f = log(n1->value.f);
2456 r->type = V_FLOAT;
2457
2458 return check_float(r);
2459 }
2460
2461
2462 static int
ar_log10(Number n1,Number r)2463 ar_log10(Number n1, Number r)
2464 { if ( !promoteToFloatNumber(n1) )
2465 return FALSE;
2466 AR_UNDEFINED_IF("log10", 1, n1->value.f <= 0.0, r);
2467 r->value.f = log10(n1->value.f);
2468 r->type = V_FLOAT;
2469
2470 return check_float(r);
2471 }
2472
2473
2474 static int
ar_lgamma(Number n1,Number r)2475 ar_lgamma(Number n1, Number r) // custom lgamma() to ensure positive inf
2476 { if ( !promoteToFloatNumber(n1) )
2477 return FALSE;
2478 r->value.f = (n1->value.f <= 0.0) ? INFINITY : lgamma(n1->value.f);
2479 r->type = V_FLOAT;
2480
2481 return check_float(r);
2482 }
2483
2484 /* IntExpr1 // IntExpr2
2485
2486 Integer division. Defined by ISO core standard as rnd(X,Y), where the
2487 direction of the rounding is conform the flag integer_rounding_function,
2488 which is one of =toward_zero= or =down=.
2489
2490 The implementation below rounds according to the C-compiler. This is not
2491 desirable, but I understand that as of C99, this is towards zero and
2492 this is precisely what we want to make this different from div/2. As we
2493 need C99 for the wide-character support anyway, we should be fairly
2494 safe.
2495 */
2496
2497 static int
ar_tdiv(Number n1,Number n2,Number r)2498 ar_tdiv(Number n1, Number n2, Number r)
2499 { if ( !toIntegerNumber(n1, 0) )
2500 return PL_error("//", 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
2501 if ( !toIntegerNumber(n2, 0) )
2502 return PL_error("//", 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
2503
2504 #ifdef O_GMP
2505 if ( n1->type == V_INTEGER && n2->type == V_INTEGER )
2506 #endif
2507 { if ( n2->value.i == 0 )
2508 return PL_error("//", 2, NULL, ERR_DIV_BY_ZERO);
2509
2510 if ( !(n2->value.i == -1 && n1->value.i == PLMININT) )
2511 { r->value.i = n1->value.i / n2->value.i;
2512 r->type = V_INTEGER;
2513
2514 succeed;
2515 }
2516 }
2517
2518 #ifdef O_GMP
2519 promoteToMPZNumber(n1);
2520 promoteToMPZNumber(n2);
2521
2522 if ( mpz_sgn(n2->value.mpz) == 0 )
2523 return PL_error("//", 2, NULL, ERR_DIV_BY_ZERO);
2524
2525 r->type = V_MPZ;
2526 mpz_init(r->value.mpz);
2527 if ( (-3 / 2) == -1 )
2528 mpz_tdiv_q(r->value.mpz, n1->value.mpz, n2->value.mpz);
2529 else
2530 mpz_fdiv_q(r->value.mpz, n1->value.mpz, n2->value.mpz);
2531
2532 succeed;
2533 #else
2534 return PL_error("//", 2, NULL, ERR_EVALUATION, ATOM_int_overflow);
2535 #endif
2536 }
2537
2538
2539 /** div(IntExpr1, IntExpr2)
2540
2541 Result is rnd_i(IntExpr1/IntExpr2), rounded towards -infinity
2542 */
2543
2544 static int
ar_div(Number n1,Number n2,Number r)2545 ar_div(Number n1, Number n2, Number r)
2546 { if ( !toIntegerNumber(n1, 0) )
2547 return PL_error("div", 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
2548 if ( !toIntegerNumber(n2, 0) )
2549 return PL_error("div", 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
2550
2551 #ifdef O_GMP
2552 if ( n1->type == V_INTEGER && n2->type == V_INTEGER )
2553 #endif
2554 { if ( n2->value.i == 0 )
2555 return PL_error("div", 2, NULL, ERR_DIV_BY_ZERO);
2556
2557 if ( !(n2->value.i == -1 && n1->value.i == PLMININT) )
2558 { r->value.i = n1->value.i / n2->value.i;
2559 if ((n1->value.i > 0) != (n2->value.i > 0) &&
2560 n1->value.i % n2->value.i != 0)
2561 --r->value.i;
2562 r->type = V_INTEGER;
2563
2564 succeed;
2565 }
2566 }
2567
2568 #ifdef O_GMP
2569 promoteToMPZNumber(n1);
2570 promoteToMPZNumber(n2);
2571
2572 if ( mpz_sgn(n2->value.mpz) == 0 )
2573 return PL_error("div", 2, NULL, ERR_DIV_BY_ZERO);
2574
2575 r->type = V_MPZ;
2576 mpz_init(r->value.mpz);
2577 mpz_fdiv_q(r->value.mpz, n1->value.mpz, n2->value.mpz);
2578
2579 succeed;
2580 #else
2581 return PL_error("div", 2, NULL, ERR_EVALUATION, ATOM_int_overflow);
2582 #endif
2583 }
2584
2585 /* Broken, at least on SunOS 5.11, gcc 4.8. No clue under what conditions.
2586 The results of configure and final linking differ. Anyway, just doing
2587 our own is most likely the safe solution.
2588 */
2589 #ifdef __sun
2590 #undef HAVE_SIGNBIT
2591 #endif
2592
2593 #ifndef HAVE_SIGNBIT /* differs for -0.0 */
2594 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2595 signbit() and copysign() are part of C99. These should be provided by
2596 most C compilers, but Microsoft decided not to adopt C99 (it is now
2597 2012).
2598
2599 Note that there is no autoconf support to verify that floats conform to
2600 the IEE754 representation, but they typically do these days. See
2601 http://www.gnu.org/software/autoconf/manual/autoconf-2.67/html_node/Floating-Point-Portability.html
2602 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2603 #define IEEE754 1
2604
2605 #ifdef IEEE754
2606 static inline int
signbit(double f)2607 signbit(double f)
2608 { union
2609 { double f;
2610 int64_t i;
2611 } v;
2612
2613 v.f = f;
2614 return v.i < 0;
2615 }
2616
2617 #ifndef copysign
2618 double
copysign(double x,double y)2619 copysign(double x, double y)
2620 { union { double f; uint64_t i; } ux, uy;
2621 const uint64_t smask = (uint64_t)1<<(sizeof(uint64_t)*8-1);
2622
2623 ux.f = x;
2624 uy.f = y;
2625 ux.i &= ~smask;
2626 ux.i |= (uy.i & smask);
2627
2628 return ux.f;
2629 }
2630 #endif
2631 #else
2632 #error "Don't know how to support signbit() and copysign()"
2633 #endif
2634 #endif
2635
2636 int
ar_sign_i(Number n1)2637 ar_sign_i(Number n1)
2638 { switch(n1->type)
2639 { case V_INTEGER:
2640 return (n1->value.i < 0 ? -1 : n1->value.i > 0 ? 1 : 0);
2641 #ifdef O_GMP
2642 case V_MPZ:
2643 return mpz_sgn(n1->value.mpz);
2644 case V_MPQ:
2645 return mpq_sgn(n1->value.mpq);
2646 #endif
2647 default:
2648 assert(0);
2649 fail;
2650 }
2651 }
2652
2653 static int
ar_sign(Number n1,Number r)2654 ar_sign(Number n1, Number r)
2655 { if ( n1->type == V_FLOAT )
2656 { r->value.f = isnan(n1->value.f) ? const_nan :
2657 (n1->value.f < 0 ? -1.0 : n1->value.f > 0.0 ? 1.0 : 0.0);
2658 r->type = V_FLOAT;
2659 } else
2660 { r->value.i = ar_sign_i(n1);
2661 r->type = V_INTEGER;
2662 }
2663
2664 succeed;
2665 }
2666
2667
2668 int
ar_signbit(Number n)2669 ar_signbit(Number n)
2670 { switch(n->type)
2671 { case V_INTEGER:
2672 return n->value.i < 0 ? -1 : 1;
2673 #ifdef O_GMP
2674 case V_MPZ:
2675 { int i = mpz_sgn(n->value.mpz);
2676 return i < 0 ? -1 : 1;
2677 }
2678 case V_MPQ:
2679 { int i = mpq_sgn(n->value.mpq);
2680 return i < 0 ? -1 : 1;
2681 }
2682 #endif
2683 case V_FLOAT:
2684 return signbit(n->value.f) ? -1 : 1;
2685 default:
2686 assert(0);
2687 return 0;
2688 }
2689 }
2690
2691
2692 static int
ar_copysign(Number n1,Number n2,Number r)2693 ar_copysign(Number n1, Number n2, Number r)
2694 {
2695 if ( n1->type == V_FLOAT && n2->type == V_FLOAT )
2696 { if ( isnan(n1->value.f) )
2697 r->value.f = const_nan;
2698 else
2699 r->value.f = copysign(n1->value.f, n2->value.f);
2700 r->type = V_FLOAT;
2701 } else
2702 { if ( ar_signbit(n1) != ar_signbit(n2) )
2703 return ar_u_minus(n1, r);
2704 else
2705 cpNumber(r, n1);
2706 }
2707
2708 return TRUE;
2709 }
2710
2711
2712 static int
ar_nexttoward(Number n1,Number n2,Number r)2713 ar_nexttoward(Number n1, Number n2, Number r)
2714 { if ( promoteToFloatNumber(n1) &&
2715 promoteToFloatNumber(n2) )
2716 { if ( n1->type == V_FLOAT && n2->type == V_FLOAT )
2717 { r->value.f = nexttoward(n1->value.f, n2->value.f);
2718 r->type = V_FLOAT;
2719
2720 return check_float(r);
2721 }
2722 }
2723
2724 return FALSE;
2725 }
2726
2727 static int
set_roundtoward(Word p,Number old ARG_LD)2728 set_roundtoward(Word p, Number old ARG_LD)
2729 { deRef(p);
2730
2731 old->type = V_INTEGER;
2732 old->value.i = fegetround();
2733
2734 if ( *p == ATOM_to_nearest )
2735 fesetround(FE_TONEAREST);
2736 else if ( *p == ATOM_to_positive )
2737 fesetround(FE_UPWARD);
2738 else if ( *p == ATOM_to_negative )
2739 fesetround(FE_DOWNWARD);
2740 else if ( *p == ATOM_to_zero )
2741 fesetround(FE_TOWARDZERO);
2742 else if ( isAtom(*p) )
2743 return PL_error(NULL, 0, NULL, ERR_PTR_DOMAIN, ATOM_round, p);
2744 else
2745 return PL_error(NULL, 0, NULL, ERR_PTR_TYPE, ATOM_atom, p);
2746
2747 return TRUE;
2748 }
2749
2750 static int
ar_roundtoward(Number n1,Number n2,Number r)2751 ar_roundtoward(Number n1, Number n2, Number r)
2752 { cpNumber(r, n1);
2753
2754 assert(n2->type == V_INTEGER);
2755 fesetround(n2->value.i);
2756
2757 return TRUE;
2758 }
2759
2760
2761 static int
ar_rem(Number n1,Number n2,Number r)2762 ar_rem(Number n1, Number n2, Number r)
2763 { if ( !toIntegerNumber(n1, 0) )
2764 return PL_error("rem", 2, NULL, ERR_AR_TYPE, ATOM_integer, n1);
2765 if ( !toIntegerNumber(n2, 0) )
2766 return PL_error("rem", 2, NULL, ERR_AR_TYPE, ATOM_integer, n2);
2767
2768 if ( !same_type_numbers(n1, n2) )
2769 return FALSE;
2770 switch(n1->type)
2771 { case V_INTEGER:
2772 if ( n2->value.i == 0 )
2773 return PL_error("rem", 2, NULL, ERR_DIV_BY_ZERO);
2774
2775 if ( n2->value.i != -1 || n1->value.i != INT64_MIN )
2776 r->value.i = n1->value.i % n2->value.i;
2777 else
2778 r->value.i = 0;
2779 r->type = V_INTEGER;
2780
2781 break;
2782 #ifdef O_GMP
2783 case V_MPZ:
2784 { if ( mpz_sgn(n2->value.mpz) == 0 )
2785 return PL_error("rem", 2, NULL, ERR_DIV_BY_ZERO);
2786
2787 r->type = V_MPZ;
2788 mpz_init(r->value.mpz);
2789 mpz_tdiv_r(r->value.mpz, n1->value.mpz, n2->value.mpz);
2790 break;
2791 }
2792 #endif
2793
2794 default:
2795 assert(0);
2796 fail;
2797 }
2798
2799 succeed;
2800 }
2801
2802
2803 #ifdef O_GMP
2804 static int
ar_rational(Number n1,Number r)2805 ar_rational(Number n1, Number r)
2806 { cpNumber(r, n1);
2807 return promoteToMPQNumber(r);
2808 }
2809
2810 static int
ar_numerator(Number n1,Number r)2811 ar_numerator(Number n1, Number r)
2812 { if ( intNumber(n1) )
2813 { cpNumber(r, n1);
2814 return TRUE;
2815 }
2816 if ( n1->type == V_MPQ )
2817 { r->type = V_MPZ;
2818 mpz_init(r->value.mpz);
2819 mpz_set(r->value.mpz, mpq_numref(n1->value.mpq));
2820 return TRUE;
2821 }
2822
2823 return PL_error("numerator", 1, NULL, ERR_AR_TYPE, ATOM_rational, n1);
2824 }
2825
2826
2827 static int
ar_denominator(Number n1,Number r)2828 ar_denominator(Number n1, Number r)
2829 { if ( intNumber(n1) )
2830 { r->type = V_INTEGER;
2831 r->value.i = 1;
2832 return TRUE;
2833 }
2834 if ( n1->type == V_MPQ )
2835 { r->type = V_MPZ;
2836 mpz_init(r->value.mpz);
2837 mpz_set(r->value.mpz, mpq_denref(n1->value.mpq));
2838 return TRUE;
2839 }
2840
2841 return PL_error("denominator", 1, NULL, ERR_AR_TYPE, ATOM_rational, n1);
2842 }
2843
2844
2845 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
2846 A is rationalize(Float)
2847
2848 Introduced on the suggestion of Richard O'Keefe after the Common Lisp
2849 standard. The algorithm is taken from figure 3 in ``A Rational Rotation
2850 Method for Robust Geometric Algorithms'' by John Canny, Bruce Donald and
2851 Eugene K. Ressler. Found at
2852
2853 http://www.cs.dartmouth.edu/~brd/papers/rotations-scg92.pdf
2854
2855 (*) Comment by Keri Harris:
2856
2857 The result of p1/q1 is retained in a FP stack register at a higher
2858 precision (80 bits); it is not stored in a variable. This extra
2859 precision skews the results when preforming the subtraction, as one
2860 operand contains extra precision:
2861
2862 (extended double precision) (double precision)
2863 d = p1/q1 - n1->value.f;
2864
2865 Forcing the result of p1/q1 to be stored in a variable produces expected
2866 results with rationalize/1:
2867
2868 volatile double p1_q1 = p1/q1;
2869 d = p1_q1 - n1->value.f;
2870 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2871
2872 static int
ar_rationalize(Number n1,Number r)2873 ar_rationalize(Number n1, Number r)
2874 { switch(n1->type)
2875 { case V_INTEGER:
2876 case V_MPZ:
2877 case V_MPQ:
2878 cpNumber(r, n1);
2879 return TRUE;
2880 case V_FLOAT:
2881 switch(fpclassify(n1->value.f))
2882 { case FP_NAN:
2883 return PL_error(NULL, 0, NULL, ERR_AR_UNDEF);
2884 case FP_INFINITE:
2885 return PL_error(NULL, 0, NULL, ERR_AR_RAT_OVERFLOW);
2886 }
2887
2888 mpq_init(r->value.mpq);
2889 mpq_set_double(r->value.mpq, n1->value.f);
2890 r->type = V_MPQ;
2891 return check_mpq(r);
2892 }
2893
2894 assert(0);
2895 fail;
2896 }
2897
2898
2899 int
ar_rdiv_mpz(Number n1,Number n2,Number r)2900 ar_rdiv_mpz(Number n1, Number n2, Number r)
2901 { if ( mpz_divisible_p(n1->value.mpz, n2->value.mpz) )
2902 { mpz_init(r->value.mpz);
2903 r->type = V_MPZ;
2904 mpz_divexact(r->value.mpz, n1->value.mpz, n2->value.mpz);
2905 } else
2906 { r->type = V_MPQ;
2907 mpq_init(r->value.mpq);
2908 mpz_set(mpq_numref(r->value.mpq), n1->value.mpz);
2909 mpz_set(mpq_denref(r->value.mpq), n2->value.mpz);
2910 mpq_canonicalize(r->value.mpq);
2911 return check_mpq(r);
2912 }
2913
2914 return TRUE;
2915 }
2916
2917
2918 static int
ar_rdiv(Number n1,Number n2,Number r)2919 ar_rdiv(Number n1, Number n2, Number r)
2920 { if ( toIntegerNumber(n1, 0) &&
2921 toIntegerNumber(n2, 0) )
2922 { promoteToMPZNumber(n1);
2923 promoteToMPZNumber(n2);
2924
2925 if ( mpz_sgn(n2->value.mpz) == 0 )
2926 return PL_error("/", 2, NULL, ERR_DIV_BY_ZERO);
2927
2928 return ar_rdiv_mpz(n1, n2, r);
2929 } else if ( ratNumber(n1) && ratNumber(n2) )
2930 { promoteToMPQNumber(n1);
2931 promoteToMPQNumber(n2);
2932
2933 if ( mpz_sgn(mpq_numref(n2->value.mpq)) == 0 )
2934 return PL_error("/", 2, NULL, ERR_DIV_BY_ZERO);
2935
2936 r->type = V_MPQ;
2937 mpq_init(r->value.mpq);
2938 mpq_div(r->value.mpq, n1->value.mpq, n2->value.mpq);
2939 return check_mpq(r);
2940 } else if ( !ratNumber(n1) )
2941 { return PL_error("rdiv", 2, NULL, ERR_AR_TYPE, ATOM_rational, n1);
2942 } else
2943 { return PL_error("rdiv", 2, NULL, ERR_AR_TYPE, ATOM_rational, n2);
2944 }
2945
2946 return TRUE;
2947 }
2948 #endif /*O_GMP*/
2949
2950
2951 static int
ar_divide(Number n1,Number n2,Number r)2952 ar_divide(Number n1, Number n2, Number r)
2953 { GET_LD
2954
2955 if ( (n2->type == V_FLOAT) && isinf(n2->value.f) ) // X/inf
2956 { if (n1->type == V_FLOAT) // float --> signed 0.0 or NaN
2957 { r->type = V_FLOAT;
2958 r->value.f = isfinite(n1->value.f) ? 0.0*sign_f(n1->value.f)*sign_f(n2->value.f) : const_nan;
2959 return check_float(r);
2960 } else // non-float --> 0
2961 { r->type = V_INTEGER;
2962 r->value.i = 0;
2963 succeed;
2964 }
2965 }
2966
2967 if ( !truePrologFlag(PLFLAG_ISO) )
2968 { if ( !same_type_numbers(n1, n2) )
2969 return FALSE;
2970
2971 switch(n1->type)
2972 { case V_INTEGER:
2973 if ( n2->value.i == LL(0) )
2974 return check_zero_div(ar_sign_i(n1), r, "/", 2);
2975 if ( n1->value.i % n2->value.i == 0 )
2976 { r->value.i = n1->value.i / n2->value.i;
2977 r->type = V_INTEGER;
2978 succeed;
2979 }
2980 #ifdef O_GMP
2981 if ( truePrologFlag(PLFLAG_RATIONAL) )
2982 { promoteToMPZNumber(n1);
2983 promoteToMPZNumber(n2);
2984 return ar_rdiv_mpz(n1, n2, r);
2985 }
2986 #endif
2987 break;
2988 #ifdef O_GMP
2989 case V_MPZ:
2990 if ( mpz_sgn(n2->value.mpz) == 0 )
2991 return check_zero_div(ar_sign_i(n1), r, "/", 2);
2992 if ( mpz_divisible_p(n1->value.mpz, n2->value.mpz) )
2993 { mpz_init(r->value.mpz);
2994 r->type = V_MPZ;
2995 mpz_divexact(r->value.mpz, n1->value.mpz, n2->value.mpz);
2996 succeed;
2997 }
2998 if ( truePrologFlag(PLFLAG_RATIONAL) )
2999 return ar_rdiv_mpz(n1, n2, r);
3000 break;
3001 case V_MPQ:
3002 if ( mpq_sgn(n2->value.mpq) == 0 )
3003 return check_zero_div(ar_sign_i(n1), r, "/", 2);
3004 mpq_init(r->value.mpq);
3005 r->type = V_MPQ;
3006 mpq_div(r->value.mpq, n1->value.mpq, n2->value.mpq);
3007 return check_mpq(r);
3008 #endif
3009 case V_FLOAT:
3010 break;
3011 }
3012 } // ! PLAG_ISO
3013
3014 /* TBD: How to handle Q? */
3015 if ( !promoteToFloatNumber(n1) ||
3016 !promoteToFloatNumber(n2) )
3017 return FALSE;
3018
3019 /* separate zero-div case from general overflow, Note: sign_f(nan)=0 */
3020 if ( (n2->value.f == 0.0) && (sign_f(n1->value.f) != 0) )
3021 return check_zero_div((signbit(n1->value.f)==signbit(n2->value.f)) ? 1 : -1,
3022 r, "/", 2);
3023
3024 r->value.f = n1->value.f / n2->value.f;
3025 r->type = V_FLOAT;
3026
3027 return check_float(r);
3028 }
3029
3030
3031 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3032 mul64(int64_t x, int64_t y, int64_t *r)
3033 *r = x*y. Returns TRUE if there is no overflow, FALSE on overflow.
3034 This is pretty complicated. Bart Demoen pointed me at "Revisiting
3035 Overflow in Integer Multiplication" by Ayeas Qawasmeh and Ahmed
3036 Dalalah. They prove nor claim their simple tests are complete
3037 (notably it is not clear whether they may falsily signal overflow).
3038 Their Multiply_using_splitting() looks promising, but is flawed
3039 as the results r2 and r3 must be shifted and split.
3040
3041 They do suggest to multiply and then divide to check the result.
3042 They claim this is not correct as the behaviour of C is undefined
3043 on overflow, but as far as I can tell, it is defined as the truncated
3044 result for the multiplication of _unsigned_ integers. Hence, we do
3045 unsigned multiplication, change back to signed and check using
3046 division.
3047
3048 As division is pretty expensive, we make a quick test to see whether
3049 we are in the danger zone.
3050
3051 Finally, we must avoid INT64_MIN/-1 :-(
3052 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
3053
3054 #define MU64_SAFE_MAX (LL(1)<<30)
3055 #ifndef INT64_MIN
3056 #define INT64_MIN (LL(1)<<63)
3057 #endif
3058
3059 static int
mul64(int64_t x,int64_t y,int64_t * r)3060 mul64(int64_t x, int64_t y, int64_t *r)
3061 { if ( x == LL(0) || y == LL(0) )
3062 { *r = LL(0);
3063 return TRUE;
3064 } else
3065 { int sign;
3066 uint64_t ax, ay;
3067 int64_t prod;
3068
3069 if ( x > LL(0) )
3070 { ax = x;
3071 if ( y > LL(0) )
3072 { ay = y;
3073 sign = 1;
3074 } else
3075 { ay = -(uint64_t)y;
3076 sign = -1;
3077 }
3078 } else
3079 { ax = -(uint64_t)x;
3080 if ( y > LL(0) )
3081 { ay = y;
3082 sign = -1;
3083 } else
3084 { ay = -y;
3085 sign = 1;
3086 }
3087 }
3088
3089 prod = (int64_t)(ax*ay);
3090 if ( sign < 0 )
3091 prod = -(uint64_t)prod;
3092 if ( (ax < MU64_SAFE_MAX && ay < MU64_SAFE_MAX) )
3093 { *r = prod;
3094 return TRUE;
3095 }
3096 if ( !(y==LL(-1) && prod == INT64_MIN) && prod/y == x )
3097 { *r = prod;
3098 return TRUE;
3099 }
3100
3101 return FALSE;
3102 }
3103 }
3104
3105
3106 int
ar_mul(Number n1,Number n2,Number r)3107 ar_mul(Number n1, Number n2, Number r)
3108 { if ( !same_type_numbers(n1, n2) )
3109 return FALSE;
3110
3111 switch(n1->type)
3112 { case V_INTEGER:
3113 if ( mul64(n1->value.i, n2->value.i, &r->value.i) )
3114 { r->type = V_INTEGER;
3115 succeed;
3116 }
3117 /*FALLTHROUGH*/
3118 #ifdef O_GMP
3119 promoteToMPZNumber(n1);
3120 promoteToMPZNumber(n2);
3121
3122 case V_MPZ:
3123 mpz_init(r->value.mpz);
3124 r->type = V_MPZ;
3125 mpz_mul(r->value.mpz, n1->value.mpz, n2->value.mpz);
3126 succeed;
3127 case V_MPQ:
3128 r->type = V_MPQ;
3129 mpq_init(r->value.mpq);
3130 mpq_mul(r->value.mpq, n1->value.mpq, n2->value.mpq);
3131 return check_mpq(r);
3132 #else
3133 return PL_error("*", 2, NULL, ERR_EVALUATION, ATOM_int_overflow);
3134 #endif
3135 case V_FLOAT:
3136 r->value.f = n1->value.f * n2->value.f;
3137 r->type = V_FLOAT;
3138
3139 return check_float(r);
3140 }
3141
3142 assert(0);
3143 fail;
3144 }
3145
3146
3147 /* min/2 and max/2 have two special cases. If one of the arguments is
3148 * NaN we must select the other and for these functions -0.0 < 0.0,
3149 * while they compare == for normal float comparison.
3150 */
3151
3152 static int
ar_max(Number n1,Number n2,Number r)3153 ar_max(Number n1, Number n2, Number r)
3154 { int diff = cmpNumbers(n1, n2);
3155
3156 if ( diff == CMP_NOTEQ ) /* one or both nan */
3157 { if ( n1->type == V_FLOAT && isnan(n1->value.f) )
3158 cpNumber(r, n2);
3159 else
3160 cpNumber(r, n1);
3161 } else if ( diff == CMP_EQUAL &&
3162 n1->type == V_FLOAT && /* is n1 -0.0 */
3163 n1->value.f == 0.0 &&
3164 signbit(n1->value.f) )
3165 { cpNumber(r, n2);
3166 } else
3167 { if ( diff >= 0 )
3168 cpNumber(r, n1);
3169 else
3170 cpNumber(r, n2);
3171 }
3172
3173 return TRUE;
3174 }
3175
3176
3177 static int
ar_min(Number n1,Number n2,Number r)3178 ar_min(Number n1, Number n2, Number r)
3179 { int diff = cmpNumbers(n1, n2);
3180
3181 if ( diff == CMP_NOTEQ ) /* if one or both nan's */
3182 { if (n1->type == V_FLOAT && isnan(n1->value.f))
3183 cpNumber(r, n2);
3184 else
3185 cpNumber(r, n1);
3186 } else if ( diff == CMP_EQUAL &&
3187 n2->type == V_FLOAT && /* is n2 -0.0 */
3188 n2->value.f == 0.0 &&
3189 signbit(n2->value.f) )
3190 { cpNumber(r, n2);
3191 } else
3192 { if ( diff <= 0 )
3193 cpNumber(r, n1);
3194 else
3195 cpNumber(r, n2);
3196 }
3197
3198 return TRUE;
3199 }
3200
3201
3202 static int
ar_negation(Number n1,Number r)3203 ar_negation(Number n1, Number r)
3204 { if ( !toIntegerNumber(n1, 0) )
3205 return PL_error("\\", 1, NULL, ERR_AR_TYPE, ATOM_integer, n1);
3206
3207 switch(n1->type)
3208 { case V_INTEGER:
3209 r->value.i = ~n1->value.i;
3210 r->type = V_INTEGER;
3211 break;
3212 #ifdef O_GMP
3213 case V_MPZ:
3214 r->type = V_MPZ;
3215 mpz_init(r->value.mpz);
3216 mpz_com(r->value.mpz, n1->value.mpz);
3217 break;
3218 #endif
3219 default:
3220 assert(0);
3221 fail;
3222 }
3223 succeed;
3224 }
3225
3226
3227 static int
notLessThanZero(const char * f,int a,Number n)3228 notLessThanZero(const char *f, int a, Number n)
3229 { return PL_error(f, a, NULL, ERR_AR_DOMAIN, ATOM_not_less_than_zero, n);
3230 }
3231
3232
3233 static int
mustBePositive(const char * f,int a,Number n)3234 mustBePositive(const char *f, int a, Number n)
3235 { return PL_error(f, a, NULL, ERR_AR_DOMAIN, ATOM_not_less_than_one, n);
3236 }
3237
3238
3239 static int
ar_msb(Number n1,Number r)3240 ar_msb(Number n1, Number r)
3241 { if ( !toIntegerNumber(n1, 0) )
3242 return PL_error("msb", 1, NULL, ERR_AR_TYPE, ATOM_integer, n1);
3243
3244 switch(n1->type)
3245 { case V_INTEGER:
3246 if ( n1->value.i <= 0 )
3247 return mustBePositive("msb", 1, n1);
3248
3249 r->value.i = MSB64(n1->value.i);
3250 r->type = V_INTEGER;
3251 succeed;
3252 #ifdef O_GMP
3253 case V_MPZ:
3254 if ( mpz_sgn(n1->value.mpz) <= 0 )
3255 return mustBePositive("msb", 1, n1);
3256 if ( mpz_sgn(n1->value.mpz) == 0 )
3257 { r->value.i = 0;
3258 } else /* is binary print-size the best we can do?? */
3259 { r->value.i = mpz_sizeinbase(n1->value.mpz, 2)-1;
3260 }
3261 r->type = V_INTEGER;
3262 succeed;
3263 #endif
3264 default:
3265 assert(0);
3266 fail;
3267 }
3268 }
3269
3270
3271 static int
lsb64(int64_t i)3272 lsb64(int64_t i)
3273 { int j = 0;
3274
3275 if ( i == 0 )
3276 return 0;
3277
3278 if (!(i & LL(0xffffffff))) {i >>= 32; j += 32;}
3279 if (!(i & LL(0xffff))) {i >>= 16; j += 16;}
3280 if (!(i & LL(0xff))) {i >>= 8; j += 8;}
3281 if (!(i & LL(0xf))) {i >>= 4; j += 4;}
3282 if (!(i & LL(0x3))) {i >>= 2; j += 2;}
3283 if (!(i & LL(0x1))) j++;
3284
3285 return j;
3286 }
3287
3288
3289 static int
ar_lsb(Number n1,Number r)3290 ar_lsb(Number n1, Number r)
3291 { if ( !toIntegerNumber(n1, 0) )
3292 return PL_error("lsb", 1, NULL, ERR_AR_TYPE, ATOM_integer, n1);
3293
3294 switch(n1->type)
3295 { case V_INTEGER:
3296 if ( n1->value.i <= 0 )
3297 return mustBePositive("lsb", 1, n1);
3298
3299 r->value.i = lsb64(n1->value.i);
3300 r->type = V_INTEGER;
3301 succeed;
3302 #ifdef O_GMP
3303 case V_MPZ:
3304 if ( mpz_sgn(n1->value.mpz) <= 0 )
3305 return mustBePositive("lsb", 1, n1);
3306 r->value.i = mpz_scan1(n1->value.mpz, 0);
3307 r->type = V_INTEGER;
3308 succeed;
3309 #endif
3310 default:
3311 assert(0);
3312 fail;
3313 }
3314 }
3315
3316
3317 static int
my_popcount64(int64_t i)3318 my_popcount64(int64_t i) /* my_: avoid NetBSD name conflict */
3319 {
3320 #ifdef HAVE__BUILTIN_POPCOUNT
3321 return __builtin_popcountll(i);
3322 #else
3323 int c;
3324 size_t j;
3325 int64_t m = LL(1);
3326
3327 for(j=0,c=0; j<sizeof(i)*8; j++, m<<=1)
3328 { if ( i&m )
3329 c++;
3330 }
3331
3332 return c;
3333 #endif
3334 }
3335
3336
3337 static int
ar_popcount(Number n1,Number r)3338 ar_popcount(Number n1, Number r)
3339 { if ( !toIntegerNumber(n1, 0) )
3340 return PL_error("popcount", 1, NULL, ERR_AR_TYPE, ATOM_integer, n1);
3341
3342 switch(n1->type)
3343 { case V_INTEGER:
3344 if ( n1->value.i < 0 )
3345 return notLessThanZero("popcount", 1, n1);
3346
3347 r->value.i = my_popcount64(n1->value.i);
3348 r->type = V_INTEGER;
3349 succeed;
3350 #ifdef O_GMP
3351 case V_MPZ:
3352 if ( mpz_sgn(n1->value.mpz) < 0 )
3353 return notLessThanZero("popcount", 1, n1);
3354 r->value.i = mpz_popcount(n1->value.mpz);
3355 r->type = V_INTEGER;
3356 succeed;
3357 #endif
3358 default:
3359 assert(0);
3360 fail;
3361 }
3362 }
3363
3364 /* bit(I,K) is the K-th bit of I
3365 */
3366
3367 #ifndef HAVE_MP_BITCNT_T
3368 typedef unsigned long mp_bitcnt_t;
3369 #endif
3370
3371 #define MP_BITCNT_T_MIN 0
3372 #define MP_BITCNT_T_MAX (~(mp_bitcnt_t)0)
3373
3374 static int
ar_getbit(Number I,Number K,Number r)3375 ar_getbit(Number I, Number K, Number r)
3376 { mp_bitcnt_t bit;
3377
3378 if ( !toIntegerNumber(I, 0) )
3379 return PL_error("bit", 2, NULL, ERR_AR_TYPE, ATOM_integer, I);
3380 if ( !toIntegerNumber(K, 0) )
3381 return PL_error("bit", 2, NULL, ERR_AR_TYPE, ATOM_integer, K);
3382
3383 switch(K->type)
3384 { case V_INTEGER:
3385 if ( K->value.i < 0 )
3386 return notLessThanZero("bit", 2, K);
3387 if ( sizeof(mp_bitcnt_t) < 8 &&
3388 K->value.i > MP_BITCNT_T_MAX )
3389 { too_large:
3390 r->value.i = 0;
3391 r->type = V_INTEGER;
3392 return TRUE;
3393 }
3394 bit = K->value.i;
3395 break;
3396 #ifdef O_GMP
3397 case V_MPZ:
3398 if ( mpz_sgn(K->value.mpz) < 0 )
3399 return notLessThanZero("bit", 2, K);
3400 if ( mpz_cmp_ui(K->value.mpz, MP_BITCNT_T_MAX) > 0 )
3401 goto too_large;
3402 bit = mpz_get_ui(K->value.mpz);
3403 break;
3404 #endif
3405 default:
3406 bit = 0;
3407 assert(0);
3408 }
3409
3410 switch(I->type)
3411 { case V_INTEGER:
3412 if ( I->value.i < 0 )
3413 return notLessThanZero("bit", 2, I);
3414
3415 if ( bit >= 8*sizeof(I->value.i) )
3416 goto too_large;
3417 r->value.i = (I->value.i & ((int64_t)1<<bit)) != 0;
3418 r->type = V_INTEGER;
3419 return TRUE;
3420 #ifdef O_GMP
3421 case V_MPZ:
3422 if ( mpz_sgn(I->value.mpz) < 0 )
3423 return notLessThanZero("bit", 2, I);
3424
3425 r->value.i = mpz_tstbit(I->value.mpz, bit);
3426 r->type = V_INTEGER;
3427 return TRUE;
3428 #endif
3429 default:
3430 assert(0);
3431 fail;
3432 }
3433 }
3434
3435
3436 static int
ar_u_minus(Number n1,Number r)3437 ar_u_minus(Number n1, Number r)
3438 { r->type = n1->type;
3439
3440 switch(n1->type)
3441 { case V_INTEGER:
3442 if ( n1->value.i == PLMININT )
3443 {
3444 #ifdef O_GMP
3445 promoteToMPZNumber(n1);
3446 r->type = V_MPZ;
3447 #else
3448 if ( !promoteToFloatNumber(n1) )
3449 return FALSE;
3450 r->type = V_FLOAT;
3451 #endif
3452 /*FALLTHROUGH*/
3453 } else
3454 { r->value.i = -n1->value.i;
3455 break;
3456 }
3457 #ifdef O_GMP
3458 case V_MPZ:
3459 mpz_init(r->value.mpz);
3460 mpz_neg(r->value.mpz, n1->value.mpz);
3461 break;
3462 case V_MPQ:
3463 mpq_init(r->value.mpq);
3464 mpq_neg(r->value.mpq, n1->value.mpq);
3465 return check_mpq(r);
3466 #endif
3467 case V_FLOAT:
3468 r->value.f = isnan(n1->value.f) ? n1->value.f : -n1->value.f;
3469 r->type = V_FLOAT;
3470 break;
3471 }
3472
3473 succeed;
3474 }
3475
3476
3477 static int
ar_u_plus(Number n1,Number r)3478 ar_u_plus(Number n1, Number r)
3479 { cpNumber(r, n1);
3480
3481 succeed;
3482 }
3483
3484
3485 static int
ar_eval(Number n1,Number r)3486 ar_eval(Number n1, Number r)
3487 { cpNumber(r, n1);
3488
3489 succeed;
3490 }
3491
3492
3493 static int
ar_abs(Number n1,Number r)3494 ar_abs(Number n1, Number r)
3495 { switch(n1->type)
3496 { case V_INTEGER:
3497 if ( n1->value.i == PLMININT )
3498 {
3499 #ifdef O_GMP
3500 promoteToMPZNumber(n1);
3501 r->type = V_MPZ;
3502 #else
3503 if ( !promoteToFloatNumber(n1) )
3504 return FALSE;
3505 r->type = V_FLOAT;
3506 #endif
3507 /*FALLTHROUGH*/
3508 } else
3509 { r->value.i = llabs(n1->value.i);
3510 r->type = V_INTEGER;
3511 break;
3512 }
3513 #ifdef O_GMP
3514 case V_MPZ:
3515 r->type = V_MPZ;
3516 mpz_init(r->value.mpz);
3517 mpz_abs(r->value.mpz, n1->value.mpz);
3518 break;
3519 case V_MPQ:
3520 r->type = V_MPQ;
3521 mpq_init(r->value.mpq);
3522 mpq_abs(r->value.mpq, n1->value.mpq);
3523 break;
3524 #endif
3525 case V_FLOAT:
3526 { if ( signbit(n1->value.f) )
3527 r->value.f = -n1->value.f;
3528 else
3529 r->value.f = n1->value.f;
3530 r->type = V_FLOAT;
3531 break;
3532 }
3533 }
3534
3535 succeed;
3536 }
3537
3538
3539 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3540 Translate argument to rounded integer. If the double is outside the
3541 PLMININT/PLMAXINT range it is integer anyway, so we do not have to
3542 consider rounding for conversion to MPZ.
3543 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
3544
3545 static int
ar_integer(Number n1,Number r)3546 ar_integer(Number n1, Number r)
3547 { switch(n1->type)
3548 { case V_INTEGER:
3549 #ifdef O_GMP
3550 case V_MPZ:
3551 #endif
3552 cpNumber(r, n1);
3553 return TRUE;
3554 #ifdef O_GMP
3555 case V_MPQ:
3556 { mpq_t q;
3557 mpq_t half;
3558
3559 mpq_init(q);
3560 mpq_init(half);
3561 mpq_set_ui(half, 1, 2); /* 1/2 */
3562 if ( mpq_sgn(n1->value.mpq) > 0 )
3563 mpq_add(q, n1->value.mpq, half);
3564 else
3565 mpq_sub(q, n1->value.mpq, half);
3566
3567 r->type = V_MPZ;
3568 mpz_init(r->value.mpz);
3569 mpz_set_q(r->value.mpz, q);
3570 mpq_clear(q);
3571 mpq_clear(half);
3572 return TRUE;
3573 }
3574 #endif
3575 case V_FLOAT:
3576 { if ( isnan(n1->value.f) || isinf(n1->value.f) )
3577 { cpNumber(r, n1);
3578 return TRUE;
3579 }
3580
3581 if ( n1->value.f <= (float)PLMAXINT && n1->value.f >= (float)PLMININT )
3582 { if ( n1->value.f > 0 )
3583 { r->value.i = (int64_t)(n1->value.f + 0.5);
3584 if ( r->value.i < 0 ) /* Why can this happen? */
3585 r->value.i = PLMAXINT;
3586 } else
3587 { r->value.i = (int64_t)(n1->value.f - 0.5);
3588 if ( r->value.i > 0 )
3589 r->value.i = PLMININT;
3590 }
3591
3592 r->type = V_INTEGER;
3593 return TRUE;
3594 }
3595 #ifdef O_GMP
3596 r->type = V_MPZ;
3597 mpz_init_set_d(r->value.mpz, n1->value.f);
3598 return TRUE;
3599 #else
3600 #ifdef HAVE_RINT
3601 r->value.f = rint(n1->value.f);
3602 r->type = V_FLOAT;
3603 return TRUE;
3604 #else
3605 return PL_error("integer", 1, NULL, ERR_EVALUATION, ATOM_int_overflow);
3606 #endif
3607 #endif
3608 }
3609 }
3610
3611 assert(0);
3612 return FALSE;
3613 }
3614
3615
3616 static int
ar_float(Number n1,Number r)3617 ar_float(Number n1, Number r)
3618 { cpNumber(r, n1);
3619
3620 return promoteToFloatNumber(r);
3621 }
3622
3623
3624 static int /* ISO Prolog: R --> Z */
ar_floor(Number n1,Number r)3625 ar_floor(Number n1, Number r)
3626 { switch(n1->type)
3627 { case V_INTEGER:
3628 cpNumber(r, n1);
3629 succeed;
3630 #ifdef O_GMP
3631 case V_MPZ:
3632 cpNumber(r, n1);
3633 succeed;
3634 case V_MPQ:
3635 r->type = V_MPZ;
3636 mpz_init(r->value.mpz);
3637 mpz_set_q(r->value.mpz, n1->value.mpq);
3638 if ( mpq_sgn(n1->value.mpq) < 0 &&
3639 mpz_cmp_si(mpq_denref(n1->value.mpq), 1L) != 0 )
3640 mpz_sub_ui(r->value.mpz, r->value.mpz, 1L);
3641 succeed;
3642 #endif
3643 case V_FLOAT:
3644 { if ( isnan(n1->value.f) || isinf(n1->value.f) )
3645 { cpNumber(r, n1);
3646 return TRUE;
3647 }
3648 #ifdef HAVE_FLOOR
3649 r->type = V_FLOAT;
3650 r->value.f = floor(n1->value.f);
3651 if ( !toIntegerNumber(r, TOINT_CONVERT_FLOAT|TOINT_TRUNCATE) )
3652 { return PL_error("floor", 1, NULL, ERR_EVALUATION, ATOM_int_overflow);
3653 }
3654 #else /*HAVE_FLOOR*/
3655 if ( n1->value.f > (double)PLMININT && n1->value.f < (double)PLMAXINT )
3656 { r->value.i = (int64_t)n1->value.f;
3657 if ( n1->value.f < 0 && (double)r->value.i > n1->value.f )
3658 r->value.i--;
3659 r->type = V_INTEGER;
3660 } else
3661 {
3662 #ifdef O_GMP
3663 r->type = V_MPZ;
3664 mpz_init_set_d(r->value.mpz, n1->value.f);
3665 if ( n1->value.f < 0 &&
3666 mpX_round(mpz_get_d(r->value.mpz)) > n1->value.f )
3667 mpz_sub_ui(r->value.mpz, r->value.mpz, 1L);
3668 #else
3669 return PL_error("floor", 1, NULL, ERR_EVALUATION, ATOM_int_overflow);
3670 #endif
3671 }
3672 #endif /*HAVE_FLOOR*/
3673 }
3674 }
3675
3676 succeed;
3677 }
3678
3679
3680 static int /* ISO Prolog: R --> Z */
ar_ceil(Number n1,Number r)3681 ar_ceil(Number n1, Number r)
3682 { switch(n1->type)
3683 { case V_INTEGER:
3684 cpNumber(r, n1);
3685 succeed;
3686 #ifdef O_GMP
3687 case V_MPZ:
3688 cpNumber(r, n1);
3689 succeed;
3690 case V_MPQ:
3691 r->type = V_MPZ;
3692 mpz_init(r->value.mpz);
3693 mpz_set_q(r->value.mpz, n1->value.mpq);
3694 if ( mpq_sgn(n1->value.mpq) > 0 &&
3695 mpz_cmp_si(mpq_denref(n1->value.mpq), 1L) != 0 )
3696 mpz_add_ui(r->value.mpz, r->value.mpz, 1L);
3697 succeed;
3698 #endif
3699 case V_FLOAT:
3700 { if ( isnan(n1->value.f) || isinf(n1->value.f) )
3701 { cpNumber(r, n1);
3702 return TRUE;
3703 }
3704 #ifdef HAVE_CEIL
3705 r->type = V_FLOAT;
3706 r->value.f = ceil(n1->value.f);
3707 if ( !toIntegerNumber(r, TOINT_CONVERT_FLOAT|TOINT_TRUNCATE) )
3708 { return PL_error("ceil", 1, NULL, ERR_EVALUATION, ATOM_int_overflow);
3709 }
3710 #else /*HAVE_CEIL*/
3711 if ( n1->value.f > (double)PLMININT && n1->value.f < (double)PLMAXINT )
3712 { r->value.i = (int64_t)n1->value.f;
3713 if ( (double)r->value.i < n1->value.f )
3714 r->value.i++;
3715 r->type = V_INTEGER;
3716 } else
3717 {
3718 #ifdef O_GMP
3719 r->type = V_MPZ;
3720 mpz_init_set_d(r->value.mpz, n1->value.f);
3721 if ( mpX_round(mpz_get_d(r->value.mpz)) < n1->value.f )
3722 mpz_add_ui(r->value.mpz, r->value.mpz, 1L);
3723 #else
3724 return PL_error("ceil", 1, NULL, ERR_EVALUATION, ATOM_int_overflow);
3725 #endif
3726 }
3727 #endif /*HAVE_CEIL*/
3728 }
3729 }
3730
3731 succeed;
3732 }
3733
3734
3735 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3736 X is float_integer_part(X) + float_fractional_part(X)
3737
3738 If X < 0, both float_integer_part(X) and float_integer_part(X) are <= 0
3739 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
3740
3741 static int
ar_float_fractional_part(Number n1,Number r)3742 ar_float_fractional_part(Number n1, Number r)
3743 { switch(n1->type)
3744 { case V_INTEGER:
3745 #ifdef O_GMP
3746 case V_MPZ:
3747 #endif
3748 r->value.i = 0;
3749 r->type = V_INTEGER;
3750 succeed;
3751 #ifdef O_GMP
3752 case V_MPQ:
3753 r->type = V_MPQ;
3754 mpq_init(r->value.mpq);
3755 mpz_tdiv_q(mpq_numref(r->value.mpq),
3756 mpq_numref(n1->value.mpq),
3757 mpq_denref(n1->value.mpq));
3758 mpz_set_ui(mpq_denref(r->value.mpq), 1);
3759 mpq_sub(r->value.mpq, n1->value.mpq, r->value.mpq);
3760 succeed;
3761 #endif
3762 case V_FLOAT:
3763 { double ip;
3764
3765 r->value.f = modf(n1->value.f, &ip);
3766 r->type = V_FLOAT;
3767 return check_float(r);
3768 }
3769 default:
3770 assert(0);
3771 return FALSE;
3772 }
3773 }
3774
3775
3776 static int
ar_float_integer_part(Number n1,Number r)3777 ar_float_integer_part(Number n1, Number r)
3778 { switch(n1->type)
3779 { case V_INTEGER:
3780 #ifdef O_GMP
3781 case V_MPZ:
3782 #endif
3783 cpNumber(r, n1);
3784 succeed;
3785 #ifdef O_GMP
3786 case V_MPQ:
3787 r->type = V_MPZ;
3788 mpz_init(r->value.mpz);
3789 mpz_tdiv_q(r->value.mpz,
3790 mpq_numref(n1->value.mpq),
3791 mpq_denref(n1->value.mpq));
3792 succeed;
3793 #endif
3794 case V_FLOAT:
3795 { double ip;
3796
3797 (void)modf(n1->value.f, &ip);
3798 r->value.f = ip;
3799 r->type = V_FLOAT;
3800 return check_float(r);
3801 }
3802 }
3803
3804 assert(0);
3805 fail;
3806 }
3807
3808
3809 static int
ar_truncate(Number n1,Number r)3810 ar_truncate(Number n1, Number r)
3811 { switch(n1->type)
3812 {
3813 #ifdef O_GMP
3814 case V_MPQ:
3815 if ( mpq_sgn(n1->value.mpq) >= 0 )
3816 return ar_floor(n1, r);
3817 else
3818 return ar_ceil(n1, r);
3819 #endif
3820 case V_FLOAT:
3821 if ( isnan(n1->value.f) || isinf(n1->value.f) )
3822 { cpNumber(r, n1);
3823 return TRUE;
3824 }
3825 if ( n1->value.f >= 0.0 )
3826 return ar_floor(n1, r);
3827 else
3828 return ar_ceil(n1, r);
3829 default:
3830 cpNumber(r, n1);
3831 return TRUE;
3832 }
3833 }
3834
3835
3836 #ifdef O_GMP
3837 #ifdef HAVE_SYS_TYPES_H
3838 #include <sys/types.h>
3839 #endif
3840 #ifdef HAVE_SYS_STAT_H
3841 #include <sys/stat.h>
3842 #endif
3843 #include <fcntl.h>
3844
3845 #define RAND_SEED_LEN 128
3846 #define MIN_RAND_SEED_LEN 16
3847
3848 static int
seed_from_dev(const char * dev ARG_LD)3849 seed_from_dev(const char *dev ARG_LD)
3850 { int done = FALSE;
3851 #if defined(S_ISCHR) && !defined(__WINDOWS__)
3852 int fd;
3853
3854 if ( (fd=open(dev, O_RDONLY)) != -1 )
3855 { struct stat buf;
3856
3857 if ( fstat(fd, &buf) == 0 && S_ISCHR(buf.st_mode) )
3858 { char seedarray[RAND_SEED_LEN];
3859 mpz_t seed;
3860 size_t rd = 0;
3861 ssize_t n;
3862
3863 while ( rd < MIN_RAND_SEED_LEN )
3864 { if ( (n=read(fd, seedarray+rd, sizeof(seedarray)-rd)) > 0 )
3865 rd += n;
3866 else
3867 break;
3868 }
3869
3870 if ( rd >= MIN_RAND_SEED_LEN )
3871 { DEBUG(1, Sdprintf("Seed random using %ld bytes from %s\n",
3872 (long)rd, dev));
3873
3874 LD->gmp.persistent++;
3875 mpz_init(seed);
3876 mpz_import(seed, rd, 1, sizeof(char), 0, 0, seedarray);
3877 gmp_randseed(LD->arith.random.state, seed);
3878 mpz_clear(seed);
3879 LD->gmp.persistent--;
3880
3881 done = TRUE;
3882 }
3883 }
3884
3885 close(fd);
3886 }
3887 #endif /*S_ISCHR*/
3888
3889 return done;
3890 }
3891
3892
3893
3894 static int
seed_from_crypt_context(ARG1_LD)3895 seed_from_crypt_context(ARG1_LD)
3896 {
3897 #ifdef __WINDOWS__
3898 HCRYPTPROV hCryptProv;
3899 char *user_name = "seed_random";
3900 BYTE seedarray[RAND_SEED_LEN];
3901 mpz_t seed;
3902
3903
3904 if ( CryptAcquireContext(&hCryptProv, user_name, NULL, PROV_RSA_FULL, 0) )
3905 { CryptGenRandom(hCryptProv, sizeof(seedarray), seedarray);
3906 } else if ( (GetLastError() == NTE_BAD_KEYSET) &&
3907 CryptAcquireContext(&hCryptProv, user_name, NULL,
3908 PROV_RSA_FULL, CRYPT_NEWKEYSET) )
3909 { CryptGenRandom(hCryptProv, sizeof(seedarray), seedarray);
3910 } else
3911 { return FALSE;
3912 }
3913
3914 LD->gmp.persistent++;
3915 mpz_init(seed);
3916 mpz_import(seed, RAND_SEED_LEN, 1, sizeof(BYTE), 0, 0, seedarray);
3917 gmp_randseed(LD->arith.random.state, seed);
3918 mpz_clear(seed);
3919 LD->gmp.persistent--;
3920
3921 return TRUE;
3922 #else
3923 #ifdef O_PLMT
3924 (void)__PL_ld;
3925 #endif
3926 return FALSE;
3927 #endif
3928 }
3929
3930
3931 static void
seed_random(ARG1_LD)3932 seed_random(ARG1_LD)
3933 { if ( !seed_from_dev("/dev/urandom" PASS_LD) &&
3934 !seed_from_dev("/dev/random" PASS_LD) &&
3935 !seed_from_crypt_context(PASS_LD1) )
3936 { union
3937 { double t;
3938 unsigned long l[sizeof(double)/sizeof(long)];
3939 } u;
3940 unsigned long key = 0;
3941 int i;
3942
3943 u.t = WallTime();
3944
3945 for(i=0; i<sizeof(double)/sizeof(long); i++)
3946 key ^= u.l[i];
3947
3948 LD->gmp.persistent++;
3949 gmp_randseed_ui(LD->arith.random.state, key);
3950 LD->gmp.persistent--;
3951 }
3952 }
3953
3954 #else /* O_GMP */
3955
3956 static void
seed_random(ARG1_LD)3957 seed_random(ARG1_LD)
3958 { setRandom(NULL);
3959 }
3960
3961 #endif /*O_GMP*/
3962
3963 static void
init_random(ARG1_LD)3964 init_random(ARG1_LD)
3965 {
3966 #ifdef O_GMP
3967 if ( !LD->arith.random.initialised )
3968 { LD->gmp.persistent++;
3969 #ifdef HAVE_GMP_RANDINIT_MT
3970 #define O_RANDOM_STATE 1
3971 gmp_randinit_mt(LD->arith.random.state);
3972 #else
3973 gmp_randinit_default(LD->arith.random.state);
3974 #endif
3975 LD->arith.random.initialised = TRUE;
3976 seed_random(PASS_LD1);
3977 LD->gmp.persistent--;
3978 }
3979 #endif
3980 }
3981
3982
3983 static
3984 PRED_IMPL("set_random", 1, set_random, 0)
3985 { PRED_LD
3986 atom_t name;
3987 size_t arity;
3988
3989 init_random(PASS_LD1);
3990
3991 if ( PL_get_name_arity(A1, &name, &arity) && arity == 1 )
3992 { term_t arg = PL_new_term_ref();
3993
3994 _PL_get_arg(1, A1, arg);
3995 if ( name == ATOM_seed )
3996 { atom_t a;
3997
3998 if ( PL_get_atom(arg, &a) && a == ATOM_random )
3999 { seed_random(PASS_LD1);
4000 return TRUE;
4001 } else
4002 { number n;
4003
4004 if ( !PL_get_number(arg, &n) )
4005 return PL_error(NULL, 0, "integer or 'random'",
4006 ERR_TYPE, ATOM_seed, arg);
4007 switch(n.type)
4008 {
4009 #ifdef O_GMP
4010 case V_INTEGER:
4011 gmp_randseed_ui(LD->arith.random.state,
4012 (unsigned long)n.value.i);
4013 return TRUE;
4014 case V_MPZ:
4015 gmp_randseed(LD->arith.random.state, n.value.mpz);
4016 return TRUE;
4017 #else
4018 case V_INTEGER:
4019 { unsigned int seed = (unsigned int)n.value.i;
4020 setRandom(&seed);
4021 return TRUE;
4022 }
4023 #endif
4024 default:
4025 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_seed, arg);
4026 }
4027 }
4028 #ifdef O_RANDOM_STATE
4029 } else if ( name == ATOM_state )
4030 { number n;
4031
4032 if ( !PL_get_number(arg, &n) ||
4033 n.type != V_MPZ )
4034 return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_state, arg);
4035
4036 mpz_set(LD->arith.random.state[0]._mp_seed, n.value.mpz);
4037 clearNumber(&n);
4038
4039 return TRUE;
4040 #endif /*O_GMP*/
4041 } else
4042 { return PL_error(NULL, 0, NULL, ERR_DOMAIN, ATOM_random_option, A1);
4043 }
4044 } else
4045 { return PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_random_option, A1);
4046 }
4047 }
4048
4049
4050 #ifdef O_RANDOM_STATE
4051 static
4052 PRED_IMPL("random_property", 1, random_property, 0)
4053 { PRED_LD
4054 atom_t name;
4055 size_t arity;
4056
4057 init_random(PASS_LD1);
4058
4059 if ( PL_get_name_arity(A1, &name, &arity) && arity == 1 )
4060 { term_t arg = PL_new_term_ref();
4061
4062 _PL_get_arg(1, A1, arg);
4063 if ( name == ATOM_state )
4064 { int rc;
4065 number seed;
4066
4067 seed.type = V_MPZ;
4068 mpz_init(seed.value.mpz);
4069 LD->arith.random.state[0]._mp_seed[0]._mp_size =
4070 LD->arith.random.state[0]._mp_seed[0]._mp_alloc; mpz_set(seed.value.mpz, LD->arith.random.state[0]._mp_seed);
4071 rc = PL_unify_number(arg, &seed);
4072 clearNumber(&seed);
4073
4074 return rc;
4075 }
4076 }
4077
4078 return FALSE;
4079 }
4080 #endif
4081
4082
4083 static int
ar_random(Number n1,Number r)4084 ar_random(Number n1, Number r)
4085 { GET_LD
4086
4087 if ( !toIntegerNumber(n1, TOINT_CONVERT_FLOAT) )
4088 return PL_error("random", 1, NULL, ERR_AR_TYPE, ATOM_integer, n1);
4089 if ( ar_sign_i(n1) <= 0 )
4090 return mustBePositive("random", 1, n1);
4091
4092 init_random(PASS_LD1);
4093
4094 switch(n1->type)
4095 {
4096 #ifdef O_GMP
4097 case V_INTEGER:
4098 promoteToMPZNumber(n1);
4099 assert(n1->type == V_MPZ);
4100 /*FALLTHROUGH*/
4101 case V_MPZ:
4102 { r->type = V_MPZ;
4103 mpz_init(r->value.mpz);
4104 mpz_urandomm(r->value.mpz, LD->arith.random.state, n1->value.mpz);
4105
4106 succeed;
4107 }
4108 #else
4109 case V_INTEGER:
4110 if ( n1->value.i < 1 )
4111 return mustBePositive("random", 1, n1);
4112 r->value.i = _PL_Random() % (uint64_t)n1->value.i;
4113 r->type = V_INTEGER;
4114
4115 succeed;
4116 #endif
4117 default:
4118 assert(0);
4119 fail;
4120 }
4121 }
4122
4123 #ifndef UINT64_MAX
4124 #define UINT64_MAX (~(uint64_t)0)
4125 #endif
4126
4127 static int
ar_random_float(Number r)4128 ar_random_float(Number r)
4129 { GET_LD
4130
4131 init_random(PASS_LD1);
4132
4133 do
4134 {
4135 #ifdef O_GMP
4136 mpf_t rop;
4137 mpf_init2(rop, sizeof(double)*8);
4138 mpf_urandomb(rop, LD->arith.random.state, sizeof(double)*8);
4139 r->value.f = mpX_round(mpf_get_d(rop));
4140 mpf_clear(rop);
4141 #else
4142 r->value.f = _PL_Random()/(float)UINT64_MAX;
4143 #endif
4144 } while (r->value.f == 0.0);
4145
4146 r->type = V_FLOAT;
4147 succeed;
4148 }
4149
4150
4151 static int
ar_pi(Number r)4152 ar_pi(Number r)
4153 { r->value.f = (fegetround() == FE_UPWARD) ? nexttoward(M_PI,INFINITY) : M_PI;
4154
4155 r->type = V_FLOAT;
4156 succeed;
4157 }
4158
4159
4160 static int
ar_e(Number r)4161 ar_e(Number r)
4162 { r->value.f = (fegetround() == FE_UPWARD) ? nexttoward(M_E,INFINITY) : M_E;
4163
4164 r->type = V_FLOAT;
4165 succeed;
4166 }
4167
4168
4169 static int
ar_epsilon(Number r)4170 ar_epsilon(Number r)
4171 { r->value.f = DBL_EPSILON;
4172
4173 r->type = V_FLOAT;
4174 succeed;
4175 }
4176
4177
4178 static int
ar_inf(Number r)4179 ar_inf(Number r)
4180 { static number n = {0};
4181
4182 if ( n.type != V_FLOAT )
4183 { n.value.f = const_inf;
4184 n.type = V_FLOAT;
4185 }
4186
4187 *r = n;
4188
4189 succeed;
4190 }
4191
4192
4193 static int
ar_nan(Number r)4194 ar_nan(Number r)
4195 { static number n = {0};
4196
4197 if ( n.type != V_FLOAT )
4198 { n.value.f = const_nan;
4199 n.type = V_FLOAT;
4200 }
4201
4202 *r = n;
4203
4204 succeed;
4205 }
4206
4207
4208 static int
ar_cputime(Number r)4209 ar_cputime(Number r)
4210 { r->value.f = CpuTime(CPU_USER);
4211
4212 r->type = V_FLOAT;
4213 succeed;
4214 }
4215
4216
4217 /********************************
4218 * PROLOG CONNECTION *
4219 *********************************/
4220
4221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4222 (*) valueExpression() cannot use GC, but can return a number whose value
4223 is a GPM number pointing at the global stack. If this is the case,
4224 PL_unify_number() may not invoke GC, so we must check that we have room
4225 for the required attribute wakeup and trailing before we start.
4226
4227 is/2 is the only victim of this issue, as the other arithmetic
4228 predicates (>/2, etc.) only use their arguments as inputs.
4229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
4230
4231 static
4232 PRED_IMPL("is", 2, is, PL_FA_ISO) /* -Value is +Expr */
4233 { PRED_LD
4234 AR_CTX
4235 number arg;
4236 int rc;
4237
4238 if ( !hasGlobalSpace(0) ) /* see (*) */
4239 { if ( (rc=ensureGlobalSpace(0, ALLOW_GC)) != TRUE )
4240 return raiseStackOverflow(rc);
4241 }
4242
4243 AR_BEGIN();
4244 if ( (rc=valueExpression(A2, &arg PASS_LD)) )
4245 { rc = PL_unify_number(A1, &arg);
4246 clearNumber(&arg);
4247 AR_END();
4248 } else
4249 { AR_CLEANUP();
4250 }
4251
4252 return rc;
4253 }
4254
4255
4256 /** current_arithmetic_function(?Term) is nondet.
4257
4258 True if Term is evaluable.
4259 */
4260
4261 static
4262 PRED_IMPL("current_arithmetic_function", 1, current_arithmetic_function,
4263 PL_FA_NONDETERMINISTIC)
4264 { PRED_LD
4265 unsigned int i;
4266 term_t head = A1;
4267
4268 switch( CTX_CNTRL )
4269 { case FRG_FIRST_CALL:
4270 { functor_t fd;
4271
4272 if ( PL_is_variable(head) )
4273 { i = 0;
4274 break;
4275 } else if ( PL_get_functor(head, &fd) )
4276 { return isCurrentArithFunction(fd) ? TRUE : FALSE;
4277 } else
4278 return PL_error(NULL, 0, NULL,
4279 ERR_TYPE, ATOM_callable, head);
4280 }
4281 case FRG_REDO:
4282 i = (int)CTX_INT;
4283 break;
4284 case FRG_CUTTED:
4285 default:
4286 succeed;
4287 }
4288
4289 for(; i<GD->arith.functions_allocated; i++)
4290 { if ( GD->arith.functions[i] )
4291 { functor_t f = functorArithFunction(i);
4292
4293 if ( PL_unify_functor(head, f) )
4294 ForeignRedoInt(i+1);
4295 }
4296 }
4297
4298 fail;
4299 }
4300
4301
4302 typedef struct
4303 { functor_t functor;
4304 ArithF function;
4305 } ar_funcdef;
4306
4307 #define F_ISO 0x1
4308 #define ADD(functor, func, flags) { functor, func }
4309
4310 static const ar_funcdef ar_funcdefs[] = {
4311 ADD(FUNCTOR_plus2, pl_ar_add, F_ISO),
4312 ADD(FUNCTOR_minus2, ar_minus, F_ISO),
4313 ADD(FUNCTOR_star2, ar_mul, F_ISO),
4314 ADD(FUNCTOR_divide2, ar_divide, F_ISO),
4315 #ifdef O_GMP
4316 ADD(FUNCTOR_rational1, ar_rational, 0),
4317 ADD(FUNCTOR_numerator1, ar_numerator, 0),
4318 ADD(FUNCTOR_denominator1, ar_denominator, 0),
4319 ADD(FUNCTOR_rationalize1, ar_rationalize, 0),
4320 ADD(FUNCTOR_rdiv2, ar_rdiv, 0),
4321 #endif
4322 ADD(FUNCTOR_minus1, ar_u_minus, F_ISO),
4323 ADD(FUNCTOR_plus1, ar_u_plus, F_ISO),
4324 ADD(FUNCTOR_abs1, ar_abs, F_ISO),
4325 ADD(FUNCTOR_max2, ar_max, F_ISO),
4326 ADD(FUNCTOR_min2, ar_min, F_ISO),
4327
4328 ADD(FUNCTOR_mod2, ar_mod, F_ISO),
4329 ADD(FUNCTOR_rem2, ar_rem, F_ISO),
4330 ADD(FUNCTOR_div2, ar_div, F_ISO), /* div/2 */
4331 ADD(FUNCTOR_gdiv2, ar_tdiv, 0), /* (//)/2 */
4332 ADD(FUNCTOR_gcd2, ar_gcd, 0),
4333 ADD(FUNCTOR_lcm2, ar_lcm, 0),
4334 ADD(FUNCTOR_sign1, ar_sign, F_ISO),
4335
4336 ADD(FUNCTOR_and2, ar_conjunct, F_ISO),
4337 ADD(FUNCTOR_bitor2, ar_disjunct, F_ISO),
4338 ADD(FUNCTOR_rshift2, ar_shift_right, F_ISO),
4339 ADD(FUNCTOR_lshift2, ar_shift_left, F_ISO),
4340 ADD(FUNCTOR_xor2, ar_xor, F_ISO),
4341 ADD(FUNCTOR_backslash1, ar_negation, F_ISO),
4342
4343 ADD(FUNCTOR_random1, ar_random, 0),
4344 ADD(FUNCTOR_random_float0, ar_random_float, 0),
4345
4346 ADD(FUNCTOR_integer1, ar_integer, F_ISO),
4347 ADD(FUNCTOR_round1, ar_integer, F_ISO),
4348 ADD(FUNCTOR_truncate1, ar_truncate, F_ISO),
4349 ADD(FUNCTOR_float1, ar_float, F_ISO),
4350 ADD(FUNCTOR_floor1, ar_floor, F_ISO),
4351 ADD(FUNCTOR_ceil1, ar_ceil, F_ISO),
4352 ADD(FUNCTOR_ceiling1, ar_ceil, F_ISO),
4353 ADD(FUNCTOR_float_fractional_part1, ar_float_fractional_part, F_ISO),
4354 ADD(FUNCTOR_float_integer_part1, ar_float_integer_part, F_ISO),
4355 ADD(FUNCTOR_copysign2, ar_copysign, 0),
4356 ADD(FUNCTOR_nexttoward2, ar_nexttoward, 0),
4357 ADD(FUNCTOR_roundtoward2, ar_roundtoward, 0),
4358
4359 ADD(FUNCTOR_sqrt1, ar_sqrt, F_ISO),
4360 ADD(FUNCTOR_sin1, ar_sin, F_ISO),
4361 ADD(FUNCTOR_cos1, ar_cos, F_ISO),
4362 ADD(FUNCTOR_tan1, ar_tan, F_ISO),
4363 ADD(FUNCTOR_asin1, ar_asin, F_ISO),
4364 ADD(FUNCTOR_acos1, ar_acos, F_ISO),
4365 ADD(FUNCTOR_atan1, ar_atan, F_ISO),
4366 ADD(FUNCTOR_atan2, ar_atan2, 0),
4367 ADD(FUNCTOR_atan22, ar_atan2, F_ISO),
4368 ADD(FUNCTOR_sinh1, ar_sinh, 0),
4369 ADD(FUNCTOR_cosh1, ar_cosh, 0),
4370 ADD(FUNCTOR_tanh1, ar_tanh, 0),
4371 ADD(FUNCTOR_asinh1, ar_asinh, 0),
4372 ADD(FUNCTOR_acosh1, ar_acosh, 0),
4373 ADD(FUNCTOR_atanh1, ar_atanh, 0),
4374 ADD(FUNCTOR_lgamma1, ar_lgamma, 0),
4375 ADD(FUNCTOR_log1, ar_log, F_ISO),
4376 ADD(FUNCTOR_erf1, ar_erf, 0),
4377 ADD(FUNCTOR_erfc1, ar_erfc, 0),
4378 ADD(FUNCTOR_exp1, ar_exp, F_ISO),
4379 ADD(FUNCTOR_log101, ar_log10, 0),
4380 ADD(FUNCTOR_hat2, ar_pow, F_ISO),
4381 ADD(FUNCTOR_doublestar2, ar_pow, F_ISO),
4382 ADD(FUNCTOR_pi0, ar_pi, F_ISO),
4383 ADD(FUNCTOR_e0, ar_e, 0),
4384 ADD(FUNCTOR_epsilon0, ar_epsilon, 0),
4385 ADD(FUNCTOR_inf0, ar_inf, 0),
4386 ADD(FUNCTOR_nan0, ar_nan, 0),
4387
4388 ADD(FUNCTOR_cputime0, ar_cputime, 0),
4389 ADD(FUNCTOR_msb1, ar_msb, 0),
4390 ADD(FUNCTOR_lsb1, ar_lsb, 0),
4391 ADD(FUNCTOR_popcount1, ar_popcount, 0),
4392 ADD(FUNCTOR_getbit2, ar_getbit, 0),
4393 ADD(FUNCTOR_powm3, ar_powm, 0),
4394
4395 ADD(FUNCTOR_eval1, ar_eval, 0)
4396 };
4397
4398 #undef ADD
4399
4400 static size_t
registerFunction(functor_t f,ArithF func)4401 registerFunction(functor_t f, ArithF func)
4402 { size_t index = indexFunctor(f);
4403
4404 DEBUG(1, Sdprintf("Register functor %ld\n", (long)index));
4405
4406 while ( index >= GD->arith.functions_allocated )
4407 { if ( GD->arith.functions_allocated == 0 )
4408 { size_t size = 512;
4409
4410 GD->arith.functions = allocHeapOrHalt(size*sizeof(ArithF));
4411 memset(GD->arith.functions, 0, size*sizeof(ArithF));
4412 GD->arith.functions_allocated = size;
4413 } else
4414 { size_t size = GD->arith.functions_allocated*2;
4415 ArithF *new = allocHeapOrHalt(size*sizeof(ArithF));
4416 size_t half = GD->arith.functions_allocated*sizeof(ArithF);
4417 ArithF *old = GD->arith.functions;
4418
4419 DEBUG(0, Sdprintf("Re-sized function-table to %ld\n", (long)size));
4420
4421 memcpy(new, old, half);
4422 memset(addPointer(new,half), 0, half);
4423 GD->arith.functions = new;
4424 GD->arith.functions_allocated = size;
4425 freeHeap(old, half);
4426 }
4427 }
4428
4429 GD->arith.functions[index] = func;
4430
4431 return index;
4432 }
4433
4434
4435 static void
registerBuiltinFunctions()4436 registerBuiltinFunctions()
4437 { int n, size = sizeof(ar_funcdefs)/sizeof(ar_funcdef);
4438 const ar_funcdef *d;
4439
4440 for(d = ar_funcdefs, n=0; n<size; n++, d++)
4441 registerFunction(d->functor, d->function);
4442 }
4443
4444
4445 static atom_t float_rounding_names[5] = {0};
4446
4447 static double
nan15(void)4448 nan15(void)
4449 { number n;
4450 unsigned char *end;
4451
4452 if ( str_number((const unsigned char *)"1.5NaN", &end, &n, 0) == NUM_OK )
4453 { assert(isnan(n.value.f));
4454 return n.value.f;
4455 } else
4456 { assert(0);
4457 return NAN;
4458 }
4459 }
4460
4461
4462 double
PL_nan(void)4463 PL_nan(void)
4464 { return const_nan;
4465 }
4466
4467
4468 void
initArith(void)4469 initArith(void)
4470 { GET_LD
4471 registerBuiltinFunctions();
4472
4473 #ifdef O_INHIBIT_FP_SIGNALS
4474 fpsetmask(fpgetmask() & ~(FP_X_DZ|FP_X_INV|FP_X_OFL));
4475 #endif
4476
4477 const_nan = nan15();
4478 const_inf = HUGE_VAL;
4479 const_neg_inf = -HUGE_VAL;
4480
4481 #ifdef O_GMP
4482 LD->arith.rat.max_rational_size = (size_t)-1;
4483 LD->arith.rat.max_rational_size_action = ATOM_error;
4484
4485 setPrologFlag("max_rational_size", FT_INTEGER, -1);
4486 setPrologFlag("max_rational_size_action", FT_ATOM, "error");
4487 #endif
4488
4489 LD->arith.f.flags = FLT_ROUND_NEAREST|FLT_UNDERFLOW;
4490 setPrologFlag("float_overflow", FT_ATOM, "error");
4491 setPrologFlag("float_zero_div", FT_ATOM, "error");
4492 setPrologFlag("float_undefined", FT_ATOM, "error");
4493 setPrologFlag("float_underflow", FT_ATOM, "ignore");
4494 setPrologFlag("float_rounding", FT_ATOM, "to_nearest");
4495 float_rounding_names[FLT_ROUND_NEAREST] = ATOM_to_nearest;
4496 float_rounding_names[FLT_ROUND_TO_POS] = ATOM_to_positive;
4497 float_rounding_names[FLT_ROUND_TO_NEG] = ATOM_to_negative;
4498 float_rounding_names[FLT_ROUND_TO_ZERO] = ATOM_to_zero;
4499 setPrologFlag("float_min", FT_FLOAT|FF_READONLY, DBL_MIN);
4500 setPrologFlag("float_max", FT_FLOAT|FF_READONLY, DBL_MAX);
4501 setPrologFlag("float_max_integer", FT_FLOAT|FF_READONLY, 9007199254740992.0);
4502 }
4503
4504
4505 void
cleanupArith(void)4506 cleanupArith(void)
4507 {
4508 #ifdef O_INHIBIT_FP_SIGNALS
4509 fpresetsticky(FP_X_DZ|FP_X_INV|FP_X_OFL);
4510 fpsetmask(FP_X_DZ|FP_X_INV|FP_X_OFL);
4511 #endif
4512
4513 if ( GD->arith.functions )
4514 { freeHeap(GD->arith.functions, GD->arith.functions_allocated*sizeof(ArithF));
4515 GD->arith.functions = 0;
4516 GD->arith.functions_allocated = 0;
4517 }
4518 }
4519
4520 int
is_arith_flag(atom_t k)4521 is_arith_flag(atom_t k)
4522 { return (
4523 #ifdef O_GMP
4524 k == ATOM_max_rational_size ||
4525 k == ATOM_max_rational_size_action ||
4526 #endif
4527 k == ATOM_float_overflow ||
4528 k == ATOM_float_zero_div ||
4529 k == ATOM_float_undefined ||
4530 k == ATOM_float_underflow ||
4531 k == ATOM_float_rounding );
4532 }
4533
4534 int
get_arith_flag(term_t val,atom_t k ARG_LD)4535 get_arith_flag(term_t val, atom_t k ARG_LD)
4536 { atom_t a;
4537 #ifdef O_GMP
4538 size_t sz;
4539
4540 if ( k == ATOM_max_rational_size &&
4541 (sz=LD->arith.rat.max_rational_size) != (size_t)-1 )
4542 return PL_unify_uint64(val, sz);
4543 if ( k == ATOM_max_rational_size_action )
4544 return PL_unify_atom(val, LD->arith.rat.max_rational_size_action);
4545 #endif
4546 if ( k == ATOM_float_overflow )
4547 a = LD->arith.f.flags & FLT_OVERFLOW ? ATOM_infinity : ATOM_error;
4548 else if ( k == ATOM_float_zero_div )
4549 a = LD->arith.f.flags & FLT_ZERO_DIV ? ATOM_infinity : ATOM_error;
4550 else if ( k == ATOM_float_undefined )
4551 a = LD->arith.f.flags & FLT_UNDEFINED ? ATOM_nan : ATOM_error;
4552 else if ( k == ATOM_float_underflow )
4553 a = LD->arith.f.flags & FLT_UNDERFLOW ? ATOM_ignore : ATOM_error;
4554 else if ( k == ATOM_float_rounding )
4555 a = float_rounding_names[LD->arith.f.flags&FLT_ROUND_MASK];
4556 else
4557 return FALSE;
4558
4559 return PL_unify_atom(val, a);
4560 }
4561
4562 #ifdef O_GMP
4563 static int
set_restraint(term_t t,size_t * valp)4564 set_restraint(term_t t, size_t *valp)
4565 { GET_LD
4566 atom_t inf;
4567
4568 if ( PL_get_atom(t, &inf) && inf == ATOM_infinite )
4569 { *valp = (size_t)-1;
4570 return TRUE;
4571 }
4572 return PL_get_size_ex(t, valp);
4573 }
4574
4575 static int
set_restraint_action(term_t t,atom_t key,atom_t * valp ARG_LD)4576 set_restraint_action(term_t t, atom_t key, atom_t *valp ARG_LD)
4577 { atom_t act;
4578
4579 if ( PL_get_atom_ex(t, &act) )
4580 { if ( act == ATOM_error )
4581 { ok:
4582 *valp = act;
4583 return TRUE;
4584 }
4585
4586 if ( key == ATOM_max_rational_size_action &&
4587 ( act == ATOM_float ) )
4588 goto ok;
4589
4590 return PL_domain_error("max_rational_size_action", t);
4591 }
4592
4593 return FALSE;
4594 }
4595 #endif
4596
4597 static int rounding_mode[5] =
4598 {0, FE_TONEAREST, FE_UPWARD, FE_DOWNWARD, FE_TOWARDZERO};
4599
4600 int
atom_to_rounding(atom_t a,int * m)4601 atom_to_rounding(atom_t a, int *m)
4602 { int i;
4603
4604 for(i=1; i<=FLT_ROUND_TO_ZERO; i++)
4605 { if ( float_rounding_names[i] == a )
4606 { *m = i;
4607 return TRUE;
4608 }
4609 }
4610
4611 return FALSE;
4612 }
4613
4614
4615 atom_t
float_rounding_name(int i)4616 float_rounding_name(int i)
4617 { return float_rounding_names[i];
4618 }
4619
4620
4621 void
set_rounding(int mode)4622 set_rounding(int mode)
4623 { fesetround(rounding_mode[mode]);
4624 }
4625
4626
4627 int
set_arith_flag(term_t val,atom_t key ARG_LD)4628 set_arith_flag(term_t val, atom_t key ARG_LD)
4629 { atom_t a;
4630
4631 #ifdef O_GMP
4632 if ( key == ATOM_max_rational_size )
4633 return set_restraint(val, &LD->arith.rat.max_rational_size);
4634 if ( key == ATOM_max_rational_size_action )
4635 return set_restraint_action(
4636 val, key,
4637 &LD->arith.rat.max_rational_size_action PASS_LD);
4638 #endif
4639
4640 if ( PL_get_atom_ex(val, &a) )
4641 { if ( key == ATOM_float_overflow )
4642 { if ( a == ATOM_error ) clear(&(LD->arith.f), FLT_OVERFLOW);
4643 else if ( a == ATOM_infinity ) set(&(LD->arith.f), FLT_OVERFLOW);
4644 else goto dom;
4645 } else if ( key == ATOM_float_zero_div )
4646 { if ( a == ATOM_error ) clear(&(LD->arith.f), FLT_ZERO_DIV);
4647 else if ( a == ATOM_infinity ) set(&(LD->arith.f), FLT_ZERO_DIV);
4648 else goto dom;
4649 } else if ( key == ATOM_float_undefined )
4650 { if ( a == ATOM_error ) clear(&(LD->arith.f), FLT_UNDEFINED);
4651 else if ( a == ATOM_nan ) set(&(LD->arith.f), FLT_UNDEFINED);
4652 else goto dom;
4653 } else if ( key == ATOM_float_underflow )
4654 { if ( a == ATOM_error ) clear(&(LD->arith.f), FLT_UNDERFLOW);
4655 else if ( a == ATOM_ignore ) set(&(LD->arith.f), FLT_UNDERFLOW);
4656 else goto dom;
4657 } else if ( key == ATOM_float_rounding )
4658 { int i;
4659
4660 if ( atom_to_rounding(a, &i) )
4661 { clear(&(LD->arith.f), FLT_ROUND_MASK);
4662 LD->arith.f.flags |= i;
4663 set_rounding(i);
4664 return TRUE;
4665 }
4666 goto dom;
4667 } else
4668 return FALSE;
4669
4670 return TRUE;
4671
4672 dom:
4673 return PL_domain_error("flag_value", val);
4674 }
4675
4676 return FALSE;
4677 }
4678
4679
4680 static
4681 PRED_IMPL("float_class", 2, float_class, 0)
4682 { PRED_LD
4683 Word p = valTermRef(A1);
4684
4685 deRef(p);
4686 if ( isFloat(*p) )
4687 { double f = valFloat(*p);
4688 atom_t a;
4689
4690 switch(fpclassify(f))
4691 { case FP_NAN: a = ATOM_nan; break;
4692 case FP_INFINITE: a = ATOM_infinite; break;
4693 case FP_ZERO: a = ATOM_zero; break;
4694 case FP_SUBNORMAL: a = ATOM_subnormal; break;
4695 case FP_NORMAL: a = ATOM_normal; break;
4696 default: assert(0); a = ATOM_false;
4697 }
4698
4699 return PL_unify_atom(A2, a);
4700 }
4701
4702 return PL_type_error("float", A1);
4703 }
4704
4705 #if O_COMPILE_ARITH
4706
4707 /********************************
4708 * VIRTUAL MACHINE SUPPORT *
4709 *********************************/
4710
4711 int
indexArithFunction(functor_t f)4712 indexArithFunction(functor_t f)
4713 { size_t index = indexFunctor(f);
4714
4715 if ( index < GD->arith.functions_allocated )
4716 { if ( GD->arith.functions[index] )
4717 return (int)index;
4718 }
4719
4720 return -1;
4721 }
4722
4723
4724 functor_t
functorArithFunction(unsigned int i)4725 functorArithFunction(unsigned int i)
4726 { FunctorDef fd = fetchFunctorArray(i);
4727
4728 return fd->functor;
4729 }
4730
4731
4732 static ArithF
FunctionFromIndex(int index)4733 FunctionFromIndex(int index)
4734 { return GD->arith.functions[index];
4735 }
4736
4737
4738 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
4739 ar_func_n(code, argc) is executed by the A_FUNC* virtual machine
4740 instructions. It invalidates all numbers it pops from the stack using
4741 clearNumber()
4742 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
4743
4744 bool
ar_func_n(int findex,int argc ARG_LD)4745 ar_func_n(int findex, int argc ARG_LD)
4746 { number result;
4747 int rval;
4748 ArithF f = FunctionFromIndex(findex);
4749 Number argv = argvArithStack(argc PASS_LD);
4750
4751 DEBUG(0, if ( !f )
4752 fatalError("No function at index %d", findex));
4753
4754 switch(argc)
4755 { case 0:
4756 rval = (*f)(&result);
4757 break;
4758 case 1:
4759 rval = (*f)(argv, &result);
4760 break;
4761 case 2:
4762 rval = (*f)(argv+1, argv, &result);
4763 break;
4764 case 3:
4765 rval = (*f)(argv+2, argv+1, argv, &result);
4766 break;
4767 default:
4768 rval = FALSE;
4769 sysError("Too many arguments to arithmetic function");
4770 }
4771
4772 popArgvArithStack(argc PASS_LD);
4773
4774 if ( rval )
4775 pushArithStack(&result PASS_LD);
4776
4777 return rval;
4778 }
4779
4780 #endif /* O_COMPILE_ARITH */
4781
4782
4783 /*******************************
4784 * MISC INTERFACE *
4785 *******************************/
4786
4787 /* Evaluate a term to a 64-bit integer. Term is of type
4788 */
4789
4790 int
PL_eval_expression_to_int64_ex(term_t t,int64_t * val)4791 PL_eval_expression_to_int64_ex(term_t t, int64_t *val)
4792 { GET_LD
4793 number n;
4794 int rval;
4795
4796 if ( valueExpression(t, &n PASS_LD) )
4797 { if ( toIntegerNumber(&n, 0) )
4798 { switch(n.type)
4799 { case V_INTEGER:
4800 *val = n.value.i;
4801 rval = TRUE;
4802 break;
4803 #ifdef O_GMP
4804 case V_MPZ:
4805 { if ( !(rval=mpz_to_int64(n.value.mpz, val)) )
4806 rval = PL_error(NULL, 0, NULL, ERR_EVALUATION, ATOM_int_overflow);
4807 break;
4808 }
4809 #endif
4810 default:
4811 assert(0);
4812 return FALSE;
4813 }
4814 } else
4815 { rval = PL_error(NULL, 0, NULL, ERR_TYPE, ATOM_integer_expression, t);
4816 }
4817
4818 clearNumber(&n);
4819 } else
4820 { rval = FALSE;
4821 }
4822
4823 return rval;
4824 }
4825
4826
4827 /*******************************
4828 * PUBLISH PREDICATES *
4829 *******************************/
4830
4831 BeginPredDefs(arith)
4832 PRED_DEF("is", 2, is, PL_FA_ISO)
4833 PRED_DEF("<", 2, lt, PL_FA_ISO)
4834 PRED_DEF(">", 2, gt, PL_FA_ISO)
4835 PRED_DEF("=<", 2, leq, PL_FA_ISO)
4836 PRED_DEF(">=", 2, geq, PL_FA_ISO)
4837 PRED_DEF("=\\=", 2, neq, PL_FA_ISO)
4838 PRED_DEF("=:=", 2, eq, PL_FA_ISO)
4839 PRED_DEF("succ", 2, succ, 0)
4840 PRED_DEF("plus", 3, plus, 0)
4841 PRED_DEF("between", 3, between, PL_FA_NONDETERMINISTIC)
4842 PRED_DEF("bounded_number", 3, bounded_number, 0)
4843 PRED_DEF("float_class", 2, float_class, 0)
4844 PRED_DEF("float_parts", 4, float_parts, 0)
4845
4846 PRED_DEF("current_arithmetic_function", 1, current_arithmetic_function,
4847 PL_FA_NONDETERMINISTIC)
4848
4849 #ifdef O_GMP
4850 PRED_DEF("divmod", 4, divmod, 0)
4851 PRED_DEF("nth_integer_root_and_remainder", 4,
4852 nth_integer_root_and_remainder, 0)
4853 PRED_DEF("rational", 3, rational, 0)
4854 #endif
4855 PRED_DEF("set_random", 1, set_random, 0)
4856 #ifdef O_RANDOM_STATE
4857 PRED_DEF("random_property", 1, random_property, 0)
4858 #endif
4859 EndPredDefs
4860