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) 2005-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 /*#define O_DEBUG 1*/
38 #include <math.h>
39 #include <fenv.h>
40 #include <float.h>
41 #include "pl-incl.h"
42 #include "pl-arith.h"
43 #include "pl-inline.h"
44 #undef LD
45 #define LD LOCAL_LD
46
47 #ifdef O_GMP
48
49 static mpz_t MPZ_MIN_TAGGED; /* Prolog tagged integers */
50 static mpz_t MPZ_MAX_TAGGED;
51 static mpz_t MPZ_MIN_PLINT; /* Prolog int64_t integers */
52 static mpz_t MPZ_MAX_PLINT;
53 static mpz_t MPZ_MAX_UINT64;
54 #if SIZEOF_LONG < SIZEOF_VOIDP
55 static mpz_t MPZ_MIN_LONG; /* Prolog int64_t integers */
56 static mpz_t MPZ_MAX_LONG;
57 #endif
58
59
60 /*******************************
61 * MEMORY MANAGEMENT *
62 *******************************/
63
64 #if O_MY_GMP_ALLOC
65
66 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67 GMP doesn't (yet) allow for handling memory overflows. You can redefine
68 the allocation handles, but you are not allowed to return NULL or abort
69 the execution using longjmp(). As our normal GMP numbers live on the
70 global stack, we however can cleanup the temporary numbers that are
71 created during the Prolog function evaluation and use longjmp() through
72 STACK_OVERFLOW_THROW. Patrick Pelissier acknowledged this should work.
73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74
75 static int
gmp_too_big()76 gmp_too_big()
77 { GET_LD
78
79 DEBUG(1, Sdprintf("Signalling GMP overflow\n"));
80
81 return (int)outOfStack((Stack)&LD->stacks.global, STACK_OVERFLOW_THROW);
82 }
83
84 #define TOO_BIG_GMP(n) ((n) > 1000 && (n) > (size_t)globalStackLimit())
85
86 static void *
mp_alloc(size_t bytes)87 mp_alloc(size_t bytes)
88 { GET_LD
89 mp_mem_header *mem;
90
91 if ( LD->gmp.persistent )
92 return malloc(bytes);
93
94 if ( TOO_BIG_GMP(bytes) ||
95 !(mem = malloc(sizeof(mp_mem_header)+bytes)) )
96 { gmp_too_big();
97 abortProlog();
98 PL_rethrow();
99 return NULL; /* make compiler happy */
100 }
101
102 GMP_LEAK_CHECK(LD->gmp.allocated += bytes);
103
104 mem->next = NULL;
105 mem->context = LD->gmp.context;
106 if ( LD->gmp.tail )
107 { mem->prev = LD->gmp.tail;
108 LD->gmp.tail->next = mem;
109 LD->gmp.tail = mem;
110 } else
111 { mem->prev = NULL;
112 LD->gmp.head = LD->gmp.tail = mem;
113 }
114 DEBUG(9, Sdprintf("GMP: alloc %ld@%p\n", bytes, &mem[1]));
115
116 return &mem[1];
117 }
118
119
120 static void *
mp_realloc(void * ptr,size_t oldsize,size_t newsize)121 mp_realloc(void *ptr, size_t oldsize, size_t newsize)
122 { GET_LD
123 mp_mem_header *oldmem, *newmem;
124
125 if ( LD->gmp.persistent )
126 return realloc(ptr, newsize);
127
128 oldmem = ((mp_mem_header*)ptr)-1;
129 if ( TOO_BIG_GMP(newsize) ||
130 !(newmem = realloc(oldmem, sizeof(mp_mem_header)+newsize)) )
131 { gmp_too_big();
132 abortProlog();
133 PL_rethrow();
134 return NULL; /* make compiler happy */
135 }
136
137 if ( oldmem != newmem ) /* re-link if moved */
138 { if ( newmem->prev )
139 newmem->prev->next = newmem;
140 else
141 LD->gmp.head = newmem;
142
143 if ( newmem->next )
144 newmem->next->prev = newmem;
145 else
146 LD->gmp.tail = newmem;
147 }
148
149 GMP_LEAK_CHECK(LD->gmp.allocated -= oldsize;
150 LD->gmp.allocated += newsize);
151 DEBUG(9, Sdprintf("GMP: realloc %ld@%p --> %ld@%p\n", oldsize, ptr, newsize, &newmem[1]));
152
153 return &newmem[1];
154 }
155
156
157 static void
mp_free(void * ptr,size_t size)158 mp_free(void *ptr, size_t size)
159 { GET_LD
160 mp_mem_header *mem;
161
162 if ( LD->gmp.persistent )
163 { free(ptr);
164 return;
165 }
166
167 mem = ((mp_mem_header*)ptr)-1;
168
169 if ( mem == LD->gmp.head )
170 { LD->gmp.head = LD->gmp.head->next;
171 if ( LD->gmp.head )
172 LD->gmp.head->prev = NULL;
173 else
174 LD->gmp.tail = NULL;
175 } else if ( mem == LD->gmp.tail )
176 { LD->gmp.tail = LD->gmp.tail->prev;
177 LD->gmp.tail->next = NULL;
178 } else
179 { mem->prev->next = mem->next;
180 mem->next->prev = mem->prev;
181 }
182
183 free(mem);
184 DEBUG(9, Sdprintf("GMP: free: %ld@%p\n", size, ptr));
185 GMP_LEAK_CHECK(LD->gmp.allocated -= size);
186 }
187
188
189 void
mp_cleanup(ar_context * ctx)190 mp_cleanup(ar_context *ctx)
191 { GET_LD
192 mp_mem_header *mem, *next;
193
194 if ( LD->gmp.context )
195 { for(mem=LD->gmp.head; mem; mem=next)
196 { next = mem->next;
197 if ( mem->context == LD->gmp.context )
198 { DEBUG(9, Sdprintf("GMP: cleanup of %p\n", &mem[1]));
199 mp_free(&mem[1], 0);
200 }
201 }
202 }
203
204 LD->gmp.context = ctx->parent;
205 }
206 #endif
207
208
209 #ifdef __WINDOWS__
210 #undef isascii
211 int
isascii(int c)212 isascii(int c) /* missing from gmp.lib */
213 { return c >= 0 && c < 128;
214 }
215 #endif
216
217 /*******************************
218 * STACK MANAGEMENT *
219 *******************************/
220
221 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222 globalMPZ() pushes an mpz type GMP integer onto the local stack. The
223 saved version is the _mp_size field, followed by the limps.
224 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225
226 static size_t
mpz_wsize(mpz_t mpz,size_t * s)227 mpz_wsize(mpz_t mpz, size_t *s)
228 { DEBUG(0, assert(sizeof(mpz->_mp_size) == sizeof(int)));
229 size_t size = sizeof(mp_limb_t)*abs(mpz->_mp_size);
230 size_t wsz = (size+sizeof(word)-1)/sizeof(word);
231
232 if ( s )
233 *s = size;
234
235 return wsz;
236 }
237
238
239 static int
globalMPZ(Word at,mpz_t mpz,int flags ARG_LD)240 globalMPZ(Word at, mpz_t mpz, int flags ARG_LD)
241 { DEBUG(CHK_SECURE, assert(!onStackArea(global, at) && !onStackArea(local, at)));
242
243 if ( mpz->_mp_alloc )
244 { size_t size, wsz;
245 Word p;
246 word m;
247
248 copy:
249 wsz = mpz_wsize(mpz, &size);
250 m = mkIndHdr(wsz+1, TAG_INTEGER);
251
252 if ( wsizeofInd(m) != wsz+1 )
253 { PL_error(NULL, 0, NULL, ERR_REPRESENTATION, ATOM_integer);
254 return 0;
255 }
256
257 if ( !hasGlobalSpace(wsz+3) )
258 { int rc = ensureGlobalSpace(wsz+3, flags);
259
260 if ( rc != TRUE )
261 return rc;
262 }
263 p = gTop;
264 gTop += wsz+3;
265
266 *at = consPtr(p, TAG_INTEGER|STG_GLOBAL);
267
268 *p++ = m;
269 p[wsz] = 0L; /* pad out */
270 p[wsz+1] = m;
271 *p++ = mpz_size_stack(mpz->_mp_size);
272 memcpy(p, mpz->_mp_d, size);
273 } else /* already on the stack */
274 { Word p = (Word)mpz->_mp_d - 2;
275 if ( !onStack(global, p) )
276 goto copy;
277 #ifndef NDEBUG
278 size_t size;
279 size_t wsz = mpz_wsize(mpz, &size);
280 assert(p[0] == mkIndHdr(wsz+1, TAG_INTEGER));
281 #endif
282 *at = consPtr(p, TAG_INTEGER|STG_GLOBAL);
283 }
284
285 return TRUE;
286 }
287
288
289 static int
globalMPQ(Word at,mpq_t mpq,int flags ARG_LD)290 globalMPQ(Word at, mpq_t mpq, int flags ARG_LD)
291 { mpz_t num, den; /* num/den */
292
293 num[0] = *mpq_numref(mpq);
294 den[0] = *mpq_denref(mpq);
295
296 if ( num->_mp_alloc || den->_mp_alloc )
297 { size_t num_size, den_size, num_wsz, den_wsz;
298 Word p;
299 word m;
300
301 copy:
302 num_wsz = mpz_wsize(num, &num_size);
303 den_wsz = mpz_wsize(den, &den_size);
304 m = mkIndHdr(num_wsz+den_wsz+2, TAG_INTEGER);
305
306 if ( wsizeofInd(m) != num_wsz+den_wsz+2 )
307 { PL_error(NULL, 0, NULL, ERR_REPRESENTATION, ATOM_rational);
308 return 0;
309 }
310
311 if ( !hasGlobalSpace(num_wsz+den_wsz+4) )
312 { int rc = ensureGlobalSpace(num_wsz+den_wsz+4, flags);
313
314 if ( rc != TRUE )
315 return rc;
316 }
317 p = gTop;
318 gTop += num_wsz+den_wsz+4;
319
320 *at = consPtr(p, TAG_INTEGER|STG_GLOBAL);
321 *p++ = m;
322 *p++ = mpq_size_stack(num->_mp_size);
323 *p++ = mpq_size_stack(den->_mp_size);
324 p[num_wsz-1] = 0L; /* pad out */
325 memcpy(p, num->_mp_d, num_size);
326 p += num_wsz;
327 p[den_wsz-1] = 0L; /* pad out */
328 memcpy(p, den->_mp_d, den_size);
329 p += den_wsz;
330 *p = m;
331 } else /* already on the stack */
332 { Word p = (Word)num->_mp_d - 3;
333 if ( !onStack(global, p) )
334 goto copy;
335 DEBUG(CHK_SECURE,
336 { size_t num_size;
337 size_t den_size;
338 size_t num_wsz = mpz_wsize(num, &num_size);
339 size_t den_wsz = mpz_wsize(den, &den_size);
340 assert(p[0] == mkIndHdr(num_wsz+den_wsz+2, TAG_INTEGER));
341 assert((Word)den->_mp_d == (Word)num->_mp_d + num_wsz);
342 });
343 *at = consPtr(p, TAG_INTEGER|STG_GLOBAL);
344 }
345
346 return TRUE;
347 }
348
349
350 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351 get_integer() fetches the value of a Prolog term known to be an integer
352 into a number structure. If the value is a MPZ number, it must be
353 handled as read-only and it only be used as intptr_t as no calls are made
354 that may force a relocation or garbage collection on the global stack.
355
356 The version without O_GMP is a macro defined in pl-gmp.h
357 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
358
359 void
get_integer(word w,Number n)360 get_integer(word w, Number n)
361 { if ( storage(w) == STG_INLINE )
362 { n->type = V_INTEGER,
363 n->value.i = valInt(w);
364 } else
365 { GET_LD
366 Word p = addressIndirect(w);
367 size_t wsize = wsizeofInd(*p);
368
369 p++;
370 if ( wsize == WORDS_PER_INT64 )
371 { n->type = V_INTEGER;
372 memcpy(&n->value.i, p, sizeof(int64_t));
373 } else
374 { n->type = V_MPZ;
375
376 n->value.mpz->_mp_size = mpz_stack_size(*p++);
377 n->value.mpz->_mp_alloc = 0;
378 n->value.mpz->_mp_d = (mp_limb_t*) p;
379 }
380 }
381 }
382
383
384 void
get_rational(word w,Number n)385 get_rational(word w, Number n)
386 { if ( storage(w) == STG_INLINE )
387 { n->type = V_INTEGER,
388 n->value.i = valInt(w);
389 } else
390 { GET_LD
391 Word p = addressIndirect(w);
392 size_t wsize = wsizeofInd(*p);
393
394 p++;
395 if ( wsize == WORDS_PER_INT64 )
396 { n->type = V_INTEGER;
397 memcpy(&n->value.i, p, sizeof(int64_t));
398 } else if ( (*p&MP_RAT_MASK) )
399 { mpz_t num, den;
400 size_t num_size;
401
402 n->type = V_MPQ;
403 num->_mp_size = mpz_stack_size(*p++);
404 num->_mp_alloc = 0;
405 num->_mp_d = (mp_limb_t*) (p+1);
406 num_size = mpz_wsize(num, NULL);
407 den->_mp_size = mpz_stack_size(*p++);
408 den->_mp_alloc = 0;
409 den->_mp_d = (mp_limb_t*) (p+num_size);
410
411 *mpq_numref(n->value.mpq) = num[0];
412 *mpq_denref(n->value.mpq) = den[0];
413 } else
414 { n->type = V_MPZ;
415
416 n->value.mpz->_mp_size = mpz_stack_size(*p++);
417 n->value.mpz->_mp_alloc = 0;
418 n->value.mpz->_mp_d = (mp_limb_t*) p;
419 }
420 }
421 }
422
423
424 Code
get_mpz_from_code(Code pc,mpz_t mpz)425 get_mpz_from_code(Code pc, mpz_t mpz)
426 { size_t wsize = wsizeofInd(*pc);
427
428 pc++;
429 mpz->_mp_size = mpz_stack_size(*pc);
430 mpz->_mp_alloc = 0;
431 mpz->_mp_d = (mp_limb_t*)(pc+1);
432
433 return pc+wsize;
434 }
435
436 Code
get_mpq_from_code(Code pc,mpq_t mpq)437 get_mpq_from_code(Code pc, mpq_t mpq)
438 { Word p = pc;
439 size_t wsize = wsizeofInd(*p);
440 p++;
441 int num_size = mpz_stack_size(*p++);
442 int den_size = mpz_stack_size(*p++);
443 size_t limpsize;
444
445 mpq_numref(mpq)->_mp_size = num_size;
446 mpq_denref(mpq)->_mp_size = den_size;
447 mpq_numref(mpq)->_mp_alloc = 0;
448 mpq_denref(mpq)->_mp_alloc = 0;
449 mpq_numref(mpq)->_mp_d = (mp_limb_t*)p;
450 limpsize = sizeof(mp_limb_t) * abs(num_size);
451 p += (limpsize+sizeof(word)-1)/sizeof(word);
452 mpq_denref(mpq)->_mp_d = (mp_limb_t*)p;
453
454 return pc+wsize+1;
455 }
456
457
458 /*******************************
459 * IMPORT/EXPORT *
460 *******************************/
461
462 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463 addMPZToBuffer(Buffer b, mpz_t mpz)
464 Add mpz in a machine independent representation to the given buffer.
465 The data is stored in limps of 1 byte, preceeded by the byte-count
466 as 4-byte big-endian number;
467
468 addMPQToBuffer(Buffer b, mpq_t mpq)
469 Similar to mpz, but first writes the two headers and then the
470 two bit patterns.
471 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
472
473 static void
add_mpz_size_buffer(Buffer b,mpz_t mpz,size_t size)474 add_mpz_size_buffer(Buffer b, mpz_t mpz, size_t size)
475 { ssize_t hdrsize;
476
477 if ( mpz_sgn(mpz) < 0 )
478 hdrsize = -(ssize_t)size;
479 else
480 hdrsize = (ssize_t)size;
481
482 *b->top++ = (char)((hdrsize>>24)&0xff);
483 *b->top++ = (char)((hdrsize>>16)&0xff);
484 *b->top++ = (char)((hdrsize>> 8)&0xff);
485 *b->top++ = (char)((hdrsize )&0xff);
486 }
487
488 static void
add_mpz_bits_buffer(Buffer b,mpz_t mpz,size_t size)489 add_mpz_bits_buffer(Buffer b, mpz_t mpz, size_t size)
490 { size_t count;
491
492 mpz_export(b->top, &count, 1, 1, 1, 0, mpz);
493 assert(count == size);
494 b->top = b->top + size;
495 }
496
497 void
addMPZToBuffer(Buffer b,mpz_t mpz)498 addMPZToBuffer(Buffer b, mpz_t mpz)
499 { size_t size = (mpz_sizeinbase(mpz, 2)+7)/8;
500
501 if ( !growBuffer(b, size+4) )
502 outOfCore();
503 add_mpz_size_buffer(b, mpz, size);
504 add_mpz_bits_buffer(b, mpz, size);
505 }
506
507 void
addMPQToBuffer(Buffer b,mpq_t mpq)508 addMPQToBuffer(Buffer b, mpq_t mpq)
509 { mpz_t num, den; /* num/den */
510 size_t num_size, den_size;
511
512 num[0] = *mpq_numref(mpq);
513 den[0] = *mpq_denref(mpq);
514 num_size = (mpz_sizeinbase(num, 2)+7)/8;
515 den_size = (mpz_sizeinbase(den, 2)+7)/8;
516
517 if ( !growBuffer(b, num_size+den_size+8) )
518 outOfCore();
519
520 add_mpz_size_buffer(b, num, num_size);
521 add_mpz_size_buffer(b, den, den_size);
522 add_mpz_bits_buffer(b, num, num_size);
523 add_mpz_bits_buffer(b, den, den_size);
524 }
525
526 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
527 loadMPZFromCharp() loads an MPZ number directly back to the global stack
528 from a char * as filled by addMPZToBuffer(). Memory allocation is
529 avoided by creating a dummy mpz that looks big enough to mpz_import().
530 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
531
532 #define SHIFTSIGN32 ((sizeof(int)-4)*8)
533
534 static char *
load_mpz_size(const char * data,int * szp)535 load_mpz_size(const char *data, int *szp)
536 { int size = 0;
537
538 size |= (data[0]&0xff)<<24;
539 size |= (data[1]&0xff)<<16;
540 size |= (data[2]&0xff)<<8;
541 size |= (data[3]&0xff);
542 size = (size << SHIFTSIGN32)>>SHIFTSIGN32; /* sign extend */
543
544 *szp = size;
545 data += 4;
546
547 return (char *)data;
548 }
549
550 static char *
load_abs_mpz_size(const char * data,int * szp,int * neg)551 load_abs_mpz_size(const char *data, int *szp, int *neg)
552 { int size;
553
554 data = load_mpz_size(data, &size);
555 if ( neg )
556 { if ( size < 0 )
557 { size = -size;
558 *neg = TRUE;
559 } else
560 { *neg = FALSE;
561 }
562 } else if ( size < 0 )
563 size = -size;
564 *szp = size;
565
566 return (char *)data;
567 }
568
569 static char *
load_mpz_bits(const char * data,size_t size,size_t limpsize,int neg,Word p)570 load_mpz_bits(const char *data, size_t size, size_t limpsize, int neg, Word p)
571 { mpz_t mpz;
572
573 mpz->_mp_size = limpsize;
574 mpz->_mp_alloc = limpsize;
575 mpz->_mp_d = (mp_limb_t*)p;
576
577 mpz_import(mpz, size, 1, 1, 1, 0, data);
578 assert((Word)mpz->_mp_d == p); /* check no (re-)allocation is done */
579
580 return (char*)data+size;
581 }
582
583 char *
loadMPZFromCharp(const char * data,Word r,Word * store)584 loadMPZFromCharp(const char *data, Word r, Word *store)
585 { GET_LD
586 int size = 0;
587 size_t limpsize;
588 size_t wsize;
589 int neg;
590 Word p;
591 word m;
592
593 data = load_abs_mpz_size(data, &size, &neg);
594
595 limpsize = (size+sizeof(mp_limb_t)-1)/sizeof(mp_limb_t);
596 wsize = (limpsize*sizeof(mp_limb_t)+sizeof(word)-1)/sizeof(word);
597 p = *store;
598 *store += (wsize+3);
599 *r = consPtr(p, TAG_INTEGER|STG_GLOBAL);
600 m = mkIndHdr(wsize+1, TAG_INTEGER);
601 *p++ = m;
602 p[wsize] = 0L; /* pad out */
603 p[wsize+1] = m;
604 *p++ = mpz_size_stack(neg ? -limpsize : limpsize);
605
606 return load_mpz_bits(data, size, limpsize, neg, p);
607 }
608
609 char *
loadMPQFromCharp(const char * data,Word r,Word * store)610 loadMPQFromCharp(const char *data, Word r, Word *store)
611 { GET_LD
612 int num_size;
613 int den_size;
614 size_t num_limpsize, num_wsize;
615 size_t den_limpsize, den_wsize;
616 int num_neg, den_neg;
617 size_t wsize;
618 Word p;
619 word m;
620
621 data = load_abs_mpz_size(data, &num_size, &num_neg);
622 data = load_abs_mpz_size(data, &den_size, &den_neg);
623
624 num_limpsize = (num_size+sizeof(mp_limb_t)-1)/sizeof(mp_limb_t);
625 num_wsize = (num_limpsize*sizeof(mp_limb_t)+sizeof(word)-1)/sizeof(word);
626 den_limpsize = (den_size+sizeof(mp_limb_t)-1)/sizeof(mp_limb_t);
627 den_wsize = (den_limpsize*sizeof(mp_limb_t)+sizeof(word)-1)/sizeof(word);
628 wsize = num_wsize+den_wsize;
629
630 p = *store;
631 *store += (wsize+4);
632 *r = consPtr(p, TAG_INTEGER|STG_GLOBAL);
633 m = mkIndHdr(wsize+2, TAG_INTEGER);
634 *p++ = m;
635 *p++ = mpq_size_stack(num_neg ? -num_limpsize : num_limpsize);
636 *p++ = mpq_size_stack(den_neg ? -den_limpsize : den_limpsize);
637 p[num_wsize-1] = 0;
638 data = load_mpz_bits(data, num_size, num_limpsize, num_neg, p);
639 p += num_wsize;
640 p[den_wsize-1] = 0;
641 data = load_mpz_bits(data, den_size, den_limpsize, den_neg, p);
642 p += den_wsize;
643 *p++ = m;
644 assert(p == *store);
645
646 return (char *)data;
647 }
648
649 char *
skipMPZOnCharp(const char * data)650 skipMPZOnCharp(const char *data)
651 { int size;
652
653 data = load_abs_mpz_size(data, &size, NULL);
654
655 return (char *)data + size;
656 }
657
658 char *
skipMPQOnCharp(const char * data)659 skipMPQOnCharp(const char *data)
660 { int num_size;
661 int den_size;
662
663 data = load_abs_mpz_size(data, &num_size, NULL);
664 data = load_abs_mpz_size(data, &den_size, NULL);
665 data += num_size;
666 data += den_size;
667
668 return (char *)data;
669 }
670
671 #undef SHIFTSIGN32
672
673
674 /*******************************
675 * CONVERSION *
676 *******************************/
677
678 #ifdef WORDS_BIGENDIAN
679 #define ORDER 1
680 #else
681 #define ORDER -1
682 #endif
683
684 void
mpz_init_set_si64(mpz_t mpz,int64_t i)685 mpz_init_set_si64(mpz_t mpz, int64_t i)
686 {
687 #if SIZEOF_LONG == 8
688 mpz_init_set_si(mpz, (long)i);
689 #else
690 DEBUG(2, Sdprintf("Converting %" PRId64 " to MPZ\n", i));
691
692 if ( i >= LONG_MIN && i <= LONG_MAX )
693 { mpz_init_set_si(mpz, (long)i);
694 } else
695 { mpz_init(mpz);
696 if ( i >= 0 )
697 { mpz_import(mpz, sizeof(i), ORDER, 1, 0, 0, &i);
698 } else
699 { i = -i;
700 mpz_import(mpz, sizeof(i), ORDER, 1, 0, 0, &i);
701 mpz_neg(mpz, mpz);
702 }
703 }
704 DEBUG(2, gmp_printf("\t--> %Zd\n", mpz));
705 #endif
706 }
707
708
709 static void
mpz_init_set_uint64(mpz_t mpz,uint64_t i)710 mpz_init_set_uint64(mpz_t mpz, uint64_t i)
711 {
712 #if SIZEOF_LONG == 8
713 mpz_init_set_ui(mpz, (unsigned long)i);
714 #else
715 mpz_init(mpz);
716 mpz_import(mpz, sizeof(i), ORDER, 1, 0, 0, &i);
717 #endif
718 }
719
720
721 static void
mpz_init_max_uint(mpz_t mpz,int bits)722 mpz_init_max_uint(mpz_t mpz, int bits)
723 { mpz_init_set_si(mpz, 1);
724 mpz_mul_2exp(mpz, mpz, bits);
725 mpz_sub_ui(mpz, mpz, 1);
726 }
727
728
729 int
promoteToMPZNumber(number * n)730 promoteToMPZNumber(number *n)
731 { switch(n->type)
732 { case V_INTEGER:
733 mpz_init_set_si64(n->value.mpz, n->value.i);
734 n->type = V_MPZ;
735 break;
736 case V_MPZ:
737 break;
738 case V_MPQ:
739 { mpz_t mpz;
740
741 mpz_init(mpz);
742 mpz_tdiv_q(mpz,
743 mpq_numref(n->value.mpq),
744 mpq_denref(n->value.mpq));
745 clearNumber(n);
746 n->type = V_MPZ;
747 n->value.mpz[0] = mpz[0];
748 break;
749 }
750 case V_FLOAT:
751 mpz_init_set_d(n->value.mpz, n->value.f);
752 n->type = V_MPZ;
753 break;
754 }
755
756 return TRUE;
757 }
758
759
760 int
promoteToMPQNumber(number * n)761 promoteToMPQNumber(number *n)
762 { switch(n->type)
763 { case V_INTEGER:
764 promoteToMPZNumber(n);
765 /*FALLTHOURGH*/
766 case V_MPZ:
767 { n->value.mpq->_mp_num = n->value.mpz[0];
768 mpz_init_set_ui(mpq_denref(n->value.mpq), 1L);
769 n->type = V_MPQ;
770 break;
771 }
772 case V_MPQ:
773 break;
774 case V_FLOAT:
775 { double v = n->value.f;
776
777 switch(fpclassify(v))
778 { case FP_NAN:
779 return PL_error(NULL, 0, NULL, ERR_AR_UNDEF);
780 case FP_INFINITE:
781 return PL_error(NULL, 0, NULL, ERR_AR_RAT_OVERFLOW);
782 }
783
784 n->type = V_MPQ;
785 mpq_init(n->value.mpq);
786 mpq_set_d(n->value.mpq, v);
787 break;
788 }
789 }
790
791 return TRUE;
792 }
793
794
795 /*******************************
796 * RW *
797 *******************************/
798
799 void
ensureWritableNumber(Number n)800 ensureWritableNumber(Number n)
801 { switch(n->type)
802 { case V_MPZ:
803 if ( !n->value.mpz->_mp_alloc )
804 { mpz_t tmp;
805
806 tmp[0] = n->value.mpz[0];
807 mpz_init_set(n->value.mpz, tmp);
808 break;
809 }
810 case V_MPQ:
811 { if ( !mpq_numref(n->value.mpq)->_mp_alloc )
812 { mpz_t tmp;
813
814 tmp[0] = mpq_numref(n->value.mpq)[0];
815 mpz_init_set(mpq_numref(n->value.mpq), tmp);
816 }
817 if ( !mpq_denref(n->value.mpq)->_mp_alloc )
818 { mpz_t tmp;
819
820 tmp[0] = mpq_denref(n->value.mpq)[0];
821 mpz_init_set(mpq_denref(n->value.mpq), tmp);
822 }
823 break;
824 }
825 default:
826 break;
827 }
828 }
829
830
831
832 /*******************************
833 * CLEAR *
834 *******************************/
835
836 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
837 Numbers may contain two type of MPZ numbers. Ones that are created by
838 the GMP library and must be cleared, and ones that have their `limbs'
839 stored somewhere in the Prolog memory. These may only be used read-only
840 and their _mp_alloc field is set to 0. clearNumber() discards MPZ
841 numbers that are created by GMP only.
842 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
843
844 void
clearGMPNumber(Number n)845 clearGMPNumber(Number n)
846 { switch(n->type)
847 { case V_MPZ:
848 if ( n->value.mpz->_mp_alloc )
849 mpz_clear(n->value.mpz);
850 break;
851 case V_MPQ:
852 if ( mpq_numref(n->value.mpq)->_mp_alloc )
853 mpz_clear(mpq_numref(n->value.mpq));
854 if ( mpq_denref(n->value.mpq)->_mp_alloc )
855 mpz_clear(mpq_denref(n->value.mpq));
856 break;
857 default:
858 break;
859 }
860 }
861
862
863 /*******************************
864 * INIT *
865 *******************************/
866
867 void
initGMP(void)868 initGMP(void)
869 { if ( !GD->gmp.initialised )
870 { GD->gmp.initialised = TRUE;
871
872 mpz_init_set_si64(MPZ_MIN_TAGGED, PLMINTAGGEDINT);
873 mpz_init_set_si64(MPZ_MAX_TAGGED, PLMAXTAGGEDINT);
874 mpz_init_set_si64(MPZ_MIN_PLINT, PLMININT);
875 mpz_init_set_si64(MPZ_MAX_PLINT, PLMAXINT);
876 mpz_init_max_uint(MPZ_MAX_UINT64, 64);
877 #if SIZEOF_LONG < SIZEOF_VOIDP
878 mpz_init_set_si64(MPZ_MIN_LONG, LONG_MIN);
879 mpz_init_set_si64(MPZ_MAX_LONG, LONG_MAX);
880 #endif
881 #ifdef O_MY_GMP_ALLOC
882 if ( !GD->gmp.keep_alloc_functions )
883 mp_set_memory_functions(mp_alloc, mp_realloc, mp_free);
884 #endif
885
886 #if __GNU_MP__ > 3 && __GNU_MP__ < 6
887 PL_license("lgplv3", "libgmp");
888 #else
889 PL_license("lgplv2+", "libgmp");
890 #endif
891 }
892 }
893
894
895 void
cleanupGMP()896 cleanupGMP()
897 { if ( GD->gmp.initialised )
898 { GD->gmp.initialised = FALSE;
899
900 #ifdef O_MY_GMP_ALLOC
901 if ( !GD->gmp.keep_alloc_functions )
902 mp_set_memory_functions(NULL, NULL, NULL);
903 #endif
904 mpz_clear(MPZ_MIN_TAGGED);
905 mpz_clear(MPZ_MAX_TAGGED);
906 mpz_clear(MPZ_MIN_PLINT);
907 mpz_clear(MPZ_MAX_PLINT);
908 #if SIZEOF_LONG < SIZEOF_VOIDP
909 mpz_clear(MPZ_MIN_LONG);
910 mpz_clear(MPZ_MAX_LONG);
911 #endif
912 }
913 }
914
915
916 /*******************************
917 * NUMBER HANDLING *
918 *******************************/
919
920 int
mpz_to_int64(mpz_t mpz,int64_t * i)921 mpz_to_int64(mpz_t mpz, int64_t *i)
922 { if ( mpz_cmp(mpz, MPZ_MIN_PLINT) >= 0 &&
923 mpz_cmp(mpz, MPZ_MAX_PLINT) <= 0 )
924 { uint64_t v;
925
926 mpz_export(&v, NULL, ORDER, sizeof(v), 0, 0, mpz);
927 DEBUG(2,
928 { char buf[256];
929 Sdprintf("Convert %s --> %I64d\n",
930 mpz_get_str(buf, 10, mpz), v);
931 });
932
933 if ( mpz_sgn(mpz) < 0 )
934 *i = -(int64_t)(v - 1) - 1;
935 else
936 *i = (int64_t)v;
937
938 return TRUE;
939 }
940
941 return FALSE;
942 }
943
944
945 /* return: <0: -1
946 >MPZ_UINT64_MAX: 1
947 (ok) 0
948 */
949
950 int
mpz_to_uint64(mpz_t mpz,uint64_t * i)951 mpz_to_uint64(mpz_t mpz, uint64_t *i)
952 { if ( mpz_sgn(mpz) < 0 )
953 return -1;
954
955 if ( mpz_cmp(mpz, MPZ_MAX_UINT64) <= 0 )
956 { uint64_t v;
957
958 mpz_export(&v, NULL, ORDER, sizeof(v), 0, 0, mpz);
959 *i = v;
960
961 return 0;
962 }
963
964 return 1;
965 }
966
967
968 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
969 put_number() transforms a number into a Prolog term. Note that this may
970 allocate on the global stack. Please note that this function uses the
971 most compact representation, which is essential to make unify() work
972 without any knowledge of the represented data.
973 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
974
975 static int
put_mpz(Word at,mpz_t mpz,int flags ARG_LD)976 put_mpz(Word at, mpz_t mpz, int flags ARG_LD)
977 { int64_t v;
978
979 DEBUG(2,
980 { char buf[256];
981 Sdprintf("put_mpz(%s)\n",
982 mpz_get_str(buf, 10, mpz));
983 });
984
985 #if SIZEOF_LONG < SIZEOF_VOIDP
986 if ( mpz_cmp(mpz, MPZ_MIN_LONG) >= 0 &&
987 mpz_cmp(mpz, MPZ_MAX_LONG) <= 0 )
988 #else
989 if ( mpz_cmp(mpz, MPZ_MIN_TAGGED) >= 0 &&
990 mpz_cmp(mpz, MPZ_MAX_TAGGED) <= 0 )
991 #endif
992 { long v = mpz_get_si(mpz);
993
994 if ( !hasGlobalSpace(0) ) /* ensure we have room for bindConst */
995 { int rc = ensureGlobalSpace(0, flags);
996
997 if ( rc != TRUE )
998 return rc;
999 }
1000
1001 *at = consInt(v);
1002 return TRUE;
1003 } else if ( mpz_to_int64(mpz, &v) )
1004 { return put_int64(at, v, flags PASS_LD);
1005 } else
1006 { return globalMPZ(at, mpz, flags PASS_LD);
1007 }
1008 }
1009
1010 #endif /*O_GMP*/
1011
1012 /* returns one of
1013
1014 TRUE: ok
1015 FALSE: some error
1016 GLOBAL_OVERFLOW: no space
1017 LOCAL_OVERFLOW: cannot represent (no GMP)
1018 */
1019
1020 int
put_uint64(Word at,uint64_t l,int flags ARG_LD)1021 put_uint64(Word at, uint64_t l, int flags ARG_LD)
1022 { if ( (int64_t)l >= 0 )
1023 { return put_int64(at, l, flags PASS_LD);
1024 } else
1025 {
1026 #ifdef O_GMP
1027 mpz_t mpz;
1028
1029 mpz_init_set_uint64(mpz, l);
1030 return put_mpz(at, mpz, flags PASS_LD);
1031 #else
1032 return LOCAL_OVERFLOW;
1033 #endif
1034 }
1035 }
1036
1037
1038 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1039 put_number() translates a number structure into its Prolog
1040 representation and ensures there is enough space for a subsequent
1041 bindConst() call. Note that `at' must point to an address that is not
1042 affected by GC/shift. The intented scenario is:
1043
1044 { word c;
1045
1046 if ( (rc=put_number(&c, n, ALLOW_GC PASS_LD)) == TRUE )
1047 bindConst(<somewhere>, c);
1048 ...
1049 }
1050 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1051
1052 int
put_number(Word at,Number n,int flags ARG_LD)1053 put_number(Word at, Number n, int flags ARG_LD)
1054 { switch(n->type)
1055 { case V_INTEGER:
1056 { word w = consInt(n->value.i);
1057
1058 if ( valInt(w) == n->value.i )
1059 { if ( !hasGlobalSpace(0) )
1060 { int rc = ensureGlobalSpace(0, flags);
1061
1062 if ( rc != TRUE )
1063 return rc;
1064 }
1065
1066 *at = w;
1067 return TRUE;
1068 }
1069
1070 return put_int64(at, n->value.i, flags PASS_LD);
1071 }
1072 #ifdef O_GMP
1073 case V_MPZ:
1074 return put_mpz(at, n->value.mpz, flags PASS_LD);
1075 case V_MPQ:
1076 { if ( mpz_cmp_ui(mpq_denref(n->value.mpq), 1L) == 0 )
1077 return put_mpz(at, mpq_numref(n->value.mpq), flags PASS_LD);
1078 else
1079 return globalMPQ(at, n->value.mpq, flags PASS_LD);
1080 }
1081 #endif
1082 case V_FLOAT:
1083 return put_double(at, n->value.f, flags PASS_LD);
1084 }
1085
1086 assert(0);
1087 return FALSE;
1088 }
1089
1090
1091 int
PL_unify_number__LD(term_t t,Number n ARG_LD)1092 PL_unify_number__LD(term_t t, Number n ARG_LD)
1093 { Word p = valTermRef(t);
1094
1095 deRef(p);
1096
1097 if ( canBind(*p) )
1098 { word w;
1099 int rc;
1100
1101 if ( (rc=put_number(&w, n, ALLOW_GC PASS_LD)) != TRUE )
1102 return raiseStackOverflow(rc);
1103
1104 p = valTermRef(t); /* put_number can shift the stacks */
1105 deRef(p);
1106
1107 bindConst(p, w);
1108 succeed;
1109 }
1110
1111 switch(n->type)
1112 { case V_INTEGER:
1113 if ( isTaggedInt(*p) )
1114 return valInt(*p) == n->value.i;
1115 /*FALLTHOURGH*/
1116 #ifdef O_GMP
1117 case V_MPZ:
1118 #endif
1119 if ( isInteger(*p) )
1120 { number n2;
1121 int rc;
1122
1123 get_integer(*p, &n2);
1124 rc = (cmpNumbers(n, &n2) == CMP_EQUAL);
1125 clearNumber(&n2);
1126
1127 return rc;
1128 }
1129 break;
1130 #ifdef O_GMP
1131 case V_MPQ:
1132 { if ( isRational(*p) )
1133 { number n2;
1134 int rc;
1135
1136 get_rational(*p, &n2);
1137 rc = (cmpNumbers(n, &n2) == CMP_EQUAL);
1138 clearNumber(&n2);
1139
1140 return rc;
1141 }
1142 break;
1143 }
1144 #endif
1145 case V_FLOAT:
1146 if ( isFloat(*p) )
1147 { Word ap = valIndirectP(*p);
1148
1149 return memcmp((char*)&n->value.f, ap, sizeof(n->value.f)) == 0;
1150 }
1151 break;
1152 }
1153
1154 fail;
1155 }
1156
1157
1158 int
PL_put_number__LD(term_t t,Number n ARG_LD)1159 PL_put_number__LD(term_t t, Number n ARG_LD)
1160 { word w;
1161 int rc;
1162
1163 if ( (rc=put_number(&w, n, ALLOW_GC PASS_LD)) != TRUE )
1164 return raiseStackOverflow(rc);
1165
1166 *valTermRef(t) = w;
1167
1168 return TRUE;
1169 }
1170
1171
1172 void
get_number(word w,Number n ARG_LD)1173 get_number(word w, Number n ARG_LD)
1174 { if ( isRational(w) )
1175 { get_rational(w, n);
1176 } else
1177 { n->type = V_FLOAT;
1178 n->value.f = valFloat(w);
1179 }
1180 }
1181
1182
1183 int
PL_get_number__LD(term_t t,Number n ARG_LD)1184 PL_get_number__LD(term_t t, Number n ARG_LD)
1185 { Word p = valTermRef(t);
1186
1187 deRef(p);
1188 if ( isRational(*p) )
1189 { get_rational(*p, n);
1190 succeed;
1191 }
1192 if ( isFloat(*p) )
1193 { n->value.f = valFloat(*p);
1194 n->type = V_FLOAT;
1195 succeed;
1196 }
1197
1198 fail;
1199 }
1200
1201 #undef PL_get_number
1202 int
PL_get_number(term_t t,Number n)1203 PL_get_number(term_t t, Number n)
1204 { GET_LD
1205
1206 return PL_get_number__LD(t, n PASS_LD);
1207 }
1208 #define PL_get_number(t, n) PL_get_number__LD(t, n PASS_LD)
1209
1210
1211 /*******************************
1212 * PROMOTION *
1213 *******************************/
1214
1215 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1216 The GMP functions which convert mpz's and mpq's to double's truncate
1217 (round to zero) if necessary ignoring the current IEEE rounding mode. To
1218 correct this incorrect rounding when the mode is `to_postive` or
1219 `to_negative` all calls to mpX_get_d() should be wrapped in this
1220 function. Note that this function has no GMP dependencies.
1221 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1222
1223 double
mpX_round(double d)1224 mpX_round(double d) {
1225 switch (fegetround()) {
1226 case FE_UPWARD : return (d>=0) ? nexttoward(d, INFINITY) : d;
1227 case FE_DOWNWARD: return (d<=0) ? nexttoward(d,-INFINITY) : d;
1228 default: return d;
1229 }
1230 }
1231
1232 int
promoteToFloatNumber(Number n)1233 promoteToFloatNumber(Number n)
1234 { switch(n->type)
1235 { case V_INTEGER:
1236 n->value.f = (double)n->value.i;
1237 n->type = V_FLOAT;
1238 break;
1239 #ifdef O_GMP
1240 case V_MPZ:
1241 { double val = mpX_round(mpz_get_d(n->value.mpz));
1242
1243 clearNumber(n);
1244 n->value.f = val;
1245 n->type = V_FLOAT;
1246 break;
1247 }
1248 case V_MPQ:
1249 { double val = mpq_to_double(n->value.mpq);
1250
1251 clearNumber(n);
1252 n->value.f = val;
1253 n->type = V_FLOAT;
1254 break;
1255 }
1256 #endif
1257 case V_FLOAT:
1258 break;
1259 }
1260
1261 return check_float(n);
1262 }
1263
1264
1265 int
promoteNumber(Number n,numtype t)1266 promoteNumber(Number n, numtype t)
1267 { switch(t)
1268 { case V_INTEGER:
1269 return TRUE;
1270 #ifdef O_GMP
1271 case V_MPZ:
1272 return promoteToMPZNumber(n);
1273 case V_MPQ:
1274 return promoteToMPQNumber(n);
1275 #endif
1276 case V_FLOAT:
1277 return promoteToFloatNumber(n);
1278 default:
1279 assert(0);
1280 return FALSE;
1281 }
1282 }
1283
1284
1285 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1286 same_type_numbers(n1, n2)
1287 Upgrade both numbers to the `highest' type of both. Number types are
1288 defined in the enum-type numtype, which is supposed to define a
1289 total ordering between the number types.
1290 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1291
1292 int
make_same_type_numbers(Number n1,Number n2)1293 make_same_type_numbers(Number n1, Number n2)
1294 { if ( (int)n1->type > (int)n2->type )
1295 return promoteNumber(n2, n1->type);
1296 else
1297 return promoteNumber(n1, n2->type);
1298 }
1299
1300
1301 /*******************************
1302 * OPERATIONS *
1303 *******************************/
1304
1305 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1306 cmpNumbers() compares two numbers. First, both numbers are promoted to
1307 the lowest (in V_* ordering) common type, after which they are compared.
1308 Note that if the common type is V_FLOAT, but not both are V_FLOAT we can
1309 run into troubles because big integers may be out of range for the
1310 double representation. We trust mpz_get_d() and mpq_to_double() return
1311 +/- float infinity and this compares correctly with the other argument,
1312 which is guaranteed to be a valid float.
1313 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1314
1315 static int
cmpFloatNumbers(Number n1,Number n2)1316 cmpFloatNumbers(Number n1, Number n2)
1317 { if ( n1->type == V_FLOAT )
1318 { double d2;
1319
1320 if ( isnan(n1->value.f) )
1321 return CMP_NOTEQ;
1322
1323 switch(n2->type)
1324 { case V_INTEGER:
1325 d2 = (double)n2->value.i;
1326 break;
1327 #ifdef O_GMP
1328 case V_MPZ:
1329 d2 = mpX_round(mpz_get_d(n2->value.mpz));
1330 break;
1331 case V_MPQ:
1332 d2 = mpq_to_double(n2->value.mpq);
1333 break;
1334 #endif
1335 default:
1336 assert(0);
1337 d2 = 0.0;
1338 }
1339
1340 return n1->value.f < d2 ? CMP_LESS :
1341 n1->value.f == d2 ? CMP_EQUAL : CMP_GREATER;
1342 } else
1343 { double d1;
1344
1345 assert(n2->type == V_FLOAT);
1346
1347 if ( isnan(n2->value.f) )
1348 return CMP_NOTEQ;
1349
1350 switch(n1->type)
1351 { case V_INTEGER:
1352 d1 = (double)n1->value.i;
1353 break;
1354 #ifdef O_GMP
1355 case V_MPZ:
1356 d1 = mpX_round(mpz_get_d(n1->value.mpz));
1357 break;
1358 case V_MPQ:
1359 d1 = mpq_to_double(n1->value.mpq);
1360 break;
1361 #endif
1362 default:
1363 assert(0);
1364 d1 = 0.0;
1365 }
1366
1367 return n2->value.f < d1 ? CMP_GREATER :
1368 n2->value.f == d1 ? CMP_EQUAL : CMP_LESS;
1369 }
1370 }
1371
1372
1373 int
cmpNumbers(Number n1,Number n2)1374 cmpNumbers(Number n1, Number n2)
1375 { if ( n1->type != n2->type )
1376 { int rc;
1377
1378 if ( n1->type == V_FLOAT || n2->type == V_FLOAT )
1379 return cmpFloatNumbers(n1, n2);
1380 rc = make_same_type_numbers(n1, n2);
1381 assert(rc != CMP_ERROR);
1382 }
1383
1384 switch(n1->type)
1385 { case V_INTEGER:
1386 return n1->value.i < n2->value.i ? CMP_LESS :
1387 n1->value.i == n2->value.i ? CMP_EQUAL : CMP_GREATER;
1388 #ifdef O_GMP
1389 case V_MPZ:
1390 { int rc = mpz_cmp(n1->value.mpz, n2->value.mpz);
1391
1392 return rc < 0 ? CMP_LESS : rc == 0 ? CMP_EQUAL : CMP_GREATER;
1393 }
1394 case V_MPQ:
1395 { int rc = mpq_cmp(n1->value.mpq, n2->value.mpq);
1396
1397 return rc < 0 ? CMP_LESS : rc == 0 ? CMP_EQUAL : CMP_GREATER;
1398 }
1399 #endif
1400 case V_FLOAT:
1401 return n1->value.f < n2->value.f ? CMP_LESS :
1402 n1->value.f == n2->value.f ? CMP_EQUAL :
1403 n1->value.f > n2->value.f ? CMP_GREATER : CMP_NOTEQ;
1404 }
1405
1406 assert(0);
1407 return CMP_EQUAL;
1408 }
1409
1410 void
cpNumber(Number to,Number from)1411 cpNumber(Number to, Number from)
1412 { to->type = from->type;
1413
1414 switch(from->type)
1415 { case V_INTEGER:
1416 to->value.i = from->value.i;
1417 break;
1418 #ifdef O_GMP
1419 case V_MPZ:
1420 mpz_init(to->value.mpz);
1421 mpz_set(to->value.mpz, from->value.mpz);
1422 break;
1423 case V_MPQ:
1424 mpq_init(to->value.mpq);
1425 mpq_set(to->value.mpq, from->value.mpq);
1426 break;
1427 #endif
1428 case V_FLOAT:
1429 to->value.f = from->value.f;
1430 }
1431 }
1432
1433 #ifdef O_GMP
1434
1435 /*******************************
1436 * FLOAT <-> RATIONAL *
1437 *******************************/
1438
1439 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1440 This code is copied from ECLiPSe 7.0_53. This code is covered by the
1441 CMPL 1.1 (Cisco-style Mozilla Public License Version 1.1), available at
1442 www.eclipse-clp.org/license.
1443 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1444
1445 /*
1446 * Divide two bignums giving a double float result. The naive solution
1447 * return mpz_to_double(num) / mpz_to_double(den);
1448 * suffers from floating point overflows when the numbers are huge
1449 * and is inefficient because it looks at unnecessarily many digits.
1450 *
1451 * IEEE double precision is 53 bits mantissa and 12 bits signed exponent.
1452 * So the largest integer representable with doubles is 1024 bits wide,
1453 * of which the first 53 are ones, i.e. it lies between 2^1023 and 2^1024.
1454 * If the dividend's MSB is more than 1024 bits higher than the divisor's,
1455 * the result will always be floating point infinity (no need to divide).
1456 * If we do divide, we first drop excess integer precision by keeping only
1457 * DBL_PRECISION_LIMBS and ignoring the lower limbs for both operands
1458 * (i.e. we effectively scale the integers down, or right-shift them).
1459 */
1460
1461 #define MIN_LIMB_DIFF (2+1024/GMP_NUMB_BITS)
1462 #define DBL_PRECISION_LIMBS (2+53/GMP_NUMB_BITS)
1463 #define MAX_ULONG_DBL ((double)~0L+1.0)
1464
1465 static double
mpz_fdiv(mpz_t num,mpz_t den)1466 mpz_fdiv(mpz_t num, mpz_t den)
1467 { mp_ptr longer_d, shorter_d;
1468 mp_size_t shorter_size, longer_size, ignored_limbs = 0;
1469 int negative, swapped;
1470 /* By declaring res volatile we make sure that the result is rounded
1471 * to double precision instead of being returned with extended precision
1472 * in a floating point register, which can have confusing consequences */
1473 volatile double res;
1474
1475 shorter_size = num->_mp_size;
1476 longer_size = den->_mp_size;
1477 negative = 0;
1478
1479 if ( shorter_size < 0 )
1480 { shorter_size = -shorter_size;
1481 negative = !negative;
1482 }
1483 if ( longer_size < 0 )
1484 { longer_size = -longer_size;
1485 negative = !negative;
1486 }
1487 if ( shorter_size > longer_size )
1488 { longer_size = shorter_size;
1489 longer_d = num->_mp_d;
1490 shorter_size = (den->_mp_size >= 0 ? den->_mp_size : -den->_mp_size);
1491 shorter_d = den->_mp_d;
1492 swapped = 1; /* abs(res) > 1 */
1493 } else
1494 { longer_d = den->_mp_d;
1495 shorter_d = num->_mp_d;
1496 swapped = 0; /* abs(res) < 1 */
1497 }
1498
1499 if ( longer_size - shorter_size > MIN_LIMB_DIFF )
1500 { res = swapped ? HUGE_VAL : 0.0;
1501 } else
1502 { double l,s;
1503 long int le, se;
1504 mpz_t li, si;
1505 int r_mode;
1506
1507 /* we ignore limbs that are not significant for the result */
1508 if ( longer_size > MIN_LIMB_DIFF ) /* more can't be represented */
1509 { ignored_limbs = longer_size - MIN_LIMB_DIFF;
1510 longer_size -= ignored_limbs;
1511 shorter_size -= ignored_limbs;
1512 }
1513 if ( shorter_size > DBL_PRECISION_LIMBS ) /* more exceeds the precision */
1514 { ignored_limbs += shorter_size - DBL_PRECISION_LIMBS;
1515 longer_size -= shorter_size - DBL_PRECISION_LIMBS;
1516 shorter_size = DBL_PRECISION_LIMBS;
1517 }
1518 longer_d += ignored_limbs;
1519 shorter_d += ignored_limbs;
1520 li->_mp_alloc = li->_mp_size = longer_size; li->_mp_d = longer_d;
1521 si->_mp_alloc = si->_mp_size = shorter_size; si->_mp_d = shorter_d;
1522
1523 l = mpz_get_d_2exp(&le, li);
1524 s = mpz_get_d_2exp(&se, si);
1525
1526 /* if result negative, some rounding modes must be swapped;
1527 avoid if unnecessary */
1528 if ( negative )
1529 { r_mode = fegetround();
1530
1531 switch (r_mode)
1532 { case FE_UPWARD: fesetround(FE_DOWNWARD); break;
1533 case FE_DOWNWARD: fesetround(FE_UPWARD);
1534 }
1535 }
1536
1537 if ( swapped )
1538 res = (l/s) * pow(2.0, le-se);
1539 else
1540 res = (s/l) * pow(2.0, se-le);
1541
1542 if ( negative )
1543 fesetround(r_mode);
1544 }
1545
1546 return negative ? -res : res;
1547 }
1548
1549 double
mpq_to_double(mpq_t q)1550 mpq_to_double(mpq_t q)
1551 { return mpz_fdiv(mpq_numref(q), mpq_denref(q));
1552 }
1553
1554
1555 /*
1556 * Try to compute a "nice" rational from a float, using continued fractions.
1557 * Stop when the rational converts back into the original float exactly.
1558 * If the process doesn't converge (due to numeric problems), or produces
1559 * a result longer than the fast bitwise conversion from the float mantissa,
1560 * then fall back to the bitwise conversion.
1561 */
1562
1563 void
mpq_set_double(mpq_t q,double f)1564 mpq_set_double(mpq_t q, double f) /* float -> nice rational */
1565 { double fabs = (f < 0.0) ? -f : f; /* get rid of the sign */
1566 double x = fabs;
1567 mpq_t b, c;
1568 MP_INT *pna = mpq_numref(q); /* use output q for a directly */
1569 MP_INT *pda = mpq_denref(q);
1570 MP_INT *pnb = mpq_numref(b);
1571 MP_INT *pdb = mpq_denref(b);
1572 MP_INT *pnc = mpq_numref(c);
1573 MP_INT *pdc = mpq_denref(c);
1574 mpz_t big_xi;
1575 int half_exp; /* predict bitwise conversion size */
1576 double fr = frexp(fabs, &half_exp);
1577 int bitwise_denominator_size = DBL_MANT_DIG-(half_exp-1);
1578
1579 (void)fr;
1580 if ( bitwise_denominator_size < 1 )
1581 bitwise_denominator_size = 1;
1582
1583 mpz_set_ui(pna, 1L); /* a = q = 1/0 */
1584 mpz_set_ui(pda, 0L);
1585 mpq_init(b); /* b = 0/1 */
1586 mpq_init(c); /* auxiliary */
1587 mpz_init(big_xi); /* auxiliary */
1588
1589 while ((mpz_fdiv(pna, pda)) != fabs)
1590 { /* infinite x indicates failure to converge */
1591 if ( !isfinite(x) )
1592 goto _bitwise_conversion_;
1593
1594 double xi = floor(x);
1595 double xf = x - xi;
1596
1597 /* compute a = a*xi + b for both numerator and denominator */
1598 mpq_swap(q, b);
1599 if ( x < (double)LONG_MAX )
1600 { unsigned long int_xi = (unsigned long) xi;
1601 mpz_mul_ui(pnc, pnb, int_xi);
1602 mpz_mul_ui(pdc, pdb, int_xi);
1603 } else
1604 { mpz_set_d(big_xi, xi);
1605 mpz_mul(pnc, pnb, big_xi);
1606 mpz_mul(pdc, pdb, big_xi);
1607 }
1608 mpz_add(pna, pna, pnc);
1609 mpz_add(pda, pda, pdc);
1610
1611 /* if it gets too long, fall back to bitwise conversion */
1612 if (mpz_sizeinbase(pda, 2) > bitwise_denominator_size)
1613 goto _bitwise_conversion_;
1614
1615 x = 1.0/xf;
1616 }
1617
1618 if ( f < 0.0 )
1619 mpq_neg(q, q);
1620 mpq_canonicalize(q); /* normally not be necessary */
1621
1622 goto _cleanup_;
1623
1624 _bitwise_conversion_:
1625 mpq_set_d(q, f); /* bitwise conversion */
1626
1627 _cleanup_:
1628 mpz_clear(big_xi);
1629 mpq_clear(c);
1630 mpq_clear(b);
1631 }
1632
1633
1634 /*******************************
1635 * PUBLIC INTERFACE *
1636 *******************************/
1637
1638 int
PL_get_mpz(term_t t,mpz_t mpz)1639 PL_get_mpz(term_t t, mpz_t mpz)
1640 { GET_LD
1641 Word p = valTermRef(t);
1642
1643 deRef(p);
1644 if ( isInteger(*p) )
1645 { number n;
1646
1647 get_integer(*p, &n);
1648 switch(n.type)
1649 { case V_INTEGER:
1650 promoteToMPZNumber(&n);
1651 mpz_set(mpz, n.value.mpz);
1652 clearNumber(&n);
1653 break;
1654 case V_MPZ:
1655 mpz_set(mpz, n.value.mpz);
1656 break;
1657 default:
1658 assert(0);
1659 }
1660
1661 return TRUE;
1662 }
1663
1664 return FALSE;
1665 }
1666
1667
1668 int
PL_get_mpq(term_t t,mpq_t mpq)1669 PL_get_mpq(term_t t, mpq_t mpq)
1670 { if ( PL_is_rational(t) )
1671 { GET_LD
1672 number n;
1673
1674 if ( valueExpression(t, &n PASS_LD) )
1675 { switch(n.type)
1676 { case V_INTEGER:
1677 if ( n.value.i >= LONG_MIN && n.value.i <= LONG_MAX )
1678 { mpq_set_si(mpq, (long)n.value.i, 1L);
1679 return TRUE;
1680 }
1681 promoteToMPZNumber(&n);
1682 /*FALLTHROUGH*/
1683 case V_MPZ:
1684 mpq_set_z(mpq, n.value.mpz);
1685 clearNumber(&n);
1686 return TRUE;
1687 case V_MPQ:
1688 mpq_set(mpq, n.value.mpq);
1689 clearNumber(&n);
1690 return TRUE;
1691 default:
1692 ;
1693 }
1694 clearNumber(&n);
1695 }
1696 }
1697
1698 return FALSE;
1699 }
1700
1701
1702 int
PL_unify_mpz(term_t t,mpz_t mpz)1703 PL_unify_mpz(term_t t, mpz_t mpz)
1704 { GET_LD
1705 number n;
1706 int rc;
1707
1708 n.type = V_MPZ;
1709 mpz_init(n.value.mpz);
1710 mpz_set(n.value.mpz, mpz);
1711
1712 rc = PL_unify_number(t, &n);
1713 clearNumber(&n);
1714
1715 return rc;
1716 }
1717
1718
1719 int
PL_unify_mpq(term_t t,mpq_t mpq)1720 PL_unify_mpq(term_t t, mpq_t mpq)
1721 { GET_LD
1722 number n;
1723 int rc;
1724
1725 n.type = V_MPQ;
1726 mpq_init(n.value.mpq);
1727 mpq_set(n.value.mpq, mpq);
1728
1729 rc = PL_unify_number(t, &n);
1730 clearNumber(&n);
1731
1732 return rc;
1733 }
1734
1735 /*******************************
1736 * WIN64 *
1737 *******************************/
1738
1739
1740 #if defined(WIN64) && _MSC_VER <= 1400 && !defined(_CRT_ASSEMBLY_VERSION)
1741 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1742 "#if (_MSC_VER <= 1400)" should suffice, but although both the VS2005
1743 (VC8) and the Microsoft Server 2003 R2 (VC8 SDK) define _MSC_VER=1400,
1744 VC8 SDK does not define the below functions, while VC8 does... The macro
1745 below distinguishes the two.
1746 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1747
1748 size_t
strnlen(const char * s,size_t maxlen)1749 strnlen(const char *s, size_t maxlen)
1750 { size_t len = 0;
1751
1752 while(*s++ && maxlen-- > 0)
1753 len++;
1754
1755 return len;
1756 }
1757
1758 void
__GSHandlerCheck()1759 __GSHandlerCheck()
1760 {
1761 }
1762 #endif
1763
1764 #endif /*O_GMP*/
1765