CoCalc -- integer.c

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

GAP 4.8.9 installation with standard packages -- copy to your CoCalc project to get it
Project: cocalc-sagemath-dev-slelievre
Views: ⁴¹⁸³⁴⁶
1
/****************************************************************************
2
**
3
*W  integer.c                   GAP source                   Martin Schönert
4
**                                                           & Alice Niemeyer
5
**                                                           & Werner  Nickel
6
**
7
**
8
*Y  Copyright (C)  1996,  Lehrstuhl D für Mathematik,  RWTH Aachen,  Germany
9
*Y  (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland
10
*Y  Copyright (C) 2002 The GAP Group
11
**
12
**  This file implements the  functions  handling  arbitrary  size  integers.
13
**
14
**  There are three integer types in GAP: 'T_INT', 'T_INTPOS' and 'T_INTNEG'.
15
**  Each integer has a unique representation, e.g., an integer  that  can  be
16
**  represented as 'T_INT' is never  represented as 'T_INTPOS' or 'T_INTNEG'.
17
**
18
**  'T_INT' is the type of those integers small enough to fit into  29  bits.
19
**  Therefore the value range of this small integers is $-2^{28}...2^{28}-1$.
20
**  This range contains about 99\% of all integers that usually occur in GAP.
21
**  (I just made up this number, obviously it depends on the application  :-)
22
**  Only these small integers can be used as index expression into sequences.
23
**
24
**  Small integers are represented by an immediate integer handle, containing
25
**  the value instead of pointing  to  it,  which  has  the  following  form:
26
**
27
**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
28
**      | guard | sign  | bit   | bit   |       | bit   | tag   | tag   |
29
**      | bit   | bit   | 27    | 26    |       | 0     | 0     | 1     |
30
**      +-------+-------+-------+-------+- - - -+-------+-------+-------+
31
**
32
**  Immediate integers handles carry the tag 'T_INT', i.e. the last bit is 1.
33
**  This distinguishes immediate integers from other handles which  point  to
34
**  structures aligned on 4 byte boundaries and therefore have last bit zero.
35
**  (The second bit is reserved as tag to allow extensions of  this  scheme.)
36
**  Using immediates as pointers and dereferencing them gives address errors.
37
**
38
**  To aid overflow check the most significant two bits must always be equal,
39
**  that is to say that the sign bit of immediate integers has a  guard  bit.
40
**
41
**  The macros 'INTOBJ_INT' and 'INT_INTOBJ' should be used to convert between
42
**  a small integer value and its representation as immediate integer handle.
43
**
44
**  'T_INTPOS' and 'T_INTNEG' are the types of positive (respectively, negative)
45
**  integer values  that  can  not  be  represented  by  immediate  integers.
46
**
47
**  This large integers values are represented in signed base 65536 notation.
48
**  That means that the bag of  a  large  integer  has  the  following  form:
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 65536 + d2 65536^2 + ... + d_n 65536^n$,
56
**  respectively the negative of this if the type of this object is 'T_INTNEG'.
57
**
58
**  Each digit is  of  course  stored  as  a  16  bit  wide  unsigned  short.
59
**  Note that base 65536 allows us to multiply 2 digits and add a carry digit
60
**  without overflow in 32 bit long arithmetic, available on most processors.
61
**
62
**  The number of digits in every  large  integer  is  a  multiple  of  four.
63
**  Therefore the leading three digits of some values will actually be  zero.
64
**  Note that the uniqueness of representation implies that not four or more
65
**  leading digits may be zero, since |d0|d1|d2|d3| and |d0|d1|d2|d3|0|0|0|0|
66
**  have the same value only one, the first, can be a  legal  representation.
67
**
68
**  Because of this it is possible to do a  little  bit  of  loop  unrolling.
69
**  Thus instead of looping <n> times, handling one digit in each  iteration,
70
**  we can loop <n>/4 times, handling  four  digits  during  each  iteration.
71
**  This reduces the overhead of the loop by a factor of  approximately four.
72
**
73
**  Using base 65536 representation has advantages over  using  other  bases.
74
**  Integers in base 65536 representation can be packed  dense and therefore
75
**  use roughly 20\% less space than integers in base  10000  representation.
76
**  'SumInt' is 20\% and 'ProdInt' is 40\% faster for 65536 than  for  10000,
77
**  as their runtime is linear respectively quadratic in the number of digits.
78
**  Dividing by 65536 and computing the remainder mod 65536 can be done  fast
79
**  by shifting 16 bit to  the  right  and  by  taking  the  lower  16  bits.
80
**  Larger bases are difficult because the product of two digits will not fit
81
**  into 32 bit, which is the word size  of  most  modern  micro  processors.
82
**  Base 10000 would have the advantage that printing is  very  much  easier,
83
**  but 'PrInt' keeps a terminal at 9600 baud busy for almost  all  integers.
84
*/
85

86
#include        "system.h"              /* Ints, UInts                     */
87

88
#ifndef USE_GMP /* use this file, otherwise ignore what follows */
89

90
#include        "gasman.h"              /* garbage collector               */
91
#include        "objects.h"             /* objects                         */
92
#include        "scanner.h"             /* scanner                         */
93

94
#include        "gvars.h"               /* global variables                */
95

96
#include        "calls.h"               /* generic call mechanism          */
97
#include        "opers.h"               /* generic operations              */
98

99
#include        "ariths.h"              /* basic arithmetic                */
100

101
#include        "bool.h"                /* booleans                        */
102

103
#include        "integer.h"             /* integers                        */
104

105
#include        "gap.h"                 /* error handling, initialisation  */
106

107
#include        "records.h"             /* generic records                 */
108
#include        "precord.h"             /* plain records                   */
109

110
#include        "lists.h"               /* generic lists                   */
111
#include        "string.h"              /* strings                         */
112

113
#include        "saveload.h"            /* saving and loading              */
114

115
#include        "intfuncs.h"
116

117
#include	"code.h"		/* coder                           */
118
#include	"thread.h"		/* threads			   */
119
#include	"tls.h"			/* thread-local storage		   */
120

121
#include <stdio.h>
122

123
/* for fallbacks to library */
124
Obj String;
125

126

127
/****************************************************************************
128
**
129
*F  TypeInt(<int>)  . . . . . . . . . . . . . . . . . . . . . type of integer
130
**
131
**  'TypeInt' returns the type of the integer <int>.
132
**
133
**  'TypeInt' is the function in 'TypeObjFuncs' for integers.
134
*/
135
Obj             TYPE_INT_SMALL_ZERO;
136
Obj             TYPE_INT_SMALL_POS;
137
Obj             TYPE_INT_SMALL_NEG;
138
Obj             TYPE_INT_LARGE_POS;
139
Obj             TYPE_INT_LARGE_NEG;
140

141
Obj             TypeIntSmall (
142
    Obj                 val )
143
{
144
    if ( 0 == INT_INTOBJ(val) ) {
145
        return TYPE_INT_SMALL_ZERO;
146
    }
147
    else if ( 0 < INT_INTOBJ(val) ) {
148
        return TYPE_INT_SMALL_POS;
149
    }
150
    else /* if ( 0 > INT_INTOBJ(val) ) */ {
151
        return TYPE_INT_SMALL_NEG;
152
    }
153
}
154

155
Obj             TypeIntLargePos (
156
    Obj                 val )
157
{
158
    return TYPE_INT_LARGE_POS;
159
}
160

161
Obj             TypeIntLargeNeg (
162
    Obj                 val )
163
{
164
    return TYPE_INT_LARGE_NEG;
165
}
166

167
/**************************************************************************
168
** The following two functions convert a C Int or UInt respectively into
169
** a GAP integer, either an immediate, small integer if possible or
170
** otherwise a new GAP bag with TNUM T_INTPOS or T_INTNEG.
171
**
172
*F ObjInt_Int(Int i)
173
*F ObjInt_UInt(UInt i)
174
**
175
****************************************************************************/
176

177
Obj ObjInt_Int(Int i)
178
{
179
    Obj n;
180
    Int bound = 1L << NR_SMALL_INT_BITS;
181
    if (i >= bound) {
182
        /* We have to make a big integer */
183
        n = NewBag(T_INTPOS,4*sizeof(TypDigit));
184
        ADDR_INT(n)[0] = (TypDigit) (i & ((Int) INTBASE - 1L));
185
        ADDR_INT(n)[1] = (TypDigit) (i >> NR_DIGIT_BITS);
186
        ADDR_INT(n)[2] = 0;
187
        ADDR_INT(n)[3] = 0;
188
        return n;
189
    } else if (-i > bound) {
190
        n = NewBag(T_INTNEG,4*sizeof(TypDigit));
191
        ADDR_INT(n)[0] = (TypDigit) ((-i) & ((Int) INTBASE - 1L));
192
        ADDR_INT(n)[1] = (TypDigit) ((-i) >> NR_DIGIT_BITS);
193
        ADDR_INT(n)[2] = 0;
194
        ADDR_INT(n)[3] = 0;
195
        return n;
196
    } else {
197
        return INTOBJ_INT(i);
198
    }
199
}
200

201
Obj ObjInt_UInt(UInt i)
202
{
203
    Obj n;
204
    UInt bound = 1UL << NR_SMALL_INT_BITS;
205
    if (i >= bound) {
206
        /* We have to make a big integer */
207
        n = NewBag(T_INTPOS,4*sizeof(TypDigit));
208
        ADDR_INT(n)[0] = (TypDigit) (i & ((UInt) INTBASE - 1L));
209
        ADDR_INT(n)[1] = (TypDigit) (i >> NR_DIGIT_BITS);
210
        ADDR_INT(n)[2] = 0;
211
        ADDR_INT(n)[3] = 0;
212
        return n;
213
    } else {
214
        return INTOBJ_INT(i);
215
    }
216
}
217

218

219

220
/****************************************************************************
221
**
222
*F  PrintInt( <int> ) . . . . . . . . . . . . . . . print an integer constant
223
**
224
**  'PrintInt'  prints  the integer  <int>   in the  usual  decimal notation.
225
**  'PrintInt' handles objects of type 'T_INT', 'T_INTPOS' and 'T_INTNEG'.
226
**
227
**  Large integers are first converted into  base  10000  and  then  printed.
228
**  The time for a conversion depends quadratically on the number of  digits.
229
**  For 2000 decimal digit integers, a screenfull,  it  is  reasonable  fast.
230
**
231
**  The number  of digits  needed in PrIntD[] is the ceiling of the logarithm
232
**  with respect to base PRINT_BASE of
233
**
234
**           ( (1<<NR_DIGIT_BITS) )^1000 - 1.
235
**
236
**  The latter is the largest number that can be represented with 1000 digits
237
**  of type TypDigit.
238
**
239
**  If NR_DIGIT_BITS is 16, we get 1205.
240
**  If NR_DIGIT_BITS is 32, we get 1071.
241
**
242
**  The subsidiary function IntToPrintBase converts an integer into base
243
**  PRINT_BASE, leaving the result in base PrIntD. It returns the index of the
244
**  most significant digits. It assumes that the argument is a large
245
**  integer small enough to fit.
246
*/
247

248
TypDigit        PrIntC [1000];          /* copy of integer to be printed   */
249

250
#ifdef SYS_IS_64_BIT
251

252
#define PRINT_BASE 1000000000L          /* 10^9                            */
253
#define PRINT_FORMAT "%09d"             /* print 9 decimals at a time      */
254
#define CHARS_PER_PRINT_BASE 9
255
TypDigit        PrIntD [1071];          /* integer converted to base 10^9  */
256
#define NR_HEX_DIGITS 8
257

258
#else
259

260
#define PRINT_BASE 10000
261
#define PRINT_FORMAT "%04d"             /* print 4 decimals at a time      */
262
#define CHARS_PER_PRINT_BASE 4
263
TypDigit        PrIntD [1205];          /* integer converted to base 10000 */
264
#define NR_HEX_DIGITS 4
265

266
#endif
267
/****************************************************************************
268
**
269
*F  FuncHexStringInt( <self>, <int> ) . . . . . . . . hex string for integer
270
*F  FuncIntHexString( <self>, <string> ) . . . . . .  integer from hex string
271
**
272
**  The  function  `FuncHexStringInt'  constructs from  an  integer  the
273
**  corresponding string in  hexadecimal notation. It has  a leading '-'
274
**  for negative numbers and the digits 10..15 are written as A..F.
275
**
276
**  The  function `FuncIntHexString'  does  the converse,  but here  the
277
**  letters a..f are also allowed in <string> instead of A..F.
278
**
279
*/
280
Obj FuncHexStringInt( Obj self, Obj integer )
281
{
282
     Int len, i, j, n;
283
     UInt nf;
284
     TypDigit d, f;
285
     UInt1 *p, a;
286
     Obj res;
287

288
     /* immediate integers */
289
     if (IS_INTOBJ(integer)) {
290
         n = INT_INTOBJ(integer);
291
         /* 0 is special */
292
         if (n == 0) {
293
             res = NEW_STRING(1);
294
             CHARS_STRING(res)[0] = '0';
295
             return res;
296
         }
297

298
         /* else we create a string big enough for any immediate integer */
299
         res = NEW_STRING(2 * NR_HEX_DIGITS + 1);
300
         p = CHARS_STRING(res);
301
         /* handle sign */
302
         if (n<0) {
303
            p[0] = '-';
304
            n = -n;
305
            p++;
306
         }
307
         else
308
            SET_LEN_STRING(res, GET_LEN_STRING(res)-1);
309
         /* collect digits, skipping leading zeros */
310
         j = 0;
311
         nf = ((UInt)15) << (4*(2*NR_HEX_DIGITS-1));
312
         for (i = 2*NR_HEX_DIGITS; i; i-- ) {
313
             a = ((UInt)n & nf) >> (4*(i-1));
314
             if (j==0 && a==0) SET_LEN_STRING(res, GET_LEN_STRING(res)-1);
315
             else if (a<10) p[j++] = a + '0';
316
             else p[j++] = a - 10 + 'A';
317
             nf = nf >> 4;
318
         }
319
         /* final null character */
320
         p[j] = 0;
321
         return res;
322
     }
323
     else if (TNUM_OBJ(integer) == T_INTNEG || TNUM_OBJ(integer) == T_INTPOS) {
324
         /* nr of digits */
325
         len = SIZE_INT(integer);
326
         for (; ADDR_INT(integer)[len-1] == 0; len--);
327

328
         /* result string and sign */
329
         if (TNUM_OBJ(integer) == T_INTNEG) {
330
             res = NEW_STRING(len * NR_HEX_DIGITS + 1);
331
             p = CHARS_STRING(res);
332
             p[0] = '-';
333
             p++;
334
         }
335
         else {
336
             res = NEW_STRING(len * NR_HEX_DIGITS);
337
             p = CHARS_STRING(res);
338
         }
339
         /* collect digits */
340
         j = 0;
341
         for (; len; len--) {
342
             d = ADDR_INT(integer)[len-1];
343
             f = 15L << (4*(NR_HEX_DIGITS-1));
344
             for (i = NR_HEX_DIGITS; i; i-- ) {
345
                 a = (d & f) >> (4*(i-1));
346
                 if (j==0 && a==0) SET_LEN_STRING(res, GET_LEN_STRING(res)-1);
347
                 else if (a<10) p[j++] = a + '0';
348
                 else p[j++] = a - 10 + 'A';
349
                 f = f >> 4;
350
             }
351
         }
352
         /* final null character */
353
         p[j] = 0;
354
         return res;
355
     }
356
     else
357
         ErrorReturnObj("HexStringInt: argument must be integer, (not a %s)",
358
           (Int)TNAM_OBJ(integer), 0L,
359
           "");
360
     return (Obj) 0L; /* please picky cc */
361
}
362

363
Obj  FuncIntHexString( Obj self,  Obj str )
364
{
365
    Obj res;
366
    Int  i, j, s, ii, len, sign, nd;
367
    UInt n;
368
    UInt1 *p, a;
369
    TypDigit d;
370

371
    if (! IsStringConv(str))
372
        ErrorReturnObj("IntHexString: argument must be string (not a %s)",
373
          (Int)TNAM_OBJ(str), 0L,
374
          "");
375

376
    /* number of hex digits and sign */
377
    len = GET_LEN_STRING(str);
378
    if (len == 0) {
379
       res = INTOBJ_INT(0);
380
       return res;
381
    }
382
    if (*(CHARS_STRING(str)) == '-') {
383
       sign = -1;
384
       i = 1;
385
    }
386
    else {
387
       sign = 1;
388
       i = 0;
389
    }
390
    
391
    /* skip leading zeros */
392
    while ((CHARS_STRING(str))[i] == '0' && i < len)
393
        i++;
394
    
395
    /* small int case */
396
    if ((len-i)*4 <= NR_SMALL_INT_BITS) {
397
       n = 0;
398
       p = CHARS_STRING(str);
399
       for (; i<len; i++) {
400
          a = p[i];
401
          if (a>96)
402
             a -= 87;
403
          else if (a>64)
404
             a -= 55;
405
          else
406
             a -= 48;
407
          if (a > 15)
408
              ErrorReturnObj("IntHexString: non-valid character in hex-string",
409
                0L, 0L, "");
410
          n = (n << 4) + a;
411
       }
412
       res = INTOBJ_INT(sign * n);
413
       return res;
414
    }
415
    else {
416
       /* number of Digits */
417
       nd = (len-i)/NR_HEX_DIGITS;
418
       if (nd * NR_HEX_DIGITS < (len-i)) nd++;
419
       nd += ((3*nd) % 4);
420
       if (sign == 1)
421
          res = NewBag( T_INTPOS, nd*sizeof(TypDigit) );
422
       else
423
          res = NewBag( T_INTNEG, nd*sizeof(TypDigit) );
424
       /* collect digits, easiest to start from the end */
425
       p = CHARS_STRING(str);
426
       for (j=0; j < nd; j++) {
427
          d = 0;
428
          for (s=0, ii=len-j*NR_HEX_DIGITS-1;
429
               ii>=i && ii>len-(j+1)*NR_HEX_DIGITS-1;
430
               s+=4, ii--) {
431
             a = p[ii];
432
             if (a>96)
433
                a -= 87;
434
             else if (a>64)
435
                a -= 55;
436
             else
437
                a -= 48;
438
             if (a > 15)
439
               ErrorReturnObj("IntHexString: non-valid character in hex-string",
440
                   0L, 0L, "");
441

442
             d += (a<<s);
443
          }
444
          ADDR_INT(res)[j] = d;
445
       }
446
       return res;
447
    }
448
}
449

450

451

452
Int IntToPrintBase ( Obj op )
453
{
454
    UInt                 i, k;           /* loop counter                    */
455
    TypDigit *          p;              /* loop pointer                    */
456
    UInt                c;              /* carry in division step          */
457

458
    i = 0;
459
    for ( k = 0; k < SIZE_INT(op); k++ )
460
      PrIntC[k] = ADDR_INT(op)[k];
461
    while ( k > 0 && PrIntC[k-1] == 0 )  k--;
462
    while ( k > 0 ) {
463
      for ( c = 0, p = PrIntC+k-1; p >= PrIntC; p-- ) {
464
	c  = (c<<NR_DIGIT_BITS) + *p;
465
	*p = (TypDigit)(c / PRINT_BASE);
466
	c  = c - PRINT_BASE * *p;
467
      }
468
      PrIntD[i++] = (TypDigit)c;
469
      while ( k > 0 && PrIntC[k-1] == 0 )  k--;
470
    }
471
    return i-1;
472

473
}
474

475
void            PrintInt (
476
    Obj                 op )
477
{
478
    Int                 i;           /* loop counter                    */
479
    Obj               str;           /* fallback to lib for large ints  */
480

481
    /* print a small integer                                               */
482
    if ( IS_INTOBJ(op) ) {
483
        Pr( "%>%d%<", INT_INTOBJ(op), 0L );
484
    }
485

486
    /* print a large integer                                               */
487
    else if ( SIZE_INT(op) < 1000 ) {
488

489
        /* start printing, %> means insert '\' before a linebreak          */
490
        Pr("%>",0L,0L);
491

492
        if ( TNUM_OBJ(op) == T_INTNEG )
493
            Pr("-",0L,0L);
494

495
        /* convert the integer into base PRINT_BASE                        */
496
	i = IntToPrintBase(op);
497

498
        /* print the base PRINT_BASE digits                                 */
499
        Pr( "%d", (Int)PrIntD[i], 0L );
500
        while ( i > 0 )
501
            Pr( PRINT_FORMAT, (Int)PrIntD[--i], 0L );
502
        Pr("%<",0L,0L);
503

504
    }
505

506
    else {
507
        str = CALL_1ARGS( String, op );
508
        Pr("%>%s%<",(Int)(CHARS_STRING(str)), 0);
509
        /* for a long time Print of large ints did not follow the general idea
510
         * that Print should produce something that can be read back into GAP:
511
           Pr("<<an integer too large to be printed>>",0L,0L); */
512
    }
513
}
514

515
/****************************************************************************
516
**  
517
**  Implementation of Log2Int for C integers.
518
*/
519

520
Int CLog2Int(Int a)
521
{
522
  Int res, mask;
523
  if (a < 0) a = -a;
524
  if (a < 1) return -1;
525
  if (a < 65536) {
526
    for(mask = 2, res = 0; ;mask *= 2, res += 1) {
527
      if(a < mask) return res;
528
    }
529
  }
530
  for(mask = 65536, res = 15; ;mask *= 2, res += 1) {
531
    if(a < mask) return res;
532
  }
533
}
534

535
/****************************************************************************
536
**
537
*F  FuncLog2Int( <self>, <int> ) . . . . . . . . . nr of bits of integer - 1
538
**
539
**  Given to GAP-Level as "Log2Int".
540
*/
541
Obj FuncLog2Int( Obj self, Obj integer)
542
{
543
  Int d;
544
  Int a, len;
545
  TypDigit dmask;
546

547
  /* case of small ints */
548
  if (IS_INTOBJ(integer)) {
549
    return INTOBJ_INT(CLog2Int(INT_INTOBJ(integer)));
550
  }
551

552
  /* case of long ints */
553
  if (TNUM_OBJ(integer) == T_INTNEG || TNUM_OBJ(integer) == T_INTPOS) {
554
    for (len = SIZE_INT(integer); ADDR_INT(integer)[len-1] == 0; len--);
555
    /* Instead of computing 
556
          res = len * NR_DIGIT_BITS - d;
557
       we keep len and d separate, because on 32 bit systems res may
558
       not fit into an Int (and not into an immediate integer).            */
559
    d = 1;
560
    a = (TypDigit)(ADDR_INT(integer)[len-1]);
561
    for(dmask = (TypDigit)1 << (NR_DIGIT_BITS - 1);
562
        (dmask & a) == 0 && dmask != (TypDigit)0;
563
        dmask = dmask >> 1, d++);
564
    return DiffInt(ProdInt(INTOBJ_INT(len), INTOBJ_INT(NR_DIGIT_BITS)), 
565
                   INTOBJ_INT(d));
566
  }
567
  else {
568
    ErrorReturnObj("Log2Int: argument must be integer, (not a %s)",
569
           (Int)TNAM_OBJ(integer), 0L,
570
           "");
571
    return (Obj) 0L; /* please picky cc */
572
  }
573
}
574

575
/****************************************************************************
576
**
577
*F  FuncSTRING_INT( <self>, <int> ) . . . . .  convert an integer to a string
578
**
579
**  `FuncSTRING_INT' returns an immutable string representing the integer
580
**  <int>
581
**
582
*/
583

584
Obj FuncSTRING_INT( Obj self, Obj integer )
585
{
586
  Int x;
587
  Obj str;
588
  Int len;
589
  Int i;
590
  Char c;
591
  Int j,top, chunk, neg;
592

593
  /* handle a small integer                                               */
594
  if ( IS_INTOBJ(integer) ) {
595
    x = INT_INTOBJ(integer);
596
    str = NEW_STRING( (NR_SMALL_INT_BITS+5)/3 );
597
    RetypeBag(str, T_STRING+IMMUTABLE);
598
    len = 0;
599
    /* Case of zero */
600
    if (x == 0)
601
      {
602
	CHARS_STRING(str)[0] = '0';
603
	CHARS_STRING(str)[1] = '\0';
604
	ResizeBag(str, SIZEBAG_STRINGLEN(1));
605
	SET_LEN_STRING(str, 1);
606

607
	return str;
608
      }
609
    /* Negative numbers */
610
    if (x < 0)
611
      {
612
	CHARS_STRING(str)[len++] = '-';
613
	x = -x;
614
	neg = 1;
615
      }
616
    else
617
      neg = 0;
618

619
    /* Now the main case */
620
    while (x != 0)
621
      {
622
	CHARS_STRING(str)[len++] = '0'+ x % 10;
623
	x /= 10;
624
      }
625
    CHARS_STRING(str)[len] = '\0';
626

627
    /* finally, reverse the digits in place */
628
    for (i = neg; i < (neg+len)/2; i++)
629
      {
630
	c = CHARS_STRING(str)[neg+len-1-i];
631
	CHARS_STRING(str)[neg+len-1-i] = CHARS_STRING(str)[i];
632
	CHARS_STRING(str)[i] = c;
633
      }
634

635
    ResizeBag(str, SIZEBAG_STRINGLEN(len));
636
    SET_LEN_STRING(str, len);
637
    return str;
638
  }
639

640
  /* handle a large integer                                               */
641
  else if ( SIZE_INT(integer) < 1000 ) {
642

643
    /* convert the integer into base PRINT_BASE                        */
644
    len = IntToPrintBase(integer);
645
    str = NEW_STRING(CHARS_PER_PRINT_BASE*(len+1)+2);
646
    RetypeBag(str, T_STRING+IMMUTABLE);
647

648
    /* sort out the length of the top group */
649
    j = 1;
650
    top = (Int)PrIntD[len];
651
    while ( top >= j)
652
      {
653
	j *= 10;
654
      }
655

656
    /* Start filling in the string */
657
    i = 0;
658
    if ( TNUM_OBJ(integer) == T_INTNEG ) {
659
      CHARS_STRING(str)[i++] = '-';
660
    }
661

662
    while (j > 1)
663
      {
664
	j /= 10;
665
        CHARS_STRING(str)[i++] = '0' + (top / j) % 10;
666
      }
667

668
    /* Now the rest of the base PRINT_BASE digits are easy */
669
    while( len > 0)
670
      {
671
	chunk = (Int)PrIntD[--len];
672
	j = PRINT_BASE/10;
673
	while (j > 0)
674
	  {
675
	    CHARS_STRING(str)[i++] = '0' + (chunk / j) % 10;
676
	    j /= 10;
677
	  }
678
      }
679

680
    CHARS_STRING(str)[i] = '\0';
681
    ResizeBag(str, SIZEBAG_STRINGLEN(i));
682
    SET_LEN_STRING(str, i);
683
    return str;
684
  }
685
  else {
686

687
      /* Very large integer, fall back on the GAP function */
688
      return CALL_1ARGS( String, integer);
689
  }
690
}
691

692

693

694
/****************************************************************************
695
**
696
*F  EqInt( <intL>, <intR> ) . . . . . . . . .  test if two integers are equal
697
**
698
**  'EqInt' returns 1  if  the two integer   arguments <intL> and  <intR> are
699
**  equal and 0 otherwise.
700
*/
701
Int             EqInt (
702
    Obj                 opL,
703
    Obj                 opR )
704
{
705
    Int                 k;              /* loop counter                    */
706
    TypDigit *          l;              /* pointer into the left operand   */
707
    TypDigit *          r;              /* pointer into the right operand  */
708

709
    /* compare two small integers                                          */
710
    if ( ARE_INTOBJS( opL, opR ) ) {
711
        if ( INT_INTOBJ(opL) == INT_INTOBJ(opR) )  return 1L;
712
        else                                       return 0L;
713
    }
714

715
    /* compare a small and a large integer                                 */
716
    else if ( IS_INTOBJ(opL) ) {
717
        return 0L;
718
    }
719
    else if ( IS_INTOBJ(opR) ) {
720
        return 0L;
721
    }
722

723
    /* compare two large integers                                          */
724
    else {
725

726
        /* compare the sign and size                                       */
727
        if ( TNUM_OBJ(opL) != TNUM_OBJ(opR)
728
          || SIZE_INT(opL) != SIZE_INT(opR) )
729
            return 0L;
730

731
        /* set up the pointers                                             */
732
        l = ADDR_INT(opL);
733
        r = ADDR_INT(opR);
734

735
        /* run through the digits, four at a time                          */
736
        for ( k = SIZE_INT(opL)/4-1; k >= 0; k-- ) {
737
            if ( *l++ != *r++ )  return 0L;
738
            if ( *l++ != *r++ )  return 0L;
739
            if ( *l++ != *r++ )  return 0L;
740
            if ( *l++ != *r++ )  return 0L;
741
        }
742

743
        /* no differences found, so they must be equal                     */
744
        return 1L;
745

746
    }
747
}
748

749

750
/****************************************************************************
751
**
752
*F  LtInt( <intL>, <intR> ) . . . . . test if an integer is less than another
753
**
754
**  'LtInt' returns 1 if the integer <intL> is strictly less than the integer
755
**  <intR> and 0 otherwise.
756
*/
757
Int             LtInt (
758
    Obj                 opL,
759
    Obj                 opR )
760
{
761
    Int                 k;              /* loop counter                    */
762
    TypDigit *          l;              /* pointer into the left operand   */
763
    TypDigit *          r;              /* pointer into the right operand  */
764

765
    /* compare two small integers                                          */
766
    if ( ARE_INTOBJS( opL, opR ) ) {
767
        if ( INT_INTOBJ(opL) <  INT_INTOBJ(opR) )  return 1L;
768
        else                                       return 0L;
769
    }
770

771
    /* compare a small and a large integer                                 */
772
    else if ( IS_INTOBJ(opL) ) {
773
        if ( TNUM_OBJ(opR) == T_INTPOS )  return 1L;
774
        else                              return 0L;
775
    }
776
    else if ( IS_INTOBJ(opR) ) {
777
        if ( TNUM_OBJ(opL) == T_INTPOS )  return 0L;
778
        else                              return 1L;
779
    }
780

781
    /* compare two large integers                                          */
782
    else {
783

784
        /* compare the sign and size                                       */
785
        if (      TNUM_OBJ(opL) == T_INTNEG
786
               && TNUM_OBJ(opR) == T_INTPOS )
787
            return 1L;
788
        else if ( TNUM_OBJ(opL) == T_INTPOS
789
               && TNUM_OBJ(opR) == T_INTNEG )
790
            return 0L;
791
        else if ( (TNUM_OBJ(opL) == T_INTPOS
792
                && SIZE_INT(opL) < SIZE_INT(opR))
793
               || (TNUM_OBJ(opL) == T_INTNEG
794
                && SIZE_INT(opL) > SIZE_INT(opR)) )
795
            return 1L;
796
        else if ( (TNUM_OBJ(opL) == T_INTPOS
797
                && SIZE_INT(opL) > SIZE_INT(opR))
798
               || (TNUM_OBJ(opL) == T_INTNEG
799
                && SIZE_INT(opL) < SIZE_INT(opR)) )
800
            return 0L;
801

802
        /* set up the pointers                                             */
803
        l = ADDR_INT(opL);
804
        r = ADDR_INT(opR);
805

806
        /* run through the digits, from the end downwards                  */
807
        for ( k = SIZE_INT(opL)-1; k >= 0; k-- ) {
808
            if ( l[k] != r[k] ) {
809
                if ( (TNUM_OBJ(opL) == T_INTPOS
810
                   && l[k] < r[k])
811
                  || (TNUM_OBJ(opL) == T_INTNEG
812
                   && l[k] > r[k]) )
813
                    return 1L;
814
                else
815
                    return 0L;
816
            }
817
        }
818

819
        /* no differences found, so they must be equal                     */
820
        return 0L;
821

822
    }
823
}
824

825

826
/****************************************************************************
827
**
828
*F  SumInt( <intL>, <intR> )  . . . . . . . . . . . . . . sum of two integers
829
**
830
**  'SumInt' returns the sum of the two integer arguments <intL> and  <intR>.
831
**  'SumInt' handles operands of type 'T_INT', 'T_INTPOS' and 'T_INTNEG'.
832
**
833
**  It can also be used in the cases that both operands  are  small  integers
834
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
835
**  This case is usually already handled in 'EvalSum' for a better  efficiency.
836
**
837
**  Is called from the 'EvalSum'  binop so both operands are already evaluated.
838
**
839
**  'SumInt' is a little bit difficult since there are 16  different cases to
840
**  handle, each operand can be positive or negative, small or large integer.
841
**  If the operands have opposite sign 'SumInt' calls 'DiffInt',  this  helps
842
**  reduce the total amount of code by a factor of two.
843
*/
844
Obj             SumInt (
845
    Obj                 opL,
846
    Obj                 opR )
847
{
848
    Int                 i;              /* loop variable                   */
849
    Int                 k;              /* loop variable                   */
850
    Int                 cs;             /* sum of two smalls */
851
    UInt                 c;              /* sum of two digits               */
852
    TypDigit *          l;              /* pointer into the left operand   */
853
    TypDigit *          r;              /* pointer into the right operand  */
854
    TypDigit *          s;              /* pointer into the sum            */
855
    UInt *              l2;             /* pointer to get 2 digits at once */
856
    UInt *              s2;             /* pointer to put 2 digits at once */
857
    Obj                 sum;            /* handle of the result bag        */
858

859
    /* adding two small integers                                           */
860
    if ( ARE_INTOBJS( opL, opR ) ) {
861

862
        /* add two small integers with a small sum                         */
863
        /* add and compare top two bits to check that no overflow occured  */
864
        if ( SUM_INTOBJS( sum, opL, opR ) ) {
865
            return sum;
866
        }
867

868
        /* add two small integers with a large sum                         */
869
        cs = INT_INTOBJ(opL) + INT_INTOBJ(opR);
870
        if ( 0 < cs ) {
871
            sum = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
872
            ADDR_INT(sum)[0] = (TypDigit)cs;
873
            ADDR_INT(sum)[1] = (TypDigit)(cs >> NR_DIGIT_BITS);
874
        }
875
        else {
876
            sum = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
877
            ADDR_INT(sum)[0] = (TypDigit)(-cs);
878
            ADDR_INT(sum)[1] = (TypDigit)((-cs) >> NR_DIGIT_BITS);
879
        }
880

881
    }
882

883
    /* adding one large integer and one small integer                      */
884
    else if ( IS_INTOBJ(opL) || IS_INTOBJ(opR) ) {
885

886
        /* make the right operand the small one                            */
887
        if ( IS_INTOBJ(opL) ) {
888
            sum = opL;  opL = opR;  opR = sum;
889
        }
890

891
        /* if the integers have different sign, let 'DiffInt' do the work  */
892
        if ( (TNUM_OBJ(opL) == T_INTNEG && 0 <= INT_INTOBJ(opR))
893
          || (TNUM_OBJ(opL) == T_INTPOS && INT_INTOBJ(opR) <  0) ) {
894
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
895
            else                              RetypeBag( opL, T_INTPOS );
896
            sum = DiffInt( opR, opL );
897
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
898
            else                              RetypeBag( opL, T_INTPOS );
899
            return sum;
900
        }
901

902
        /* allocate the result bag and set up the pointers                 */
903
        if ( TNUM_OBJ(opL) == T_INTPOS ) {
904
            i   = INT_INTOBJ(opR);
905
            sum = NewBag( T_INTPOS, (SIZE_INT(opL)+4)*sizeof(TypDigit) );
906
        }
907
        else {
908
            i   = -INT_INTOBJ(opR);
909
            sum = NewBag( T_INTNEG, (SIZE_INT(opL)+4)*sizeof(TypDigit) );
910
        }
911
        l = ADDR_INT(opL);
912
        s = ADDR_INT(sum);
913

914
        /* add the first four digits,the right operand has only two digits */
915
        c = (UInt)*l++ + (TypDigit)i;                             *s++ = (TypDigit)c;
916
        c = (UInt)*l++ + (i>>NR_DIGIT_BITS) + (c>>NR_DIGIT_BITS); *s++ = (TypDigit)c;
917
        c = (UInt)*l++                      + (c>>NR_DIGIT_BITS); *s++ = (TypDigit)c;
918
        c = (UInt)*l++                      + (c>>NR_DIGIT_BITS); *s++ = (TypDigit)c;
919

920
        /* propagate the carry, this loop is almost never executed         */
921
        for ( k = SIZE_INT(opL)/4-1; k != 0 && (c>>NR_DIGIT_BITS) != 0; k-- ) {
922
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
923
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
924
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
925
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
926
        }
927

928
        /* just copy the remaining digits, do it two digits at once        */
929
        for ( l2 = (UInt*)l, s2 = (UInt*)s; k != 0; k-- ) {
930
            *s2++ = *l2++;
931
            *s2++ = *l2++;
932
        }
933

934
        /* if there is a carry, enter it, otherwise shrink the sum         */
935
        if ( (c>>NR_DIGIT_BITS) != 0 )
936
            *s++ = (TypDigit)(c>>NR_DIGIT_BITS);
937
        else
938
            ResizeBag( sum, (SIZE_INT(sum)-4)*sizeof(TypDigit) );
939

940
    }
941

942
    /* add two large integers                                              */
943
    else {
944

945
        /* if the integers have different sign, let 'DiffInt' do the work  */
946
        if ( (TNUM_OBJ(opL) == T_INTPOS && TNUM_OBJ(opR) == T_INTNEG)
947
          || (TNUM_OBJ(opL) == T_INTNEG && TNUM_OBJ(opR) == T_INTPOS) ) {
948
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
949
            else                              RetypeBag( opL, T_INTPOS );
950
            sum = DiffInt( opR, opL );
951
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
952
            else                              RetypeBag( opL, T_INTPOS );
953
            return sum;
954
        }
955

956
        /* make the right operand the smaller one                          */
957
        if ( SIZE_INT(opL) < SIZE_INT(opR) ) {
958
            sum = opL;  opL = opR;  opR = sum;
959
        }
960

961
        /* allocate the result bag and set up the pointers                 */
962
        if ( TNUM_OBJ(opL) == T_INTPOS ) {
963
            sum = NewBag( T_INTPOS, (SIZE_INT(opL)+4)*sizeof(TypDigit) );
964
        }
965
        else {
966
            sum = NewBag( T_INTNEG, (SIZE_INT(opL)+4)*sizeof(TypDigit) );
967
        }
968
        l = ADDR_INT(opL);
969
        r = ADDR_INT(opR);
970
        s = ADDR_INT(sum);
971

972
        /* add the digits, convert to UInt to get maximum precision         */
973
        c = 0;
974
        for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
975
            c = (UInt)*l++ + (UInt)*r++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
976
            c = (UInt)*l++ + (UInt)*r++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
977
            c = (UInt)*l++ + (UInt)*r++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
978
            c = (UInt)*l++ + (UInt)*r++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
979
        }
980

981
        /* propagate the carry, this loop is almost never executed         */
982
        for ( k=(SIZE_INT(opL)-SIZE_INT(opR))/4;
983
             k!=0 && (c>>NR_DIGIT_BITS)!=0; k-- ) {
984
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
985
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
986
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
987
            c = (UInt)*l++ + (c>>NR_DIGIT_BITS);  *s++ = (TypDigit)c;
988
        }
989

990
        /* just copy the remaining digits, do it two digits at once        */
991
        for ( l2 = (UInt*)l, s2 = (UInt*)s; k != 0; k-- ) {
992
            *s2++ = *l2++;
993
            *s2++ = *l2++;
994
        }
995

996
        /* if there is a carry, enter it, otherwise shrink the sum         */
997
        if ( (c>>NR_DIGIT_BITS) != 0 )
998
            *s++ = (TypDigit)(c>>NR_DIGIT_BITS);
999
        else
1000
            ResizeBag( sum, (SIZE_INT(sum)-4)*sizeof(TypDigit) );
1001

1002
    }
1003

1004
    /* return the sum                                                      */
1005
    return sum;
1006
}
1007

1008

1009
/****************************************************************************
1010
**
1011
*F  ZeroInt(<int>)  . . . . . . . . . . . . . . . . . . . .  zero of integers
1012
*/
1013
Obj         ZeroInt (
1014
    Obj                 op )
1015
{
1016
    return INTOBJ_INT(0);
1017
}
1018

1019

1020
/****************************************************************************
1021
**
1022
*F  AInvInt(<int>)  . . . . . . . . . . . . .  additive inverse of an integer
1023
*/
1024
Obj         AInvInt (
1025
    Obj                 op )
1026
{
1027
    Obj                 inv;
1028
    UInt                i;
1029

1030
    /* handle small integer                                                */
1031
    if ( IS_INTOBJ( op ) ) {
1032

1033
        /* special case (ugh)                                              */
1034
        if ( op == INTOBJ_INT( -(1L<<NR_SMALL_INT_BITS) ) ) {
1035
            inv = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
1036
            ADDR_INT(inv)[0] = 0;
1037
            ADDR_INT(inv)[1] = (TypDigit)(1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS));
1038
        }
1039

1040
        /* general case                                                    */
1041
        else {
1042
            inv = INTOBJ_INT( - INT_INTOBJ( op ) );
1043
        }
1044

1045
    }
1046

1047
    /* invert a large integer                                              */
1048
    else {
1049

1050
        /* special case (ugh)                                              */
1051
        if ( TNUM_OBJ(op) == T_INTPOS && SIZE_INT(op) == 4
1052
          && ADDR_INT(op)[3] == 0
1053
          && ADDR_INT(op)[2] == 0
1054
          && ADDR_INT(op)[1] == (1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS))
1055
          && ADDR_INT(op)[0] == 0 ) {
1056
            inv = INTOBJ_INT( -(1L<<NR_SMALL_INT_BITS) );
1057
        }
1058

1059
        /* general case                                                    */
1060
        else {
1061
            if ( TNUM_OBJ(op) == T_INTPOS ) {
1062
                inv = NewBag( T_INTNEG, SIZE_OBJ(op) );
1063
            }
1064
            else {
1065
                inv = NewBag( T_INTPOS, SIZE_OBJ(op) );
1066
            }
1067
            for ( i = 0; i < SIZE_INT(op); i++ ) {
1068
                ADDR_INT(inv)[i] = ADDR_INT(op)[i];
1069
            }
1070
        }
1071

1072
    }
1073

1074
    /* return the inverse                                                  */
1075
    return inv;
1076
}
1077

1078

1079
/****************************************************************************
1080
**
1081
*F  DiffInt( <intL>, <intR> ) . . . . . . . . . .  difference of two integers
1082
**
1083
**  'DiffInt' returns the difference of the two integer arguments <intL>  and
1084
**  <intR>.  'DiffInt' handles  operands  of  type  'T_INT',  'T_INTPOS'  and
1085
**  'T_INTNEG'.
1086
**
1087
**  It can also be used in the cases that both operands  are  small  integers
1088
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
1089
**  This case is usually already handled in 'EvalDiff' for a better efficiency.
1090
**
1091
**  Is called from the 'EvalDiff' binop so both operands are already evaluated.
1092
**
1093
**  'DiffInt' is a little bit difficult since there are 16 different cases to
1094
**  handle, each operand can be positive or negative, small or large integer.
1095
**  If the operands have opposite sign 'DiffInt' calls 'SumInt',  this  helps
1096
**  reduce the total amount of code by a factor of two.
1097
*/
1098
Obj             DiffInt (
1099
    Obj                 opL,
1100
    Obj                 opR )
1101
{
1102
    Int                 i;              /* loop variable                   */
1103
    Int                 k;              /* loop variable                   */
1104
    Int                 c;              /* difference of two digits        */
1105
    TypDigit *          l;              /* pointer into the left operand   */
1106
    TypDigit *          r;              /* pointer into the right operand  */
1107
    TypDigit *          d;              /* pointer into the difference     */
1108
    UInt *              l2;             /* pointer to get 2 digits at once */
1109
    UInt *              d2;             /* pointer to put 2 digits at once */
1110
    Obj                 dif;            /* handle of the result bag        */
1111

1112
    /* subtracting two small integers                                      */
1113
    if ( ARE_INTOBJS( opL, opR ) ) {
1114

1115
        /* subtract two small integers with a small difference             */
1116
        /* sub and compare top two bits to check that no overflow occured  */
1117
        if ( DIFF_INTOBJS( dif, opL, opR ) ) {
1118
            return dif;
1119
        }
1120

1121
        /* subtract two small integers with a large difference             */
1122
        c = INT_INTOBJ(opL) - INT_INTOBJ(opR);
1123
        if ( 0 < c ) {
1124
            dif = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
1125
            ADDR_INT(dif)[0] = (TypDigit)c;
1126
            ADDR_INT(dif)[1] = (TypDigit)(c >> NR_DIGIT_BITS);
1127
        }
1128
        else {
1129
            dif = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
1130
            ADDR_INT(dif)[0] = (TypDigit)(-c);
1131
            ADDR_INT(dif)[1] = (TypDigit)((-c) >> NR_DIGIT_BITS);
1132
        }
1133

1134
    }
1135

1136
    /* subtracting one small integer and one large integer                 */
1137
    else if ( IS_INTOBJ( opL ) || IS_INTOBJ( opR ) ) {
1138

1139
        /* make the right operand the small one                            */
1140
        if ( IS_INTOBJ( opL ) ) {
1141
            dif = opL;  opL = opR;  opR = dif;
1142
            c = -1;
1143
        }
1144
        else {
1145
            c =  1;
1146
        }
1147

1148
        /* if the integers have different sign, let 'SumInt' do the work   */
1149
        if ( (TNUM_OBJ(opL) == T_INTNEG && 0 <= INT_INTOBJ(opR))
1150
          || (TNUM_OBJ(opL) == T_INTPOS && INT_INTOBJ(opR) < 0)  ) {
1151
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
1152
            else                              RetypeBag( opL, T_INTPOS );
1153
            dif = SumInt( opL, opR );
1154
            if ( TNUM_OBJ(opL) == T_INTPOS )  RetypeBag( opL, T_INTNEG );
1155
            else                              RetypeBag( opL, T_INTPOS );
1156
            if ( c == 1 ) {
1157
                if ( TNUM_OBJ(dif) == T_INTPOS )  RetypeBag( dif, T_INTNEG );
1158
                else                              RetypeBag( dif, T_INTPOS );
1159
            }
1160
            return dif;
1161
        }
1162

1163
        /* allocate the result bag and set up the pointers                 */
1164
        if ( TNUM_OBJ(opL) == T_INTPOS ) {
1165
            i   = INT_INTOBJ(opR);
1166
            if ( c == 1 )  dif = NewBag( T_INTPOS, SIZE_OBJ(opL) );
1167
            else           dif = NewBag( T_INTNEG, SIZE_OBJ(opL) );
1168
        }
1169
        else {
1170
            i   = - INT_INTOBJ(opR);
1171
            if ( c == 1 )  dif = NewBag( T_INTNEG, SIZE_OBJ(opL) );
1172
            else           dif = NewBag( T_INTPOS, SIZE_OBJ(opL) );
1173
        }
1174
        l = ADDR_INT(opL);
1175
        d = ADDR_INT(dif);
1176

1177
        /* sub the first four digit, note the left operand has only two    */
1178
        /*N (c>>16<) need not work, replace by (c<0?-1:0)                   */
1179
        c = (Int)*l++ - (TypDigit)i;                              *d++ = (TypDigit)c;
1180
        c = (Int)*l++ - (TypDigit)(i/(1L<<NR_DIGIT_BITS)) + (c<0?-1:0);  *d++ = (TypDigit)c;
1181
        c = (Int)*l++                      + (c<0?-1:0);  *d++ = (TypDigit)c;
1182
        c = (Int)*l++                      + (c<0?-1:0);  *d++ = (TypDigit)c;
1183

1184
        /* propagate the carry, this loop is almost never executed         */
1185
        for ( k = SIZE_INT(opL)/4-1; k != 0 && c < 0; k-- ) {
1186
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1187
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1188
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1189
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1190
        }
1191

1192
        /* just copy the remaining digits, do it two digits at once        */
1193
        for ( l2 = (UInt*)l, d2 = (UInt*)d; k != 0; k-- ) {
1194
            *d2++ = *l2++;
1195
            *d2++ = *l2++;
1196
        }
1197

1198
        /* no underflow since we subtracted a small int from a large one   */
1199
        /* but there may be leading zeroes in the result, get rid of them  */
1200
        /* occurs almost never, so it doesn't matter that it is expensive  */
1201
        if ( ((UInt*)d == d2
1202
          && d[-4] == 0 && d[-3] == 0 && d[-2] == 0 && d[-1] == 0)
1203
          || (SIZE_INT(dif) == 4 && d[-2] == 0 && d[-1] == 0) ) {
1204

1205
            /* find the number of significant digits                       */
1206
            d = ADDR_INT(dif);
1207
            for ( k = SIZE_INT(dif); k != 0; k-- ) {
1208
                if ( d[k-1] != 0 )
1209
                    break;
1210
            }
1211

1212
            /* reduce to small integer if possible, otherwise shrink bag   */
1213
            if ( k <= 2 && TNUM_OBJ(dif) == T_INTPOS
1214
              && (UInt)(INTBASE*d[1]+d[0])<(1L<<NR_SMALL_INT_BITS) )
1215
                dif = INTOBJ_INT( INTBASE*d[1]+d[0] );
1216
            else if ( k <= 2 && TNUM_OBJ(dif) == T_INTNEG
1217
              && (UInt)(INTBASE*d[1]+d[0])<=(1L<<NR_SMALL_INT_BITS) )
1218
                dif = INTOBJ_INT( -(Int)(INTBASE*d[1]+d[0]) );
1219
            else
1220
                ResizeBag( dif, (((k + 3) / 4) * 4) * sizeof(TypDigit) );
1221
        }
1222

1223
    }
1224

1225
    /* subtracting two large integers                                      */
1226
    else {
1227

1228
        /* if the integers have different sign, let 'SumInt' do the work   */
1229
        if ( (TNUM_OBJ(opL) == T_INTPOS && TNUM_OBJ(opR) == T_INTNEG)
1230
          || (TNUM_OBJ(opL) == T_INTNEG && TNUM_OBJ(opR) == T_INTPOS) ) {
1231
            if ( TNUM_OBJ(opR) == T_INTPOS )  RetypeBag( opR, T_INTNEG );
1232
            else                              RetypeBag( opR, T_INTPOS );
1233
            dif = SumInt( opL, opR );
1234
            if ( TNUM_OBJ(opR) == T_INTPOS )  RetypeBag( opR, T_INTNEG );
1235
            else                              RetypeBag( opR, T_INTPOS );
1236
            return dif;
1237
        }
1238

1239
        /* make the right operand the smaller one                          */
1240
        if ( SIZE_INT(opL) <  SIZE_INT(opR)
1241
          || (TNUM_OBJ(opL) == T_INTPOS && LtInt(opL,opR) )
1242
          || (TNUM_OBJ(opL) == T_INTNEG && LtInt(opR,opL) ) ) {
1243
            dif = opL;  opL = opR;  opR = dif;  c = -1;
1244
        }
1245
        else {
1246
            c = 1;
1247
        }
1248

1249
        /* allocate the result bag and set up the pointers                 */
1250
        if ( (TNUM_OBJ(opL) == T_INTPOS && c ==  1)
1251
          || (TNUM_OBJ(opL) == T_INTNEG && c == -1) )
1252
            dif = NewBag( T_INTPOS, SIZE_OBJ(opL) );
1253
        else
1254
            dif = NewBag( T_INTNEG, SIZE_OBJ(opL) );
1255
        l = ADDR_INT(opL);
1256
        r = ADDR_INT(opR);
1257
        d = ADDR_INT(dif);
1258

1259
        /* subtract the digits                                             */
1260
        c = 0;
1261
        for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
1262
            c = (Int)*l++ - (Int)*r++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1263
            c = (Int)*l++ - (Int)*r++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1264
            c = (Int)*l++ - (Int)*r++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1265
            c = (Int)*l++ - (Int)*r++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1266
        }
1267

1268
        /* propagate the carry, this loop is almost never executed         */
1269
        for ( k=(SIZE_INT(opL)-SIZE_INT(opR))/4;
1270
             k!=0 && c < 0; k-- ) {
1271
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1272
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1273
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1274
            c = (Int)*l++ + (c < 0 ? -1 : 0);  *d++ = (TypDigit)c;
1275
        }
1276

1277
        /* just copy the remaining digits, do it two digits at once        */
1278
        for ( d2 = (UInt*)d, l2 = (UInt*)l; k != 0; k-- ) {
1279
            *d2++ = *l2++;
1280
            *d2++ = *l2++;
1281
        }
1282

1283
        /* no underflow since we subtracted a small int from a large one   */
1284
        /* but there may be leading zeroes in the result, get rid of them  */
1285
        /* occurs almost never, so it doesn't matter that it is expensive  */
1286
        if ( ((UInt*)d == d2
1287
          && d[-4] == 0 && d[-3] == 0 && d[-2] == 0 && d[-1] == 0)
1288
          || (SIZE_INT(dif) == 4 && d[-2] == 0 && d[-1] == 0) ) {
1289

1290
            /* find the number of significant digits                       */
1291
            d = ADDR_INT(dif);
1292
            for ( k = SIZE_INT(dif); k != 0; k-- ) {
1293
                if ( d[k-1] != 0 )
1294
                    break;
1295
            }
1296

1297
            /* reduce to small integer if possible, otherwise shrink bag   */
1298
            if ( k <= 2 && TNUM_OBJ(dif) == T_INTPOS
1299
              && (UInt)(INTBASE*d[1]+d[0]) < (1L<<NR_SMALL_INT_BITS) )
1300
                dif = INTOBJ_INT( INTBASE*d[1]+d[0] );
1301
            else if ( k <= 2 && TNUM_OBJ(dif) == T_INTNEG
1302
              && (UInt)(INTBASE*d[1]+d[0])<=(1L<<NR_SMALL_INT_BITS))
1303
                dif = INTOBJ_INT( -(Int)(INTBASE*d[1]+d[0]) );
1304
            else
1305
                ResizeBag( dif, (((k + 3) / 4) * 4) * sizeof(TypDigit) );
1306

1307
        }
1308

1309
    }
1310

1311
    /* return the difference                                               */
1312
    return dif;
1313
}
1314

1315

1316
/****************************************************************************
1317
**
1318
*F  ProdInt( <intL>, <intR> ) . . . . . . . . . . . . product of two integers
1319
**
1320
**  'ProdInt' returns the product of the two  integer  arguments  <intL>  and
1321
**  <intR>.  'ProdInt' handles  operands  of  type  'T_INT',  'T_INTPOS'  and
1322
**  'T_INTNEG'.
1323
**
1324
**  It can also be used in the cases that both operands  are  small  integers
1325
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
1326
**  This case is usually already handled in 'EvalProd' for a better efficiency.
1327
**
1328
**  Is called from the 'EvalProd' binop so both operands are already evaluated.
1329
**
1330
**  The only difficulty about this function is the fact that is has to handle
1331
**  3 different situations, depending on how many arguments  are  small  ints.
1332
*/
1333
Obj             ProdInt (
1334
    Obj                 opL,
1335
    Obj                 opR )
1336
{
1337
    Int                 i;              /* loop count, value for small int */
1338
    Int                 k;              /* loop count, value for small int */
1339
    UInt                c;              /* product of two digits           */
1340
    TypDigit            l;              /* one digit of the left operand   */
1341
    TypDigit *          r;              /* pointer into the right operand  */
1342
    TypDigit *          p;              /* pointer into the product        */
1343
    Obj                 prd;            /* handle of the result bag        */
1344

1345
    /* multiplying two small integers                                      */
1346
    if ( ARE_INTOBJS( opL, opR ) ) {
1347

1348
        /* multiply two small integers with a small product                */
1349
        /* multiply and divide back to check that no overflow occured      */
1350
        if ( PROD_INTOBJS( prd, opL, opR ) ) {
1351
            return prd;
1352
        }
1353

1354
        /* get the integer values                                          */
1355
        i = INT_INTOBJ(opL);
1356
        k = INT_INTOBJ(opR);
1357

1358
        /* allocate the product bag                                        */
1359
        if ( (0 < i && 0 < k) || (i < 0 && k < 0) )
1360
            prd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
1361
        else
1362
            prd = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
1363
        p = ADDR_INT(prd);
1364

1365
        /* make both operands positive                                     */
1366
        if ( i < 0 )  i = -i;
1367
        if ( k < 0 )  k = -k;
1368

1369
        /* multiply digitwise                                              */
1370
        c = (UInt)(TypDigit)i * (TypDigit)k;            p[0] = (TypDigit)c;
1371
        c = (UInt)(TypDigit)i * (((UInt)k)>>NR_DIGIT_BITS)
1372
          + (c>>NR_DIGIT_BITS);                        p[1] = (TypDigit)c;
1373
        p[2] = c>>NR_DIGIT_BITS;
1374

1375
        c = (UInt)(TypDigit)(((UInt)i)>>NR_DIGIT_BITS) * (TypDigit)k
1376
          + p[1];                                      p[1] = (TypDigit)c;
1377
        c = (UInt)(TypDigit)(((UInt)i)>>NR_DIGIT_BITS) * (TypDigit)(((UInt)k)>>NR_DIGIT_BITS)
1378
          + p[2] + (c>>NR_DIGIT_BITS);                 p[2] = (TypDigit)c;
1379
        p[3] = (TypDigit)(c>>NR_DIGIT_BITS);
1380

1381
    }
1382

1383
    /* multiply a small and a large integer                                */
1384
    else if ( IS_INTOBJ(opL) || IS_INTOBJ(opR) ) {
1385

1386
        /* make the left operand the small one                             */
1387
        if ( IS_INTOBJ(opR) ) {
1388
            i = INT_INTOBJ(opR);  opR = opL;
1389
        }
1390
        else {
1391
            i = INT_INTOBJ(opL);
1392
        }
1393

1394
        /* handle trivial cases first                                      */
1395
        if ( i == 0 )
1396
            return INTOBJ_INT(0);
1397
        if ( i == 1 )
1398
            return opR;
1399

1400
        /* the large integer 1<<28 times -1 is the small integer -(1<<28)  */
1401
        if ( i == -1
1402
          && TNUM_OBJ(opR) == T_INTPOS && SIZE_INT(opR) == 4
1403
          && ADDR_INT(opR)[3] == 0
1404
          && ADDR_INT(opR)[2] == 0
1405
          && ADDR_INT(opR)[1] == (1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS))
1406
          && ADDR_INT(opR)[0] == 0 )
1407
            return INTOBJ_INT( -(Int)(1L<<NR_SMALL_INT_BITS) );
1408

1409
        /* multiplication by -1 is easy, just switch the sign and copy     */
1410
        if ( i == -1 ) {
1411
            if ( TNUM_OBJ(opR) == T_INTPOS )
1412
                prd = NewBag( T_INTNEG, SIZE_OBJ(opR) );
1413
            else
1414
                prd = NewBag( T_INTPOS, SIZE_OBJ(opR) );
1415
            r = ADDR_INT(opR);
1416
            p = ADDR_INT(prd);
1417
            for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
1418
                /*N should be: *p2++=*r2++;  *p2++=*r2++;                  */
1419
                *p++ = *r++;  *p++ = *r++;  *p++ = *r++;  *p++ = *r++;
1420
            }
1421
            return prd;
1422
        }
1423

1424
        /* allocate a bag for the result                                   */
1425
        if ( (0 < i && TNUM_OBJ(opR) == T_INTPOS)
1426
          || (i < 0 && TNUM_OBJ(opR) == T_INTNEG) )
1427
            prd = NewBag( T_INTPOS, (SIZE_INT(opR)+4)*sizeof(TypDigit) );
1428
        else
1429
            prd = NewBag( T_INTNEG, (SIZE_INT(opR)+4)*sizeof(TypDigit) );
1430
        if ( i < 0 )  i = -i;
1431

1432
        /* multiply with the lower digit of the left operand               */
1433
        l = (TypDigit)i;
1434
        if ( l != 0 ) {
1435

1436
            r = ADDR_INT(opR);
1437
            p = ADDR_INT(prd);
1438
            c = 0;
1439

1440
            /* multiply the right with this digit and store in the product */
1441
            for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
1442
                c = (UInt)l * (UInt)*r++ + (c>>NR_DIGIT_BITS);  *p++ = (TypDigit)c;
1443
                c = (UInt)l * (UInt)*r++ + (c>>NR_DIGIT_BITS);  *p++ = (TypDigit)c;
1444
                c = (UInt)l * (UInt)*r++ + (c>>NR_DIGIT_BITS);  *p++ = (TypDigit)c;
1445
                c = (UInt)l * (UInt)*r++ + (c>>NR_DIGIT_BITS);  *p++ = (TypDigit)c;
1446
            }
1447
            *p = (TypDigit)(c>>NR_DIGIT_BITS);
1448
        }
1449

1450
        /* multiply with the larger digit of the left operand              */
1451
        l = ((UInt)i) >> NR_DIGIT_BITS;
1452
        if ( l != 0 ) {
1453

1454
            r = ADDR_INT(opR);
1455
            p = ADDR_INT(prd) + 1;
1456
            c = 0;
1457

1458
            /* multiply the right with this digit and add into the product */
1459
            for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
1460
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1461
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1462
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1463
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1464
            }
1465
            *p = (TypDigit)(c>>NR_DIGIT_BITS);
1466
        }
1467

1468
        /* remove the leading zeroes, note that there can't be more than 6 */
1469
        p = ADDR_INT(prd) + SIZE_INT(prd);
1470
        if ( p[-4] == 0 && p[-3] == 0 && p[-2] == 0 && p[-1] == 0 ) {
1471
            ResizeBag( prd, (SIZE_INT(prd)-4)*sizeof(TypDigit) );
1472
        }
1473

1474
    }
1475

1476
    /* multiply two large integers                                         */
1477
    else {
1478

1479
        /* make the left operand the smaller one, for performance          */
1480
        if ( SIZE_INT(opL) > SIZE_INT(opR) ) {
1481
            prd = opR;  opR = opL;  opL = prd;
1482
        }
1483

1484
        /* allocate a bag for the result                                   */
1485
        if ( TNUM_OBJ(opL) == TNUM_OBJ(opR) )
1486
            prd = NewBag( T_INTPOS, SIZE_OBJ(opL)+SIZE_OBJ(opR) );
1487
        else
1488
            prd = NewBag( T_INTNEG, SIZE_OBJ(opL)+SIZE_OBJ(opR) );
1489

1490
        /* run through the digits of the left operand                      */
1491
        for ( i = 0; i < (Int)SIZE_INT(opL); i++ ) {
1492

1493
            /* set up pointer for one loop iteration                       */
1494
            l = ADDR_INT(opL)[i];
1495
            if ( l == 0 )  continue;
1496
            r = ADDR_INT(opR);
1497
            p = ADDR_INT(prd) + i;
1498
            c = 0;
1499

1500
            /* multiply the right with this digit and add into the product */
1501
            for ( k = SIZE_INT(opR)/4; k != 0; k-- ) {
1502
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1503
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1504
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1505
                c = (UInt)l * (UInt)*r++ + (UInt)*p + (c>>NR_DIGIT_BITS); *p++ = (TypDigit)c;
1506
            }
1507
            *p = (TypDigit)(c>>NR_DIGIT_BITS);
1508
        }
1509

1510
        /* remove the leading zeroes, note that there can't be more than 7 */
1511
        p = ADDR_INT(prd) + SIZE_INT(prd);
1512
        if ( p[-4] == 0 && p[-3] == 0 && p[-2] == 0 && p[-1] == 0 ) {
1513
            ResizeBag( prd, (SIZE_INT(prd)-4)*sizeof(TypDigit) );
1514
        }
1515

1516
    }
1517

1518
    /* return the product                                                  */
1519
    return prd;
1520
}
1521

1522

1523
/****************************************************************************
1524
**
1525
*F  ProdIntObj(<n>,<op>)  . . . . . . . . product of an integer and an object
1526
*/
1527
Obj             ProdIntObj (
1528
    Obj                 n,
1529
    Obj                 op )
1530
{
1531
    Obj                 res = 0;        /* result                          */
1532
    UInt                i, k, l;        /* loop variables                  */
1533

1534
    /* if the integer is zero, return the neutral element of the operand   */
1535
    if      ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) ==  0 ) {
1536
        res = ZERO( op );
1537
    }
1538

1539
    /* if the integer is one, return the object if immutable
1540
       if mutable, add the object to its ZeroSameMutability to
1541
       ensure correct mutability propagation */
1542
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) ==  1 ) {
1543
      if (IS_MUTABLE_OBJ(op))
1544
        res = SUM(ZERO(op),op);
1545
      else
1546
	res = op;
1547
    }
1548

1549
    /* if the integer is minus one, return the inverse of the operand      */
1550
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) == -1 ) {
1551
        res = AINV( op );
1552
    }
1553

1554
    /* if the integer is negative, invert the operand and the integer      */
1555
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) <  -1 ) {
1556
        res = AINV( op );
1557
        if ( res == Fail ) {
1558
            return ErrorReturnObj(
1559
                "Operations: <obj> must have an additive inverse",
1560
                0L, 0L,
1561
                "you can supply an inverse <inv> for <obj> via 'return <inv>;'" );
1562
        }
1563
        res = PROD( AINV( n ), res );
1564
    }
1565

1566
    /* if the integer is negative, invert the operand and the integer      */
1567
    else if ( TNUM_OBJ(n) == T_INTNEG ) {
1568
        res = AINV( op );
1569
        if ( res == Fail ) {
1570
            return ErrorReturnObj(
1571
                "Operations: <obj> must have an additive inverse",
1572
                0L, 0L,
1573
                "you can supply an inverse <inv> for <obj> via 'return <inv>;'" );
1574
        }
1575
        res = PROD( AINV( n ), res );
1576
    }
1577

1578
    /* if the integer is small, compute the product by repeated doubling   */
1579
    /* the loop invariant is <result> = <k>*<res> + <l>*<op>, <l> < <k>    */
1580
    /* <res> = 0 means that <res> is the neutral element                   */
1581
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) >   1 ) {
1582
        res = 0;
1583
        k = 1L << (NR_SMALL_INT_BITS+1);
1584
        l = INT_INTOBJ(n);
1585
        while ( 1 < k ) {
1586
            res = (res == 0 ? res : SUM( res, res ));
1587
            k = k / 2;
1588
            if ( k <= l ) {
1589
                res = (res == 0 ? op : SUM( res, op ));
1590
                l = l - k;
1591
            }
1592
        }
1593
    }
1594

1595
    /* if the integer is large, compute the product by repeated doubling   */
1596
    else if ( TNUM_OBJ(n) == T_INTPOS ) {
1597
        res = 0;
1598
        for ( i = SIZE_OBJ(n)/sizeof(TypDigit); 0 < i; i-- ) {
1599
            k = 1L << (8*sizeof(TypDigit));
1600
            l = ((TypDigit*) ADDR_OBJ(n))[i-1];
1601
            while ( 1 < k ) {
1602
                res = (res == 0 ? res : SUM( res, res ));
1603
                k = k / 2;
1604
                if ( k <= l ) {
1605
                    res = (res == 0 ? op : SUM( res, op ));
1606
                    l = l - k;
1607
                }
1608
            }
1609
        }
1610
    }
1611

1612
    /* return the result                                                   */
1613
    return res;
1614
}
1615

1616
Obj             ProdIntObjFunc;
1617

1618
Obj             FuncPROD_INT_OBJ (
1619
    Obj                 self,
1620
    Obj                 opL,
1621
    Obj                 opR )
1622
{
1623
    return ProdIntObj( opL, opR );
1624
}
1625

1626

1627
/****************************************************************************
1628
**
1629
*F  OneInt(<int>) . . . . . . . . . . . . . . . . . . . . . one of an integer
1630
*/
1631
Obj             OneInt (
1632
    Obj                 op )
1633
{
1634
    return INTOBJ_INT( 1L );
1635
}
1636

1637

1638
/****************************************************************************
1639
**
1640
*F  PowInt( <intL>, <intR> )  . . . . . . . . . . . . . . power of an integer
1641
**
1642
**  'PowInt' returns the <intR>-th (an integer) power of the integer  <intL>.
1643
**  'PowInt' handles operands of type 'T_INT', 'T_INTPOS' and 'T_INTNEG'.
1644
**
1645
**  It can also be used in the cases that both operands  are  small  integers
1646
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
1647
**  This case is usually already handled in 'EvalPow' for a better  efficiency.
1648
**
1649
**  Is called from the 'EvalPow'  binop so both operands are already evaluated.
1650
*/
1651
Obj             PowInt (
1652
    Obj                 opL,
1653
    Obj                 opR )
1654
{
1655
    Int                 i;
1656
    Obj                 pow;
1657

1658
    /* power with a large exponent                                         */
1659
    ReadGuard(opL);
1660
    ReadGuard(opR);
1661
    if ( ! IS_INTOBJ(opR) ) {
1662
        if ( opL == INTOBJ_INT(0) )
1663
            pow = INTOBJ_INT(0);
1664
        else if ( opL == INTOBJ_INT(1) )
1665
            pow = INTOBJ_INT(1);
1666
        else if ( opL == INTOBJ_INT(-1) && ADDR_INT(opR)[0] % 2 == 0 )
1667
            pow = INTOBJ_INT(1);
1668
        else if ( opL == INTOBJ_INT(-1) && ADDR_INT(opR)[0] % 2 != 0 )
1669
            pow = INTOBJ_INT(-1);
1670
        else {
1671
            opR = ErrorReturnObj(
1672
                "Integer operands: <exponent> is too large",
1673
                0L, 0L,
1674
                "you can replace the integer <exponent> via 'return <exponent>;'" );
1675
            return POW( opL, opR );
1676
        }
1677
    }
1678

1679
    /* power with a negative exponent                                      */
1680
    else if ( INT_INTOBJ(opR) < 0 ) {
1681
        if ( opL == INTOBJ_INT(0) ) {
1682
            opL = ErrorReturnObj(
1683
                "Integer operands: <base> must not be zero",
1684
                0L, 0L,
1685
                "you can replace the integer <base> via 'return <base>;'" );
1686
            return POW( opL, opR );
1687
        }
1688
        else if ( opL == INTOBJ_INT(1) )
1689
            pow = INTOBJ_INT(1);
1690
        else if ( opL == INTOBJ_INT(-1) && INT_INTOBJ(opR) % 2 == 0 )
1691
            pow = INTOBJ_INT(1);
1692
        else if ( opL == INTOBJ_INT(-1) && INT_INTOBJ(opR) % 2 != 0 )
1693
            pow = INTOBJ_INT(-1);
1694
        else
1695
            pow = QUO( INTOBJ_INT(1),
1696
                       PowInt( opL, INTOBJ_INT( -INT_INTOBJ(opR)) ) );
1697
    }
1698

1699
    /* power with a small positive exponent, do it by a repeated squaring  */
1700
    else {
1701
        pow = INTOBJ_INT(1);
1702
        i = INT_INTOBJ(opR);
1703
        while ( i != 0 ) {
1704
            if ( i % 2 == 1 )  pow = ProdInt( pow, opL );
1705
            if ( i     >  1 )  opL = ProdInt( opL, opL );
1706
            i = i / 2;
1707
        }
1708
    }
1709

1710
    /* return the power                                                    */
1711
    return pow;
1712
}
1713

1714

1715
/****************************************************************************
1716
**
1717
*F  PowObjInt(<op>,<n>) . . . . . . . . . . power of an object and an integer
1718
*/
1719
static Obj OneAttr;
1720

1721
Obj             PowObjInt (
1722
    Obj                 op,
1723
    Obj                 n )
1724
{
1725
    Obj                 res = 0;        /* result                          */
1726
    UInt                i, k, l;        /* loop variables                  */
1727

1728
    /* if the integer is zero, return the neutral element of the operand   */
1729
    if      ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) ==  0 ) {
1730
      return ONE_MUT( op );
1731
    }
1732

1733
    /* if the integer is one, return a copy of the operand                 */
1734
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) ==  1 ) {
1735
        res = CopyObj( op, 1 );
1736
    }
1737

1738
    /* if the integer is minus one, return the inverse of the operand      */
1739
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) == -1 ) {
1740
        res = INV_MUT( op );
1741
    }
1742

1743
    /* if the integer is negative, invert the operand and the integer      */
1744
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) <   0 ) {
1745
        res = INV_MUT( op );
1746
        if ( res == Fail ) {
1747
            return ErrorReturnObj(
1748
                "Operations: <obj> must have an inverse",
1749
                0L, 0L,
1750
                "you can supply an inverse <inv> for <obj> via 'return <inv>;'" );
1751
        }
1752
        res = POW( res, AINV( n ) );
1753
    }
1754

1755
    /* if the integer is negative, invert the operand and the integer      */
1756
    else if ( TNUM_OBJ(n) == T_INTNEG ) {
1757
        res = INV_MUT( op );
1758
        if ( res == Fail ) {
1759
            return ErrorReturnObj(
1760
                "Operations: <obj> must have an inverse",
1761
                0L, 0L,
1762
                "you can supply an inverse <inv> for <obj> via 'return <inv>;'" );
1763
        }
1764
        res = POW( res, AINV( n ) );
1765
    }
1766

1767
    /* if the integer is small, compute the power by repeated squaring     */
1768
    /* the loop invariant is <result> = <res>^<k> * <op>^<l>, <l> < <k>    */
1769
    /* <res> = 0 means that <res> is the neutral element                   */
1770
    else if ( TNUM_OBJ(n) == T_INT && INT_INTOBJ(n) >   0 ) {
1771
        res = 0;
1772
        k = 1L << (NR_SMALL_INT_BITS+1);
1773
        l = INT_INTOBJ(n);
1774
        while ( 1 < k ) {
1775
            res = (res == 0 ? res : PROD( res, res ));
1776
            k = k / 2;
1777
            if ( k <= l ) {
1778
                res = (res == 0 ? op : PROD( res, op ));
1779
                l = l - k;
1780
            }
1781
        }
1782
    }
1783

1784
    /* if the integer is large, compute the power by repeated squaring     */
1785
    else if ( TNUM_OBJ(n) == T_INTPOS ) {
1786
        res = 0;
1787
        for ( i = SIZE_OBJ(n)/sizeof(TypDigit); 0 < i; i-- ) {
1788
            k = 1L << (8*sizeof(TypDigit));
1789
            l = ((TypDigit*) ADDR_OBJ(n))[i-1];
1790
            while ( 1 < k ) {
1791
                res = (res == 0 ? res : PROD( res, res ));
1792
                k = k / 2;
1793
                if ( k <= l ) {
1794
                    res = (res == 0 ? op : PROD( res, op ));
1795
                    l = l - k;
1796
                }
1797
            }
1798
        }
1799
    }
1800

1801
    /* return the result                                                   */
1802
    return res;
1803
}
1804

1805
Obj             PowObjIntFunc;
1806

1807
Obj             FuncPOW_OBJ_INT (
1808
    Obj                 self,
1809
    Obj                 opL,
1810
    Obj                 opR )
1811
{
1812
    return PowObjInt( opL, opR );
1813
}
1814

1815

1816
/****************************************************************************
1817
**
1818
*F  ModInt( <intL>, <intR> )  . representative of residue class of an integer
1819
**
1820
**  'ModInt' returns the smallest positive representant of the residue  class
1821
**  of the  integer  <intL>  modulo  the  integer  <intR>.  'ModInt'  handles
1822
**  operands of type 'T_INT', 'T_INTPOS', 'T_INTNEG'.
1823
**
1824
**  It can also be used in the cases that both operands  are  small  integers
1825
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
1826
**  This case is usually already handled in 'EvalMod' for a better efficiency.
1827
p**
1828
**  Is called from the 'EvalMod'  binop so both operands are already evaluated.
1829
*/
1830
Obj             ModInt (
1831
    Obj                 opL,
1832
    Obj                 opR )
1833
{
1834
    Int                 i;              /* loop count, value for small int */
1835
    Int                 k;              /* loop count, value for small int */
1836
    UInt                c;              /* product of two digits           */
1837
    TypDigit            d;              /* carry into the next digit       */
1838
    TypDigit *          l;              /* pointer into the left operand   */
1839
    TypDigit *          r;              /* pointer into the right operand  */
1840
    TypDigit            r1;             /* leading digit of the right oper */
1841
    TypDigit            r2;             /* next digit of the right operand */
1842
    UInt                rs;             /* size of the right operand       */
1843
    UInt                e;              /* we mult r by 2^e so r1 >= 32768 */
1844
    Obj                 mod;            /* handle of the remainder bag     */
1845
    TypDigit *          m;              /* pointer into the remainder      */
1846
    UInt                m01;            /* leading two digits of the rem.  */
1847
    TypDigit            m2;             /* next digit of the remainder     */
1848
    TypDigit            qi;             /* guessed digit of the quotient   */
1849

1850
    /* compute the remainder of two small integers                         */
1851
    if ( ARE_INTOBJS( opL, opR ) ) {
1852

1853
        /* pathological case first                                         */
1854
        if ( opR == INTOBJ_INT(0) ) {
1855
            opR = ErrorReturnObj(
1856
                "Integer operations: <divisor> must be nonzero",
1857
                0L, 0L,
1858
                "you can replace the integer <divisor> via 'return <divisor>;'" );
1859
            return MOD( opL, opR );
1860
        }
1861

1862
        /* get the integer values                                          */
1863
        i = INT_INTOBJ(opL);
1864
        k = INT_INTOBJ(opR);
1865

1866
        /* compute the remainder, make sure we divide only positive numbers*/
1867
        if (      0 <= i && 0 <= k )  i =       (  i %  k );
1868
        else if ( 0 <= i && k <  0 )  i =       (  i % -k );
1869
        else if ( i < 0  && 0 <= k )  i = ( k - ( -i %  k )) % k;
1870
        else if ( i < 0  && k <  0 )  i = (-k - ( -i % -k )) % k;
1871
        mod = INTOBJ_INT( i );
1872

1873
    }
1874

1875
    /* compute the remainder of a small integer by a large integer         */
1876
    else if ( IS_INTOBJ(opL) ) {
1877

1878
        /* the small int -(1<<28) mod the large int (1<<28) is 0           */
1879
        if ( opL == INTOBJ_INT((UInt)-(Int)(1L<<NR_SMALL_INT_BITS))
1880
          && TNUM_OBJ(opR) == T_INTPOS && SIZE_INT(opR) == 4
1881
          && ADDR_INT(opR)[3] == 0
1882
          && ADDR_INT(opR)[2] == 0
1883
          && ADDR_INT(opR)[1] == (NR_SMALL_INT_BITS-NR_DIGIT_BITS)
1884
          && ADDR_INT(opR)[0] == 0 )
1885
            mod = INTOBJ_INT(0);
1886

1887
        /* in all other cases the remainder is equal the left operand      */
1888
        else if ( 0 <= INT_INTOBJ(opL) )
1889
            mod = opL;
1890
        else if ( TNUM_OBJ(opR) == T_INTPOS )
1891
            mod = SumInt( opL, opR );
1892
        else
1893
            mod = DiffInt( opL, opR );
1894

1895
    }
1896

1897
    /* compute the remainder of a large integer by a small integer         */
1898
    else if ( IS_INTOBJ(opR)
1899
           && INT_INTOBJ(opR) < INTBASE
1900
           && -(Int)INTBASE <= INT_INTOBJ(opR) ) {
1901

1902
        /* pathological case first                                         */
1903
        if ( opR == INTOBJ_INT(0) ) {
1904
            opR = ErrorReturnObj(
1905
                "Integer operations: <divisor> must be nonzero",
1906
                0L, 0L,
1907
                "you can replace the integer <divisor> via 'return <divisor>;'" );
1908
            return MOD( opL, opR );
1909
        }
1910

1911
        /* get the integer value, make positive                            */
1912
        i = INT_INTOBJ(opR);  if ( i < 0 )  i = -i;
1913

1914
        /* maybe its trivial                                               */
1915
        if ( INTBASE % i == 0 ) {
1916
            c = ADDR_INT(opL)[0] % i;
1917
        }
1918

1919
        /* otherwise run through the left operand and divide digitwise     */
1920
        else {
1921
            l = ADDR_INT(opL) + SIZE_INT(opL) - 1;
1922
            c = 0;
1923
            for ( ; l >= ADDR_INT(opL); l-- ) {
1924
                c  = (c<<NR_DIGIT_BITS) + (Int)*l;
1925
                c  = c % i;
1926
            }
1927
        }
1928

1929
        /* now c is the result, it has the same sign as the left operand   */
1930
        if ( TNUM_OBJ(opL) == T_INTPOS )
1931
            mod = INTOBJ_INT( c );
1932
        else if ( c == 0 )
1933
            mod = INTOBJ_INT( c );
1934
        else if ( 0 <= INT_INTOBJ(opR) )
1935
            mod = SumInt( INTOBJ_INT( -(Int)c ), opR );
1936
        else
1937
            mod = DiffInt( INTOBJ_INT( -(Int)c ), opR );
1938

1939
    }
1940

1941
    /* compute the remainder of a large integer modulo a large integer     */
1942
    else {
1943
        /* a small divisor larger than one digit isn't handled above       */
1944
        if ( IS_INTOBJ(opR) ) {
1945
            if ( 0 < INT_INTOBJ(opR) ) {
1946
                mod = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
1947
                ADDR_INT(mod)[0] = (TypDigit)(INT_INTOBJ(opR));
1948
                ADDR_INT(mod)[1] = (TypDigit)(INT_INTOBJ(opR)>>NR_DIGIT_BITS);
1949
                opR = mod;
1950
            }
1951
            else {
1952
                mod = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
1953
                ADDR_INT(mod)[0] = (TypDigit)(-INT_INTOBJ(opR));
1954
                ADDR_INT(mod)[1] = (TypDigit)((-INT_INTOBJ(opR))>>NR_DIGIT_BITS);
1955
                opR = mod;
1956
            }
1957
        }
1958

1959
        /* trivial case first                                              */
1960
        if ( SIZE_INT(opL) < SIZE_INT(opR) ) {
1961
            if ( TNUM_OBJ(opL) == T_INTPOS )
1962
                return opL;
1963
            else if ( TNUM_OBJ(opR) == T_INTPOS )
1964
                return SumInt( opL, opR );
1965
            else
1966
                return DiffInt( opL, opR );
1967
        }
1968

1969
        /* copy the left operand into a new bag, this holds the remainder  */
1970
        mod = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+4)*sizeof(TypDigit) );
1971
        l = ADDR_INT(opL);
1972
        m = ADDR_INT(mod);
1973
        for ( k = SIZE_INT(opL)-1; k >= 0; k-- )
1974
            *m++ = *l++;
1975

1976
        /* get the size of the right operand, and get the leading 2 digits */
1977
        rs = SIZE_INT(opR);
1978
        r  = ADDR_INT(opR);
1979
        while ( r[rs-1] == 0 )  rs--;
1980
        for ( e = 0;
1981
             ((Int)r[rs-1]<<e) + (e ?  r[rs-2]>>(NR_DIGIT_BITS-e) : 0)
1982
                               < INTBASE/2; e++ ) ;
1983

1984
        r1 = ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e) : 0);
1985
        r2 = ((Int)r[rs-2]<<e) + ((rs>=3 && e)? r[rs-3]>>(NR_DIGIT_BITS-e) : 0);
1986

1987
        /* run through the digits in the quotient                          */
1988
        for ( i = SIZE_INT(mod)-SIZE_INT(opR)-1; i >= 0; i-- ) {
1989

1990
            /* guess the factor                                            */
1991
            m = ADDR_INT(mod) + rs + i;
1992
            m01 = ((INTBASE*m[0]+m[-1])<<e)
1993
                           + (e ? m[-2]>>(NR_DIGIT_BITS-e) : 0);
1994
            if ( m01 == 0 )  continue;
1995
            m2  = ((Int)m[-2]<<e) + ((e && rs+i>=3) ? m[-3]>>(NR_DIGIT_BITS-e) : 0);
1996
            if ( ((Int)m[0]<<e)+(e ? m[-1]>>(NR_DIGIT_BITS-e) : 0) < r1 )
1997
                    qi = m01 / r1;
1998
            else    qi = INTBASE - 1;
1999
            while ( m01-(Int)qi*r1 < (Int)INTBASE
2000
                    && INTBASE*(m01-(Int)qi*r1)+m2 < (Int)qi*r2 )
2001
                    qi--;
2002

2003
            /* m = m - qi * r;                                             */
2004
            d = 0;
2005
            m = ADDR_INT(mod) + i;
2006
            r = ADDR_INT(opR);
2007
            for ( k = 0; k < (Int)rs; ++k, ++m, ++r ) {
2008
                c = (Int)*m - (Int)qi * *r - (Int)d;
2009
		*m = (TypDigit)c;
2010
		d = -(TypDigit)(c>>NR_DIGIT_BITS);
2011
            }
2012
            c = (Int)*m - d;  *m = (TypDigit)c;  d = -(TypDigit)(c>>NR_DIGIT_BITS);
2013

2014
            /* if we have a borrow then add back                           */
2015
            if ( d != 0 ) {
2016
                d = 0;
2017
                m = ADDR_INT(mod) + i;
2018
                r = ADDR_INT(opR);
2019
                for ( k = 0; k < (Int)rs; ++k, ++m, ++r ) {
2020
                    c = (Int)*m + (Int)*r + (Int)d;
2021
                    *m = (TypDigit)c;
2022
                    d = (TypDigit)(c>>NR_DIGIT_BITS);
2023
                }
2024
                c = (Int)*m + d;  *m = (TypDigit)c;  d = (TypDigit)(c>>NR_DIGIT_BITS);
2025
                qi--;
2026
            }
2027

2028
        }
2029

2030
        /* remove the leading zeroes                                       */
2031
        m = ADDR_INT(mod) + SIZE_INT(mod);
2032
        if ( (m[-4] == 0 && m[-3] == 0 && m[-2] == 0 && m[-1] == 0)
2033
          || (SIZE_INT(mod) == 4 && m[-2] == 0 && m[-1] == 0) ) {
2034

2035
            /* find the number of significant digits                       */
2036
            m = ADDR_INT(mod);
2037
            for ( k = SIZE_INT(mod); k != 0; k-- ) {
2038
                if ( m[k-1] != 0 )
2039
                    break;
2040
            }
2041

2042
            /* reduce to small integer if possible, otherwise shrink bag   */
2043

2044
            if ( k <= 2 && TNUM_OBJ(mod) == T_INTPOS
2045
              && (UInt)(INTBASE*m[1]+m[0])<(1L<<NR_SMALL_INT_BITS) )
2046
                mod = INTOBJ_INT( INTBASE*m[1]+m[0] );
2047
            else if ( k <= 2 && TNUM_OBJ(mod) == T_INTNEG
2048
              && (UInt)(INTBASE*m[1]+m[0])<=(1L<<NR_SMALL_INT_BITS) )
2049
                mod = INTOBJ_INT( -(Int)(INTBASE*m[1]+m[0]) );
2050
            else
2051
                ResizeBag( mod, (((k + 3) / 4) * 4) * sizeof(TypDigit) );
2052
        }
2053

2054
        /* make the representative positive                                  */
2055
        if ( (TNUM_OBJ(mod) == T_INT && INT_INTOBJ(mod) < 0)
2056
          || TNUM_OBJ(mod) == T_INTNEG ) {
2057
            if ( TNUM_OBJ(opR) == T_INTPOS )
2058
                mod = SumInt( mod, opR );
2059
            else
2060
                mod = DiffInt( mod, opR );
2061
        }
2062

2063
    }
2064

2065
    /* return the result                                                   */
2066
    return mod;
2067
}
2068

2069

2070
/****************************************************************************
2071
**
2072
*F  FuncIS_INT( <self>, <val> ) . . . . . . . . . . internal function 'IsInt'
2073
**
2074
**  'FuncIS_INT' implements the internal filter 'IsInt'.
2075
**
2076
**  'IsInt( <val> )'
2077
**
2078
**  'IsInt'  returns 'true'  if the  value  <val>  is an integer  and 'false'
2079
**  otherwise.
2080
*/
2081
Obj IsIntFilt;
2082

2083
Obj FuncIS_INT (
2084
    Obj                 self,
2085
    Obj                 val )
2086
{
2087
    /* return 'true' if <obj> is an integer and 'false' otherwise          */
2088
    if ( TNUM_OBJ(val) == T_INT
2089
      || TNUM_OBJ(val) == T_INTPOS
2090
      || TNUM_OBJ(val) == T_INTNEG ) {
2091
        return True;
2092
    }
2093
    else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
2094
        return False;
2095
    }
2096
    else {
2097
        return DoFilter( self, val );
2098
    }
2099
}
2100

2101

2102
/****************************************************************************
2103
**
2104
*F  QuoInt( <intL>, <intR> )  . . . . . . . . . . . quotient of two integers
2105
**
2106
**  'QuoInt' returns the integer part of the two integers <intL> and  <intR>.
2107
**  'QuoInt' handles operands of type  'T_INT',  'T_INTPOS'  and  'T_INTNEG'.
2108
**
2109
**  It can also be used in the cases that both operands  are  small  integers
2110
**  and the result is a small integer too,  i.e., that  no  overflow  occurs.
2111
**
2112
**  Note that this routine is not called from 'EvalQuo', the  division  of  two
2113
**  integers yields  a  rational  and  is  therefor  performed  in  'QuoRat'.
2114
**  This operation is however available through the internal function 'Quo'.
2115
*/
2116
Obj             QuoInt (
2117
    Obj                 opL,
2118
    Obj                 opR )
2119
{
2120
    Int                 i;              /* loop count, value for small int */
2121
    Int                 k;              /* loop count, value for small int */
2122
    UInt                c;              /* product of two digits           */
2123
    TypDigit            d;              /* carry into the next digit       */
2124
    TypDigit *          l;              /* pointer into the left operand   */
2125
    UInt                l01;            /* leading two digits of the left  */
2126
    TypDigit            l2;             /* next digit of the left operand  */
2127
    TypDigit *          r;              /* pointer into the right operand  */
2128
    TypDigit            r1;             /* leading digit of the right oper */
2129
    TypDigit            r2;             /* next digit of the right operand */
2130
    UInt                rs;             /* size of the right operand       */
2131
    UInt                e;              /* we mult r by 2^e so r1 >= 32768 */
2132
    TypDigit *          q;              /* pointer into the quotient       */
2133
    Obj                 quo;            /* handle of the result bag        */
2134
    TypDigit            qi;             /* guessed digit of the quotient   */
2135

2136
    /* divide to small integers                                            */
2137
    if ( ARE_INTOBJS( opL, opR ) ) {
2138

2139
        /* pathological case first                                         */
2140
        if ( opR == INTOBJ_INT(0) ) {
2141
            opR = ErrorReturnObj(
2142
                "Integer operations: <divisor> must be nonzero",
2143
                0L, 0L,
2144
                "you can replace the integer <divisor> via 'return <divisor>;'" );
2145
            return QUO( opL, opR );
2146
        }
2147

2148
        /* the small int -(1<<28) divided by -1 is the large int (1<<28)   */
2149
        if ( opL == INTOBJ_INT(-(Int)(1L<<NR_SMALL_INT_BITS))
2150
            && opR == INTOBJ_INT(-1) ) {
2151
            quo = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2152
            ADDR_INT(quo)[1] = 1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS);
2153
            ADDR_INT(quo)[0] = 0;
2154
            return quo;
2155
        }
2156

2157
        /* get the integer values                                          */
2158
        i = INT_INTOBJ(opL);
2159
        k = INT_INTOBJ(opR);
2160

2161
        /* divide, make sure we divide only positive numbers               */
2162
        if (      0 <= i && 0 <= k )  i =    (  i /  k );
2163
        else if ( 0 <= i && k <  0 )  i =  - (  i / -k );
2164
        else if ( i < 0  && 0 <= k )  i =  - ( -i /  k );
2165
        else if ( i < 0  && k <  0 )  i =    ( -i / -k );
2166
        quo = INTOBJ_INT( i );
2167

2168
    }
2169

2170
    /* divide a small integer by a large one                               */
2171
    else if ( IS_INTOBJ(opL) ) {
2172

2173
        /* the small int -(1<<28) divided by the large int (1<<28) is -1   */
2174

2175
        if ( opL == INTOBJ_INT(-(Int)(1L<<NR_SMALL_INT_BITS))
2176
          && TNUM_OBJ(opR) == T_INTPOS && SIZE_INT(opR) == 4
2177
          && ADDR_INT(opR)[3] == 0
2178
          && ADDR_INT(opR)[2] == 0
2179
          && ADDR_INT(opR)[1] == 1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS)
2180
          && ADDR_INT(opR)[0] == 0 )
2181
            quo = INTOBJ_INT(-1);
2182

2183
        /* in all other cases the quotient is of course zero               */
2184
        else
2185
            quo = INTOBJ_INT(0);
2186

2187
    }
2188

2189
    /* divide a large integer by a small integer                           */
2190
    else if ( IS_INTOBJ(opR)
2191
           && INT_INTOBJ(opR) < INTBASE
2192
           && -(Int)INTBASE  <= INT_INTOBJ(opR) ) {
2193

2194
        /* pathological case first                                         */
2195
        if ( opR == INTOBJ_INT(0) ) {
2196
            opR = ErrorReturnObj(
2197
                "Integer operations: <divisor> must be nonzero",
2198
                0L, 0L,
2199
                "you can replace the integer <divisor> via 'return <divisor>;'" );
2200
            return QUO( opL, opR );
2201
        }
2202

2203
        /* get the integer value, make positive                            */
2204
        i = INT_INTOBJ(opR);  if ( i < 0 )  i = -i;
2205

2206
        /* allocate a bag for the result and set up the pointers           */
2207
        if ( (TNUM_OBJ(opL)==T_INTPOS && 0 < INT_INTOBJ(opR))
2208
          || (TNUM_OBJ(opL)==T_INTNEG && INT_INTOBJ(opR) < 0) )
2209
            quo = NewBag( T_INTPOS, SIZE_OBJ(opL) );
2210
        else
2211
            quo = NewBag( T_INTNEG, SIZE_OBJ(opL) );
2212
        l = ADDR_INT(opL) + SIZE_INT(opL) - 1;
2213
        q = ADDR_INT(quo) + SIZE_INT(quo) - 1;
2214

2215
        /* run through the left operand and divide digitwise               */
2216
        c = 0;
2217
        for ( ; l >= ADDR_INT(opL); l--, q-- ) {
2218
            c  = (c<<NR_DIGIT_BITS) + (Int)*l;
2219
            *q = (TypDigit)(c / i);
2220
            c  = c - i * *q;
2221
            /*N clever compilers may prefer:  c  = c % i;                  */
2222
        }
2223

2224
        /* remove the leading zeroes, note that there can't be more than 5 */
2225
        q = ADDR_INT(quo) + SIZE_INT(quo);
2226
        if ( q[-4] == 0 && q[-3] == 0 && q[-2] == 0 && q[-1] == 0 ) {
2227
            ResizeBag( quo, (SIZE_INT(quo)-4)*sizeof(TypDigit) );
2228
        }
2229

2230
        /* reduce to small integer if possible                             */
2231
        q = ADDR_INT(quo) + SIZE_INT(quo);
2232
        if ( SIZE_INT(quo) == 4 && q[-2] == 0 && q[-1] == 0 ) {
2233
            if ( TNUM_OBJ(quo) == T_INTPOS
2234
              &&(UInt)(INTBASE*q[-3]+q[-4])
2235
                < (1L<<NR_SMALL_INT_BITS))
2236
                quo = INTOBJ_INT( INTBASE*q[-3]+q[-4] );
2237
            else if ( TNUM_OBJ(quo) == T_INTNEG
2238
              && (UInt)(INTBASE*q[-3]+q[-4])
2239
                     <= (1L<<NR_SMALL_INT_BITS) )
2240
                quo = INTOBJ_INT( -(Int)(INTBASE*q[-3]+q[-4]) );
2241
        }
2242

2243
    }
2244

2245
    /* divide a large integer by a large integer                           */
2246
    else {
2247

2248
        /* a small divisor larger than one digit isn't handled above       */
2249
        if ( IS_INTOBJ(opR) ) {
2250
            if ( 0 < INT_INTOBJ(opR) ) {
2251
                quo = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2252
                ADDR_INT(quo)[0] = (TypDigit)(INT_INTOBJ(opR));
2253
                ADDR_INT(quo)[1] = (TypDigit)(INT_INTOBJ(opR)>>NR_DIGIT_BITS);
2254
                opR = quo;
2255
            }
2256
            else {
2257
                quo = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
2258
                ADDR_INT(quo)[0] = (TypDigit)(-INT_INTOBJ(opR));
2259
                ADDR_INT(quo)[1] = (TypDigit)((-INT_INTOBJ(opR))>>NR_DIGIT_BITS);
2260
                opR = quo;
2261
            }
2262
        }
2263

2264
        /* trivial case first                                              */
2265
        if ( SIZE_INT(opL) < SIZE_INT(opR) )
2266
            return INTOBJ_INT(0);
2267

2268
        /* copy the left operand into a new bag, this holds the remainder  */
2269
        quo = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+4)*sizeof(TypDigit) );
2270
        l = ADDR_INT(opL);
2271
        q = ADDR_INT(quo);
2272
        for ( k = SIZE_INT(opL)-1; k >= 0; k-- )
2273
            *q++ = *l++;
2274
        opL = quo;
2275

2276
        /* get the size of the right operand, and get the leading 2 digits */
2277
        rs = SIZE_INT(opR);
2278
        r  = ADDR_INT(opR);
2279
        while ( r[rs-1] == 0 )  rs--;
2280
        for ( e = 0;
2281
             ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e) : 0)
2282
                               < INTBASE/2 ; e++ ) ;
2283

2284
        r1 = ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e) : 0);
2285
        r2 = ((Int)r[rs-2]<<e) + ((e && rs>=3) ? r[rs-3]>>(NR_DIGIT_BITS-e) : 0);
2286

2287
        /* allocate a bag for the quotient                                 */
2288
        if ( TNUM_OBJ(opL) == TNUM_OBJ(opR) )
2289
            quo = NewBag( T_INTPOS, SIZE_OBJ(opL)-SIZE_OBJ(opR) );
2290
        else
2291
            quo = NewBag( T_INTNEG, SIZE_OBJ(opL)-SIZE_OBJ(opR) );
2292

2293
        /* run through the digits in the quotient                          */
2294
        for ( i = SIZE_INT(opL)-SIZE_INT(opR)-1; i >= 0; i-- ) {
2295

2296
            /* guess the factor                                            */
2297
            l = ADDR_INT(opL) + rs + i;
2298
            l01 = ((INTBASE*l[0]+l[-1])<<e)
2299
              + (e ? l[-2]>>(NR_DIGIT_BITS-e):0);
2300

2301
            if ( l01 == 0 )  continue;
2302
            l2  = ((Int)l[-2]<<e) + (( e && rs+i>=3) ? l[-3]>>(NR_DIGIT_BITS-e) : 0);
2303
            if ( ((Int)l[0]<<e)+(e ? l[-1]>>(NR_DIGIT_BITS-e):0) < r1 )
2304
                     qi = l01 / r1;
2305
            else     qi = INTBASE - 1;
2306
            while ( l01-(Int)qi*r1 < (Int)INTBASE
2307
                    && INTBASE*(l01-(UInt)qi*r1)+l2 < (UInt)qi*r2 )
2308
                qi--;
2309

2310
            /* l = l - qi * r;                                             */
2311
            d = 0;
2312
            l = ADDR_INT(opL) + i;
2313
            r = ADDR_INT(opR);
2314
            for ( k = 0; k < (Int)rs; ++k, ++l, ++r ) {
2315
                c = *l - (Int)qi * *r - d;  *l = c;  d = -(TypDigit)(c>>NR_DIGIT_BITS);
2316
            }
2317
            c = (Int)*l - d; d = -(TypDigit)(c>>NR_DIGIT_BITS);
2318

2319
            /* if we have a borrow then add back                           */
2320
            if ( d != 0 ) {
2321
                d = 0;
2322
                l = ADDR_INT(opL) + i;
2323
                r = ADDR_INT(opR);
2324
                for ( k = 0; k < (Int)rs; ++k, ++l, ++r ) {
2325
                    c = (Int)*l + (Int)*r + (Int)d;
2326
                    *l = (TypDigit)c;
2327
                    d = (TypDigit)(c>>NR_DIGIT_BITS);
2328
                }
2329
                c = *l + d; d = (TypDigit)(c>>NR_DIGIT_BITS);
2330
                qi--;
2331
            }
2332

2333
            /* store the digit in the quotient                             */
2334
            ADDR_INT(quo)[i] = qi;
2335

2336
        }
2337

2338
        /* remove the leading zeroes, note that there can't be more than 7 */
2339
        q = ADDR_INT(quo) + SIZE_INT(quo);
2340
        if ( SIZE_INT(quo) > 4
2341
          && q[-4] == 0 && q[-3] == 0 && q[-2] == 0 && q[-1] == 0 ) {
2342
            ResizeBag( quo, (SIZE_INT(quo)-4)*sizeof(TypDigit) );
2343
        }
2344

2345
        /* reduce to small integer if possible                             */
2346
        q = ADDR_INT(quo) + SIZE_INT(quo);
2347
        if ( SIZE_INT(quo) == 4 && q[-2] == 0 && q[-1] == 0 ) {
2348
            if ( TNUM_OBJ(quo) == T_INTPOS
2349
              && (UInt)(INTBASE*q[-3]+q[-4]) < (1L<<NR_SMALL_INT_BITS) )
2350
                quo = INTOBJ_INT( INTBASE*q[-3]+q[-4] );
2351
            else if ( TNUM_OBJ(quo) == T_INTNEG
2352
              && (UInt)(INTBASE*q[-3]+q[-4]) <= (1L<<NR_SMALL_INT_BITS) )
2353
                quo = INTOBJ_INT( -(Int)(INTBASE*q[-3]+q[-4]) );
2354
        }
2355

2356
    }
2357

2358
    /* return the result                                                   */
2359
    return quo;
2360
}
2361

2362

2363
/****************************************************************************
2364
**
2365
*F  FuncQUO_INT(<self>,<opL>,<opR>) . . . . . . .  internal function 'QuoInt'
2366
**
2367
**  'FuncQUO_INT' implements the internal function 'QuoInt'.
2368
**
2369
**  'QuoInt( <i>, <k> )'
2370
**
2371
**  'Quo' returns the  integer part of the quotient  of its integer operands.
2372
**  If <i>  and <k> are  positive 'Quo( <i>,  <k> )' is  the largest positive
2373
**  integer <q>  such that '<q> * <k>  \<= <i>'.  If  <i> or  <k> or both are
2374
**  negative we define 'Abs( Quo(<i>,<k>) ) = Quo( Abs(<i>), Abs(<k>) )'  and
2375
**  'Sign( Quo(<i>,<k>) ) = Sign(<i>) * Sign(<k>)'.  Dividing by 0  causes an
2376
**  error.  'Rem' (see "Rem") can be used to compute the remainder.
2377
*/
2378
Obj FuncQUO_INT (
2379
    Obj                 self,
2380
    Obj                 opL,
2381
    Obj                 opR )
2382
{
2383
    /* check the arguments                                                 */
2384
    while ( TNUM_OBJ(opL) != T_INT
2385
         && TNUM_OBJ(opL) != T_INTPOS
2386
         && TNUM_OBJ(opL) != T_INTNEG ) {
2387
        opL = ErrorReturnObj(
2388
            "QuoInt: <left> must be an integer (not a %s)",
2389
            (Int)TNAM_OBJ(opL), 0L,
2390
            "you can replace <left> via 'return <left>;'" );
2391
    }
2392
    while ( TNUM_OBJ(opR) != T_INT
2393
         && TNUM_OBJ(opR) != T_INTPOS
2394
         && TNUM_OBJ(opR) != T_INTNEG ) {
2395
        opR = ErrorReturnObj(
2396
            "QuoInt: <right> must be an integer (not a %s)",
2397
            (Int)TNAM_OBJ(opR), 0L,
2398
            "you can replace <right> via 'return <right>;'" );
2399
    }
2400

2401
    /* return the quotient                                                 */
2402
    return QuoInt( opL, opR );
2403
}
2404

2405

2406
/****************************************************************************
2407
**
2408
*F  RemInt( <intL>, <intR> )  . . . . . . . . . . . remainder of two integers
2409
**
2410
**  'RemInt' returns the remainder of the quotient  of  the  integers  <intL>
2411
**  and <intR>.  'RemInt' handles operands of type  'T_INT',  'T_INTPOS'  and
2412
**  'T_INTNEG'.
2413
**
2414
**  Note that the remainder is different from the value returned by the 'mod'
2415
**  operator which is always positive.
2416
**
2417
**  'RemInt' is called from 'FunRemInt'.
2418
*/
2419
Obj             RemInt (
2420
    Obj                 opL,
2421
    Obj                 opR )
2422
{
2423
    Int                 i;              /* loop count, value for small int */
2424
    Int                 k;              /* loop count, value for small int */
2425
    UInt                c;              /* product of two digits           */
2426
    TypDigit            d;              /* carry into the next digit       */
2427
    TypDigit *          l;              /* pointer into the left operand   */
2428
    TypDigit *          r;              /* pointer into the right operand  */
2429
    TypDigit            r1;             /* leading digit of the right oper */
2430
    TypDigit            r2;             /* next digit of the right operand */
2431
    UInt                rs;             /* size of the right operand       */
2432
    UInt                e;              /* we mult r by 2^e so r1 >= 32768 */
2433
    Obj                 rem;            /* handle of the remainder bag     */
2434
    TypDigit *          m;              /* pointer into the remainder      */
2435
    UInt                m01;            /* leading two digits of the rem.  */
2436
    TypDigit            m2;             /* next digit of the remainder     */
2437
    TypDigit            qi;             /* guessed digit of the quotient   */
2438

2439
    /* compute the remainder of two small integers                         */
2440
    if ( ARE_INTOBJS( opL, opR ) ) {
2441

2442
        /* pathological case first                                         */
2443
        if ( opR == INTOBJ_INT(0) ) {
2444
            opR = ErrorReturnObj(
2445
                "Integer operations: <divisor> must be nonzero",
2446
                0L, 0L,
2447
                "you can replace the integer <divisor> via 'return <divisor>;'" );
2448
            return QUO( opL, opR );
2449
        }
2450

2451
        /* get the integer values                                          */
2452
        i = INT_INTOBJ(opL);
2453
        k = INT_INTOBJ(opR);
2454

2455
        /* compute the remainder, make sure we divide only positive numbers*/
2456
        if (      0 <= i && 0 <= k )  i =    (  i %  k );
2457
        else if ( 0 <= i && k <  0 )  i =    (  i % -k );
2458
        else if ( i < 0  && 0 <= k )  i =  - ( -i %  k );
2459
        else if ( i < 0  && k <  0 )  i =  - ( -i % -k );
2460
        rem = INTOBJ_INT( i );
2461

2462
    }
2463

2464
    /* compute the remainder of a small integer by a large integer         */
2465
    else if ( IS_INTOBJ(opL) ) {
2466

2467
        /* the small int -(1<<28) rem the large int (1<<28) is 0           */
2468
        if ( opL == INTOBJ_INT(-(Int)(1L<<NR_SMALL_INT_BITS))
2469
          && TNUM_OBJ(opR) == T_INTPOS && SIZE_INT(opR) == 4
2470
          && ADDR_INT(opR)[3] == 0
2471
          && ADDR_INT(opR)[2] == 0
2472
          && ADDR_INT(opR)[1] == 1L << (NR_SMALL_INT_BITS-NR_DIGIT_BITS)
2473
          && ADDR_INT(opR)[0] == 0 )
2474
            rem = INTOBJ_INT(0);
2475

2476
        /* in all other cases the remainder is equal the left operand      */
2477
        else
2478
            rem = opL;
2479

2480
    }
2481

2482
    /* compute the remainder of a large integer by a small integer         */
2483
    else if ( IS_INTOBJ(opR)
2484
           && INT_INTOBJ(opR) < INTBASE
2485
           && -(Int)INTBASE <= INT_INTOBJ(opR) ) {
2486

2487
        /* pathological case first                                         */
2488
        if ( opR == INTOBJ_INT(0) ) {
2489
            opR = ErrorReturnObj(
2490
                "Integer operations: <divisor> must be nonzero",
2491
                0L, 0L,
2492
                "you can replace the integer <divisor> via 'return <divisor>;'" );
2493
            return QUO( opL, opR );
2494
        }
2495

2496
        /* get the integer value, make positive                            */
2497
        i = INT_INTOBJ(opR);  if ( i < 0 )  i = -i;
2498

2499
        /* maybe its trivial                                               */
2500
        if ( INTBASE % i == 0 ) {
2501
            c = ADDR_INT(opL)[0] % i;
2502
        }
2503

2504
        /* otherwise run through the left operand and divide digitwise     */
2505
        else {
2506
            l = ADDR_INT(opL) + SIZE_INT(opL) - 1;
2507
            c = 0;
2508
            for ( ; l >= ADDR_INT(opL); l-- ) {
2509
                c  = (c<<NR_DIGIT_BITS) + (Int)*l;
2510
                c  = c % i;
2511
            }
2512
        }
2513

2514
        /* now c is the result, it has the same sign as the left operand   */
2515
        if ( TNUM_OBJ(opL) == T_INTPOS )
2516
            rem = INTOBJ_INT(  c );
2517
        else
2518
            rem = INTOBJ_INT( -(Int)c );
2519

2520
    }
2521

2522
    /* compute the remainder of a large integer modulo a large integer     */
2523
    else {
2524

2525
        /* a small divisor larger than one digit isn't handled above       */
2526
        if ( IS_INTOBJ(opR) ) {
2527
            if ( 0 < INT_INTOBJ(opR) ) {
2528
                rem = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2529
                ADDR_INT(rem)[0] = (TypDigit)(INT_INTOBJ(opR));
2530
                ADDR_INT(rem)[1] = (TypDigit)(INT_INTOBJ(opR)>>NR_DIGIT_BITS);
2531
                opR = rem;
2532
            }
2533
            else {
2534
                rem = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
2535
                ADDR_INT(rem)[0] = (TypDigit)(-INT_INTOBJ(opR));
2536
                ADDR_INT(rem)[1] = (TypDigit)((-INT_INTOBJ(opR))>>NR_DIGIT_BITS);
2537
                opR = rem;
2538
            }
2539
        }
2540

2541
        /* trivial case first                                              */
2542
        if ( SIZE_INT(opL) < SIZE_INT(opR) )
2543
            return opL;
2544

2545
        /* copy the left operand into a new bag, this holds the remainder  */
2546
        rem = NewBag( TNUM_OBJ(opL), (SIZE_INT(opL)+4)*sizeof(TypDigit) );
2547
        l = ADDR_INT(opL);
2548
        m = ADDR_INT(rem);
2549
        for ( k = SIZE_INT(opL)-1; k >= 0; k-- )
2550
            *m++ = *l++;
2551

2552
        /* get the size of the right operand, and get the leading 2 digits */
2553
        rs = SIZE_INT(opR);
2554
        r  = ADDR_INT(opR);
2555
        while ( r[rs-1] == 0 )  rs--;
2556
        for ( e = 0;
2557
             ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e): 0)
2558
                               < INTBASE/2; e++ ) ;
2559

2560
        r1 = ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e):0);
2561
        r2 = ((Int)r[rs-2]<<e) + ((e && rs>=3) ? r[rs-3]>>(NR_DIGIT_BITS-e) : 0);
2562

2563
        /* run through the digits in the quotient                          */
2564
        for ( i = SIZE_INT(rem)-SIZE_INT(opR)-1; i >= 0; i-- ) {
2565

2566
            /* guess the factor                                            */
2567
            m = ADDR_INT(rem) + rs + i;
2568
            m01 = ((INTBASE*m[0]+m[-1])<<e) + (e ? m[-2]>>(NR_DIGIT_BITS-e):0);
2569
            if ( m01 == 0 )  continue;
2570
            m2  = ((Int)m[-2]<<e) + ((e && rs+i>=3) ? m[-3]>>(NR_DIGIT_BITS-e) : 0);
2571
            if ( ((Int)m[0]<<e)+(e ? m[-1]>>(NR_DIGIT_BITS-e): 0) < r1 )
2572
                    qi = m01 / r1;
2573
            else    qi = INTBASE - 1;
2574
            while ( m01-(Int)qi*r1 < (Int)INTBASE
2575
                   && INTBASE*(m01-(Int)qi*r1)+m2 < (Int)qi*r2 )
2576
                qi--;
2577

2578
            /* m = m - qi * r;                                             */
2579
            d = 0;
2580
            m = ADDR_INT(rem) + i;
2581
            r = ADDR_INT(opR);
2582
            for ( k = 0; k < (Int)rs; ++k, ++m, ++r ) {
2583
                c = (Int)*m - (Int)qi * *r - (Int)d;
2584
                *m = (TypDigit)c;
2585
                d = -(TypDigit)(c>>NR_DIGIT_BITS);
2586
            }
2587
            c = (Int)*m - d;  *m = (TypDigit)c;  d = -(TypDigit)(c>>NR_DIGIT_BITS);
2588

2589
            /* if we have a borrow then add back                           */
2590
            if ( d != 0 ) {
2591
                d = 0;
2592
                m = ADDR_INT(rem) + i;
2593
                r = ADDR_INT(opR);
2594
                for ( k = 0; k < (Int)rs; ++k, ++m, ++r ) {
2595
                    c = (Int)*m + (Int)*r + (Int)d;
2596
                    *m = (TypDigit)c;
2597
                    d = (TypDigit)(c>>NR_DIGIT_BITS);
2598
                }
2599
                c = (Int)*m + d;  *m = (TypDigit)c;  d = (TypDigit)(c>>NR_DIGIT_BITS);
2600
                qi--;
2601
            }
2602

2603
        }
2604

2605
        /* remove the leading zeroes                                       */
2606
        m = ADDR_INT(rem) + SIZE_INT(rem);
2607
        if ( (m[-4] == 0 && m[-3] == 0 && m[-2] == 0 && m[-1] == 0)
2608
          || (SIZE_INT(rem) == 4 && m[-2] == 0 && m[-1] == 0) ) {
2609

2610
            /* find the number of significant digits                       */
2611
            m = ADDR_INT(rem);
2612
            for ( k = SIZE_INT(rem); k != 0; k-- ) {
2613
                if ( m[k-1] != 0 )
2614
                    break;
2615
            }
2616

2617
            /* reduce to small integer if possible, otherwise shrink bag   */
2618
            if ( k <= 2 && TNUM_OBJ(rem) == T_INTPOS
2619
              && (UInt)(INTBASE*m[1]+m[0]) < (1L<<NR_SMALL_INT_BITS) )
2620
                rem = INTOBJ_INT( INTBASE*m[1]+m[0] );
2621
            else if ( k <= 2 && TNUM_OBJ(rem) == T_INTNEG
2622
              && (UInt)(INTBASE*m[1]+m[0]) <= (1L<<NR_SMALL_INT_BITS) )
2623
                rem = INTOBJ_INT( -(Int)(INTBASE*m[1]+m[0]) );
2624
            else
2625
                ResizeBag( rem, (((k + 3) / 4) * 4) * sizeof(TypDigit) );
2626
        }
2627

2628
    }
2629

2630
    /* return the result                                                   */
2631
    return rem;
2632
}
2633

2634

2635
/****************************************************************************
2636
**
2637
*F  FuncREM_INT(<self>,<opL>,<opR>)  . . . . . . .  internal function 'RemInt'
2638
**
2639
**  'FuncREM_INT' implements the internal function 'RemInt'.
2640
**
2641
**  'RemInt( <i>, <k> )'
2642
**
2643
**  'Rem' returns the remainder of its two integer operands,  i.e., if <k> is
2644
**  not equal to zero 'Rem( <i>, <k> ) = <i> - <k> *  Quo( <i>, <k> )'.  Note
2645
**  that the rules given  for 'Quo' (see "Quo") imply  that 'Rem( <i>, <k> )'
2646
**  has the same sign as <i> and its absolute value is strictly less than the
2647
**  absolute value of <k>.  Dividing by 0 causes an error.
2648
*/
2649
Obj             FuncREM_INT (
2650
    Obj                 self,
2651
    Obj                 opL,
2652
    Obj                 opR )
2653
{
2654
    /* check the arguments                                                 */
2655
    while ( TNUM_OBJ(opL) != T_INT
2656
         && TNUM_OBJ(opL) != T_INTPOS
2657
         && TNUM_OBJ(opL) != T_INTNEG ) {
2658
        opL = ErrorReturnObj(
2659
            "RemInt: <left> must be an integer (not a %s)",
2660
            (Int)TNAM_OBJ(opL), 0L,
2661
            "you can replace <left> via 'return <left>;'" );
2662
    }
2663
    while ( TNUM_OBJ(opR) != T_INT
2664
         && TNUM_OBJ(opR) != T_INTPOS
2665
         && TNUM_OBJ(opR) != T_INTNEG ) {
2666
        opR = ErrorReturnObj(
2667
            "RemInt: <right> must be an integer (not a %s)",
2668
            (Int)TNAM_OBJ(opR), 0L,
2669
            "you can replace <right> via 'return <right>;'" );
2670
    }
2671

2672
    /* return the remainder                                                */
2673
    return RemInt( opL, opR );
2674
}
2675

2676

2677
/****************************************************************************
2678
**
2679
*F  GcdInt( <opL>, <opR> )  . . . . . . . . . . . . . . . gcd of two integers
2680
**
2681
**  'GcdInt' returns the gcd of the two integers <opL> and <opR>.
2682
**
2683
**  It is called from 'FunGcdInt' and the rational package.
2684
*/
2685
Obj             GcdInt (
2686
    Obj                 opL,
2687
    Obj                 opR )
2688
{
2689
    Int                 i;              /* loop count, value for small int */
2690
    Int                 k;              /* loop count, value for small int */
2691
    UInt                c;              /* product of two digits           */
2692
    TypDigit            d;              /* carry into the next digit       */
2693
    TypDigit *          r;              /* pointer into the right operand  */
2694
    TypDigit            r1;             /* leading digit of the right oper */
2695
    TypDigit            r2;             /* next digit of the right operand */
2696
    UInt                rs;             /* size of the right operand       */
2697
    UInt                e;              /* we mult r by 2^e so r1 >= 32768 */
2698
    TypDigit *          l;              /* pointer into the left operand   */
2699
    UInt                l01;            /* leading two digits of the rem.  */
2700
    TypDigit            l2;             /* next digit of the remainder     */
2701
    UInt                ls;             /* size of the left operand        */
2702
    TypDigit            qi;             /* guessed digit of the quotient   */
2703
    Obj                 gcd;            /* handle of the result            */
2704

2705
    /* compute the gcd of two small integers                               */
2706
    if ( ARE_INTOBJS( opL, opR ) ) {
2707

2708
        /* get the integer values, make them positive                      */
2709
        i = INT_INTOBJ(opL);  if ( i < 0 )  i = -i;
2710
        k = INT_INTOBJ(opR);  if ( k < 0 )  k = -k;
2711

2712
        /* compute the gcd using Euclids algorithm                         */
2713
        while ( k != 0 ) {
2714
            c = k;
2715
            k = i % k;
2716
            i = c;
2717
        }
2718

2719
        /* now i is the result                                             */
2720
        if ( i == (1L<<NR_SMALL_INT_BITS) ) {
2721
            gcd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2722
            ADDR_INT(gcd)[0] = (TypDigit)i;
2723
            ADDR_INT(gcd)[1] = (TypDigit)(i>>NR_DIGIT_BITS);
2724
        }
2725
        else {
2726
            gcd = INTOBJ_INT( i );
2727
        }
2728

2729
    }
2730

2731
    /* compute the gcd of a small and a large integer                      */
2732
    else if ( (IS_INTOBJ(opL)
2733
               && INT_INTOBJ(opL) < INTBASE && -(Int)INTBASE <= INT_INTOBJ(opL))
2734
           || (IS_INTOBJ(opR)
2735
               && INT_INTOBJ(opR) < INTBASE && -(Int)INTBASE <= INT_INTOBJ(opR)) ) {
2736

2737
        /* make the right operand the small one                            */
2738
        if ( IS_INTOBJ(opL) ) {
2739
            gcd = opL;  opL = opR;  opR = gcd;
2740
        }
2741

2742
        /* maybe it's trivial                                              */
2743
        if ( opR == INTOBJ_INT(0) ) {
2744
            if( TNUM_OBJ( opL ) == T_INTNEG ) {
2745
                /* If opL is negative, change the sign.  We do this by    */
2746
                /* copying opL into a bag of type T_INTPOS.  Note that    */
2747
                /* opL is a large negative number, so it cannot be the    */
2748
                /* the negative of 1 << NR_SMALL_INT_BITS.                */
2749
                gcd = NewBag( T_INTPOS, SIZE_OBJ(opL) );
2750
                l = ADDR_INT(opL); r = ADDR_INT(gcd);
2751
                for ( k = SIZE_INT(opL); k != 0; k-- ) *r++ = *l++;
2752

2753
                return gcd;
2754
            }
2755
            else return opL;
2756
        }
2757

2758
        /* get the right operand value, make it positive                   */
2759
        i = INT_INTOBJ(opR);  if ( i < 0 )  i = -i;
2760

2761
        /* do one remainder operation                                      */
2762
        l = ADDR_INT(opL) + SIZE_INT(opL) - 1;
2763
        c = 0;
2764
        for ( ; l >= ADDR_INT(opL); l-- ) {
2765
            c  = (c<<NR_DIGIT_BITS) + (Int)*l;
2766
            c  = c % i;
2767
        }
2768
        k = c;
2769

2770
        /* compute the gcd using Euclids algorithm                         */
2771
        while ( k != 0 ) {
2772
            c = k;
2773
            k = i % k;
2774
            i = c;
2775
        }
2776

2777
        /* now i is the result                                             */
2778
        if ( i == (1L<<NR_SMALL_INT_BITS) ) {
2779
            gcd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2780
            ADDR_INT(gcd)[0] = 0;
2781
            ADDR_INT(gcd)[1] = 1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS);
2782
        }
2783
        else {
2784
            gcd = INTOBJ_INT( i );
2785
        }
2786

2787
    }
2788

2789
    /* compute the gcd of two large integers                               */
2790
    else {
2791

2792
        /* a small divisor larger than one digit isn't handled above       */
2793
        if ( IS_INTOBJ(opL) ) {
2794
            if ( 0 < INT_INTOBJ(opL) ) {
2795
                gcd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2796
                ADDR_INT(gcd)[0] = (TypDigit)(INT_INTOBJ(opL));
2797
                ADDR_INT(gcd)[1] = (TypDigit)(INT_INTOBJ(opL)>>NR_DIGIT_BITS);
2798
                opL = gcd;
2799
            }
2800
            else {
2801
                gcd = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
2802
                ADDR_INT(gcd)[0] = (TypDigit)(-INT_INTOBJ(opL));
2803
                ADDR_INT(gcd)[1] = (TypDigit)((-INT_INTOBJ(opL))>>NR_DIGIT_BITS);
2804
                opL = gcd;
2805
            }
2806
        }
2807

2808
        /* a small dividend larger than one digit isn't handled above       */
2809
        if ( IS_INTOBJ(opR) ) {
2810
            if ( 0 < INT_INTOBJ(opR) ) {
2811
                gcd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2812
                ADDR_INT(gcd)[0] = (TypDigit)(INT_INTOBJ(opR));
2813
                ADDR_INT(gcd)[1] = (TypDigit)(INT_INTOBJ(opR)>>NR_DIGIT_BITS);
2814
                opR = gcd;
2815
            }
2816
            else {
2817
                gcd = NewBag( T_INTNEG, 4*sizeof(TypDigit) );
2818
                ADDR_INT(gcd)[0] = (TypDigit)(-INT_INTOBJ(opR));
2819
                ADDR_INT(gcd)[1] = (TypDigit)((-INT_INTOBJ(opR))>>NR_DIGIT_BITS);
2820
                opR = gcd;
2821
            }
2822
        }
2823

2824
        /* copy the left operand into a new bag                            */
2825
        gcd = NewBag( T_INTPOS, (SIZE_INT(opL)+4)*sizeof(TypDigit) );
2826
        l = ADDR_INT(opL);
2827
        r = ADDR_INT(gcd);
2828
        for ( k = SIZE_INT(opL)-1; k >= 0; k-- )
2829
            *r++ = *l++;
2830
        opL = gcd;
2831

2832
        /* get the size of the left operand                                */
2833
        ls = SIZE_INT(opL);
2834
        l  = ADDR_INT(opL);
2835
        while ( ls >= 1 && l[ls-1] == 0 )  ls--;
2836

2837
        /* copy the right operand into a new bag                           */
2838
        gcd = NewBag( T_INTPOS, (SIZE_INT(opR)+4)*sizeof(TypDigit) );
2839
        r = ADDR_INT(opR);
2840
        l = ADDR_INT(gcd);
2841
        for ( k = SIZE_INT(opR)-1; k >= 0; k-- )
2842
            *l++ = *r++;
2843
        opR = gcd;
2844

2845
        /* get the size of the right operand                               */
2846
        rs = SIZE_INT(opR);
2847
        r  = ADDR_INT(opR);
2848
        while ( rs >= 1 && r[rs-1] == 0 )  rs--;
2849

2850
        /* repeat while the right operand is large                         */
2851
        while ( rs >= 2 ) {
2852

2853
            /* get the leading two digits                                  */
2854
            for ( e = 0;
2855
                 ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e) : 0)
2856
                                   < INTBASE/2; e++ ) ;
2857
            r1 = ((Int)r[rs-1]<<e) + (e ? r[rs-2]>>(NR_DIGIT_BITS-e): 0);
2858
            r2 = ((Int)r[rs-2]<<e) + ((e && rs>=3) ? r[rs-3]>>(NR_DIGIT_BITS-e) : 0);
2859

2860
            /* run through the digits in the quotient                      */
2861
            for ( i = ls - rs; i >= 0; i-- ) {
2862

2863
                /* guess the factor                                        */
2864
                l = ADDR_INT(opL) + rs + i;
2865
                l01 = ((INTBASE*l[0]+l[-1])<<e) + (e ? l[-2]>>(NR_DIGIT_BITS-e):0);
2866
                if ( l01 == 0 )  continue;
2867
                l2  = ((Int)l[-2]<<e) + ((e && rs+i>=3) ? l[-3]>>(NR_DIGIT_BITS-e):0);
2868
                if ( ((Int)l[0]<<e)+(e ? l[-1]>>(NR_DIGIT_BITS-e):0) < r1 )
2869
                        qi = l01 / r1;
2870
                else    qi = INTBASE - 1;
2871
                while ( l01-(Int)qi*r1 < (Int)INTBASE
2872
                        && (INTBASE*(l01-(Int)qi*r1)+l2) < ((Int)qi*r2) )
2873
                    qi--;
2874

2875
                /* l = l - qi * r;                                         */
2876
                d = 0;
2877
                l = ADDR_INT(opL) + i;
2878
                r = ADDR_INT(opR);
2879
                for ( k = 0; k < (Int)rs; ++k, ++l, ++r ) {
2880
                    c = (Int)*l - (Int)qi * *r - (Int)d;
2881
                    *l = (TypDigit)c;
2882
                    d = -(TypDigit)(c>>NR_DIGIT_BITS);
2883
                }
2884
                c = (Int)*l - d;  *l = (TypDigit)c;  d = -(TypDigit)(c>>NR_DIGIT_BITS);
2885

2886
                /* if we have a borrow then add back                       */
2887
                if ( d != 0 ) {
2888
                    d = 0;
2889
                    l = ADDR_INT(opL) + i;
2890
                    r = ADDR_INT(opR);
2891
                    for ( k = 0; k < (Int)rs; ++k, ++l, ++r ) {
2892
                        c = (Int)*l + (Int)*r + (Int)d;
2893
                        *l = (TypDigit)c;
2894
                        d = (TypDigit)(c>>NR_DIGIT_BITS);
2895
                    }
2896
                    c = *l + d;  *l = (TypDigit)c;  d = (TypDigit)(c>>NR_DIGIT_BITS);
2897
                    qi--;
2898
                }
2899
            }
2900

2901
            /* exchange the two operands                                   */
2902
            gcd = opL;  opL = opR;  opR = gcd;
2903
            ls = rs;
2904

2905
            /* get the size of the right operand                           */
2906
            rs = SIZE_INT(opR);
2907
            r  = ADDR_INT(opR);
2908
            while ( rs >= 1 && r[rs-1] == 0 )  rs--;
2909

2910
        }
2911

2912
        /* if the right operand is zero now, the left is the gcd           */
2913
        if ( rs == 0 ) {
2914

2915
            /* remove the leading zeroes                                   */
2916
            l = ADDR_INT(opL) + SIZE_INT(opL);
2917
            if ( (l[-4] == 0 && l[-3] == 0 && l[-2] == 0 && l[-1] == 0)
2918
              || (SIZE_INT(opL) == 4 && l[-2] == 0 && l[-1] == 0) ) {
2919

2920
                /* find the number of significant digits                   */
2921
                l = ADDR_INT(opL);
2922
                for ( k = SIZE_INT(opL); k != 0; k-- ) {
2923
                    if ( l[k-1] != 0 )
2924
                        break;
2925
                }
2926

2927
                /* reduce to small integer if possible, otherwise shrink b */
2928
                if ( k <= 2 && TNUM_OBJ(opL) == T_INTPOS
2929
                  && (UInt)(INTBASE*l[1]+l[0]) < (1L<<NR_SMALL_INT_BITS) )
2930
                    opL = INTOBJ_INT( INTBASE*l[1]+l[0] );
2931
                else if ( k <= 2 && TNUM_OBJ(opL) == T_INTNEG
2932
                  && (UInt)(INTBASE*l[1]+l[0]) <= (1L<<NR_SMALL_INT_BITS) )
2933
                    opL = INTOBJ_INT( -(Int)(INTBASE*l[1]+l[0]) );
2934
                else
2935
                    ResizeBag( opL, (((k + 3) / 4) * 4) * sizeof(TypDigit) );
2936
            }
2937

2938
            gcd = opL;
2939

2940
        }
2941

2942
        /* otherwise handle one large and one small integer as above       */
2943
        else {
2944

2945
            /* get the right operand value, make it positive               */
2946
            i = r[0];
2947

2948
            /* do one remainder operation                                  */
2949
            l = ADDR_INT(opL) + SIZE_INT(opL) - 1;
2950
            c = 0;
2951
            for ( ; l >= ADDR_INT(opL); l-- ) {
2952
                c  = (c<<NR_DIGIT_BITS) + (Int)*l;
2953
                c  = c % i;
2954
            }
2955
            k = c;
2956

2957
            /* compute the gcd using Euclids algorithm                     */
2958
            while ( k != 0 ) {
2959
                c = k;
2960
                k = i % k;
2961
                i = c;
2962
            }
2963

2964
            /* now i is the result                                         */
2965
            if ( i == (1L<<NR_SMALL_INT_BITS) ) {
2966
                gcd = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
2967
                ADDR_INT(gcd)[0] = 0;
2968
                ADDR_INT(gcd)[1] = 1L<<(NR_SMALL_INT_BITS-NR_DIGIT_BITS);
2969
            }
2970
            else {
2971
                gcd = INTOBJ_INT( i );
2972
            }
2973

2974
        }
2975

2976
    }
2977

2978
    /* return the result                                                   */
2979
    return gcd;
2980
}
2981

2982

2983
/****************************************************************************
2984
**
2985
*F  FuncGCD_INT(<self>,<opL>,<opR>)  . . . . . . .  internal function 'GcdInt'
2986
**
2987
**  'FuncGCD_INT' implements the internal function 'GcdInt'.
2988
**
2989
**  'GcdInt( <i>, <k> )'
2990
**
2991
**  'Gcd'  returns the greatest common divisor   of the two  integers <m> and
2992
**  <n>, i.e.,  the  greatest integer that  divides  both <m>  and  <n>.  The
2993
**  greatest common divisor is never negative, even if the arguments are.  We
2994
**  define $gcd( m, 0 ) = gcd( 0, m ) = abs( m )$ and $gcd( 0, 0 ) = 0$.
2995
*/
2996
Obj             FuncGCD_INT (
2997
    Obj                 self,
2998
    Obj                 opL,
2999
    Obj                 opR )
3000
{
3001
    /* check the arguments                                                 */
3002
    while ( TNUM_OBJ(opL) != T_INT
3003
         && TNUM_OBJ(opL) != T_INTPOS
3004
         && TNUM_OBJ(opL) != T_INTNEG ) {
3005
        opL = ErrorReturnObj(
3006
            "GcdInt: <left> must be an integer (not a %s)",
3007
            (Int)TNAM_OBJ(opL), 0L,
3008
            "you can replace <left> via 'return <left>;'" );
3009
    }
3010
    while ( TNUM_OBJ(opR) != T_INT
3011
         && TNUM_OBJ(opR) != T_INTPOS
3012
         && TNUM_OBJ(opR) != T_INTNEG ) {
3013
        opR = ErrorReturnObj(
3014
            "GcdInt: <right> must be an integer (not a %s)",
3015
            (Int)TNAM_OBJ(opR), 0L,
3016
            "you can replace <right> via 'return <right>;'" );
3017
    }
3018

3019
    /* return the gcd                                                      */
3020
    return GcdInt( opL, opR );
3021
}
3022

3023

3024

3025

3026

3027

3028
/****************************************************************************
3029
**
3030
*F  SaveInt( <int> )
3031
**
3032
**  Since the type is saved, we don't need to worry about sign
3033
*/
3034

3035
void SaveInt( Obj bigint)
3036
{
3037
  TypDigit *ptr;
3038
  UInt i;
3039
  ptr = (TypDigit *)ADDR_OBJ(bigint);
3040
  for (i = 0; i < SIZE_INT(bigint); i++)
3041
#ifdef SYS_IS_64_BIT
3042
    SaveUInt4(*ptr++);
3043
#else
3044
    SaveUInt2(*ptr++);
3045
#endif
3046
  return;
3047
}
3048

3049
/****************************************************************************
3050
**
3051
*F  LoadInt( <int> )
3052
**
3053
**  Since the type is loaded, we don't need to worry about sign
3054
*/
3055

3056
void LoadInt( Obj bigint)
3057
{
3058
  TypDigit *ptr;
3059
  UInt i;
3060
  ptr = (TypDigit *)ADDR_OBJ(bigint);
3061
  for (i = 0; i < SIZE_INT(bigint); i++)
3062
#ifdef SYS_IS_64_BIT
3063
    *ptr++ = LoadUInt4();
3064
#else
3065
    *ptr++ = LoadUInt2();
3066
#endif
3067
  return;
3068
}
3069

3070

3071
/****************************************************************************
3072
**
3073
** * * * * * * * "Mersenne twister" random numbers  * * * * * * * * * * * * *
3074
**
3075
**  Part of this code for fast generation of 32 bit pseudo random numbers with
3076
**  a period of length 2^19937-1 and a 623-dimensional equidistribution is
3077
**  taken from:
3078
**          http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
3079
**  (Also look in Wikipedia for "Mersenne twister".)
3080
**  We use the file mt19937ar.c, version 2002/1/26.
3081
*/
3082

3083

3084

3085
/****************************************************************************
3086
**
3087
*F  RandomIntegerMT( <mtstr>, <nrbits> )
3088
**
3089
**  Returns an integer with at most <nrbits> bits in uniform distribution.
3090
**  <nrbits> must be a small integer. <mtstr> is a string as returned by
3091
**  InitRandomMT.
3092
**
3093
**  Implementation details are a bit tricky to obtain the same random
3094
**  integers on 32 bit and 64 bit machines (which have different long
3095
**  integer digit lengths and different ranges of small integers).
3096
**
3097
*/
3098
/* for comparison in case result is small int */
3099
Obj SMALLEST_INTPOS = NULL;
3100

3101
Obj FuncRandomIntegerMT(Obj self, Obj mtstr, Obj nrbits)
3102
{
3103
  Obj res;
3104
  Int i, n, q, r, qoff, len, nlen;
3105
  UInt4 *mt, rand;
3106
  TypDigit *pt;
3107
  while (! IsStringConv(mtstr)) {
3108
     mtstr = ErrorReturnObj(
3109
         "<mtstr> must be a string, not a %s)",
3110
         (Int)TNAM_OBJ(mtstr), 0L,
3111
         "you can replace <mtstr> via 'return <mtstr>;'" );
3112
  }
3113
  while ((! IsStringConv(mtstr)) || GET_LEN_STRING(mtstr) < 2500) {
3114
     mtstr = ErrorReturnObj(
3115
         "<mtstr> must be a string with at least 2500 characters, ",
3116
         0L, 0L,
3117
         "you can replace <mtstr> via 'return <mtstr>;'" );
3118
  }
3119
  while ((! IS_INTOBJ(nrbits)) || INT_INTOBJ(nrbits) < 0) {
3120
     nrbits = ErrorReturnObj(
3121
         "<nrbits> must be a small non-negative integer, not a %s)",
3122
         (Int)TNAM_OBJ(nrbits), 0L,
3123
         "you can replace <mtstr> via 'return <mtstr>;'" );
3124
  }
3125
  n = INT_INTOBJ(nrbits);
3126

3127
  /* small int case */
3128
  if (n <= NR_SMALL_INT_BITS) {
3129
     mt = (UInt4*) CHARS_STRING(mtstr);
3130
#ifdef SYS_IS_64_BIT
3131
     if (n <= 32) {
3132
       res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4) -1L >> (32-n))));
3133
     }
3134
     else {
3135
       unsigned long  rd;
3136
       rd = nextrandMT_int32(mt);
3137
       rd += (unsigned long) ((UInt4) nextrandMT_int32(mt) &
3138
                              ((UInt4) -1L >> (64-n))) << 32;
3139
       res = INTOBJ_INT((Int)rd);
3140
     }
3141
#else
3142
     res = INTOBJ_INT((Int)(nextrandMT_int32(mt) & ((UInt4) -1L >> (32-n))));
3143
#endif
3144
  }
3145
  else {
3146
     /* number of Digits */
3147
     q = n / NR_DIGIT_BITS;
3148
     r = n - q*NR_DIGIT_BITS;
3149
     qoff = q + (r==0 ? 0:1);
3150
     len = qoff;
3151
     len = 4*((len+3) / 4);
3152
     res = NewBag( T_INTPOS, len*sizeof(TypDigit) );
3153
     pt = ADDR_INT(res);
3154
     mt = (UInt4*) CHARS_STRING(mtstr);
3155
#ifdef SYS_IS_64_BIT
3156
     for (i = 0; i < qoff; i++, pt++) {
3157
       rand = (TypDigit) nextrandMT_int32(mt);
3158
       *pt = rand;
3159
     }
3160
#else
3161
     for (i = 0; i < qoff/2; i++, pt++) {
3162
       rand = nextrandMT_int32(mt);
3163
       *pt = (TypDigit) (rand & 0xFFFFL);
3164
       pt++;
3165
       *pt = (TypDigit) ((rand & 0xFFFF0000L) >> 16);
3166
     }
3167
     if (2*i != qoff) {
3168
       rand = nextrandMT_int32(mt);
3169
       *pt = (TypDigit) (rand & 0xFFFFL);
3170
     }
3171
#endif
3172
     if (r != 0) {
3173
       ADDR_INT(res)[qoff-1] = ADDR_INT(res)[qoff-1] & ((TypDigit)(-1)
3174
                                                      >> (NR_DIGIT_BITS-r));
3175
     }
3176
     /* shrink bag if necessary */
3177
     for (nlen = len, i=0; nlen >0 && ADDR_INT(res)[nlen-1] == 0L;
3178
                           nlen--, i++);
3179
     if (i/4 != 0) {
3180
        ResizeBag(res, 4*((nlen+3) / 4)*sizeof(TypDigit));
3181
     }
3182
     /* convert result if small int */
3183
     if (LtInt(res, SMALLEST_INTPOS)) {
3184
       res = INTOBJ_INT((Int)(ADDR_INT(res)[0] +
3185
                              ((UInt)ADDR_INT(res)[1] << NR_DIGIT_BITS)));
3186
     }
3187
  }
3188

3189
  return res;
3190
}
3191

3192

3193

3194
/****************************************************************************
3195
**
3196
*F * * * * * * * * * * * * * initialize package * * * * * * * * * * * * * * *
3197
*/
3198

3199
/****************************************************************************
3200
**
3201
*V  GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export
3202
*/
3203
static StructGVarFilt GVarFilts [] = {
3204

3205
    { "IS_INT", "obj", &IsIntFilt,
3206
      FuncIS_INT, "src/integer.c:IS_INT" },
3207

3208
    { 0 }
3209

3210
};
3211

3212

3213
/****************************************************************************
3214
**
3215
*V  GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
3216
*/
3217
static StructGVarFunc GVarFuncs [] = {
3218

3219
    { "QUO_INT", 2, "int1, int2",
3220
      FuncQUO_INT, "src/integer.c:QUO_INT" },
3221

3222
    { "REM_INT", 2, "int1, int2",
3223
      FuncREM_INT, "src/integer.c:REM_INT" },
3224

3225
    { "GCD_INT", 2, "int1, int2",
3226
      FuncGCD_INT, "src/integer.c:GCD_INT" },
3227

3228
    { "PROD_INT_OBJ", 2, "int, obj",
3229
      FuncPROD_INT_OBJ, "src/integer.c:PROD_INT_OBJ" },
3230

3231
    { "POW_OBJ_INT", 2, "obj, int",
3232
      FuncPOW_OBJ_INT, "src/integer.c:POW_OBJ_INT" },
3233

3234
    { "HexStringInt", 1, "integer",
3235
      FuncHexStringInt, "src/integer.c:HexStringInt" },
3236

3237
    { "IntHexString", 1, "string",
3238
      FuncIntHexString, "src/integer.c:IntHexString" },
3239

3240
    { "Log2Int", 1, "int",
3241
      FuncLog2Int, "src/integer.c:Log2Int" },
3242

3243
    { "STRING_INT", 1, "int",
3244
      FuncSTRING_INT, "src/integer.c:STRING_INT" },
3245

3246

3247
    { "RandomIntegerMT", 2, "mtstr, nrbits",
3248
      FuncRandomIntegerMT, "src/integer.c:RandomIntegerMT" },
3249

3250

3251
    { 0 }
3252

3253
};
3254

3255

3256
/****************************************************************************
3257
**
3258
*F  InitKernel( <module> )  . . . . . . . . initialise kernel data structures
3259
*/
3260
static Int InitKernel (
3261
    StructInitInfo *    module )
3262
{
3263
    UInt                t1,  t2;
3264

3265
    /* init filters and functions                                          */
3266
    InitHdlrFiltsFromTable( GVarFilts );
3267
    InitHdlrFuncsFromTable( GVarFuncs );
3268

3269
    /* install the marking functions                                       */
3270
    InfoBags[         T_INT    ].name = "integer";
3271
#ifdef SYS_IS_64_BIT
3272
    InfoBags[         T_INTPOS ].name = "integer (>= 2^60)";
3273
    InfoBags[         T_INTNEG ].name = "integer (< -2^60)";
3274
#else
3275
    InfoBags[         T_INTPOS ].name = "integer (>= 2^28)";
3276
    InfoBags[         T_INTNEG ].name = "integer (< -2^28)";
3277
#endif
3278
    InitMarkFuncBags( T_INTPOS, MarkNoSubBags );
3279
    InitMarkFuncBags( T_INTNEG, MarkNoSubBags );
3280

3281
    /* Install the saving methods */
3282
    SaveObjFuncs [ T_INTPOS ] = SaveInt;
3283
    SaveObjFuncs [ T_INTNEG ] = SaveInt;
3284
    LoadObjFuncs [ T_INTPOS ] = LoadInt;
3285
    LoadObjFuncs [ T_INTNEG ] = LoadInt;
3286

3287
    /* install the printing function                                       */
3288
    PrintObjFuncs[ T_INT    ] = PrintInt;
3289
    PrintObjFuncs[ T_INTPOS ] = PrintInt;
3290
    PrintObjFuncs[ T_INTNEG ] = PrintInt;
3291

3292
    /* install the comparison methods                                      */
3293
    for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
3294
        for ( t2 = T_INT; t2 <= T_INTNEG; t2++ ) {
3295
            EqFuncs  [ t1 ][ t2 ] = EqInt;
3296
            LtFuncs  [ t1 ][ t2 ] = LtInt;
3297
        }
3298
    }
3299

3300
    /* install the unary arithmetic methods                                */
3301
    for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
3302
        ZeroFuncs[ t1 ] = ZeroInt;
3303
        ZeroMutFuncs[ t1 ] = ZeroInt;
3304
        AInvFuncs[ t1 ] = AInvInt;
3305
        AInvMutFuncs[ t1 ] = AInvInt;
3306
        OneFuncs [ t1 ] = OneInt;
3307
        OneMutFuncs [ t1 ] = OneInt;
3308
    }
3309

3310
    /* install the default product and power methods                       */
3311
    for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
3312
        for ( t2 = FIRST_CONSTANT_TNUM;  t2 <= LAST_CONSTANT_TNUM;  t2++ ) {
3313
            ProdFuncs[ t1 ][ t2 ] = ProdIntObj;
3314
            PowFuncs [ t2 ][ t1 ] = PowObjInt;
3315
        }
3316
        for ( t2 = FIRST_RECORD_TNUM;  t2 <= LAST_RECORD_TNUM;  t2++ ) {
3317
            ProdFuncs[ t1 ][ t2 ] = ProdIntObj;
3318
            PowFuncs [ t2 ][ t1 ] = PowObjInt;
3319
        }
3320
        for ( t2 = FIRST_LIST_TNUM;    t2 <= LAST_LIST_TNUM;    t2++ ) {
3321
            ProdFuncs[ t1 ][ t2 ] = ProdIntObj;
3322
            PowFuncs [ t2 ][ t1 ] = PowObjInt;
3323
        }
3324
    }
3325

3326
    /* install the binary arithmetic methods                               */
3327
    for ( t1 = T_INT; t1 <= T_INTNEG; t1++ ) {
3328
        for ( t2 = T_INT; t2 <= T_INTNEG; t2++ ) {
3329
            EqFuncs  [ t1 ][ t2 ] = EqInt;
3330
            LtFuncs  [ t1 ][ t2 ] = LtInt;
3331
            SumFuncs [ t1 ][ t2 ] = SumInt;
3332
            DiffFuncs[ t1 ][ t2 ] = DiffInt;
3333
            ProdFuncs[ t1 ][ t2 ] = ProdInt;
3334
            PowFuncs [ t1 ][ t2 ] = PowInt;
3335
            ModFuncs [ t1 ][ t2 ] = ModInt;
3336
        }
3337
    }
3338

3339
    /* gvars to import from the library                                    */
3340
    ImportGVarFromLibrary( "TYPE_INT_SMALL_ZERO", &TYPE_INT_SMALL_ZERO );
3341
    ImportGVarFromLibrary( "TYPE_INT_SMALL_POS",  &TYPE_INT_SMALL_POS );
3342
    ImportGVarFromLibrary( "TYPE_INT_SMALL_NEG",  &TYPE_INT_SMALL_NEG );
3343
    ImportGVarFromLibrary( "TYPE_INT_LARGE_POS",  &TYPE_INT_LARGE_POS );
3344
    ImportGVarFromLibrary( "TYPE_INT_LARGE_NEG",  &TYPE_INT_LARGE_NEG );
3345
    ImportGVarFromLibrary( "SMALLEST_INTPOS", &SMALLEST_INTPOS );
3346

3347
    ImportFuncFromLibrary( "String", &String );
3348
    ImportFuncFromLibrary( "One", &OneAttr);
3349

3350
    /* install the type functions                                          */
3351
    TypeObjFuncs[ T_INT    ] = TypeIntSmall;
3352
    TypeObjFuncs[ T_INTPOS ] = TypeIntLargePos;
3353
    TypeObjFuncs[ T_INTNEG ] = TypeIntLargeNeg;
3354

3355
    MakeBagTypePublic( T_INTPOS );
3356
    MakeBagTypePublic( T_INTNEG );
3357

3358
    /* return success                                                      */
3359
    return 0;
3360
}
3361

3362

3363
/****************************************************************************
3364
**
3365
*F  InitLibrary( <module> ) . . . . . . .  initialise library data structures
3366
*/
3367
static Int InitLibrary (
3368
    StructInitInfo *    module )
3369
{
3370
    UInt gvar;
3371
    /* init filters and functions                                          */
3372
    InitGVarFiltsFromTable( GVarFilts );
3373
    InitGVarFuncsFromTable( GVarFuncs );
3374

3375
    /* hold smallest large integer */
3376
    SMALLEST_INTPOS = NewBag( T_INTPOS, 4*sizeof(TypDigit) );
3377
    ADDR_INT(SMALLEST_INTPOS)[1] =
3378
                 (TypDigit) (1L << (NR_SMALL_INT_BITS - NR_DIGIT_BITS));
3379
    gvar = GVarName("SMALLEST_INTPOS");
3380
    MakeReadWriteGVar( gvar);
3381
    AssGVar( gvar, SMALLEST_INTPOS );
3382
    MakeReadOnlyGVar(gvar);
3383

3384
    /* return success                                                      */
3385
    return 0;
3386
}
3387

3388

3389
/****************************************************************************
3390
**
3391
*F  InitInfoInt() . . . . . . . . . . . . . . . . . . table of init functions
3392
*/
3393
static StructInitInfo module = {
3394
    MODULE_BUILTIN,                     /* type                           */
3395
    "integer",                          /* name                           */
3396
    0,                                  /* revision entry of c file       */
3397
    0,                                  /* revision entry of h file       */
3398
    0,                                  /* version                        */
3399
    0,                                  /* crc                            */
3400
    InitKernel,                         /* initKernel                     */
3401
    InitLibrary,                        /* initLibrary                    */
3402
    0,                                  /* checkInit                      */
3403
    0,                                  /* preSave                        */
3404
    0,                                  /* postSave                       */
3405
    0                                   /* postRestore                    */
3406
};
3407

3408
StructInitInfo * InitInfoInt ( void )
3409
{
3410
    return &module;
3411
}
3412

3413
/* matches the USE_GMP test at the top of the file */
3414
#endif
3415

3416
/****************************************************************************
3417
**
3418
*E  integer.c . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here
3419
*/
3420

3421
Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.

Product

Resources

Company

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more, all in one place.

Real-time collaboration for Jupyter Notebooks, Linux Terminals, LaTeX, VS Code, R IDE, and more,
all in one place.
