1 /****************************************************************************
2 **
3 ** This file is part of GAP, a system for computational discrete algebra.
4 **
5 ** Copyright of GAP belongs to its developers, whose names are too numerous
6 ** to list here. Please refer to the COPYRIGHT file for details.
7 **
8 ** SPDX-License-Identifier: GPL-2.0-or-later
9 **
10 ** This file implements the functions handling GMP integers.
11 **
12 ** There are three integer types in GAP: 'T_INT', 'T_INTPOS' and 'T_INTNEG'.
13 ** Each integer has a unique representation, e.g., an integer that can be
14 ** represented as 'T_INT' is never represented as 'T_INTPOS' or 'T_INTNEG'.
15 **
16 ** In the following, let 'N' be the number of bits in an UInt (so 32 or
17 ** 64, depending on the system). 'T_INT' is the type of those integers small
18 ** enough to fit into N-3 bits. Therefore the value range of this small
19 ** integers is $-2^{N-4}...2^{N-4}-1$. Only these small integers can be used
20 ** as index expression into sequences.
21 **
22 ** Small integers are represented by an immediate integer handle, containing
23 ** the value instead of pointing to it, which has the following form:
24 **
25 ** +-------+-------+-------+-------+- - - -+-------+-------+-------+
26 ** | guard | sign | bit | bit | | bit | tag | tag |
27 ** | bit | bit | N-5 | N-6 | | 0 | = 0 | = 1 |
28 ** +-------+-------+-------+-------+- - - -+-------+-------+-------+
29 **
30 ** Immediate integers handles carry the tag 'T_INT', i.e. the last bit is 1.
31 ** This distinguishes immediate integers from other handles which point to
32 ** structures aligned on even boundaries and therefore have last bit zero.
33 ** (The second bit is reserved as tag to allow extensions of this scheme.)
34 ** Using immediates as pointers and dereferencing them gives address errors.
35 **
36 ** To aid overflow check the most significant two bits must always be equal,
37 ** that is to say that the sign bit of immediate integers has a guard bit.
38 **
39 ** The macros 'INTOBJ_INT' and 'INT_INTOBJ' should be used to convert between
40 ** a small integer value and its representation as immediate integer handle.
41 **
42 ** 'T_INTPOS' and 'T_INTNEG' are the types of positive respectively negative
43 ** integer values that can not be represented by immediate integers.
44 **
45 ** These large integers values are represented as low-level GMP integer
46 ** objects, that is, in base 2^N. That means that the bag of a large integer
47 ** has the following form, where the "digits" are also more commonly referred
48 ** to as "limbs".:
49 **
50 ** +-------+-------+-------+-------+- - - -+-------+-------+-------+
51 ** | digit | digit | digit | digit | | digit | digit | digit |
52 ** | 0 | 1 | 2 | 3 | | <n>-2 | <n>-1 | <n> |
53 ** +-------+-------+-------+-------+- - - -+-------+-------+-------+
54 **
55 ** The value of this is: $d0 + d1 2^N + d2 (2^N)^2 + ... + d_n (2^N)^n$,
56 ** respectively the negative of this if the object type is 'T_INTNEG'.
57 **
58 ** Each digit resp. limb is stored as a N bit wide unsigned integer.
59 **
60 ** Note that we require that all large integers be normalized (that is, they
61 ** must not contain leading zero limbs) and reduced (they do not fit into a
62 ** small integer). Internally, a large integer may temporarily not be
63 ** normalized or not be reduced, but all kernel functions must make sure
64 ** that they eventually return normalized and reduced values. The function
65 ** GMP_NORMALIZE and GMP_REDUCE can be used to ensure this.
66 */
67
68 #include "integer.h"
69
70 #include "ariths.h"
71 #include "bool.h"
72 #include "calls.h"
73 #include "error.h"
74 #include "intfuncs.h"
75 #include "io.h"
76 #include "modules.h"
77 #include "opers.h"
78 #include "saveload.h"
79 #include "stats.h"
80 #include "stringobj.h"
81
82
83 /* TODO: Remove after Ward2 */
84 #ifndef WARD_ENABLED
85
86 #include <gmp.h>
87
88 #if GMP_NAIL_BITS != 0
89 #error Aborting compile: GAP does not support non-zero GMP nail size
90 #endif
91 #if !defined(__GNU_MP_RELEASE)
92 #if __GMP_MP_RELEASE < 50002
93 #error Aborting compile: GAP requires GMP 5.0.2 or newer
94 #endif
95 #endif
96
97 GAP_STATIC_ASSERT(GMP_LIMB_BITS == 8 * sizeof(UInt),
98 "GMP_LIMB_BITS != 8 * sizeof(UInt)");
99 GAP_STATIC_ASSERT(sizeof(mp_limb_t) == sizeof(UInt),
100 "sizeof(mp_limb_t) != sizeof(UInt)");
101
102 static Obj ObjInt_UIntInv( UInt i );
103
104
105 /* debugging */
106 #ifndef DEBUG_GMP
107 #define DEBUG_GMP 0
108 #endif
109
110 #if DEBUG_GMP
111
112 #if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901)
113 # define CURRENT_FUNCTION __func__
114 #elif defined(_MSC_VER)
115 # define CURRENT_FUNCTION __FUNCTION__
116 #else
117 # define CURRENT_FUNCTION "<unknown>"
118 #endif
119
120 #define CHECK_INT(op) IS_NORMALIZED_AND_REDUCED(op, CURRENT_FUNCTION, __LINE__)
121 #else
122 #define CHECK_INT(op) do { } while(0);
123 #endif
124
125
126 /* macros to save typing later :) */
127 #define VAL_LIMB0(obj) (*CONST_ADDR_INT(obj))
128 #define SET_VAL_LIMB0(obj,val) do { *ADDR_INT(obj) = val; } while(0)
129 #define IS_INTPOS(obj) (TNUM_OBJ(obj) == T_INTPOS)
130 #define IS_INTNEG(obj) (TNUM_OBJ(obj) == T_INTNEG)
131
132 #define SIZE_INT_OR_INTOBJ(obj) (IS_INTOBJ(obj) ? 1 : SIZE_INT(obj))
133
134 #define RequireNonzero(funcname, op, argname) \
135 do { \
136 if (op == INTOBJ_INT(0) ) { \
137 ErrorMayQuit(funcname ": <" argname "> must be nonzero", 0, 0); \
138 } \
139 } while (0)
140
141
142 GAP_STATIC_ASSERT( sizeof(mp_limb_t) == sizeof(UInt), "gmp limb size incompatible with GAP word size");
143
144
145 /* This ensures that all memory underlying a bag is actually committed
146 ** to physical memory and can be written to.
147 ** This is a workaround to a bug specific to Cygwin 64-bit and bad
148 ** interaction with GMP, so this is only needed specifically for new
149 ** bags created in this module to hold the outputs of GMP routines.
150 **
151 ** Thus, any time NewBag is called, it is also necessary to call
152 ** ENSURE_BAG(bag) on the newly created bag if some GMP function will be
153 ** the first place that bag's data is written to.
154 **
155 ** To give a counter-example, ENSURE_BAG is *not* needed in ObjInt_Int,
156 ** because it just creates a bag to hold a single mp_limb_t, and
157 ** immediately assigns it a value.
158 **
159 ** The bug this works around is explained more in
160 ** https://github.com/gap-system/gap/issues/3434
161 */
ENSURE_BAG(Bag bag)162 static inline void ENSURE_BAG(Bag bag)
163 {
164 // Note: This workaround is only required with the original GMP and not with
165 // MPIR
166 #if defined(SYS_IS_CYGWIN32) && defined(SYS_IS_64_BIT) && \
167 !defined(__MPIR_VERSION)
168 memset(PTR_BAG(bag), 0, SIZE_BAG(bag));
169 #endif
170 }
171
172
173 /* for fallbacks to library */
174 static Obj String;
175 static Obj OneAttr;
176 static Obj IsIntFilt;
177
178
179 /****************************************************************************
180 **
181 *F TypeInt(<val>) . . . . . . . . . . . . . . . . . . . . . type of integer
182 **
183 ** 'TypeInt' returns the type of the integer <val>.
184 **
185 ** 'TypeInt' is the function in 'TypeObjFuncs' for integers.
186 */
187 static Obj TYPE_INT_SMALL_ZERO;
188 static Obj TYPE_INT_SMALL_POS;
189 static Obj TYPE_INT_SMALL_NEG;
190 static Obj TYPE_INT_LARGE_POS;
191 static Obj TYPE_INT_LARGE_NEG;
192
TypeIntSmall(Obj val)193 static Obj TypeIntSmall(Obj val)
194 {
195 if ( 0 == INT_INTOBJ(val) ) {
196 return TYPE_INT_SMALL_ZERO;
197 }
198 else if ( 0 < INT_INTOBJ(val) ) {
199 return TYPE_INT_SMALL_POS;
200 }
201 else {
202 return TYPE_INT_SMALL_NEG;
203 }
204 }
205
TypeIntLargePos(Obj val)206 static Obj TypeIntLargePos(Obj val)
207 {
208 return TYPE_INT_LARGE_POS;
209 }
210
TypeIntLargeNeg(Obj val)211 static Obj TypeIntLargeNeg(Obj val)
212 {
213 return TYPE_INT_LARGE_NEG;
214 }
215
216
217 /****************************************************************************
218 **
219 *F FiltIS_INT( <self>, <val> ) . . . . . . . . . . internal function 'IsInt'
220 **
221 ** 'FiltIS_INT' implements the internal filter 'IsInt'.
222 **
223 ** 'IsInt( <val> )'
224 **
225 ** 'IsInt' returns 'true' if the value <val> is a small integer or a
226 ** large int, and 'false' otherwise.
227 */
FiltIS_INT(Obj self,Obj val)228 static Obj FiltIS_INT(Obj self, Obj val)
229 {
230 if ( IS_INT(val) ) {
231 return True;
232 }
233 else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
234 return False;
235 }
236 else {
237 return DoFilter( self, val );
238 }
239 }
240
241
242 /****************************************************************************
243 **
244 *F SaveInt( <op> )
245 **
246 **
247 */
SaveInt(Obj op)248 static void SaveInt(Obj op)
249 {
250 const UInt * ptr = CONST_ADDR_INT(op);
251 for (UInt i = 0; i < SIZE_INT(op); i++)
252 SaveUInt(*ptr++);
253 return;
254 }
255
256
257 /****************************************************************************
258 **
259 *F LoadInt( <op> )
260 **
261 **
262 */
LoadInt(Obj op)263 static void LoadInt(Obj op)
264 {
265 UInt * ptr = ADDR_INT(op);
266 for (UInt i = 0; i < SIZE_INT(op); i++)
267 *ptr++ = LoadUInt();
268 return;
269 }
270
271
272 /****************************************************************************
273 **
274 ** In order to use the high-level GMP mpz_* functions conveniently while
275 ** retaining a low overhead, we introduce fake_mpz_t, which can be thought
276 ** of as a "subclass" of mpz_t, extending it by temporary storage for
277 ** single limb of data, as well as a reference to a corresponding GAP
278 ** object.
279 **
280 ** For an example on how to correctly use this, see GcdInt.
281 */
282 typedef struct {
283 mpz_t v;
284 mp_limb_t tmp;
285 Obj obj;
286 } fake_mpz_t[1];
287
288
289 /****************************************************************************
290 **
291 *F NEW_FAKEMPZ( <fake>, <size> )
292 **
293 ** Setup a fake mpz_t object for capturing the output of a GMP mpz_ function,
294 ** with space for up to <size> limbs allocated.
295 */
NEW_FAKEMPZ(fake_mpz_t fake,UInt size)296 static void NEW_FAKEMPZ( fake_mpz_t fake, UInt size )
297 {
298 fake->v->_mp_alloc = size;
299 fake->v->_mp_size = 0;
300 if (size == 1) {
301 fake->obj = 0;
302 }
303 else {
304 fake->obj = NewBag( T_INTPOS, size * sizeof(mp_limb_t) );
305 ENSURE_BAG(fake->obj);
306 }
307 }
308
309
310 /****************************************************************************
311 **
312 *F FAKEMPZ_GMPorINTOBJ( <fake>, <op> )
313 **
314 ** Initialize <fake> to reference the content of <op>. For this, <op>
315 ** must be an integer, either small or large, but this is *not* checked.
316 ** The calling code is responsible for any verification.
317 */
FAKEMPZ_GMPorINTOBJ(fake_mpz_t fake,Obj op)318 static void FAKEMPZ_GMPorINTOBJ( fake_mpz_t fake, Obj op )
319 {
320 if (IS_INTOBJ(op)) {
321 fake->obj = 0;
322 fake->v->_mp_alloc = 1;
323 const Int i = INT_INTOBJ(op);
324 if ( i >= 0 ) {
325 fake->tmp = i;
326 fake->v->_mp_size = i ? 1 : 0;
327 }
328 else {
329 fake->tmp = -i;
330 fake->v->_mp_size = -1;
331 }
332 }
333 else {
334 fake->obj = op;
335 fake->v->_mp_alloc = SIZE_INT(op);
336 fake->v->_mp_size = IS_INTPOS(op) ? SIZE_INT(op) : -SIZE_INT(op);
337 }
338 }
339
340
341 /****************************************************************************
342 **
343 *F GMPorINTOBJ_FAKEMPZ( <fake> )
344 **
345 ** This function converts a fake mpz_t into a GAP integer object.
346 */
GMPorINTOBJ_FAKEMPZ(fake_mpz_t fake)347 static Obj GMPorINTOBJ_FAKEMPZ( fake_mpz_t fake )
348 {
349 Obj obj = fake->obj;
350 if ( fake->v->_mp_size == 0 ) {
351 obj = INTOBJ_INT(0);
352 }
353 else if ( obj != 0 ) {
354 if ( fake->v->_mp_size < 0 ) {
355 /* Warning: changing the bag type is only correct if the object was
356 not yet visible to the outside world. Thus, it is safe to use
357 with an fake_mpz_t initialized with NEW_FAKEMPZ, but not with
358 one that was setup by FAKEMPZ_GMPorINTOBJ. */
359 RetypeBag( obj, T_INTNEG );
360 }
361 obj = GMP_NORMALIZE( obj );
362 obj = GMP_REDUCE( obj );
363 }
364 else {
365 if ( fake->v->_mp_size == 1 )
366 obj = ObjInt_UInt( fake->tmp );
367 else
368 obj = ObjInt_UIntInv( fake->tmp );
369 }
370 return obj;
371 }
372
373 /****************************************************************************
374 **
375 *F GMPorINTOBJ_MPZ( <fake> )
376 **
377 ** This function converts an mpz_t into a GAP integer object.
378 */
GMPorINTOBJ_MPZ(mpz_t v)379 static Obj GMPorINTOBJ_MPZ( mpz_t v )
380 {
381 return MakeObjInt((const UInt *)v->_mp_d, v->_mp_size);
382 }
383
384
385 /****************************************************************************
386 **
387 ** MPZ_FAKEMPZ( <fake> )
388 **
389 ** This converts a fake_mpz_t into an mpz_t. As a side effect, it updates
390 ** fake->v->_mp_d. This allows us to use SWAP on fake_mpz_t objects, and
391 ** also protects against garbage collection moving around data.
392 */
393 #define MPZ_FAKEMPZ(fake) (UPDATE_FAKEMPZ(fake), fake->v)
394
395 /* UPDATE_FAKEMPZ is a helper function for the MPZ_FAKEMPZ macro */
UPDATE_FAKEMPZ(fake_mpz_t fake)396 static inline void UPDATE_FAKEMPZ( fake_mpz_t fake )
397 {
398 fake->v->_mp_d = fake->obj ? (mp_ptr)ADDR_INT(fake->obj) : &fake->tmp;
399 }
400
401 /* some extra debugging tools for FAKMPZ objects */
402 #if DEBUG_GMP
403 #define CHECK_FAKEMPZ(fake) \
404 assert( ((fake)->v->_mp_d == ((fake)->obj ? (mp_ptr)ADDR_INT((fake)->obj) : &(fake)->tmp )) \
405 && (fake->v->_mp_alloc == ((fake)->obj ? SIZE_INT((fake)->obj) : 1 )) )
406 #else
407 #define CHECK_FAKEMPZ(fake) do { } while(0);
408 #endif
409
410
411 /****************************************************************************
412 **
413 *F GMP_NORMALIZE( <op> ) . . . . . . . remove leading zeros from a GMP bag
414 **
415 ** 'GMP_NORMALIZE' removes any leading zeros from a large integer object
416 ** and returns a small int or resizes the bag if possible.
417 **
418 */
GMP_NORMALIZE(Obj op)419 Obj GMP_NORMALIZE(Obj op)
420 {
421 mp_size_t size;
422 if (IS_INTOBJ(op)) {
423 return op;
424 }
425 for (size = SIZE_INT(op); size != (mp_size_t)1; size--) {
426 if (CONST_ADDR_INT(op)[(size - 1)] != 0) {
427 break;
428 }
429 }
430 if (size < SIZE_INT(op)) {
431 ResizeBag(op, size * sizeof(mp_limb_t));
432 }
433 return op;
434 }
435
GMP_REDUCE(Obj op)436 Obj GMP_REDUCE(Obj op)
437 {
438 if (IS_INTOBJ(op)) {
439 return op;
440 }
441 if (SIZE_INT(op) == 1) {
442 if ((VAL_LIMB0(op) <= INT_INTOBJ_MAX) ||
443 (IS_INTNEG(op) && VAL_LIMB0(op) == -INT_INTOBJ_MIN)) {
444 if (IS_INTNEG(op)) {
445 return INTOBJ_INT(-(Int)VAL_LIMB0(op));
446 }
447 else {
448 return INTOBJ_INT((Int)VAL_LIMB0(op));
449 }
450 }
451 }
452 return op;
453 }
454
455 /****************************************************************************
456 **
457 ** This is a helper function for the CHECK_INT macro, which checks that
458 ** the given integer object <op> is normalized and reduced.
459 **
460 */
461 #if DEBUG_GMP
IS_NORMALIZED_AND_REDUCED(Obj op,const char * func,int line)462 static int IS_NORMALIZED_AND_REDUCED(Obj op, const char * func, int line)
463 {
464 mp_size_t size;
465 if ( IS_INTOBJ( op ) ) {
466 return 1;
467 }
468 if ( !IS_LARGEINT( op ) ) {
469 /* ignore non-integers */
470 return 0;
471 }
472 for ( size = SIZE_INT(op); size != (mp_size_t)1; size-- ) {
473 if ( CONST_ADDR_INT(op)[(size - 1)] != 0 ) {
474 break;
475 }
476 }
477 if ( size < SIZE_INT(op) ) {
478 Pr("WARNING: non-normalized gmp value (%s:%d)\n",(Int)func,line);
479 }
480 if ( SIZE_INT(op) == 1) {
481 if ( ( VAL_LIMB0(op) <= INT_INTOBJ_MAX ) ||
482 ( IS_INTNEG(op) && VAL_LIMB0(op) == -INT_INTOBJ_MIN ) ) {
483 if ( IS_INTNEG(op) ) {
484 Pr("WARNING: non-reduced negative gmp value (%s:%d)\n",(Int)func,line);
485 return 0;
486 }
487 else {
488 Pr("WARNING: non-reduced positive gmp value (%s:%d)\n",(Int)func,line);
489 return 0;
490 }
491 }
492 }
493 return 1;
494 }
495 #endif
496
497
498 /****************************************************************************
499 **
500 *F ObjInt_Int( <cint> ) . . . . . . . . . . convert c int to integer object
501 **
502 ** 'ObjInt_Int' takes the C integer <cint> and returns the equivalent
503 ** GMP obj or int obj, according to the value of <cint>.
504 **
505 */
ObjInt_Int(Int i)506 Obj ObjInt_Int( Int i )
507 {
508 Obj gmp;
509
510 if (INT_INTOBJ_MIN <= i && i <= INT_INTOBJ_MAX) {
511 return INTOBJ_INT(i);
512 }
513 else if (i < 0 ) {
514 gmp = NewBag( T_INTNEG, sizeof(mp_limb_t) );
515 i = -i;
516 }
517 else {
518 gmp = NewBag( T_INTPOS, sizeof(mp_limb_t) );
519 }
520 SET_VAL_LIMB0( gmp, i );
521 return gmp;
522 }
523
ObjInt_UInt(UInt i)524 Obj ObjInt_UInt( UInt i )
525 {
526 Obj gmp;
527 if (i <= INT_INTOBJ_MAX) {
528 return INTOBJ_INT(i);
529 }
530 else {
531 gmp = NewBag( T_INTPOS, sizeof(mp_limb_t) );
532 SET_VAL_LIMB0( gmp, i );
533 return gmp;
534 }
535 }
536
537 // initialize to -i
ObjInt_UIntInv(UInt i)538 Obj ObjInt_UIntInv( UInt i )
539 {
540 Obj gmp;
541 // we need to test INT_INTOBJ_MIN <= -i; to express this with unsigned
542 // values, we must avoid all negative terms, which leads to this equivalent
543 // check:
544 if (i <= -INT_INTOBJ_MIN) {
545 return INTOBJ_INT(-i);
546 }
547 else {
548 gmp = NewBag( T_INTNEG, sizeof(mp_limb_t) );
549 SET_VAL_LIMB0( gmp, i );
550 return gmp;
551 }
552 }
553
554
ObjInt_Int8(Int8 i)555 Obj ObjInt_Int8( Int8 i )
556 {
557 #ifdef SYS_IS_64_BIT
558 return ObjInt_Int(i);
559 #else
560 if (i == (Int4)i) {
561 return ObjInt_Int((Int4)i);
562 }
563
564 /* we need two limbs to store this integer */
565 assert( sizeof(mp_limb_t) == 4 );
566 Obj gmp;
567 if (i >= 0) {
568 gmp = NewBag( T_INTPOS, 2 * sizeof(mp_limb_t) );
569 } else {
570 gmp = NewBag( T_INTNEG, 2 * sizeof(mp_limb_t) );
571 i = -i;
572 }
573
574 UInt *ptr = ADDR_INT(gmp);
575 ptr[0] = (UInt4)i;
576 ptr[1] = ((UInt8)i) >> 32;
577 return gmp;
578 #endif
579 }
580
ObjInt_UInt8(UInt8 i)581 Obj ObjInt_UInt8( UInt8 i )
582 {
583 #ifdef SYS_IS_64_BIT
584 return ObjInt_UInt(i);
585 #else
586 if (i == (UInt4)i) {
587 return ObjInt_UInt((UInt4)i);
588 }
589
590 /* we need two limbs to store this integer */
591 assert( sizeof(mp_limb_t) == 4 );
592 Obj gmp = NewBag( T_INTPOS, 2 * sizeof(mp_limb_t) );
593 UInt *ptr = ADDR_INT(gmp);
594 ptr[0] = (UInt4)i;
595 ptr[1] = ((UInt8)i) >> 32;
596 return gmp;
597 #endif
598 }
599
600 /****************************************************************************
601 **
602 ** Convert GAP Integers to various C types -- see header file
603 */
Int_ObjInt(Obj i)604 Int Int_ObjInt(Obj i)
605 {
606 UInt sign = 0;
607 if (IS_INTOBJ(i))
608 return INT_INTOBJ(i);
609 // must be a single limb
610 if (TNUM_OBJ(i) == T_INTPOS)
611 sign = 0;
612 else if (TNUM_OBJ(i) == T_INTNEG)
613 sign = 1;
614 else
615 ErrorMayQuit("Conversion error, expecting an integer, not a %s",
616 (Int)TNAM_OBJ(i), 0);
617 if (SIZE_BAG(i) != sizeof(mp_limb_t))
618 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
619
620 // now check if val is small enough to fit in the signed Int type
621 // that has a range from -2^N to 2^N-1 so we need to check both ends
622 // Since -2^N is the same bit pattern as the UInt 2^N (N is 31 or 63)
623 // we can do it as below which avoids some compiler warnings
624 UInt val = VAL_LIMB0(i);
625 #ifdef SYS_IS_64_BIT
626 if ((!sign && (val > INT64_MAX)) || (sign && (val > (UInt)INT64_MIN)))
627 #else
628 if ((!sign && (val > INT32_MAX)) || (sign && (val > (UInt)INT32_MIN)))
629 #endif
630 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
631 return sign ? -(Int)val : (Int)val;
632 }
633
UInt_ObjInt(Obj i)634 UInt UInt_ObjInt(Obj i)
635 {
636 if (IS_NEG_INT(i))
637 ErrorMayQuit("Conversion error, cannot convert negative integer to unsigned type", 0, 0);
638 if (IS_INTOBJ(i))
639 return (UInt)INT_INTOBJ(i);
640 if (TNUM_OBJ(i) != T_INTPOS)
641 ErrorMayQuit("Conversion error, expecting an integer, not a %s",
642 (Int)TNAM_OBJ(i), 0);
643
644 // must be a single limb
645 if (SIZE_INT(i) != 1)
646 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
647 return VAL_LIMB0(i);
648 }
649
Int8_ObjInt(Obj i)650 Int8 Int8_ObjInt(Obj i)
651 {
652 #ifdef SYS_IS_64_BIT
653 // in this case Int8 is Int
654 return Int_ObjInt(i);
655 #else
656 if (IS_INTOBJ(i))
657 return (Int8)INT_INTOBJ(i);
658
659 UInt sign = 0;
660 if (TNUM_OBJ(i) == T_INTPOS)
661 sign = 0;
662 else if (TNUM_OBJ(i) == T_INTNEG)
663 sign = 1;
664 else
665 ErrorMayQuit("Conversion error, expecting an integer, not a %s",
666 (Int)TNAM_OBJ(i), 0);
667
668 // must be at most two limbs
669 if (SIZE_INT(i) > 2)
670 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
671 UInt vall = VAL_LIMB0(i);
672 UInt valh = (SIZE_INT(i) == 1) ? 0 : CONST_ADDR_INT(i)[1];
673 UInt8 val = (UInt8)vall + ((UInt8)valh << 32);
674 // now check if val is small enough to fit in the signed Int8 type
675 // that has a range from -2^63 to 2^63-1 so we need to check both ends
676 // Since -2^63 is the same bit pattern as the UInt8 2^63 we can do it
677 // this way which avoids some compiler warnings
678 if ((!sign && (val > INT64_MAX)) || (sign && (val > (UInt8)INT64_MIN)))
679 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
680 return sign ? -(Int8)val : (Int8)val;
681 #endif
682 }
683
UInt8_ObjInt(Obj i)684 UInt8 UInt8_ObjInt(Obj i)
685 {
686 #ifdef SYS_IS_64_BIT
687 // in this case UInt8 is UInt
688 return UInt_ObjInt(i);
689 #else
690 if (IS_NEG_INT(i))
691 ErrorMayQuit("Conversion error, cannot convert negative integer to unsigned type", 0, 0);
692 if (IS_INTOBJ(i))
693 return (UInt8)INT_INTOBJ(i);
694 if (TNUM_OBJ(i) != T_INTPOS)
695 ErrorMayQuit("Conversion error, expecting an integer, not a %s",
696 (Int)TNAM_OBJ(i), 0);
697 if (SIZE_INT(i) > 2)
698 ErrorMayQuit("Conversion error, integer too large", 0L, 0L);
699 UInt vall = VAL_LIMB0(i);
700 UInt valh = (SIZE_INT(i) == 1) ? 0 : CONST_ADDR_INT(i)[1];
701 return (UInt8)vall + ((UInt8)valh << 32);
702 #endif
703 }
704
705
706 /****************************************************************************
707 **
708 *F MakeObjInt(<limbs>, <size>) . . . . . . . . . create a new integer object
709 **
710 */
MakeObjInt(const UInt * limbs,int size)711 Obj MakeObjInt(const UInt * limbs, int size)
712 {
713 Obj obj;
714 if (size == 0)
715 obj = INTOBJ_INT(0);
716 else if (size == 1)
717 obj = ObjInt_UInt(limbs[0]);
718 else if (size == -1)
719 obj = ObjInt_UIntInv(limbs[0]);
720 else {
721 UInt tnum = (size > 0 ? T_INTPOS : T_INTNEG);
722 if (size < 0) size = -size;
723 obj = NewBag(tnum, size * sizeof(mp_limb_t));
724 memcpy(ADDR_INT(obj), limbs, size * sizeof(mp_limb_t));
725
726 obj = GMP_NORMALIZE(obj);
727 obj = GMP_REDUCE(obj);
728 }
729 return obj;
730 }
731
732
733 // This function returns an immediate integer, or
734 // an integer object with exactly one limb, and returns
735 // its absolute value as an unsigned integer.
AbsOfSmallInt(Obj x)736 static inline UInt AbsOfSmallInt(Obj x)
737 {
738 if (!IS_INTOBJ(x)) {
739 GAP_ASSERT(SIZE_INT(x) == 1);
740 return VAL_LIMB0(x);
741 }
742 Int val = INT_INTOBJ(x);
743 return val > 0 ? val : -val;
744 }
745
746
747 /****************************************************************************
748 **
749 *F PrintInt( <op> ) . . . . . . . . . . . . . . . . print an integer object
750 **
751 ** 'PrintInt' prints the integer <op> in the usual decimal notation.
752 */
PrintInt(Obj op)753 void PrintInt ( Obj op )
754 {
755 /* print a small integer */
756 if ( IS_INTOBJ(op) ) {
757 Pr( "%>%d%<", INT_INTOBJ(op), 0L );
758 }
759
760 /* print a large integer */
761 else if ( SIZE_INT(op) < 1000 ) {
762 CHECK_INT(op);
763
764 /* Use GMP to print the integer to a buffer. We are looking at an
765 integer with less than 1000 limbs, hence in decimal notation, it
766 will take up at most LogInt( 2^(1000 * GMP_LIMB_BITS), 10)
767 digits. Since 1000*Log(2)/Log(10) = 301.03, we get the following
768 estimate for the buffer size (the overestimate is big enough to
769 include space for a sign and a null terminator). */
770 Char buf[302 * GMP_LIMB_BITS];
771 mpz_t v;
772 v->_mp_alloc = SIZE_INT(op);
773 v->_mp_size = IS_INTPOS(op) ? v->_mp_alloc : -v->_mp_alloc;
774 v->_mp_d = (mp_ptr)ADDR_INT(op);
775 mpz_get_str(buf, 10, v);
776
777 /* print the buffer, %> means insert '\' before a linebreak */
778 Pr("%>%s%<",(Int)buf, 0);
779 }
780 else {
781 Obj str = CALL_1ARGS( String, op );
782 Pr("%>", 0, 0);
783 PrintString1(str);
784 Pr("%<", 0, 0);
785 /* for a long time Print of large ints did not follow the general idea
786 * that Print should produce something that can be read back into GAP:
787 Pr("<<an integer too large to be printed>>",0L,0L); */
788 }
789 }
790
791
792 /****************************************************************************
793 **
794 *F StringIntBase( <op>, <base> )
795 **
796 ** Convert the integer <op> to a string relative to the given base <base>.
797 ** Here, base may range from 2 to 36.
798 */
StringIntBase(Obj op,int base)799 static Obj StringIntBase(Obj op, int base)
800 {
801 int len;
802 Obj res;
803 fake_mpz_t v;
804
805 GAP_ASSERT(IS_INT(op));
806 CHECK_INT(op);
807 GAP_ASSERT(2 <= base && base <= 36);
808
809 /* 0 is special */
810 if (op == INTOBJ_INT(0)) {
811 res = NEW_STRING(1);
812 CHARS_STRING(res)[0] = '0';
813 return res;
814 }
815
816 /* convert integer to fake_mpz_t */
817 FAKEMPZ_GMPorINTOBJ(v, op);
818
819 /* allocate the result string */
820 len = mpz_sizeinbase( MPZ_FAKEMPZ(v), base ) + 2;
821 res = NEW_STRING( len );
822
823 /* ask GMP to perform the actual conversion */
824 mpz_get_str( CSTR_STRING( res ), -base, MPZ_FAKEMPZ(v) );
825
826 /* we may have to shrink the string */
827 int real_len = strlen( CONST_CSTR_STRING(res) );
828 if ( real_len != GET_LEN_STRING(res) ) {
829 SET_LEN_STRING(res, real_len);
830 }
831
832 return res;
833 }
834
835 /****************************************************************************
836 **
837 *F FuncHexStringInt( <self>, <n> ) . . . . . . . . . hex string for GAP int
838 *F FuncIntHexString( <self>, <string> ) . . . . . . GAP int from hex string
839 **
840 ** The function `FuncHexStringInt' constructs from a GAP integer the
841 ** corresponding string in hexadecimal notation. It has a leading '-'
842 ** for negative numbers and the digits 10..15 are written as A..F.
843 **
844 ** The function `FuncIntHexString' does the converse, but here the
845 ** letters a..f are also allowed in <string> instead of A..F.
846 **
847 */
FuncHexStringInt(Obj self,Obj n)848 static Obj FuncHexStringInt(Obj self, Obj n)
849 {
850 RequireInt("HexStringInt", n);
851 return StringIntBase(n, 16);
852 }
853
854
855 /****************************************************************************
856 **
857 ** This helper function for IntHexString reads <len> bytes in the string
858 ** <p> and parses them as a hexadecimal. The resulting integer is returned.
859 ** This function does not check for overflow, so make sure that len*4 does
860 ** not exceed the number of bits in mp_limb_t.
861 **/
hexstr2int(const UInt1 * p,UInt len)862 static mp_limb_t hexstr2int( const UInt1 *p, UInt len )
863 {
864 mp_limb_t n = 0;
865 UInt1 a;
866 while (len--) {
867 a = *p++;
868 if (a >= 'a')
869 a -= 'a' - 10;
870 else if ( a>= 'A')
871 a -= 'A' - 10;
872 else
873 a -= '0';
874 if (a > 15)
875 ErrorMayQuit("IntHexString: invalid character in hex-string", 0L, 0L);
876 n = (n << 4) + a;
877 }
878 return n;
879 }
880
FuncIntHexString(Obj self,Obj str)881 static Obj FuncIntHexString(Obj self, Obj str)
882 {
883 return IntHexString(str);
884 }
885
IntHexString(Obj str)886 Obj IntHexString(Obj str)
887 {
888 Obj res;
889 Int i, len, sign, nd;
890 mp_limb_t n;
891 const UInt1 *p;
892 UInt *limbs;
893
894 RequireStringRep("IntHexString", str);
895
896 len = GET_LEN_STRING(str);
897 if (len == 0) {
898 res = INTOBJ_INT(0);
899 return res;
900 }
901 p = CONST_CHARS_STRING(str);
902 if (*p == '-') {
903 sign = -1;
904 i = 1;
905 }
906 else {
907 sign = 1;
908 i = 0;
909 }
910
911 while (p[i] == '0' && i < len)
912 i++;
913 len -= i;
914
915 if (len*4 <= NR_SMALL_INT_BITS) {
916 n = hexstr2int( p + i, len );
917 res = INTOBJ_INT(sign * n);
918 return res;
919 }
920
921 else {
922 /* Each hex digit corresponds to to 4 bits, and each GMP limb has sizeof(UInt)
923 bytes, thus 2*sizeof(UInt) hex digits fit into one limb. We use this
924 to compute the number of limbs minus 1: */
925 nd = (len - 1) / (2*sizeof(UInt));
926 res = NewBag( (sign == 1) ? T_INTPOS : T_INTNEG, (nd + 1) * sizeof(mp_limb_t) );
927
928 /* update pointer, in case a garbage collection happened */
929 p = CONST_CHARS_STRING(str) + i;
930 limbs = ADDR_INT(res);
931
932 /* if len is not divisible by 2*sizeof(UInt), then take care of the extra bytes */
933 UInt diff = len - nd * (2*sizeof(UInt));
934 if ( diff ) {
935 n = hexstr2int( p, diff );
936 p += diff;
937 len -= diff;
938 limbs[nd--] = n;
939 }
940
941 /* */
942 while ( len ) {
943 n = hexstr2int( p, 2*sizeof(UInt) );
944 p += 2*sizeof(UInt);
945 len -= 2*sizeof(UInt);
946 limbs[nd--] = n;
947 }
948
949 res = GMP_NORMALIZE(res);
950 res = GMP_REDUCE(res);
951 return res;
952 }
953 }
954
955
956 /****************************************************************************
957 **
958 ** Implementation of Log2Int for C integers.
959 **
960 ** When available, we try to use GCC builtins. Otherwise, fall back to code
961 ** based on
962 ** https://graphics.stanford.edu/~seander/bithacks.html#IntegerLogLookup
963 ** On a test machine with x86 64bit, the builtins are about 4 times faster
964 ** than the generic code.
965 **
966 */
967
CLog2UInt(UInt a)968 static Int CLog2UInt(UInt a)
969 {
970 #if SIZEOF_VOID_P == SIZEOF_INT && defined(HAVE___BUILTIN_CLZ)
971 return GMP_LIMB_BITS - 1 - __builtin_clz(a);
972 #elif SIZEOF_VOID_P == SIZEOF_LONG && defined(HAVE___BUILTIN_CLZL)
973 return GMP_LIMB_BITS - 1 - __builtin_clzl(a);
974 #elif SIZEOF_VOID_P == SIZEOF_LONG_LONG && defined(HAVE___BUILTIN_CLZLL)
975 return GMP_LIMB_BITS - 1 - __builtin_clzll(a);
976 #else
977 static const char LogTable256[256] = {
978 -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
979 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
980 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
981 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
982 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
983 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
984 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
985 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
986 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
987 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
988 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
989 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
990 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
991 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
992 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
993 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7
994 };
995
996 Int res = 0;
997 UInt b;
998 b = a >> 32; if (b) { res+=32; a=b; }
999 b = a >> 16; if (b) { res+=16; a=b; }
1000 b = a >> 8; if (b) { res+= 8; a=b; }
1001 return res + LogTable256[a];
1002 #endif
1003 }
1004
CLog2Int(Int a)1005 Int CLog2Int(Int a)
1006 {
1007 if (a == 0) return -1;
1008 if (a < 0) a = -a;
1009 return CLog2UInt(a);
1010 }
1011
1012 /****************************************************************************
1013 **
1014 *F FuncLog2Int( <self>, <n> ) . . . . . . . . . . . nr of bits of a GAP int
1015 **
1016 ** Given to GAP-Level as "Log2Int".
1017 */
FuncLog2Int(Obj self,Obj n)1018 static Obj FuncLog2Int(Obj self, Obj n)
1019 {
1020 RequireInt("Log2Int", n);
1021
1022 if (IS_INTOBJ(n)) {
1023 return INTOBJ_INT(CLog2Int(INT_INTOBJ(n)));
1024 }
1025
1026 UInt len = SIZE_INT(n) - 1;
1027 UInt a = CLog2UInt(CONST_ADDR_INT(n)[len]);
1028
1029 CHECK_INT(n);
1030
1031 #ifdef SYS_IS_64_BIT
1032 return INTOBJ_INT(len * GMP_LIMB_BITS + a);
1033 #else
1034 /* The final result is len * GMP_LIMB_BITS - d, which may not
1035 fit into an immediate integer (at least on a 32bit system) */
1036 return SumInt(ProdInt(INTOBJ_INT(len), INTOBJ_INT(GMP_LIMB_BITS)),
1037 INTOBJ_INT(a));
1038 #endif
1039 }
1040
1041 /****************************************************************************
1042 **
1043 *F FuncSTRING_INT( <self>, <n> ) . . . . . . convert an integer to a string
1044 **
1045 ** `FuncSTRING_INT' returns an immutable string representing the integer <n>
1046 **
1047 */
FuncSTRING_INT(Obj self,Obj n)1048 static Obj FuncSTRING_INT(Obj self, Obj n)
1049 {
1050 RequireInt("STRING_INT", n);
1051 return StringIntBase(n, 10);
1052 }
1053
1054
IntStringInternal(Obj string,const Char * str)1055 Obj IntStringInternal(Obj string, const Char *str)
1056 {
1057 Obj val; // value = <upp> * <pow> + <low>
1058 Obj upp; // upper part
1059 Int pow; // power
1060 Int low; // lower part
1061 Int sign; // is the integer negative
1062 UInt i; // loop variable
1063
1064 // if <string> is given, then we ignore <str>
1065 if (string)
1066 str = CONST_CSTR_STRING(string);
1067
1068 // get the sign, if any
1069 sign = 1;
1070 i = 0;
1071 if (str[i] == '-') {
1072 sign = -sign;
1073 i++;
1074 }
1075
1076 // collect the digits in groups of 8, for improved performance
1077 // note that 2^26 < 10^8 < 2^27, so the intermediate
1078 // values always fit into an immediate integer
1079 low = 0;
1080 pow = 1;
1081 upp = INTOBJ_INT(0);
1082 while (str[i] != '\0') {
1083 if (str[i] < '0' || str[i] > '9') {
1084 return Fail;
1085 }
1086 low = 10 * low + str[i] - '0';
1087 pow = 10 * pow;
1088 if (pow == 100000000L) {
1089 upp = ProdInt(upp, INTOBJ_INT(pow));
1090 upp = SumInt(upp, INTOBJ_INT(sign*low));
1091 // refresh 'str', in case the arithmetic operations triggered
1092 // a garbage collection
1093 if (string)
1094 str = CONST_CSTR_STRING(string);
1095 pow = 1;
1096 low = 0;
1097 }
1098 i++;
1099 }
1100
1101 // check if 0 char does not mark the end of the string
1102 if (string && i < GET_LEN_STRING(string))
1103 return Fail;
1104
1105 // compose the integer value
1106 if (upp == INTOBJ_INT(0)) {
1107 val = INTOBJ_INT(sign * low);
1108 }
1109 else if (pow == 1) {
1110 val = upp;
1111 }
1112 else {
1113 upp = ProdInt(upp, INTOBJ_INT(pow));
1114 val = SumInt(upp, INTOBJ_INT(sign * low));
1115 }
1116
1117 // return the integer value
1118 return val;
1119 }
1120
1121 /****************************************************************************
1122 **
1123 *F FuncINT_STRING( <self>, <string> ) . . . . convert a string to an integer
1124 **
1125 ** `FuncINT_STRING' returns an integer representing the string, or
1126 ** fail if the string is not a valid integer.
1127 **
1128 */
FuncINT_STRING(Obj self,Obj string)1129 static Obj FuncINT_STRING(Obj self, Obj string)
1130 {
1131 if( !IS_STRING(string) ) {
1132 return Fail;
1133 }
1134
1135 if( !IS_STRING_REP(string) ) {
1136 string = CopyToStringRep(string);
1137 }
1138
1139 return IntStringInternal(string, 0);
1140 }
1141
1142 /****************************************************************************
1143 **
1144 *F EqInt( <opL>, <opR> ) . . . . . . . . . . test if two integers are equal
1145 **
1146 ** 'EqInt' returns 1 if the two integer arguments <opL> and <opR> are equal
1147 ** and 0 otherwise.
1148 */
EqInt(Obj opL,Obj opR)1149 Int EqInt(Obj opL, Obj opR)
1150 {
1151 CHECK_INT(opL);
1152 CHECK_INT(opR);
1153
1154 // if at least one input is a small int, a naive equality test suffices
1155 if (IS_INTOBJ(opL) || IS_INTOBJ(opR))
1156 return opL == opR;
1157
1158 // compare the sign and size
1159 // note: at this point we know that opL and opR are proper bags, so
1160 // we can use TNUM_BAG instead of TNUM_OBJ
1161 if (TNUM_BAG(opL) != TNUM_BAG(opR) || SIZE_INT(opL) != SIZE_INT(opR))
1162 return 0;
1163
1164 return !mpn_cmp((mp_srcptr)CONST_ADDR_INT(opL),
1165 (mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opL));
1166 }
1167
1168 /****************************************************************************
1169 **
1170 *F LtInt( <opL>, <opR> ) . . . . . . test if an integer is less than another
1171 **
1172 ** 'LtInt' returns 1 if the integer <opL> is strictly less than the
1173 ** integer <opR> and 0 otherwise.
1174 */
LtInt(Obj opL,Obj opR)1175 Int LtInt(Obj opL, Obj opR)
1176 {
1177 Int res;
1178
1179 CHECK_INT(opL);
1180 CHECK_INT(opR);
1181
1182 /* compare two small integers */
1183 if (ARE_INTOBJS(opL, opR))
1184 return (Int)opL < (Int)opR;
1185
1186 // a small int is always less than a positive large int,
1187 // and always more than a negative large int
1188 if (IS_INTOBJ(opL))
1189 return IS_INTPOS(opR);
1190 if (IS_INTOBJ(opR))
1191 return IS_INTNEG(opL);
1192
1193 // at this point, both inputs are large integers, and we compare their
1194 // signs first
1195 if (TNUM_OBJ(opL) != TNUM_OBJ(opR))
1196 return IS_INTNEG(opL);
1197
1198 /* signs are equal; compare sizes and absolute values */
1199 if (SIZE_INT(opL) < SIZE_INT(opR))
1200 res = -1;
1201 else if (SIZE_INT(opL) > SIZE_INT(opR))
1202 res = +1;
1203 else
1204 res = mpn_cmp((mp_srcptr)CONST_ADDR_INT(opL),
1205 (mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opL));
1206
1207 /* if both arguments are negative, flip the result */
1208 if (IS_INTNEG(opL))
1209 res = -res;
1210
1211 return res < 0;
1212 }
1213
1214
1215 /****************************************************************************
1216 **
1217 *F SumOrDiffInt( <opL>, <opR>, <sign> ) . . . . sum or diff of two integers
1218 **
1219 ** 'SumOrDiffInt' returns the sum or difference of the two integer arguments
1220 ** <opL> and <opR>, depending on whether sign is +1 or -1.
1221 **
1222 ** 'SumOrDiffInt' is a little bit tricky since there are many different
1223 ** cases to handle: each operand can be positive or negative, small or large
1224 ** integer.
1225 */
SumOrDiffInt(Obj opL,Obj opR,Int sign)1226 static Obj SumOrDiffInt(Obj opL, Obj opR, Int sign)
1227 {
1228 UInt sizeL, sizeR;
1229 fake_mpz_t mpzL, mpzR, mpzResult;
1230 Obj result;
1231
1232 CHECK_INT(opL);
1233 CHECK_INT(opR);
1234
1235 /* handle trivial cases first */
1236 if (opR == INTOBJ_INT(0))
1237 return opL;
1238 if (opL == INTOBJ_INT(0)) {
1239 if (sign == 1)
1240 return opR;
1241 else
1242 return AInvInt(opR);
1243 }
1244
1245 sizeL = SIZE_INT_OR_INTOBJ(opL);
1246 sizeR = SIZE_INT_OR_INTOBJ(opR);
1247
1248 NEW_FAKEMPZ( mpzResult, sizeL > sizeR ? sizeL+1 : sizeR+1 );
1249 FAKEMPZ_GMPorINTOBJ(mpzL, opL);
1250 FAKEMPZ_GMPorINTOBJ(mpzR, opR);
1251
1252 /* add or subtract */
1253 if (sign == 1)
1254 mpz_add( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
1255 else
1256 mpz_sub( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
1257
1258 /* convert result to GAP object and return it */
1259 CHECK_FAKEMPZ(mpzResult);
1260 CHECK_FAKEMPZ(mpzL);
1261 CHECK_FAKEMPZ(mpzR);
1262 result = GMPorINTOBJ_FAKEMPZ( mpzResult );
1263 CHECK_INT(result);
1264 return result;
1265 }
1266
1267
1268 /****************************************************************************
1269 **
1270 *F SumInt( <opL>, <opR> ) . . . . . . . . . . . . . . . sum of two integers
1271 **
1272 */
SumInt(Obj opL,Obj opR)1273 inline Obj SumInt(Obj opL, Obj opR)
1274 {
1275 Obj sum;
1276
1277 if (!ARE_INTOBJS(opL, opR) || !SUM_INTOBJS(sum, opL, opR))
1278 sum = SumOrDiffInt(opL, opR, +1);
1279
1280 CHECK_INT(sum);
1281 return sum;
1282 }
1283
1284
1285 /****************************************************************************
1286 **
1287 *F DiffInt( <opL>, <opR> ) . . . . . . . . . . . difference of two integers
1288 **
1289 */
DiffInt(Obj opL,Obj opR)1290 inline Obj DiffInt(Obj opL, Obj opR)
1291 {
1292 Obj dif;
1293
1294 if (!ARE_INTOBJS(opL, opR) || !DIFF_INTOBJS(dif, opL, opR))
1295 dif = SumOrDiffInt(opL, opR, -1);
1296
1297 CHECK_INT(dif);
1298 return dif;
1299 }
1300
1301
1302 /****************************************************************************
1303 **
1304 *F ZeroInt(<op>) . . . . . . . . . . . . . . . . . . . . . zero of integers
1305 */
ZeroInt(Obj op)1306 static Obj ZeroInt(Obj op)
1307 {
1308 return INTOBJ_INT(0);
1309 }
1310
1311
1312 /****************************************************************************
1313 **
1314 *F AInvInt( <op> ) . . . . . . . . . . . . . additive inverse of an integer
1315 */
AInvInt(Obj op)1316 Obj AInvInt(Obj op)
1317 {
1318 Obj inv;
1319
1320 CHECK_INT(op);
1321
1322 // handle small integer
1323 if (IS_INTOBJ(op)) {
1324
1325 // special case (ugh)
1326 if (op == INTOBJ_MIN) {
1327 inv = NewBag( T_INTPOS, sizeof(mp_limb_t) );
1328 SET_VAL_LIMB0( inv, -INT_INTOBJ_MIN );
1329 }
1330
1331 // general case
1332 else {
1333 inv = INTOBJ_INT(-INT_INTOBJ(op));
1334 }
1335
1336 }
1337
1338 else {
1339 if (IS_INTPOS(op)) {
1340 // special case
1341 if (SIZE_INT(op) == 1 && VAL_LIMB0(op) == -INT_INTOBJ_MIN) {
1342 return INTOBJ_MIN;
1343 }
1344 else {
1345 inv = NewBag( T_INTNEG, SIZE_OBJ(op) );
1346 }
1347 }
1348
1349 else {
1350 inv = NewBag( T_INTPOS, SIZE_OBJ(op) );
1351 }
1352
1353 memcpy( ADDR_INT(inv), CONST_ADDR_INT(op), SIZE_OBJ(op) );
1354 }
1355
1356 CHECK_INT(inv);
1357 return inv;
1358
1359 }
1360
1361 /****************************************************************************
1362 **
1363 *F AbsInt(<op>) . . . . . . . . . . . . . . . . absolute value of an integer
1364 */
AbsInt(Obj op)1365 Obj AbsInt( Obj op )
1366 {
1367 Obj a;
1368 if ( IS_INTOBJ(op) ) {
1369 if ( ((Int)op) > 0 ) /* non-negative? */
1370 return op;
1371 else if ( op == INTOBJ_MIN ) {
1372 a = NewBag( T_INTPOS, sizeof(mp_limb_t) );
1373 SET_VAL_LIMB0( a, -INT_INTOBJ_MIN );
1374 return a;
1375 } else
1376 return (Obj)( 2 - (Int)op );
1377 }
1378 CHECK_INT(op);
1379 if ( IS_INTPOS(op) ) {
1380 return op;
1381 } else if ( IS_INTNEG(op) ) {
1382 a = NewBag( T_INTPOS, SIZE_OBJ(op) );
1383 memcpy( ADDR_INT(a), CONST_ADDR_INT(op), SIZE_OBJ(op) );
1384 return a;
1385 }
1386 return Fail;
1387 }
1388
FuncABS_INT(Obj self,Obj n)1389 static Obj FuncABS_INT(Obj self, Obj n)
1390 {
1391 RequireInt("AbsInt", n);
1392 Obj res = AbsInt(n);
1393 CHECK_INT(res);
1394 return res;
1395 }
1396
1397 /****************************************************************************
1398 **
1399 *F SignInt(<op>) . . . . . . . . . . . . . . . . . . . . sign of an integer
1400 */
SignInt(Obj op)1401 Obj SignInt( Obj op )
1402 {
1403 if ( IS_INTOBJ(op) ) {
1404 if ( op == INTOBJ_INT(0) )
1405 return INTOBJ_INT(0);
1406 else if ( ((Int)op) > (Int)INTOBJ_INT(0) )
1407 return INTOBJ_INT(1);
1408 else
1409 return INTOBJ_INT(-1);
1410 }
1411 CHECK_INT(op);
1412 if ( IS_INTPOS(op) ) {
1413 return INTOBJ_INT(1);
1414 } else if ( IS_INTNEG(op) ) {
1415 return INTOBJ_INT(-1);
1416 }
1417 return Fail;
1418 }
1419
FuncSIGN_INT(Obj self,Obj n)1420 static Obj FuncSIGN_INT(Obj self, Obj n)
1421 {
1422 RequireInt("SignInt", n);
1423 Obj res = SignInt(n);
1424 CHECK_INT(res);
1425 return res;
1426 }
1427
1428
1429 /****************************************************************************
1430 **
1431 *F ProdInt( <opL>, <opR> ) . . . . . . . . . . . . . product of two integers
1432 **
1433 ** 'ProdInt' returns the product of the two integer arguments <opL> and
1434 ** <opR>.
1435 */
ProdInt(Obj opL,Obj opR)1436 Obj ProdInt(Obj opL, Obj opR)
1437 {
1438 Obj prd; /* handle of the result bag */
1439 UInt sizeL, sizeR;
1440 fake_mpz_t mpzL, mpzR, mpzResult;
1441
1442 CHECK_INT(opL);
1443 CHECK_INT(opR);
1444
1445 /* multiplying two small integers */
1446 if ( ARE_INTOBJS( opL, opR ) ) {
1447
1448 /* multiply two small integers with a small product */
1449 if ( PROD_INTOBJS( prd, opL, opR ) ) {
1450 CHECK_INT(prd);
1451 return prd;
1452 }
1453 }
1454
1455 /* handle trivial cases first */
1456 if ( opL == INTOBJ_INT(0) || opR == INTOBJ_INT(1) )
1457 return opL;
1458 if ( opR == INTOBJ_INT(0) || opL == INTOBJ_INT(1) )
1459 return opR;
1460 if ( opR == INTOBJ_INT(-1) )
1461 return AInvInt( opL );
1462 if ( opL == INTOBJ_INT(-1) )
1463 return AInvInt( opR );
1464
1465 sizeL = SIZE_INT_OR_INTOBJ(opL);
1466 sizeR = SIZE_INT_OR_INTOBJ(opR);
1467
1468 NEW_FAKEMPZ( mpzResult, sizeL + sizeR );
1469 FAKEMPZ_GMPorINTOBJ( mpzL, opL );
1470 FAKEMPZ_GMPorINTOBJ( mpzR, opR );
1471
1472 /* multiply */
1473 mpz_mul( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
1474
1475 /* convert result to GAP object and return it */
1476 CHECK_FAKEMPZ(mpzResult);
1477 CHECK_FAKEMPZ(mpzL);
1478 CHECK_FAKEMPZ(mpzR);
1479 prd = GMPorINTOBJ_FAKEMPZ( mpzResult );
1480 CHECK_INT(prd);
1481 return prd;
1482 }
1483
1484
1485 /****************************************************************************
1486 **
1487 *F ProdIntObj(<n>,<op>) . . . . . . . . product of an integer and an object
1488 */
ProdIntObj(Obj n,Obj op)1489 static Obj ProdIntObj ( Obj n, Obj op )
1490 {
1491 Obj res = 0; /* result */
1492 UInt i, k; /* loop variables */
1493 mp_limb_t l; /* loop variable */
1494
1495 CHECK_INT(n);
1496
1497 /* if the integer is zero, return the neutral element of the operand */
1498 if ( n == INTOBJ_INT(0) ) {
1499 res = ZERO( op );
1500 }
1501
1502 /* if the integer is one, return the object if immutable -
1503 if mutable, add the object to its ZeroSameMutability to
1504 ensure correct mutability propagation */
1505 else if ( n == INTOBJ_INT(1) ) {
1506 if (IS_MUTABLE_OBJ(op))
1507 res = SUM(ZERO(op),op);
1508 else
1509 res = op;
1510 }
1511
1512 /* if the integer is minus one, return the inverse of the operand */
1513 else if ( n == INTOBJ_INT(-1) ) {
1514 res = AINV( op );
1515 }
1516
1517 /* if the integer is negative, invert the operand and the integer */
1518 else if ( IS_NEG_INT(n) ) {
1519 res = AINV( op );
1520 if ( res == Fail ) {
1521 ErrorMayQuit("Operations: <obj> must have an additive inverse", 0, 0);
1522 }
1523 res = PROD( AINV( n ), res );
1524 }
1525
1526 /* if the integer is small, compute the product by repeated doubling */
1527 /* the loop invariant is <result> = <k>*<res> + <l>*<op>, <l> < <k> */
1528 /* <res> = 0 means that <res> is the neutral element */
1529 else if ( IS_INTOBJ(n) && INT_INTOBJ(n) > 1 ) {
1530 res = 0;
1531 k = 1L << NR_SMALL_INT_BITS;
1532 l = INT_INTOBJ(n);
1533 while ( 0 < k ) {
1534 res = (res == 0 ? res : SUM( res, res ));
1535 if ( k <= l ) {
1536 res = (res == 0 ? op : SUM( res, op ));
1537 l = l - k;
1538 }
1539 k = k / 2;
1540 }
1541 }
1542
1543 /* if the integer is large, compute the product by repeated doubling */
1544 else if ( IS_INTPOS(n) ) {
1545 res = 0;
1546 for ( i = SIZE_INT(n); 0 < i; i-- ) {
1547 k = 8*sizeof(mp_limb_t);
1548 l = CONST_ADDR_INT(n)[i-1];
1549 while ( 0 < k ) {
1550 res = (res == 0 ? res : SUM( res, res ));
1551 k--;
1552 if ( (l >> k) & 1 ) {
1553 res = (res == 0 ? op : SUM( res, op ));
1554 }
1555 }
1556 }
1557 }
1558
1559 /* return the result */
1560 return res;
1561 }
1562
FuncPROD_INT_OBJ(Obj self,Obj opL,Obj opR)1563 static Obj FuncPROD_INT_OBJ(Obj self, Obj opL, Obj opR)
1564 {
1565 return ProdIntObj( opL, opR );
1566 }
1567
1568
1569 /****************************************************************************
1570 **
1571 *F OneInt(<op>) . . . . . . . . . . . . . . . . . . . . . one of an integer
1572 */
OneInt(Obj op)1573 static Obj OneInt(Obj op)
1574 {
1575 return INTOBJ_INT( 1 );
1576 }
1577
1578
1579 /****************************************************************************
1580 **
1581 *F PowInt( <opL>, <opR> ) . . . . . . . . . . . . . . . power of an integer
1582 **
1583 ** 'PowInt' returns the <opR>-th (an integer) power of the integer <opL>.
1584 */
PowInt(Obj opL,Obj opR)1585 Obj PowInt ( Obj opL, Obj opR )
1586 {
1587 Int i;
1588 Obj pow;
1589
1590 CHECK_INT(opL);
1591 CHECK_INT(opR);
1592
1593 if ( opR == INTOBJ_INT(0) ) {
1594 pow = INTOBJ_INT(1);
1595 }
1596 else if ( opL == INTOBJ_INT(0) ) {
1597 if ( IS_NEG_INT( opR ) ) {
1598 ErrorMayQuit("Integer operands: <base> must not be zero", 0, 0);
1599 }
1600 pow = INTOBJ_INT(0);
1601 }
1602 else if ( opL == INTOBJ_INT(1) ) {
1603 pow = INTOBJ_INT(1);
1604 }
1605 else if ( opL == INTOBJ_INT(-1) ) {
1606 pow = IS_EVEN_INT(opR) ? INTOBJ_INT(1) : INTOBJ_INT(-1);
1607 }
1608
1609 /* power with a large exponent */
1610 else if ( ! IS_INTOBJ(opR) ) {
1611 ErrorMayQuit("Integer operands: <exponent> is too large", 0, 0);
1612 }
1613
1614 /* power with a negative exponent */
1615 else if ( INT_INTOBJ(opR) < 0 ) {
1616 pow = QUO( INTOBJ_INT(1),
1617 PowInt( opL, INTOBJ_INT( -INT_INTOBJ(opR)) ) );
1618 }
1619
1620 /* findme - can we use the gmp function mpz_n_pow_ui? */
1621
1622 /* power with a small positive exponent, do it by a repeated squaring */
1623 else {
1624 pow = INTOBJ_INT(1);
1625 i = INT_INTOBJ(opR);
1626 while ( i != 0 ) {
1627 if ( i % 2 == 1 ) pow = ProdInt( pow, opL );
1628 if ( i > 1 ) opL = ProdInt( opL, opL );
1629 TakeInterrupt();
1630 i = i / 2;
1631 }
1632 }
1633
1634 /* return the power */
1635 CHECK_INT(pow);
1636 return pow;
1637 }
1638
1639
1640 /****************************************************************************
1641 **
1642 *F PowObjInt(<op>,<n>) . . . . . . . . . . power of an object and an integer
1643 */
PowObjInt(Obj op,Obj n)1644 static Obj PowObjInt(Obj op, Obj n)
1645 {
1646 Obj res = 0; /* result */
1647 UInt i, k; /* loop variables */
1648 mp_limb_t l; /* loop variable */
1649
1650 CHECK_INT(n);
1651
1652 /* if the integer is zero, return the neutral element of the operand */
1653 if ( n == INTOBJ_INT(0) ) {
1654 return ONE_MUT( op );
1655 }
1656
1657 /* if the integer is one, return a copy of the operand */
1658 else if ( n == INTOBJ_INT(1) ) {
1659 res = CopyObj( op, 1 );
1660 }
1661
1662 /* if the integer is minus one, return the inverse of the operand */
1663 else if ( n == INTOBJ_INT(-1) ) {
1664 res = INV_MUT( op );
1665 }
1666
1667 /* if the integer is negative, invert the operand and the integer */
1668 else if ( IS_NEG_INT(n) ) {
1669 res = INV_MUT( op );
1670 if ( res == Fail ) {
1671 ErrorMayQuit("Operations: <obj> must have an inverse", 0, 0);
1672 }
1673 res = POW( res, AINV( n ) );
1674 }
1675
1676 /* if the integer is small, compute the power by repeated squaring */
1677 /* the loop invariant is <result> = <res>^<k> * <op>^<l>, <l> < <k> */
1678 /* <res> = 0 means that <res> is the neutral element */
1679 else if ( IS_INTOBJ(n) && INT_INTOBJ(n) > 0 ) {
1680 res = 0;
1681 k = 1L << NR_SMALL_INT_BITS;
1682 l = INT_INTOBJ(n);
1683 while ( 0 < k ) {
1684 res = (res == 0 ? res : PROD( res, res ));
1685 if ( k <= l ) {
1686 res = (res == 0 ? op : PROD( res, op ));
1687 l = l - k;
1688 }
1689 k = k / 2;
1690 }
1691 }
1692
1693 /* if the integer is large, compute the power by repeated squaring */
1694 else if ( IS_INTPOS(n) ) {
1695 res = 0;
1696 for ( i = SIZE_INT(n); 0 < i; i-- ) {
1697 k = 8*sizeof(mp_limb_t);
1698 l = CONST_ADDR_INT(n)[i-1];
1699 while ( 0 < k ) {
1700 res = (res == 0 ? res : PROD( res, res ));
1701 k--;
1702 if ( (l >> k) & 1 ) {
1703 res = (res == 0 ? op : PROD( res, op ));
1704 }
1705 }
1706 }
1707 }
1708
1709 /* return the result */
1710 return res;
1711 }
1712
FuncPOW_OBJ_INT(Obj self,Obj opL,Obj opR)1713 static Obj FuncPOW_OBJ_INT(Obj self, Obj opL, Obj opR)
1714 {
1715 return PowObjInt( opL, opR );
1716 }
1717
1718
1719 /****************************************************************************
1720 **
1721 *F ModInt( <opL>, <opR> ) . . representative of residue class of an integer
1722 **
1723 ** 'ModInt' returns the smallest positive representative of the residue
1724 ** class of the integer <opL> modulo the integer <opR>.
1725 */
ModInt(Obj opL,Obj opR)1726 Obj ModInt(Obj opL, Obj opR)
1727 {
1728 Int i; /* loop count, value for small int */
1729 Int k; /* loop count, value for small int */
1730 UInt c; /* product of two digits */
1731 Obj mod; /* handle of the remainder bag */
1732 Obj quo; /* handle of the quotient bag */
1733
1734 CHECK_INT(opL);
1735 CHECK_INT(opR);
1736
1737 /* pathological case first */
1738 RequireNonzero("Integer operations", opR, "divisor");
1739
1740 /* compute the remainder of two small integers */
1741 if ( ARE_INTOBJS( opL, opR ) ) {
1742
1743 /* get the integer values */
1744 i = INT_INTOBJ(opL);
1745 k = INT_INTOBJ(opR);
1746
1747 /* compute the remainder, make sure we divide only positive numbers */
1748 i %= k;
1749 if (i < 0)
1750 i += k > 0 ? k : -k;
1751 mod = INTOBJ_INT(i);
1752 }
1753
1754 /* compute the remainder of a small integer by a large integer */
1755 else if ( IS_INTOBJ(opL) ) {
1756
1757 /* the small int -(1<<28) mod the large int (1<<28) is 0 */
1758 if ( opL == INTOBJ_MIN
1759 && ( IS_INTPOS(opR) )
1760 && ( SIZE_INT(opR) == 1 )
1761 && ( VAL_LIMB0(opR) == -INT_INTOBJ_MIN ) )
1762 mod = INTOBJ_INT(0);
1763
1764 /* in all other cases the remainder is equal the left operand */
1765 else if ( 0 <= INT_INTOBJ(opL) )
1766 mod = opL;
1767 else if ( IS_INTPOS(opR) )
1768 mod = SumOrDiffInt( opL, opR, 1 );
1769 else
1770 mod = SumOrDiffInt( opL, opR, -1 );
1771 }
1772
1773 /* compute the remainder of a large integer by a small integer */
1774 else if ( IS_INTOBJ(opR) ) {
1775
1776 /* get the integer value, make positive */
1777 i = INT_INTOBJ(opR); if ( i < 0 ) i = -i;
1778
1779 /* check whether right operand is a small power of 2 */
1780 if ( !(i & (i-1)) ) {
1781 c = VAL_LIMB0(opL) & (i-1);
1782 }
1783
1784 /* otherwise use the gmp function to divide */
1785 else {
1786 c = mpn_mod_1( (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL), i );
1787 }
1788
1789 // now c is the absolute value of the actual result. Thus, if the left
1790 // operand is negative, and c is non-zero, we have to adjust it.
1791 if (IS_INTPOS(opL) || c == 0)
1792 mod = INTOBJ_INT( c );
1793 else
1794 // even if opR is INT_INTOBJ_MIN, and hence i is INT_INTOBJ_MAX+1, we
1795 // have 0 <= i-c <= INT_INTOBJ_MAX, so i-c fits into a small integer
1796 mod = INTOBJ_INT( i - (Int)c );
1797
1798 }
1799
1800 /* compute the remainder of a large integer modulo a large integer */
1801 else {
1802
1803 /* trivial case first */
1804 if ( SIZE_INT(opL) < SIZE_INT(opR) ) {
1805 if ( IS_INTPOS(opL) )
1806 return opL;
1807 else if ( IS_INTPOS(opR) )
1808 mod = SumOrDiffInt( opL, opR, 1 );
1809 else
1810 mod = SumOrDiffInt( opL, opR, -1 );
1811 #if DEBUG_GMP
1812 assert( !IS_NEG_INT(mod) );
1813 #endif
1814 CHECK_INT(mod);
1815 return mod;
1816 }
1817
1818 mod = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
1819 ENSURE_BAG(mod);
1820
1821 quo = NewBag( T_INTPOS,
1822 (SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
1823 ENSURE_BAG(quo);
1824
1825 /* and let gmp do the work */
1826 mpn_tdiv_qr( (mp_ptr)ADDR_INT(quo), (mp_ptr)ADDR_INT(mod), 0,
1827 (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
1828 (mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opR) );
1829
1830 /* reduce to small integer if possible, otherwise shrink bag */
1831 mod = GMP_NORMALIZE( mod );
1832 mod = GMP_REDUCE( mod );
1833
1834 /* make the representative positive */
1835 if ( IS_NEG_INT(mod) ) {
1836 if ( IS_INTPOS(opR) )
1837 mod = SumOrDiffInt( mod, opR, 1 );
1838 else
1839 mod = SumOrDiffInt( mod, opR, -1 );
1840 }
1841
1842 }
1843
1844 /* return the result */
1845 #if DEBUG_GMP
1846 assert( !IS_NEG_INT(mod) );
1847 #endif
1848 CHECK_INT(mod);
1849 return mod;
1850 }
1851
1852
1853 /****************************************************************************
1854 **
1855 *F QuoInt( <opL>, <opR> ) . . . . . . . . . . . . . quotient of two integers
1856 **
1857 ** 'QuoInt' returns the integer part of the two integers <opL> and <opR>.
1858 **
1859 ** Note that this routine is not called from 'EvalQuo', the division of two
1860 ** integers yields a rational and is therefor performed in 'QuoRat'. This
1861 ** operation is however available through the internal function 'Quo'.
1862 */
QuoInt(Obj opL,Obj opR)1863 Obj QuoInt(Obj opL, Obj opR)
1864 {
1865 Int i; /* loop count, value for small int */
1866 Int k; /* loop count, value for small int */
1867 Obj quo; /* handle of the result bag */
1868 Obj rem; /* handle of the remainder bag */
1869
1870 CHECK_INT(opL);
1871 CHECK_INT(opR);
1872
1873 /* pathological case first */
1874 RequireNonzero("Integer operations", opR, "divisor");
1875
1876 /* divide two small integers */
1877 if ( ARE_INTOBJS( opL, opR ) ) {
1878
1879 /* the small int -(1<<28) divided by -1 is the large int (1<<28) */
1880 if ( opL == INTOBJ_MIN && opR == INTOBJ_INT(-1) ) {
1881 quo = NewBag( T_INTPOS, sizeof(mp_limb_t) );
1882 SET_VAL_LIMB0( quo, -INT_INTOBJ_MIN );
1883 return quo;
1884 }
1885
1886 /* get the integer values */
1887 i = INT_INTOBJ(opL);
1888 k = INT_INTOBJ(opR);
1889
1890 // divide, truncated towards zero; this is also what section 6.5.5 of the
1891 // C99 standard guarantees for integer division
1892 quo = INTOBJ_INT(i / k);
1893 }
1894
1895 /* divide a small integer by a large one */
1896 else if ( IS_INTOBJ(opL) ) {
1897
1898 /* the small int -(1<<28) divided by the large int (1<<28) is -1 */
1899 if ( opL == INTOBJ_MIN
1900 && IS_INTPOS(opR) && SIZE_INT(opR) == 1
1901 && VAL_LIMB0(opR) == -INT_INTOBJ_MIN )
1902 quo = INTOBJ_INT(-1);
1903
1904 /* in all other cases the quotient is of course zero */
1905 else
1906 quo = INTOBJ_INT(0);
1907
1908 }
1909
1910 /* divide a large integer by a small integer */
1911 else if ( IS_INTOBJ(opR) ) {
1912
1913 k = INT_INTOBJ(opR);
1914
1915 /* allocate a bag for the result and set up the pointers */
1916 if ( IS_INTNEG(opL) == ( k < 0 ) )
1917 quo = NewBag( T_INTPOS, SIZE_OBJ(opL) );
1918 else
1919 quo = NewBag( T_INTNEG, SIZE_OBJ(opL) );
1920
1921 ENSURE_BAG(quo);
1922
1923 if ( k < 0 ) k = -k;
1924
1925 /* use gmp function for dividing by a 1-limb number */
1926 mpn_divrem_1( (mp_ptr)ADDR_INT(quo), 0,
1927 (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
1928 k );
1929 }
1930
1931 /* divide a large integer by a large integer */
1932 else {
1933
1934 /* trivial case first */
1935 if ( SIZE_INT(opL) < SIZE_INT(opR) )
1936 return INTOBJ_INT(0);
1937
1938 /* create a new bag for the remainder */
1939 rem = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
1940 ENSURE_BAG(rem);
1941
1942 /* allocate a bag for the quotient */
1943 if ( TNUM_OBJ(opL) == TNUM_OBJ(opR) )
1944 quo = NewBag( T_INTPOS,
1945 (SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
1946 else
1947 quo = NewBag( T_INTNEG,
1948 (SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
1949 ENSURE_BAG(quo);
1950
1951 mpn_tdiv_qr( (mp_ptr)ADDR_INT(quo), (mp_ptr)ADDR_INT(rem), 0,
1952 (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
1953 (mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opR) );
1954 }
1955
1956 /* normalize and return the result */
1957 quo = GMP_NORMALIZE(quo);
1958 quo = GMP_REDUCE( quo );
1959 return quo;
1960 }
1961
1962
1963 /****************************************************************************
1964 **
1965 *F FuncQUO_INT( <self>, <a>, <b> ) . . . . . . . internal function 'QuoInt'
1966 **
1967 ** 'FuncQUO_INT' implements the internal function 'QuoInt'.
1968 **
1969 ** 'QuoInt( <a>, <b> )'
1970 **
1971 ** 'Quo' returns the integer part of the quotient of its integer operands.
1972 ** If <a> and <b> are positive 'Quo( <a>, <b> )' is the largest positive
1973 ** integer <q> such that '<q> * <b> \<= <a>'. If <a> or <b> or both are
1974 ** negative we define 'Abs( Quo(<a>,<b>) ) = Quo( Abs(<a>), Abs(<b>) )' and
1975 ** 'Sign( Quo(<a>,<b>) ) = Sign(<a>) * Sign(<b>)'. Dividing by 0 causes an
1976 ** error. 'Rem' (see "Rem") can be used to compute the remainder.
1977 */
FuncQUO_INT(Obj self,Obj a,Obj b)1978 static Obj FuncQUO_INT(Obj self, Obj a, Obj b)
1979 {
1980 RequireInt("QuoInt", a);
1981 RequireInt("QuoInt", b);
1982 return QuoInt(a, b);
1983 }
1984
1985
1986 /****************************************************************************
1987 **
1988 *F RemInt( <opL>, <opR> ) . . . . . . . . . . . . remainder of two integers
1989 **
1990 ** 'RemInt' returns the remainder of the quotient of the integers <opL>
1991 ** and <opR>.
1992 **
1993 ** Note that the remainder is different from the value returned by the 'mod'
1994 ** operator which is always positive, while the result of 'RemInt' has
1995 ** the same sign as <opL>.
1996 */
RemInt(Obj opL,Obj opR)1997 Obj RemInt(Obj opL, Obj opR)
1998 {
1999 Int i; /* loop count, value for small int */
2000 Int k; /* loop count, value for small int */
2001 UInt c;
2002 Obj rem; /* handle of the remainder bag */
2003 Obj quo; /* handle of the quotient bag */
2004
2005 CHECK_INT(opL);
2006 CHECK_INT(opR);
2007
2008 /* pathological case first */
2009 RequireNonzero("Integer operations", opR, "divisor");
2010
2011 /* compute the remainder of two small integers */
2012 if ( ARE_INTOBJS( opL, opR ) ) {
2013
2014 /* get the integer values */
2015 i = INT_INTOBJ(opL);
2016 k = INT_INTOBJ(opR);
2017
2018 // compute the remainder with sign matching that of the dividend i;
2019 // this matches the sign of i % k as specified in section 6.5.5 of
2020 // the C99 standard, which indicates division truncates towards zero,
2021 // and the invariant i%k == i - (i/k)*k holds
2022 rem = INTOBJ_INT(i % k);
2023 }
2024
2025 /* compute the remainder of a small integer by a large integer */
2026 else if ( IS_INTOBJ(opL) ) {
2027
2028 /* the small int -(1<<28) rem the large int (1<<28) is 0 */
2029 if ( opL == INTOBJ_MIN
2030 && IS_INTPOS(opR) && SIZE_INT(opR) == 1
2031 && VAL_LIMB0(opR) == -INT_INTOBJ_MIN )
2032 rem = INTOBJ_INT(0);
2033
2034 /* in all other cases the remainder is equal the left operand */
2035 else
2036 rem = opL;
2037 }
2038
2039 /* compute the remainder of a large integer by a small integer */
2040 else if ( IS_INTOBJ(opR) ) {
2041
2042 /* get the integer value, make positive */
2043 i = INT_INTOBJ(opR); if ( i < 0 ) i = -i;
2044
2045 /* check whether right operand is a small power of 2 */
2046 if ( !(i & (i-1)) ) {
2047 c = VAL_LIMB0(opL) & (i-1);
2048 }
2049
2050 /* otherwise use the gmp function to divide */
2051 else {
2052 c = mpn_mod_1( (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL), i );
2053 }
2054
2055 // adjust c for the sign of the left operand
2056 if ( IS_INTPOS(opL) )
2057 rem = INTOBJ_INT( c );
2058 else
2059 rem = INTOBJ_INT( -(Int)c );
2060
2061 }
2062
2063 /* compute the remainder of a large integer modulo a large integer */
2064 else {
2065
2066 /* trivial case first */
2067 if ( SIZE_INT(opL) < SIZE_INT(opR) )
2068 return opL;
2069
2070 rem = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+1)*sizeof(mp_limb_t) );
2071 ENSURE_BAG(rem);
2072
2073 quo = NewBag( T_INTPOS,
2074 (SIZE_INT(opL)-SIZE_INT(opR)+1)*sizeof(mp_limb_t) );
2075 ENSURE_BAG(quo);
2076
2077 /* and let gmp do the work */
2078 mpn_tdiv_qr( (mp_ptr)ADDR_INT(quo), (mp_ptr)ADDR_INT(rem), 0,
2079 (mp_srcptr)CONST_ADDR_INT(opL), SIZE_INT(opL),
2080 (mp_srcptr)CONST_ADDR_INT(opR), SIZE_INT(opR) );
2081
2082 /* reduce to small integer if possible, otherwise shrink bag */
2083 rem = GMP_NORMALIZE( rem );
2084 rem = GMP_REDUCE( rem );
2085
2086 }
2087
2088 /* return the result */
2089 CHECK_INT(rem);
2090 return rem;
2091 }
2092
2093
2094 /****************************************************************************
2095 **
2096 *F FuncREM_INT( <self>, <a>, <b> ) . . . . . . . internal function 'RemInt'
2097 **
2098 ** 'FuncREM_INT' implements the internal function 'RemInt'.
2099 **
2100 ** 'RemInt( <a>, <b> )'
2101 **
2102 ** 'RemInt' returns the remainder of its two integer operands, i.e., if <k>
2103 ** is not equal to zero 'Rem( <i>, <k> ) = <i> - <k> * Quo( <i>, <k> )'.
2104 ** Note that the rules given for 'Quo' imply that 'Rem( <i>, <k> )' has the
2105 ** same sign as <i> and its absolute value is strictly less than the
2106 ** absolute value of <k>. Dividing by 0 causes an error.
2107 */
FuncREM_INT(Obj self,Obj a,Obj b)2108 static Obj FuncREM_INT(Obj self, Obj a, Obj b)
2109 {
2110 RequireInt("RemInt", a);
2111 RequireInt("RemInt", b);
2112 return RemInt(a, b);
2113 }
2114
2115
2116 /****************************************************************************
2117 **
2118 *F GcdInt( <opL>, <opR> ) . . . . . . . . . . . . . . . gcd of two integers
2119 **
2120 ** 'GcdInt' returns the gcd of the two integers <opL> and <opR>.
2121 **
2122 ** It is called from 'FuncGCD_INT' and from the rational package.
2123 */
GcdInt(Obj opL,Obj opR)2124 Obj GcdInt ( Obj opL, Obj opR )
2125 {
2126 UInt sizeL, sizeR;
2127 fake_mpz_t mpzL, mpzR, mpzResult;
2128 Obj result;
2129
2130 CHECK_INT(opL);
2131 CHECK_INT(opR);
2132
2133 if (opL == INTOBJ_INT(0)) return AbsInt(opR);
2134 if (opR == INTOBJ_INT(0)) return AbsInt(opL);
2135
2136 sizeL = SIZE_INT_OR_INTOBJ(opL);
2137 sizeR = SIZE_INT_OR_INTOBJ(opR);
2138
2139 // for small inputs, run Euclid directly
2140 if (sizeL == 1 || sizeR == 1) {
2141 if (sizeR != 1) {
2142 SWAP(Obj, opL, opR);
2143 }
2144 UInt r = AbsOfSmallInt(opR);
2145 FAKEMPZ_GMPorINTOBJ(mpzL, opL);
2146 r = mpz_gcd_ui(0, MPZ_FAKEMPZ(mpzL), r);
2147 CHECK_FAKEMPZ(mpzL);
2148 return ObjInt_UInt(r);
2149 }
2150
2151 NEW_FAKEMPZ( mpzResult, sizeL < sizeR ? sizeL : sizeR );
2152 FAKEMPZ_GMPorINTOBJ( mpzL, opL );
2153 FAKEMPZ_GMPorINTOBJ( mpzR, opR );
2154
2155 /* compute the gcd */
2156 mpz_gcd( MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
2157
2158 /* convert result to GAP object and return it */
2159 CHECK_FAKEMPZ(mpzResult);
2160 CHECK_FAKEMPZ(mpzL);
2161 CHECK_FAKEMPZ(mpzR);
2162 result = GMPorINTOBJ_FAKEMPZ( mpzResult );
2163 CHECK_INT(result);
2164 return result;
2165 }
2166
2167
2168 /****************************************************************************
2169 **
2170 *F FuncGCD_INT(<self>,<a>,<b>) . . . . . . . internal function 'GcdInt'
2171 **
2172 ** 'FuncGCD_INT' implements the internal function 'GcdInt'.
2173 **
2174 ** 'GcdInt( <a>, <b> )'
2175 **
2176 ** 'Gcd' returns the greatest common divisor of the two integers <a> and
2177 ** <b>, i.e., the greatest integer that divides both <a> and <b>. The
2178 ** greatest common divisor is never negative, even if the arguments are. We
2179 ** define $gcd( a, 0 ) = gcd( 0, a ) = abs( a )$ and $gcd( 0, 0 ) = 0$.
2180 */
FuncGCD_INT(Obj self,Obj a,Obj b)2181 static Obj FuncGCD_INT(Obj self, Obj a, Obj b)
2182 {
2183 RequireInt("GcdInt", a);
2184 RequireInt("GcdInt", b);
2185 return GcdInt(a, b);
2186 }
2187
2188 /****************************************************************************
2189 **
2190 *F LcmInt( <opL>, <opR> ) . . . . . . . . . . . . . . . lcm of two integers
2191 **
2192 ** 'LcmInt' returns the lcm of the two integers <opL> and <opR>.
2193 */
LcmInt(Obj opL,Obj opR)2194 Obj LcmInt(Obj opL, Obj opR)
2195 {
2196 UInt sizeL, sizeR;
2197 fake_mpz_t mpzL, mpzR, mpzResult;
2198 Obj result;
2199
2200 CHECK_INT(opL);
2201 CHECK_INT(opR);
2202
2203 if (opL == INTOBJ_INT(0) || opR == INTOBJ_INT(0))
2204 return INTOBJ_INT(0);
2205
2206 if (IS_INTOBJ(opL) || IS_INTOBJ(opR)) {
2207 if (!IS_INTOBJ(opR)) {
2208 SWAP(Obj, opL, opR);
2209 }
2210 Obj gcd = GcdInt(opL, opR);
2211 opR = QuoInt(opR, gcd);
2212 return AbsInt(ProdInt(opL, opR));
2213 }
2214
2215 sizeL = SIZE_INT_OR_INTOBJ(opL);
2216 sizeR = SIZE_INT_OR_INTOBJ(opR);
2217
2218 NEW_FAKEMPZ(mpzResult, sizeL + sizeR);
2219 FAKEMPZ_GMPorINTOBJ(mpzL, opL);
2220 FAKEMPZ_GMPorINTOBJ(mpzR, opR);
2221
2222 /* compute the gcd */
2223 mpz_lcm(MPZ_FAKEMPZ(mpzResult), MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR));
2224
2225 /* convert result to GAP object and return it */
2226 CHECK_FAKEMPZ(mpzResult);
2227 CHECK_FAKEMPZ(mpzL);
2228 CHECK_FAKEMPZ(mpzR);
2229 result = GMPorINTOBJ_FAKEMPZ(mpzResult);
2230 CHECK_INT(result);
2231 return result;
2232 }
2233
FuncLCM_INT(Obj self,Obj a,Obj b)2234 static Obj FuncLCM_INT(Obj self, Obj a, Obj b)
2235 {
2236 RequireInt("LcmInt", a);
2237 RequireInt("LcmInt", b);
2238 return LcmInt(a, b);
2239 }
2240
2241 /****************************************************************************
2242 **
2243 */
FuncFACTORIAL_INT(Obj self,Obj n)2244 static Obj FuncFACTORIAL_INT(Obj self, Obj n)
2245 {
2246 RequireNonnegativeSmallInt("Factorial", n);
2247
2248 mpz_t mpzResult;
2249 mpz_init(mpzResult);
2250 mpz_fac_ui(mpzResult, INT_INTOBJ(n));
2251
2252 // convert mpzResult into a GAP integer object.
2253 Obj result = GMPorINTOBJ_MPZ(mpzResult);
2254
2255 // free mpzResult
2256 mpz_clear(mpzResult);
2257
2258 return result;
2259 }
2260
2261
2262 /****************************************************************************
2263 **
2264 */
BinomialInt(Obj n,Obj k)2265 Obj BinomialInt(Obj n, Obj k)
2266 {
2267 Int negate_result = 0;
2268
2269 // deal with k <= 1
2270 if (k == INTOBJ_INT(0))
2271 return INTOBJ_INT(1);
2272 if (k == INTOBJ_INT(1))
2273 return n;
2274 if (IS_NEG_INT(k))
2275 return INTOBJ_INT(0);
2276
2277 // deal with n < 0
2278 if (IS_NEG_INT(n)) {
2279 // use the identity Binomial(n,k) = (-1)^k * Binomial(-n+k-1, k)
2280 negate_result = IS_ODD_INT(k);
2281 n = DiffInt(DiffInt(k,n), INTOBJ_INT(1));
2282 }
2283
2284 // deal with n <= k
2285 if (n == k)
2286 return negate_result ? INTOBJ_INT(-1) : INTOBJ_INT(1);
2287 if (LtInt(n, k))
2288 return INTOBJ_INT(0);
2289
2290 // deal with n-k < k <=> n < 2k
2291 Obj k2 = DiffInt(n, k);
2292 if (LtInt(k2, k))
2293 k = k2;
2294
2295 // From here on, we only support single limb integers for k. Anything else
2296 // would lead to output too big for storage anyway, at least on a 64 bit
2297 // system. To be specific, in that case n >= 2k and k >= 2^60. Thus, in
2298 // the *best* case, we are trying to compute the central binomial
2299 // coefficient Binomial(2k k) for k = 2^60. This value is approximately
2300 // (4^k / sqrt(pi*k)); taking the logarithm and dividing by 8 yields that
2301 // we need about k/4 = 2^58 bytes to store this value. No computer in the
2302 // foreseeable future will be able to store the result (and much less
2303 // compute it in a reasonable time).
2304 //
2305 // On 32 bit systems, the limit is k = 2^28, and then the central binomial
2306 // coefficient Binomial(2k,k) takes up about 64 MB, so that would still be
2307 // feasible. However, GMP does not support computing binomials when k is
2308 // larger than a single limb, so we'd have to implement this on our own.
2309 //
2310 // Since GAP previously was effectively unable to compute such binomial
2311 // coefficients (unless you were willing to wait for a few days or so), we
2312 // simply do not implement this on 32bit systems, and instead return Fail,
2313 // jut as we do on 64 bit. If somebody complains about this, we can still
2314 // look into implementing this (and in the meantime, tell the user to use
2315 // the old GAP version of this function).
2316
2317 if (SIZE_INT_OR_INTOBJ(k) > 1)
2318 return Fail;
2319
2320 UInt K = IS_INTOBJ(k) ? INT_INTOBJ(k) : VAL_LIMB0(k);
2321 mpz_t mpzResult;
2322 mpz_init( mpzResult );
2323
2324 if (SIZE_INT_OR_INTOBJ(n) == 1) {
2325 UInt N = IS_INTOBJ(n) ? INT_INTOBJ(n) : VAL_LIMB0(n);
2326 mpz_bin_uiui(mpzResult, N, K);
2327 } else {
2328 fake_mpz_t mpzN;
2329 FAKEMPZ_GMPorINTOBJ( mpzN, n );
2330 mpz_bin_ui(mpzResult, MPZ_FAKEMPZ(mpzN), K);
2331 }
2332
2333 // adjust sign of result
2334 if (negate_result)
2335 mpzResult->_mp_size = -mpzResult->_mp_size;
2336
2337 // convert mpzResult into a GAP integer object.
2338 Obj result = GMPorINTOBJ_MPZ(mpzResult);
2339
2340 // free mpzResult
2341 mpz_clear( mpzResult );
2342
2343 return result;
2344 }
2345
FuncBINOMIAL_INT(Obj self,Obj n,Obj k)2346 static Obj FuncBINOMIAL_INT(Obj self, Obj n, Obj k)
2347 {
2348 RequireInt("Binomial", n);
2349 RequireInt("Binomial", k);
2350 return BinomialInt(n, k);
2351 }
2352
2353
2354 /****************************************************************************
2355 **
2356 */
FuncJACOBI_INT(Obj self,Obj n,Obj m)2357 static Obj FuncJACOBI_INT(Obj self, Obj n, Obj m)
2358 {
2359 fake_mpz_t mpzL, mpzR;
2360 int result;
2361
2362 RequireInt("Jacobi", n);
2363 RequireInt("Jacobi", m);
2364
2365 CHECK_INT(n);
2366 CHECK_INT(m);
2367
2368 FAKEMPZ_GMPorINTOBJ(mpzL, n);
2369 FAKEMPZ_GMPorINTOBJ(mpzR, m);
2370
2371 result = mpz_kronecker( MPZ_FAKEMPZ(mpzL), MPZ_FAKEMPZ(mpzR) );
2372 CHECK_FAKEMPZ(mpzL);
2373 CHECK_FAKEMPZ(mpzR);
2374
2375 return INTOBJ_INT( result );
2376 }
2377
2378
2379 /****************************************************************************
2380 **
2381 */
FuncPVALUATION_INT(Obj self,Obj n,Obj p)2382 static Obj FuncPVALUATION_INT(Obj self, Obj n, Obj p)
2383 {
2384 fake_mpz_t mpzN, mpzP;
2385 mpz_t mpzResult;
2386 int k;
2387
2388 RequireInt("PValuation", n);
2389 RequireInt("PValuation", p);
2390
2391 CHECK_INT(n);
2392 CHECK_INT(p);
2393
2394 RequireNonzero("PValuation", p, "p");
2395
2396 if (SIZE_INT_OR_INTOBJ(n) == 1 && SIZE_INT_OR_INTOBJ(p) == 1) {
2397 UInt N = AbsOfSmallInt(n);
2398 UInt P = AbsOfSmallInt(p);
2399 if (N == 0 || P == 1) return INTOBJ_INT(0);
2400 k = 0;
2401 while (N % P == 0) {
2402 N /= P;
2403 k++;
2404 }
2405 return INTOBJ_INT(k);
2406 }
2407
2408 /* For certain values of p, mpz_remove replaces its "dest" argument
2409 and tries to deallocate the original mpz_t in it. This means
2410 we cannot use a fake_mpz_t for it. However, we are not really
2411 interested in it anyway. */
2412 mpz_init( mpzResult );
2413 FAKEMPZ_GMPorINTOBJ( mpzN, n );
2414 FAKEMPZ_GMPorINTOBJ( mpzP, p );
2415
2416 k = mpz_remove( mpzResult, MPZ_FAKEMPZ(mpzN), MPZ_FAKEMPZ(mpzP) );
2417 CHECK_FAKEMPZ(mpzN);
2418 CHECK_FAKEMPZ(mpzP);
2419
2420 /* throw away mpzResult -- it equals m / p^k */
2421 mpz_clear( mpzResult );
2422
2423 return INTOBJ_INT( k );
2424 }
2425
2426
2427 /****************************************************************************
2428 **
2429 */
FuncROOT_INT(Obj self,Obj n,Obj k)2430 static Obj FuncROOT_INT(Obj self, Obj n, Obj k)
2431 {
2432 fake_mpz_t n_mpz, result_mpz;
2433
2434 RequireInt("Root", n);
2435 RequireInt("Root", k);
2436
2437 if (!IS_POS_INT(k))
2438 ErrorMayQuit("Root: <k> must be a positive integer", 0, 0);
2439 if (IS_NEG_INT(n) && IS_EVEN_INT(k))
2440 ErrorMayQuit("Root: <n> is negative but <k> is even", 0, 0);
2441
2442 if (k == INTOBJ_INT(1) || n == INTOBJ_INT(0))
2443 return n;
2444
2445 if (!IS_INTOBJ(k)) {
2446 #ifdef SYS_IS_64_BIT
2447 // if k is not immediate, i.e., k >= 2^60, then the root is 1 unless
2448 // n >= 2^k >= 2^(2^60), which means storage for n must be at least
2449 // 2^60 bits, i.e. 2^57 bytes. That's more RAM than anybody has for
2450 // the near future.
2451 return IS_NEG_INT(n) ? INTOBJ_INT(-1) : INTOBJ_INT(1);
2452 #else
2453 // if k is not immediate, i.e., k >= 2^28, then the root is 1 unless
2454 // n >= 2^k >= 2^(2^28), which means storage for n must be at least
2455 // 2^28 bits, i.e., 2^25 bytes, or 2^23 words (each 32 bits).
2456 if (SIZE_INT_OR_INTOBJ(n) < (1 << 23))
2457 return IS_NEG_INT(n) ? INTOBJ_INT(-1) : INTOBJ_INT(1);
2458 else
2459 return Fail; // return fail so that high level code can handle
2460 // this
2461 #endif
2462 }
2463
2464 UInt K = INT_INTOBJ(k);
2465 UInt root_size = 1 + (SIZE_INT_OR_INTOBJ(n) - 1) / K;
2466 NEW_FAKEMPZ(result_mpz, root_size);
2467 FAKEMPZ_GMPorINTOBJ(n_mpz, n);
2468
2469 if (K == 2)
2470 mpz_sqrt(MPZ_FAKEMPZ(result_mpz), MPZ_FAKEMPZ(n_mpz));
2471 else
2472 mpz_root(MPZ_FAKEMPZ(result_mpz), MPZ_FAKEMPZ(n_mpz), K);
2473
2474 CHECK_FAKEMPZ(result_mpz);
2475 CHECK_FAKEMPZ(n_mpz);
2476
2477 return GMPorINTOBJ_FAKEMPZ(result_mpz);
2478 }
2479
2480
2481 /****************************************************************************
2482 **
2483 */
InverseModInt(Obj base,Obj mod)2484 Obj InverseModInt(Obj base, Obj mod)
2485 {
2486 fake_mpz_t base_mpz, mod_mpz, result_mpz;
2487 int success;
2488
2489 CHECK_INT(base);
2490 CHECK_INT(mod);
2491
2492 RequireNonzero("InverseModInt", mod, "mod");
2493 if (mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1))
2494 return INTOBJ_INT(0);
2495 if (base == INTOBJ_INT(0))
2496 return Fail;
2497
2498 // handle small inputs separately
2499 if (IS_INTOBJ(mod)) {
2500
2501 Int a = INT_INTOBJ(mod);
2502 if (a < 0)
2503 a = -a;
2504
2505 Int b = INT_INTOBJ(ModInt(base, mod));
2506
2507 Int aL = 0; // cofactor of a
2508 Int bL = 1; // cofactor of b
2509
2510 // extended Euclidean algorithm
2511 while (b != 0) {
2512 Int hdQ = a / b;
2513 Int c = b;
2514 Int cL = bL;
2515 b = a - hdQ * b;
2516 bL = aL - hdQ * bL;
2517 a = c;
2518 aL = cL;
2519 }
2520 if (a != 1)
2521 return Fail;
2522 return ModInt(INTOBJ_INT(aL), mod);
2523 }
2524
2525 NEW_FAKEMPZ(result_mpz, SIZE_INT_OR_INTOBJ(mod) + 1);
2526 FAKEMPZ_GMPorINTOBJ(base_mpz, base);
2527 FAKEMPZ_GMPorINTOBJ(mod_mpz, mod);
2528
2529 success = mpz_invert(MPZ_FAKEMPZ(result_mpz), MPZ_FAKEMPZ(base_mpz),
2530 MPZ_FAKEMPZ(mod_mpz));
2531
2532 if (!success)
2533 return Fail;
2534
2535 CHECK_FAKEMPZ(result_mpz);
2536 CHECK_FAKEMPZ(base_mpz);
2537 CHECK_FAKEMPZ(mod_mpz);
2538
2539 return GMPorINTOBJ_FAKEMPZ(result_mpz);
2540 }
2541
2542 /****************************************************************************
2543 **
2544 */
FuncINVMODINT(Obj self,Obj base,Obj mod)2545 static Obj FuncINVMODINT(Obj self, Obj base, Obj mod)
2546 {
2547 RequireInt("InverseModInt", base);
2548 RequireInt("InverseModInt", mod);
2549 return InverseModInt(base, mod);
2550 }
2551
2552
2553 /****************************************************************************
2554 **
2555 */
FuncPOWERMODINT(Obj self,Obj base,Obj exp,Obj mod)2556 static Obj FuncPOWERMODINT(Obj self, Obj base, Obj exp, Obj mod)
2557 {
2558 fake_mpz_t base_mpz, exp_mpz, mod_mpz, result_mpz;
2559
2560 RequireInt("PowerModInt", base);
2561 RequireInt("PowerModInt", exp);
2562 RequireInt("PowerModInt", mod);
2563
2564 CHECK_INT(base);
2565 CHECK_INT(exp);
2566 CHECK_INT(mod);
2567
2568 RequireNonzero("PowerModInt", mod, "mod");
2569 if ( mod == INTOBJ_INT(1) || mod == INTOBJ_INT(-1) )
2570 return INTOBJ_INT(0);
2571
2572 if ( IS_NEG_INT(exp) ) {
2573 base = InverseModInt( base, mod );
2574 if (base == Fail)
2575 ErrorMayQuit( "PowerModInt: negative <exp> but <base> is not invertible modulo <mod>", 0L, 0L );
2576 exp = AInvInt(exp);
2577 }
2578
2579 NEW_FAKEMPZ( result_mpz, SIZE_INT_OR_INTOBJ(mod) );
2580 FAKEMPZ_GMPorINTOBJ( base_mpz, base );
2581 FAKEMPZ_GMPorINTOBJ( exp_mpz, exp );
2582 FAKEMPZ_GMPorINTOBJ( mod_mpz, mod );
2583
2584 mpz_powm( MPZ_FAKEMPZ(result_mpz), MPZ_FAKEMPZ(base_mpz),
2585 MPZ_FAKEMPZ(exp_mpz), MPZ_FAKEMPZ(mod_mpz) );
2586
2587 CHECK_FAKEMPZ(result_mpz);
2588 CHECK_FAKEMPZ(base_mpz);
2589 CHECK_FAKEMPZ(exp_mpz);
2590 CHECK_FAKEMPZ(mod_mpz);
2591
2592 return GMPorINTOBJ_FAKEMPZ( result_mpz );
2593 }
2594
2595
2596 /****************************************************************************
2597 **
2598 */
FuncIS_PROBAB_PRIME_INT(Obj self,Obj n,Obj reps)2599 static Obj FuncIS_PROBAB_PRIME_INT(Obj self, Obj n, Obj reps)
2600 {
2601 fake_mpz_t n_mpz;
2602 Int res;
2603
2604 RequireInt("IsProbablyPrimeInt", n);
2605 UInt r = GetPositiveSmallInt("IsProbablyPrimeInt", reps);
2606
2607 CHECK_INT(n);
2608
2609 FAKEMPZ_GMPorINTOBJ( n_mpz, n );
2610
2611 res = mpz_probab_prime_p(MPZ_FAKEMPZ(n_mpz), r);
2612
2613 if (res == 2) return True; /* definitely prime */
2614 if (res == 0) return False; /* definitely not prime */
2615 return Fail; /* probably prime */
2616 }
2617
2618
2619 /****************************************************************************
2620 **
2621 ** * * * * * * * "Mersenne twister" random numbers * * * * * * * * * * * * *
2622 **
2623 ** Part of this code for fast generation of 32 bit pseudo random numbers
2624 ** with a period of length 2^19937-1 and a 623-dimensional equidistribution
2625 ** is taken from:
2626 ** http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
2627 ** (Also look in Wikipedia for "Mersenne twister".)
2628 */
2629
2630 /****************************************************************************
2631 **
2632 *F RandomIntegerMT( <mtstr>, <nrbits> )
2633 **
2634 ** Returns an integer with at most <nrbits> bits in uniform distribution.
2635 ** <nrbits> must be a small integer. <mtstr> is a string as returned by
2636 ** InitRandomMT.
2637 **
2638 ** Implementation details are a bit tricky to obtain the same random
2639 ** integers on 32 bit and 64 bit machines (which have different long
2640 ** integer digit lengths and different ranges of small integers).
2641 **
2642 */
FuncRandomIntegerMT(Obj self,Obj mtstr,Obj nrbits)2643 static Obj FuncRandomIntegerMT(Obj self, Obj mtstr, Obj nrbits)
2644 {
2645 Obj res;
2646 Int i, n, q, r, qoff, len;
2647 UInt4 *mt;
2648 UInt4 *pt;
2649 RequireStringRep("RandomIntegerMT", mtstr);
2650 if (GET_LEN_STRING(mtstr) < 2500) {
2651 ErrorMayQuit(
2652 "RandomIntegerMT: <mtstr> must be a string with at least 2500 characters",
2653 0, 0);
2654 }
2655 RequireNonnegativeSmallInt("RandomIntegerMT", nrbits);
2656 n = INT_INTOBJ(nrbits);
2657
2658 /* small int case */
2659 if (n <= NR_SMALL_INT_BITS) {
2660 mt = (UInt4 *)(ADDR_OBJ(mtstr) + 1);
2661 #ifdef SYS_IS_64_BIT
2662 if (n <= 32) {
2663 res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4) -1L >> (32-n))));
2664 }
2665 else {
2666 unsigned long rd;
2667 rd = nextrandMT_int32(mt);
2668 rd += (unsigned long) ((UInt4) nextrandMT_int32(mt) &
2669 ((UInt4) -1L >> (64-n))) << 32;
2670 res = INTOBJ_INT((Int)rd);
2671 }
2672 #else
2673 res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4) -1L >> (32-n))));
2674 #endif
2675 }
2676 else {
2677 /* large int case */
2678 q = n / 32;
2679 r = n - q * 32;
2680 /* qoff = number of 32 bit words we need */
2681 qoff = q + (r==0 ? 0:1);
2682 /* len = number of limbs we need (limbs currently are either 32 or 64 bit wide) */
2683 len = (qoff*4 + sizeof(mp_limb_t) - 1) / sizeof(mp_limb_t);
2684 res = NewBag( T_INTPOS, len*sizeof(mp_limb_t) );
2685 pt = (UInt4*) ADDR_INT(res);
2686 mt = (UInt4 *)(ADDR_OBJ(mtstr) + 1);
2687 for (i = 0; i < qoff; i++, pt++) {
2688 *pt = nextrandMT_int32(mt);
2689 }
2690 if (r != 0) {
2691 /* we generated too many random bits -- chop of the extra bits */
2692 pt = (UInt4*) ADDR_INT(res);
2693 pt[qoff-1] = pt[qoff-1] & ((UInt4)(-1) >> (32-r));
2694 }
2695 #if defined(SYS_IS_64_BIT) && defined(WORDS_BIGENDIAN)
2696 // swap the halves of the 64bit words to match the
2697 // little endian resp. 32 bit versions of this code
2698 pt = (UInt4 *)ADDR_INT(res);
2699 for (i = 0; i < qoff; i += 2, pt += 2) {
2700 SWAP(UInt4, pt[0], pt[1]);
2701 }
2702 #endif
2703 /* shrink bag if necessary */
2704 res = GMP_NORMALIZE(res);
2705 /* convert result if small int */
2706 res = GMP_REDUCE(res);
2707 }
2708
2709 return res;
2710 }
2711
2712 /****************************************************************************
2713 **
2714 ** The following functions only exist to enable use to test the conversion
2715 ** functions (Int_ObjInt, ObjInt_Int, etc.) via a regular .tst file.
2716 */
FuncINTERNAL_TEST_CONV_INT(Obj self,Obj val)2717 static Obj FuncINTERNAL_TEST_CONV_INT(Obj self, Obj val)
2718 {
2719 Int ival = Int_ObjInt(val);
2720 return ObjInt_Int(ival);
2721 }
2722
FuncINTERNAL_TEST_CONV_UINT(Obj self,Obj val)2723 static Obj FuncINTERNAL_TEST_CONV_UINT(Obj self, Obj val)
2724 {
2725 UInt ival = UInt_ObjInt(val);
2726 return ObjInt_UInt(ival);
2727 }
2728
FuncINTERNAL_TEST_CONV_UINTINV(Obj self,Obj val)2729 static Obj FuncINTERNAL_TEST_CONV_UINTINV(Obj self, Obj val)
2730 {
2731 UInt ival = UInt_ObjInt(val);
2732 return ObjInt_UIntInv(ival);
2733 }
2734
FuncINTERNAL_TEST_CONV_INT8(Obj self,Obj val)2735 static Obj FuncINTERNAL_TEST_CONV_INT8(Obj self, Obj val)
2736 {
2737 Int8 ival = Int8_ObjInt(val);
2738 return ObjInt_Int8(ival);
2739 }
2740
FuncINTERNAL_TEST_CONV_UINT8(Obj self,Obj val)2741 static Obj FuncINTERNAL_TEST_CONV_UINT8(Obj self, Obj val)
2742 {
2743 UInt8 ival = UInt8_ObjInt(val);
2744 return ObjInt_UInt8(ival);
2745 }
2746
2747
2748 /****************************************************************************
2749 **
2750 *F * * * * * * * * * * * * * initialize module * * * * * * * * * * * * * * *
2751 */
2752
2753
2754 /****************************************************************************
2755 **
2756 *V BagNames . . . . . . . . . . . . . . . . . . . . . . . list of bag names
2757 */
2758 static StructBagNames BagNames[] = {
2759 { T_INT, "integer" },
2760 { T_INTPOS, "large positive integer" },
2761 { T_INTNEG, "large negative integer" },
2762 { -1, "" }
2763 };
2764
2765
2766 /****************************************************************************
2767 **
2768 *V GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export
2769 */
2770 static StructGVarFilt GVarFilts [] = {
2771
2772 GVAR_FILT(IS_INT, "obj", &IsIntFilt),
2773 { 0, 0, 0, 0, 0 }
2774
2775 };
2776
2777
2778 /****************************************************************************
2779 **
2780 *V GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
2781 */
2782 static StructGVarFunc GVarFuncs[] = {
2783
2784 GVAR_FUNC(QUO_INT, 2, "a, b"),
2785 GVAR_FUNC(ABS_INT, 1, "n"),
2786 GVAR_FUNC(SIGN_INT, 1, "n"),
2787 GVAR_FUNC(REM_INT, 2, "a, b"),
2788 GVAR_FUNC(GCD_INT, 2, "a, b"),
2789 GVAR_FUNC(LCM_INT, 2, "a, b"),
2790 GVAR_FUNC(PROD_INT_OBJ, 2, "opL, opR"),
2791 GVAR_FUNC(POW_OBJ_INT, 2, "opL, opR"),
2792 GVAR_FUNC(JACOBI_INT, 2, "n, m"),
2793 GVAR_FUNC(FACTORIAL_INT, 1, "n"),
2794 GVAR_FUNC(BINOMIAL_INT, 2, "n, k"),
2795 GVAR_FUNC(PVALUATION_INT, 2, "n, p"),
2796 GVAR_FUNC(ROOT_INT, 2, "n, k"),
2797 GVAR_FUNC(POWERMODINT, 3, "base, exp, mod"),
2798 GVAR_FUNC(IS_PROBAB_PRIME_INT, 2, "n, reps"),
2799 GVAR_FUNC(INVMODINT, 2, "base, mod"),
2800 GVAR_FUNC(HexStringInt, 1, "n"),
2801 GVAR_FUNC(IntHexString, 1, "string"),
2802 GVAR_FUNC(Log2Int, 1, "n"),
2803 GVAR_FUNC(STRING_INT, 1, "n"),
2804 GVAR_FUNC(INT_STRING, 1, "string"),
2805 GVAR_FUNC(RandomIntegerMT, 2, "mtstr, nrbits"),
2806
2807 GVAR_FUNC(INTERNAL_TEST_CONV_INT, 1, "val"),
2808 GVAR_FUNC(INTERNAL_TEST_CONV_UINT, 1, "val"),
2809 GVAR_FUNC(INTERNAL_TEST_CONV_UINTINV, 1, "val"),
2810 GVAR_FUNC(INTERNAL_TEST_CONV_INT8, 1, "val"),
2811 GVAR_FUNC(INTERNAL_TEST_CONV_UINT8, 1, "val"),
2812
2813 { 0, 0, 0, 0, 0 }
2814
2815 };
2816
2817
2818 /****************************************************************************
2819 **
2820 *F InitKernel( <module> ) . . . . . . . . initialise kernel data structures
2821 */
InitKernel(StructInitInfo * module)2822 static Int InitKernel ( StructInitInfo * module )
2823 {
2824 UInt t1, t2;
2825
2826 if (mp_bits_per_limb != GMP_LIMB_BITS) {
2827 Panic("GMP limb size mismatch");
2828 }
2829 if (INTOBJ_MIN != INTOBJ_INT(INT_INTOBJ_MIN)) {
2830 Panic("INTOBJ_MIN mismatch");
2831 }
2832 if (INTOBJ_MAX != INTOBJ_INT(INT_INTOBJ_MAX)) {
2833 Panic("INTOBJ_MAX mismatch");
2834 }
2835
2836 /* init filters and functions */
2837 InitHdlrFiltsFromTable( GVarFilts );
2838 InitHdlrFuncsFromTable( GVarFuncs );
2839
2840 // set the bag type names (for error messages and debugging)
2841 InitBagNamesFromTable( BagNames );
2842
2843 /* install the marking functions */
2844 InitMarkFuncBags( T_INTPOS, MarkNoSubBags );
2845 InitMarkFuncBags( T_INTNEG, MarkNoSubBags );
2846
2847 /* Install the saving methods */
2848 SaveObjFuncs [ T_INTPOS ] = SaveInt;
2849 SaveObjFuncs [ T_INTNEG ] = SaveInt;
2850 LoadObjFuncs [ T_INTPOS ] = LoadInt;
2851 LoadObjFuncs [ T_INTNEG ] = LoadInt;
2852
2853 /* install the printing functions */
2854 PrintObjFuncs[ T_INT ] = PrintInt;
2855 PrintObjFuncs[ T_INTPOS ] = PrintInt;
2856 PrintObjFuncs[ T_INTNEG ] = PrintInt;
2857
2858 /* install the comparison methods */
2859 for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
2860 for ( t2 = T_INT; t2 <= T_INTNEG; t2++ ) {
2861 EqFuncs [ t1 ][ t2 ] = EqInt;
2862 LtFuncs [ t1 ][ t2 ] = LtInt;
2863 }
2864 }
2865
2866 /* install the unary arithmetic methods */
2867 for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
2868 ZeroFuncs[ t1 ] = ZeroInt;
2869 ZeroMutFuncs[ t1 ] = ZeroInt;
2870 AInvFuncs[ t1 ] = AInvInt;
2871 AInvMutFuncs[ t1 ] = AInvInt;
2872 OneFuncs [ t1 ] = OneInt;
2873 OneMutFuncs [ t1 ] = OneInt;
2874 }
2875
2876 /* install the default power methods */
2877 for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
2878 for ( t2 = FIRST_MULT_TNUM; t2 <= LAST_MULT_TNUM; t2++ ) {
2879 PowFuncs [ t2 ][ t1 ] = PowObjInt;
2880 }
2881 for ( t2 = FIRST_PLIST_TNUM; t2 <= LAST_PLIST_TNUM; t2++ ) {
2882 PowFuncs [ t2 ][ t1 ] = PowObjInt;
2883 }
2884 PowFuncs [ T_RANGE_NSORT ][ t1 ] = PowObjInt;
2885 PowFuncs [ T_RANGE_SSORT ][ t1 ] = PowObjInt;
2886 }
2887
2888 /* install the binary arithmetic methods */
2889 for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
2890 for ( t2 = T_INT; t2 <= T_INTNEG; t2++ ) {
2891 EqFuncs [ t1 ][ t2 ] = EqInt;
2892 LtFuncs [ t1 ][ t2 ] = LtInt;
2893 SumFuncs [ t1 ][ t2 ] = SumInt;
2894 DiffFuncs[ t1 ][ t2 ] = DiffInt;
2895 ProdFuncs[ t1 ][ t2 ] = ProdInt;
2896 PowFuncs [ t1 ][ t2 ] = PowInt;
2897 ModFuncs [ t1 ][ t2 ] = ModInt;
2898 }
2899 }
2900
2901 /* gvars to import from the library */
2902 ImportGVarFromLibrary( "TYPE_INT_SMALL_ZERO", &TYPE_INT_SMALL_ZERO );
2903 ImportGVarFromLibrary( "TYPE_INT_SMALL_POS", &TYPE_INT_SMALL_POS );
2904 ImportGVarFromLibrary( "TYPE_INT_SMALL_NEG", &TYPE_INT_SMALL_NEG );
2905 ImportGVarFromLibrary( "TYPE_INT_LARGE_POS", &TYPE_INT_LARGE_POS );
2906 ImportGVarFromLibrary( "TYPE_INT_LARGE_NEG", &TYPE_INT_LARGE_NEG );
2907
2908 ImportFuncFromLibrary( "String", &String );
2909 ImportFuncFromLibrary( "One", &OneAttr);
2910
2911 /* install the type functions */
2912 TypeObjFuncs[ T_INT ] = TypeIntSmall;
2913 TypeObjFuncs[ T_INTPOS ] = TypeIntLargePos;
2914 TypeObjFuncs[ T_INTNEG ] = TypeIntLargeNeg;
2915
2916 #ifdef HPCGAP
2917 MakeBagTypePublic( T_INTPOS );
2918 MakeBagTypePublic( T_INTNEG );
2919 #endif
2920
2921 /* return success */
2922 return 0;
2923 }
2924
2925
2926 /****************************************************************************
2927 **
2928 *F InitLibrary( <module> ) . . . . . . . initialise library data structures
2929 */
InitLibrary(StructInitInfo * module)2930 static Int InitLibrary ( StructInitInfo * module )
2931 {
2932 /* init filters and functions */
2933 InitGVarFiltsFromTable( GVarFilts );
2934 InitGVarFuncsFromTable( GVarFuncs );
2935
2936 /* return success */
2937 return 0;
2938 }
2939
2940
2941 /****************************************************************************
2942 **
2943 *F InitInfoInt() . . . . . . . . . . . . . . . . . . table of init functions
2944 */
2945 static StructInitInfo module = {
2946 // init struct using C99 designated initializers; for a full list of
2947 // fields, please refer to the definition of StructInitInfo
2948 .type = MODULE_BUILTIN,
2949 .name = "integer",
2950 .initKernel = InitKernel,
2951 .initLibrary = InitLibrary,
2952 };
2953
InitInfoInt(void)2954 StructInitInfo * InitInfoInt ( void )
2955 {
2956 return &module;
2957 }
2958
2959 #endif /* ! WARD_ENABLED */
2960