1 /* bignum.c -*- mode:c; coding:utf-8; -*-
2 *
3 * Copyright (c) 2010-2021 Takashi Kato <ktakashi@ymail.com>
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 *
16 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
17 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
18 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
19 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
20 * OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
21 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
22 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
23 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
24 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
25 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
26 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 *
28 * $Id: $
29 */
30
31 #include <math.h>
32 #include <string.h>
33
34 #define NO_NBITS
35 #define LIBSAGITTARIUS_BODY
36 #include "sagittarius/private/bignum.h"
37 #include "sagittarius/private/core.h"
38 #include "sagittarius/private/number.h"
39 #include "sagittarius/private/error.h"
40 #include "sagittarius/private/arith.h"
41 #include "sagittarius/private/bits.h"
42 #include "sagittarius/private/pair.h"
43 #include "sagittarius/private/port.h"
44 #include "sagittarius/private/string.h"
45 #include "sagittarius/private/vm.h"
46
47 #undef min
48 #define min(x, y) (((x) < (y))? (x) : (y))
49 #undef max
50 #define max(x, y) (((x) > (y))? (x) : (y))
51
52 /* will be used in bignum.inc
53 TODO move to bignum.inc
54 */
55 #define clear_buffer(v, size) \
56 do { \
57 long __i, __size = (size); \
58 for (__i = 0; __i < __size; __i++) (v)[__i] = 0L; \
59 } while (0)
60
61 #ifdef HAVE_ALLOCA
62 #define ALLOC_TEMP_BUFFER_REC(v, type, size) \
63 (v) = (type *)alloca(sizeof(type) * (size));
64 #else
65 #define ALLOC_TEMP_BUFFER_REC(v, type, size) \
66 (v) = SG_NEW_ATOMIC2(type *, sizeof(type) * (size));
67 #endif
68
69 #define ALLOC_TEMP_BUFFER(v, type, size) \
70 do { \
71 ALLOC_TEMP_BUFFER_REC(v, type, size) \
72 clear_buffer(v, size); \
73 } while (0)
74
75
76 /* debug utility macro */
77 #define dump_array_rec(flag, array, size) \
78 do { \
79 int __i, __size = (size); \
80 fprintf(stderr, #array " = [ "); \
81 for (__i = 0; __i < __size; __i++) { \
82 fprintf(stderr, flag, (array)[__i]); \
83 } \
84 fprintf(stderr, "]\n"); \
85 } while (0)
86
87 #define dump_sarray(array, size) dump_array_rec("%ld ", array, size)
88 #define dump_uarray(array, size) dump_array_rec("%lu ", array, size)
89 #if SIZEOF_LONG == 8
90 #define dump_xarray(array, size) dump_array_rec("%016lx ", array, size)
91 #else
92 #define dump_xarray(array, size) dump_array_rec("%08lx ", array, size)
93 #endif
94
95 #define dump_rec(flag, v) fprintf(stderr, #v" "flag, v)
96 #define dump_s(v) dump_rec("%ld\n", v)
97 #define dump_u(v) dump_rec("%lu\n", v)
98
99 #define dump_bignum_s(b) dump_sarray((b)->elements, (b)->size)
100 #define dump_bignum_u(b) dump_uarray((b)->elements, (b)->size)
101 #define dump_bignum_x(b) dump_xarray((b)->elements, (b)->size)
102
103 /* TODO maybe we want to put this in header file? */
104 #define SG_LEFT_SHIFT_SPACE(size, amount) \
105 (int)((size) + ((amount) + WORD_BITS -1)/WORD_BITS)
106 #include "bignum.inc"
107
108 static long bignum_safe_size_for_add(SgBignum *x, SgBignum *y);
109 static SgBignum *bignum_add_int(SgBignum *br, SgBignum *bx, SgBignum *by);
110
bignum_clear(SgBignum * b,long size)111 static SgBignum* bignum_clear(SgBignum *b, long size)
112 {
113 long i;
114 for (i = 0; i < size; i++) b->elements[i] = 0L;
115 return b;
116 }
117
make_bignum_rec(long size,int need_clear)118 static SgBignum* make_bignum_rec(long size, int need_clear)
119 {
120 SgBignum *b;
121 long real_size = size;
122 if (size < 0) Sg_Error(UC("[internal error] invalid bignum size: %d"), size);
123 if (size > SG_BIGNUM_MAX_DIGITS) Sg_Error(UC("too large bignum"));
124 if (real_size == 0) real_size++; /* to avoid minus allocation */
125 b = SG_NEW_ATOMIC2(SgBignum*, BIGNUM_SIZE(real_size));
126 SG_SET_CLASS(b, SG_CLASS_INTEGER);
127 if (size == 0) {
128 SG_BIGNUM_SET_ZERO(b);
129 } else {
130 SG_BIGNUM_SET_COUNT(b, size);
131 SG_BIGNUM_SET_SIGN(b, 1);
132 }
133 if (need_clear)
134 return bignum_clear(b, real_size);
135 else
136 return b;
137 }
138
Sg_AllocateBignum(int size)139 SgBignum* Sg_AllocateBignum(int size)
140 {
141 /* should we initialise to zero? */
142 return make_bignum_rec(size, TRUE);
143 }
144
145 #define make_bignum(size) make_bignum_rec(size, TRUE)
146
147 /* Should we expose this? */
148 #ifdef HAVE_ALLOCA
149 #define ALLOC_TEMP_BIGNUM_REC(var, size) \
150 do { \
151 (var) = SG_BIGNUM(alloca(BIGNUM_SIZE(size))); \
152 SG_SET_CLASS(var, SG_CLASS_INTEGER); \
153 SG_BIGNUM_SET_COUNT(var, size); \
154 SG_BIGNUM_SET_SIGN(var, 1); \
155 } while(0)
156
157 #define ALLOC_TEMP_BIGNUM(var, size) \
158 do { \
159 ALLOC_TEMP_BIGNUM_REC(var, size); \
160 bignum_clear(var, size); \
161 } while (0)
162 #else
163 #define ALLOC_TEMP_BIGNUM_REC(var, size) \
164 do { \
165 (var) = make_bignum_rec(size, FALSE); \
166 } while (0)
167
168 #define ALLOC_TEMP_BIGNUM(var, size) \
169 do { \
170 ALLOC_TEMP_BIGNUM_REC(var, size); \
171 bignum_clear(var, size); \
172 } while (0)
173 #endif
174
175
Sg_MakeBignumFromSI(long value)176 SgObject Sg_MakeBignumFromSI(long value)
177 {
178 SgBignum *b;
179 if (value == 0) return make_bignum(0);
180 else if (value == LONG_MIN) {
181 b = make_bignum(1);
182 SG_BIGNUM_SET_SIGN(b, -1);
183 b->elements[0] = (unsigned long)LONG_MAX + 1;
184 } else if (value < 0) {
185 b = make_bignum(1);
186 SG_BIGNUM_SET_SIGN(b, -1);
187 b->elements[0] = -value;
188 } else {
189 b = make_bignum(1);
190 SG_BIGNUM_SET_SIGN(b, 1);
191 b->elements[0] = value;
192 }
193 return SG_OBJ(b);
194 }
195
Sg_MakeBignumFromUI(unsigned long value)196 SgObject Sg_MakeBignumFromUI(unsigned long value)
197 {
198 SgBignum *b;
199 if (value == 0) return make_bignum(0);
200 b = make_bignum(1);
201 SG_BIGNUM_SET_SIGN(b, 1);
202 b->elements[0] = value;
203 return SG_OBJ(b);
204 }
205
Sg_MakeBignumFromU64(uint64_t value)206 SgObject Sg_MakeBignumFromU64(uint64_t value)
207 {
208 if (value) {
209 #if SIZEOF_LONG >= 8
210 SgBignum *ans = make_bignum(1);
211 SG_BIGNUM_SET_SIGN(ans, 1);
212 ans->elements[0] = value;
213 return ans;
214 #else
215 SgBignum *ans;
216 if ((value >> WORD_BITS) != 0) {
217 ans = make_bignum(2);
218 ans->elements[0] = value & SG_ULONG_MAX;
219 ans->elements[1] = value >> WORD_BITS;
220 } else {
221 ans = make_bignum(1);
222 ans->elements[0] = (unsigned long)value;
223 }
224 SG_BIGNUM_SET_SIGN(ans, 1);
225 return ans;
226 #endif
227 }
228 return make_bignum(0);
229 }
230
Sg_MakeBignumFromS64(int64_t value)231 SgObject Sg_MakeBignumFromS64(int64_t value)
232 {
233 if (value) {
234 #if SIZEOF_LONG >= 8
235 SgBignum *ans = make_bignum(1);
236 if (value > 0) {
237 SG_BIGNUM_SET_SIGN(ans, 1);
238 ans->elements[0] = value;
239 } else {
240 SG_BIGNUM_SET_SIGN(ans, -1);
241 ans->elements[0] = -value;
242 }
243 return ans;
244 #else
245 int sign;
246 SgBignum *ans;
247 if (value > 0) {
248 sign = 1;
249 } else {
250 sign = -1;
251 value = -value;
252 }
253 if ((value >> WORD_BITS) != 0) {
254 ans = make_bignum(2);
255 ans->elements[0] = value & SG_ULONG_MAX;
256 ans->elements[1] = value >> WORD_BITS;
257 } else {
258 ans = make_bignum(1);
259 ans->elements[0] = (long)value;
260 }
261 SG_BIGNUM_SET_SIGN(ans, sign);
262 return ans;
263 #endif
264 }
265 return make_bignum(0);
266 }
267
Sg_MakeBignumFromDouble(double value)268 SgObject Sg_MakeBignumFromDouble(double value)
269 {
270 int expornent, sign;
271 SgObject mantissa, b;
272 if (value >= LONG_MIN && value <= LONG_MAX) {
273 return Sg_MakeBignumFromSI((long)value);
274 }
275 mantissa = Sg_DecodeFlonum(value, &expornent, &sign);
276 if (!SG_NUMBERP(mantissa)) {
277 Sg_Error(UC("can't convert %lf to an integer"), value);
278 }
279 b = Sg_Ash(mantissa, expornent);
280 if (sign < 0) b = Sg_Negate(b);
281 if (SG_INTP(b)) {
282 return Sg_MakeBignumFromSI(SG_INT_VALUE(b));
283 } else {
284 return b;
285 }
286 }
287
bignum_copy(SgBignum * dst,SgBignum * src)288 static inline void bignum_copy(SgBignum *dst, SgBignum *src)
289 {
290 long i;
291 long size = SG_BIGNUM_GET_COUNT(src);
292 ASSERT(SG_BIGNUM_GET_COUNT(dst) >= size);
293 SG_BIGNUM_SET_SIGN(dst, SG_BIGNUM_GET_SIGN(src));
294 SG_BIGNUM_SET_COUNT(dst, size);
295 for (i = 0; i < size; i++) dst->elements[i] = src->elements[i];
296 }
297
Sg_BignumCopy(SgBignum * b)298 SgObject Sg_BignumCopy(SgBignum *b)
299 {
300 SgBignum *c = make_bignum_rec(SG_BIGNUM_GET_COUNT(b), FALSE);
301 bignum_copy(c, b);
302 return SG_OBJ(c);
303 }
304
bignum_normalize_rec(SgBignum * bn,int convertp)305 static SgObject bignum_normalize_rec(SgBignum *bn, int convertp)
306 {
307 long size = SG_BIGNUM_GET_COUNT(bn);
308 long i;
309 for (i = size - 1; i > 0; i--) {
310 if (bn->elements[i] == 0) size--;
311 else break;
312 }
313 /* if given bignum was 0, i can be -1 */
314 if (SG_BIGNUM_GET_SIGN(bn) == 0 && convertp) return SG_MAKE_INT(0);
315 if (i <= 0) {
316 if (convertp) {
317 if (SG_BIGNUM_GET_SIGN(bn) == 0 ||
318 (size == 1 && bn->elements[0] == 0)) {
319 return SG_MAKE_INT(0);
320 }
321 if (SG_BIGNUM_GET_SIGN(bn) > 0
322 && bn->elements[0] <= (unsigned long)SG_INT_MAX) {
323 return SG_MAKE_INT(bn->elements[0]);
324 }
325 if (SG_BIGNUM_GET_SIGN(bn) < 0
326 && bn->elements[0] <= (unsigned long)-SG_INT_MIN) {
327 return SG_MAKE_INT(-((signed long)bn->elements[0]));
328 }
329 } else {
330 /* make 0 bignum */
331 if (size == 1 && bn->elements[0] == 0) {
332 SG_BIGNUM_SET_ZERO(bn);
333 return bn;
334 }
335 }
336 }
337 SG_BIGNUM_SET_COUNT(bn, size);
338 return SG_OBJ(bn);
339 }
340
bignum_normalize(SgBignum * bn)341 static SgBignum* bignum_normalize(SgBignum *bn)
342 {
343 return SG_BIGNUM(bignum_normalize_rec(bn, FALSE));
344 }
345
Sg_NormalizeBignum(SgBignum * bn)346 SgObject Sg_NormalizeBignum(SgBignum *bn)
347 {
348 return bignum_normalize_rec(bn, TRUE);
349 }
350
351
352 /* assume bignums are normalized */
Sg_BignumCmp(SgBignum * lhs,SgBignum * rhs)353 int Sg_BignumCmp(SgBignum *lhs, SgBignum *rhs)
354 {
355 long lhs_count = SG_BIGNUM_GET_COUNT(lhs);
356 long rhs_count = SG_BIGNUM_GET_COUNT(rhs);
357 int lhs_sign = SG_BIGNUM_GET_SIGN(lhs);
358 int rhs_sign = SG_BIGNUM_GET_SIGN(rhs);
359 long i;
360
361 if (lhs_sign < rhs_sign) return -1;
362 if (lhs_sign > rhs_sign) return 1;
363 if (lhs_count < rhs_count) return (lhs_sign > 0) ? -1 : 1;
364 if (lhs_count > rhs_count) return (lhs_sign > 0) ? 1 : -1;
365 for (i = lhs_count - 1; i >= 0; i--) {
366 if (lhs->elements[i] < rhs->elements[i]) return (lhs_sign > 0) ? -1 : 1;
367 if (lhs->elements[i] > rhs->elements[i]) return (lhs_sign > 0) ? 1 : -1;
368 }
369 return 0;
370 }
371
Sg_BignumCmp3U(SgBignum * bx,SgBignum * off,SgBignum * by)372 int Sg_BignumCmp3U(SgBignum *bx, SgBignum *off, SgBignum *by)
373 {
374 long xsize = SG_BIGNUM_GET_COUNT(bx);
375 long ysize = SG_BIGNUM_GET_COUNT(by);
376 long osize = SG_BIGNUM_GET_COUNT(off);
377 long tsize;
378 long i;
379 SgBignum *br;
380
381 if (xsize > ysize) return 1;
382 if (xsize < ysize) {
383 if (osize < ysize && by->elements[ysize - 1] > 1) {
384 return -1;
385 }
386 if (osize == ysize) {
387 if (off->elements[osize - 1] > by->elements[ysize - 1]) return 1;
388 if (off->elements[osize - 1] < by->elements[ysize - 1] - 1) return -1;
389 }
390 } else {
391 /* xsize == ysize */
392 unsigned long w;
393 int c = 0;
394 if (osize > ysize) return 1;
395 if (bx->elements[xsize - 1] > by->elements[ysize - 1]) return 1;
396 if (osize < xsize) {
397 if (bx->elements[xsize - 1] < by->elements[ysize - 1] - 1) return -1;
398 } else {
399 /* osize == ysize */
400 unsigned long xx = bx->elements[xsize - 1], oo = off->elements[osize-1];
401 UADD(w, c, xx, oo);
402 if (c > 0 || w > by->elements[ysize - 1]) return 1;
403 if (w < by->elements[ysize - 1] - 1) return -1;
404 }
405 }
406 tsize = bignum_safe_size_for_add(bx, off);
407 ALLOC_TEMP_BIGNUM(br, tsize);
408 bignum_add_int(br, bx, off);
409
410 if (SG_BIGNUM_GET_COUNT(br) < SG_BIGNUM_GET_COUNT(by)) return -1;
411 for (i = SG_BIGNUM_GET_COUNT(br) - 1; i >= 0; i--) {
412 if (i >= SG_BIGNUM_GET_COUNT(by)) {
413 if (br->elements[i]) return 1;
414 continue;
415 }
416 if (br->elements[i] < by->elements[i]) return -1;
417 if (br->elements[i] > by->elements[i]) return 1;
418 }
419 return 0;
420 }
421
u64_to_double(uint64_t u)422 static double u64_to_double(uint64_t u)
423 {
424 union { double f64; uint64_t u64; } d;
425 d.u64 = u;
426 return d.f64;
427 }
428
lsb(SgBignum * b)429 static long lsb(SgBignum *b)
430 {
431 long count = SG_BIGNUM_GET_COUNT(b), i;
432 long l;
433 if (count == 0) return -1;
434
435 for (i = 0; i < count; i++) {
436 l = b->elements[i];
437 if (l != 0) break;
438 }
439 return (i<<SHIFT_MAGIC) + ntz(l);
440 }
441
Sg_BignumToDouble(SgBignum * b)442 double Sg_BignumToDouble(SgBignum *b)
443 {
444 long count = SG_BIGNUM_GET_COUNT(b);
445 long exponent, shift, increment;
446 uint64_t twiceSigFloor, sigFloor, sigRounded, bits;
447 SgObject o, ab;
448
449 if (count == 0) return 0.0;
450
451 exponent = Sg_BignumBitSize(b) - 1;
452 if (exponent < 64 - 1) {
453 return Sg_BignumToS64(b, SG_CLAMP_NONE, NULL);
454 } else if (exponent > 1023) {
455 return SG_BIGNUM_GET_SIGN(b) > 0
456 ? u64_to_double(0x7ff0000000000000ULL)
457 : u64_to_double(0xfff0000000000000ULL);
458 }
459 shift = exponent - 53; /* significand width */
460
461 /* TODO improve performance... */
462 ab = Sg_Abs(b);
463 o = Sg_BignumShiftRight(ab, shift);
464 if (SG_INTP(o)) {
465 twiceSigFloor = SG_INT_VALUE(o);
466 } else {
467 twiceSigFloor = Sg_BignumToU64(o, SG_CLAMP_NONE, NULL);
468 }
469 sigFloor = twiceSigFloor >> 1;
470 sigFloor &= 0x000FFFFFFFFFFFFFL;
471 increment = (twiceSigFloor & 1) != 0
472 && ((sigFloor & 1) != 0 || lsb(ab) < shift);
473 sigRounded = increment ? sigFloor + 1: sigFloor;
474 bits = (uint64_t)(exponent + 1023) << (53 - 1);
475 bits += sigRounded;
476 bits |= SG_BIGNUM_GET_SIGN(b) & 0x8000000000000000L;
477 return u64_to_double(bits);
478 }
479
bn_norm_pred(SgBignum * bn)480 static inline int bn_norm_pred(SgBignum *bn)
481 {
482 long bn_count = SG_BIGNUM_GET_COUNT(bn);
483 return (bn_count == 0) || (bn->elements[bn_count - 1] != 0);
484 }
485
Sg_BignumToInteger(SgBignum * bn)486 SgObject Sg_BignumToInteger(SgBignum *bn)
487 {
488 ASSERT(bn_norm_pred(bn));
489 ASSERT(SG_BIGNUM_GET_SIGN(bn) != 0);
490 if (SG_BIGNUM_GET_COUNT(bn) == 0) return SG_MAKE_INT(0);
491 if (SG_BIGNUM_GET_COUNT(bn) == 1) {
492 unsigned long n = bn->elements[0];
493 if (SG_BIGNUM_GET_SIGN(bn) < 0) {
494 if (n < SG_ULONG_MAX) {
495 n = -n;
496 } else {
497 /* more than long */
498 return Sg_BignumCopy(bn);
499 }
500 }
501 if ((n >= SG_INT_MIN) && (n <= SG_INT_MAX)) return SG_MAKE_INT(n);
502 }
503 return Sg_BignumCopy(bn);
504 }
505
Sg_BignumToSI(SgBignum * b,int clamp,int * oor)506 long Sg_BignumToSI(SgBignum *b, int clamp, int *oor)
507 {
508 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
509 if (SG_BIGNUM_GET_SIGN(b) >= 0) {
510 if (b->elements[0] > LONG_MAX || SG_BIGNUM_GET_COUNT(b) >= 2) {
511 if (clamp & SG_CLAMP_HI) return LONG_MAX;
512 else goto err;
513 } else {
514 return (long)b->elements[0];
515 }
516 } else {
517 if (b->elements[0] > (unsigned long)LONG_MAX + 1 ||
518 SG_BIGNUM_GET_COUNT(b) >= 2) {
519 if (clamp & SG_CLAMP_LO) return LONG_MIN;
520 else goto err;
521 } else {
522 return -(long)b->elements[0];
523 }
524 }
525 err:
526 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = TRUE;
527 else {
528 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
529 }
530 return 0;
531 }
532
Sg_BignumToUI(SgBignum * b,int clamp,int * oor)533 unsigned long Sg_BignumToUI(SgBignum *b, int clamp, int *oor)
534 {
535 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
536 if (SG_BIGNUM_GET_SIGN(b) >= 0) {
537 if (SG_BIGNUM_GET_COUNT(b) >= 2) {
538 if (clamp & SG_CLAMP_HI) return SG_ULONG_MAX;
539 else goto err;
540 } else {
541 return b->elements[0];
542 }
543 } else {
544 if (clamp & SG_CLAMP_LO) return 0;
545 else goto err;
546 }
547 err:
548 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = TRUE;
549 else {
550 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
551 }
552 return 0;
553 }
554
555 #if SIZEOF_LONG >= 8
Sg_BignumToS32(SgBignum * b,int clamp,int * oor)556 int32_t Sg_BignumToS32(SgBignum *b, int clamp, int *oor)
557 {
558 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
559 if (SG_BIGNUM_GET_SIGN(b) >= 0) {
560 if (b->elements[0] > INT32_MAX || SG_BIGNUM_GET_COUNT(b) >= 2) {
561 if (clamp & SG_CLAMP_HI) return INT32_MAX;
562 else goto err;
563 } else {
564 return (int32_t)b->elements[0];
565 }
566 } else {
567 if (b->elements[0] > (uint32_t)INT32_MAX + 1 ||
568 SG_BIGNUM_GET_COUNT(b) >= 2) {
569 if (clamp & SG_CLAMP_LO) return INT32_MIN;
570 else goto err;
571 } else {
572 return -(int32_t)b->elements[0];
573 }
574 }
575 err:
576 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = TRUE;
577 else {
578 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
579 }
580 return 0;
581 }
582
Sg_BignumToU32(SgBignum * b,int clamp,int * oor)583 uint32_t Sg_BignumToU32(SgBignum *b, int clamp, int *oor)
584 {
585 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
586 if (SG_BIGNUM_GET_SIGN(b) > 0) {
587 if (b->elements[0] <= UINT32_MAX) {
588 return (uint32_t)b->elements[0];
589 } else {
590 if (!(clamp & SG_CLAMP_HI)) goto err;
591 return UINT32_MAX;
592 }
593 } else {
594 if (clamp & SG_CLAMP_LO) return 0;
595 else goto err;
596 }
597 err:
598 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = TRUE;
599 else {
600 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
601 }
602 return 0;
603 }
604
605 #else
Sg_BignumToS64(SgBignum * b,int clamp,int * oor)606 int64_t Sg_BignumToS64(SgBignum *b, int clamp, int *oor)
607 {
608 int64_t r = 0;
609 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
610 if (SG_BIGNUM_GET_SIGN(b) > 0) {
611 if (SG_BIGNUM_GET_COUNT(b) == 1) {
612 r = b->elements[0];
613 } else if (SG_BIGNUM_GET_COUNT(b) > 2 || b->elements[1] > LONG_MAX) {
614 if (!(clamp & SG_CLAMP_HI)) goto err;
615 r = (((int64_t)LONG_MAX) << 32) + (int64_t)ULONG_MAX;
616 } else {
617 r = ((int64_t)b->elements[1] << 32) + (uint64_t)b->elements[0];
618 }
619 } else {
620 if (SG_BIGNUM_GET_COUNT(b) == 1) {
621 r = -(int64_t)b->elements[0];
622 } else if (SG_BIGNUM_GET_COUNT(b) > 2 ||
623 b->elements[1] > (unsigned long)LONG_MAX+1) {
624 if (!(clamp & SG_CLAMP_HI)) goto err;
625 r = (((int64_t)LONG_MAX + 1) << 32);
626 } else {
627 r = -(int64_t)(((int64_t)b->elements[1] << 32)+(uint64_t)b->elements[0]);
628 }
629 }
630 return r;
631 err:
632 if (clamp == SG_CLAMP_NONE && oor != NULL) {
633 *oor = TRUE;
634 } else {
635 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
636 }
637 return r;
638 }
639
Sg_BignumToU64(SgBignum * b,int clamp,int * oor)640 uint64_t Sg_BignumToU64(SgBignum *b, int clamp, int *oor)
641 {
642 uint64_t r = 0;
643 if (clamp == SG_CLAMP_NONE && oor != NULL) *oor = FALSE;
644 if (SG_BIGNUM_GET_SIGN(b) > 0) {
645 if (SG_BIGNUM_GET_COUNT(b) > 2) {
646 if (!(clamp & SG_CLAMP_HI)) goto err;
647 r = (((uint64_t)ULONG_MAX) << 32) + (uint64_t)ULONG_MAX;
648
649 } else if (SG_BIGNUM_GET_COUNT(b) == 2) {
650 r = (((uint64_t)b->elements[1]) << 32) + (uint64_t)b->elements[0];
651 } else {
652 r = (uint64_t)b->elements[0];
653 }
654 } else {
655 if (!(clamp & SG_CLAMP_LO)) goto err;
656 }
657 return r;
658 err:
659 if (clamp == SG_CLAMP_NONE && oor != NULL) {
660 *oor = TRUE;
661 } else {
662 Sg_Error(UC("argument out of range: %S"), SG_OBJ(b));
663 }
664 return r;
665 }
666
667 #endif
668
669
bignum_2scmpl(SgBignum * br)670 static inline SgBignum* bignum_2scmpl(SgBignum *br)
671 {
672 mp_2scmpl(br->elements, br->size);
673 return br;
674 }
675
Sg_BignumComplement(SgBignum * bx)676 SgObject Sg_BignumComplement(SgBignum *bx)
677 {
678 SgBignum *r = SG_BIGNUM(Sg_BignumCopy(bx));
679 return SG_OBJ(bignum_2scmpl(r));
680 }
681
682
Sg_BignumBitCount(SgBignum * b)683 long Sg_BignumBitCount(SgBignum *b)
684 {
685 unsigned long *bits;
686 SgBignum *z;
687 long size;
688 if (SG_BIGNUM_GET_SIGN(b) == 0) {
689 return 0;
690 } else if (SG_BIGNUM_GET_SIGN(b) > 0) {
691 z = b;
692 goto calc_bit_count;
693 }
694 ALLOC_TEMP_BIGNUM(z, SG_BIGNUM_GET_COUNT(b));
695 bignum_copy(z, b);
696 bignum_2scmpl(z);
697
698 calc_bit_count:
699 size = SG_BIGNUM_GET_COUNT(z) * WORD_BITS;
700 bits = z->elements;
701 if (SG_BIGNUM_GET_SIGN(b) > 0) {
702 return Sg_BitsCount1(bits, 0, size);
703 } else {
704 return ~Sg_BitsCount0(bits, 0, size);
705 }
706 }
707
Sg_BignumBitSize(SgBignum * b)708 long Sg_BignumBitSize(SgBignum *b)
709 {
710 if (SG_BIGNUM_GET_SIGN(b) == 0) return 0;
711 return mp_bit_size(b->elements, b->size);
712 }
713
Sg_BignumFirstBitSet(SgBignum * b)714 long Sg_BignumFirstBitSet(SgBignum *b)
715 {
716 long bit = 0, i, size;
717 SgBignum *z;
718 if (SG_BIGNUM_GET_SIGN(b) == 0) return 0;
719 else if (SG_BIGNUM_GET_SIGN(b) > 0) {
720 z = b;
721 goto calc_bit_set;
722 }
723 ALLOC_TEMP_BIGNUM_REC(z, SG_BIGNUM_GET_COUNT(b));
724 bignum_copy(z, b);
725 bignum_2scmpl(z);
726
727 calc_bit_set:
728 size = SG_BIGNUM_GET_COUNT(z);
729 for (i = 0; i < size; i++) {
730 unsigned long n = z->elements[i];
731 if (n == 0) { bit += WORD_BITS; continue; }
732 bit += ntz(n);
733 return bit;
734 }
735 /* Sg_Write(b, Sg_StandardErrorPort(), 0); */
736 ASSERT(FALSE);
737 return -1; /* dummy */
738 }
739
bignum_test_bit(SgBignum * b,long p)740 static inline int bignum_test_bit(SgBignum *b, long p)
741 {
742 long pos = p >> SHIFT_MAGIC;
743 ulong v;
744 if (pos >= (long)SG_BIGNUM_GET_COUNT(b)) {
745 return FALSE;
746 }
747 v = b->elements[pos];
748 if (pos) {
749 p %= WORD_BITS;
750 }
751 return (v & (1L << p)) != 0;
752 }
753
Sg_BignumBitSetP(SgBignum * b,long n)754 int Sg_BignumBitSetP(SgBignum *b, long n)
755 {
756 return bignum_test_bit(b, n);
757 }
758
Sg_BignumAbsCmp(SgBignum * bx,SgBignum * by)759 int Sg_BignumAbsCmp(SgBignum *bx, SgBignum *by)
760 {
761 long i;
762 long xsize = SG_BIGNUM_GET_COUNT(bx);
763 long ysize = SG_BIGNUM_GET_COUNT(by);
764
765 if (xsize < ysize) return -1;
766 if (xsize > ysize) return 1;
767 for (i = (int)xsize-1; i >= 0; i--) {
768 if (bx->elements[i] < by->elements[i]) return -1;
769 if (bx->elements[i] > by->elements[i]) return 1;
770 }
771 return 0;
772 }
773
774 /* returns bignum, so weird name... */
Sg_BignumAshSI(long lx,long count)775 SgObject Sg_BignumAshSI(long lx, long count)
776 {
777 SgBignum *x;
778 if (count > 0 && lx == 1) {
779 long size = SG_LEFT_SHIFT_SPACE(1, count);
780 long words = count / WORD_BITS;
781 long nbits = count % WORD_BITS;
782 ulong t = (ulong)lx << nbits;
783 x = make_bignum(size);
784 x->elements[words] = t;
785 return Sg_NormalizeBignum(x);
786 } else {
787 ALLOC_TEMP_BIGNUM(x, 1);
788 if (lx < 0) {
789 x->elements[0] = -lx;
790 x->sign = -1;
791 } else {
792 x->elements[0] = lx;
793 }
794 return Sg_BignumAsh(x, count);
795 }
796 }
797
Sg_BignumAsh(SgBignum * b,long count)798 SgObject Sg_BignumAsh(SgBignum *b, long count)
799 {
800 if (count == 0) return Sg_NormalizeBignum(b);
801 else if (count > 0) return Sg_BignumShiftLeft(b, count);
802 else return Sg_BignumShiftRight(b, -count);
803 }
804
bignum_lshift(SgBignum * br,SgBignum * bx,long amount)805 static SgBignum* bignum_lshift(SgBignum *br, SgBignum *bx, long amount)
806 {
807 mp_lshift(br->elements, br->size, bx->elements, bx->size, amount);
808 if (br != bx) {
809 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
810 }
811 return br;
812 }
813
Sg_BignumShiftLeft(SgBignum * b,long shift)814 SgObject Sg_BignumShiftLeft(SgBignum *b, long shift)
815 {
816 int rsize = SG_LEFT_SHIFT_SPACE(b->size, shift);
817 SgBignum *r = make_bignum(rsize);
818 return Sg_NormalizeBignum(bignum_lshift(r, b, shift));
819 }
820
bignum_rshift(SgBignum * br,SgBignum * bx,long amount)821 static SgBignum* bignum_rshift(SgBignum *br, SgBignum *bx, long amount)
822 {
823 ulong size = mp_rshift(br->elements, br->size,
824 bx->elements, bx->size, amount);
825
826 SG_BIGNUM_SET_COUNT(br, size);
827 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
828 return br;
829 }
830
Sg_BignumShiftRight(SgBignum * b,long shift)831 SgObject Sg_BignumShiftRight(SgBignum *b, long shift)
832 {
833 long rsize = SG_BIGNUM_GET_COUNT(b) + (-shift)/WORD_BITS;
834 if (rsize < 1) {
835 if (SG_BIGNUM_GET_SIGN(b) < 0) {
836 return SG_MAKE_INT(-1);
837 } else {
838 return SG_MAKE_INT(0);
839 }
840 } else {
841 if (SG_BIGNUM_GET_SIGN(b) < 0) {
842 SgObject r = Sg_Quotient(Sg_Add(SG_OBJ(b), SG_MAKE_INT(1)),
843 Sg_Ash(SG_MAKE_INT(1), shift),
844 NULL);
845 return Sg_Add(r, SG_MAKE_INT(-1));
846 } else {
847 SgBignum *r = make_bignum(rsize);
848 return Sg_NormalizeBignum(bignum_rshift(r, b, shift));
849 }
850 }
851 }
852
853 #define DEF_BIGNUM_LOG_OP(name, op) \
854 static inline SgBignum* name(SgBignum *z, SgBignum *x, SgBignum *y, \
855 int x2sc, int y2sc) \
856 { \
857 long i; \
858 long xs = SG_BIGNUM_GET_COUNT(x); \
859 long ys = SG_BIGNUM_GET_COUNT(y); \
860 long zs = SG_BIGNUM_GET_COUNT(z); \
861 for (i = zs-1; i >= 0; i--) { \
862 ulong lx = (i < xs) ? x->elements[i] : (x2sc ? SG_ULONG_MAX : 0); \
863 ulong ly = (i < ys) ? y->elements[i] : (y2sc ? SG_ULONG_MAX : 0); \
864 z->elements[i] = lx op ly; \
865 } \
866 return z; \
867 }
868
869 DEF_BIGNUM_LOG_OP(bignum_and, &)
870
Sg_BignumLogAnd(SgBignum * x,SgBignum * y)871 SgObject Sg_BignumLogAnd(SgBignum *x, SgBignum *y)
872 {
873 long xsize = SG_BIGNUM_GET_COUNT(x);
874 int xsign = SG_BIGNUM_GET_SIGN(x);
875 long ysize = SG_BIGNUM_GET_COUNT(y);
876 int ysign = SG_BIGNUM_GET_SIGN(y);
877 long zsize, minsize = min(xsize, ysize);
878 SgBignum *xx, *yy, *z;
879
880 /* handle obvious case first */
881 if (xsign == 0 || ysign == 0) return SG_MAKE_INT(0);
882
883 if (xsign > 0) {
884 if (ysign > 0) {
885 /* (+, +) */
886 z = bignum_and(make_bignum(minsize), x, y, FALSE, FALSE);
887 return Sg_NormalizeBignum(z);
888 } else {
889 /* (+, -) */
890 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
891 bignum_copy(yy, y);
892 bignum_2scmpl(yy);
893 z = bignum_and(make_bignum(xsize), x, yy, FALSE, TRUE);
894 return Sg_NormalizeBignum(z);
895 }
896 } else {
897 if (ysign > 0) {
898 /* (-, +) */
899 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
900 bignum_copy(xx, x);
901 bignum_2scmpl(xx);
902 z = bignum_and(make_bignum(ysize), xx, y, TRUE, FALSE);
903 return Sg_NormalizeBignum(z);
904 } else {
905 /* (-, -) */
906 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
907 bignum_copy(xx, x);
908 bignum_2scmpl(xx);
909 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
910 bignum_copy(yy, y);
911 bignum_2scmpl(yy);
912
913 zsize = max(xsize, ysize);
914 z = bignum_and(make_bignum(zsize), xx, yy, TRUE, TRUE);
915 SG_BIGNUM_SET_SIGN(z, -1);
916 bignum_2scmpl(z);
917 return Sg_NormalizeBignum(z);
918 }
919 }
920 }
921
922 DEF_BIGNUM_LOG_OP(bignum_ior, |)
923
Sg_BignumLogIor(SgBignum * x,SgBignum * y)924 SgObject Sg_BignumLogIor(SgBignum *x, SgBignum *y)
925 {
926 long xsize = SG_BIGNUM_GET_COUNT(x);
927 int xsign = SG_BIGNUM_GET_SIGN(x);
928 long ysize = SG_BIGNUM_GET_COUNT(y);
929 int ysign = SG_BIGNUM_GET_SIGN(y);
930 long zsize, minsize = min(xsize, ysize);
931 SgBignum *xx, *yy, *z;
932
933 /* handle 0 case first */
934 if (xsign == 0 || ysign == 0) {
935 if (xsign) return Sg_NormalizeBignum(x);
936 if (ysign) return Sg_NormalizeBignum(y);
937 return SG_MAKE_INT(0);
938 }
939
940 if (xsign > 0) {
941 if (ysign > 0) {
942 /* (+, +) */
943 zsize = max(xsize, ysize);
944 z = bignum_ior(make_bignum(zsize), x, y, FALSE, FALSE);
945 return Sg_NormalizeBignum(z);
946 } else {
947 /* (+, -) */
948 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
949 bignum_copy(yy, y);
950 bignum_2scmpl(yy);
951 z = bignum_ior(make_bignum(ysize), x, yy, FALSE, TRUE);
952 SG_BIGNUM_SET_SIGN(z, -1);
953 bignum_2scmpl(z);
954 return Sg_NormalizeBignum(z);
955 }
956 } else {
957 if (ysign > 0) {
958 /* (-, +) */
959 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
960 bignum_copy(xx, x);
961 bignum_2scmpl(xx);
962 z = bignum_ior(make_bignum(xsize), xx, y, TRUE, FALSE);
963 SG_BIGNUM_SET_SIGN(z, -1);
964 bignum_2scmpl(z);
965 return Sg_NormalizeBignum(z);
966 } else {
967 /* (-, -) */
968 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
969 bignum_copy(xx, x);
970 bignum_2scmpl(xx);
971 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
972 bignum_copy(yy, y);
973 bignum_2scmpl(yy);
974
975 z = bignum_ior(make_bignum(minsize), xx, yy, TRUE, TRUE);
976 SG_BIGNUM_SET_SIGN(z, -1);
977 bignum_2scmpl(z);
978 return Sg_NormalizeBignum(z);
979 }
980 }
981 }
982
983 DEF_BIGNUM_LOG_OP(bignum_xor, ^)
984
Sg_BignumLogXor(SgBignum * x,SgBignum * y)985 SgObject Sg_BignumLogXor(SgBignum *x, SgBignum *y)
986 {
987 long xsize = SG_BIGNUM_GET_COUNT(x);
988 int xsign = SG_BIGNUM_GET_SIGN(x);
989 long ysize = SG_BIGNUM_GET_COUNT(y);
990 int ysign = SG_BIGNUM_GET_SIGN(y);
991 long zsize = max(xsize, ysize);
992 SgBignum *xx, *yy, *z;
993
994 /* handle 0 case first */
995 if (xsign == 0 || ysign == 0) {
996 if (xsign) return Sg_NormalizeBignum(x);
997 if (ysign) return Sg_NormalizeBignum(y);
998 return SG_MAKE_INT(0);
999 }
1000
1001 if (xsign > 0) {
1002 if (ysign > 0) {
1003 /* (+,+) */
1004 z = bignum_xor(make_bignum(zsize), x, y, FALSE, FALSE);
1005 return Sg_NormalizeBignum(z);
1006 } else {
1007 /* (+,-) */
1008 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
1009 bignum_copy(yy, y);
1010 bignum_2scmpl(yy);
1011 z = bignum_xor(make_bignum(zsize), x, yy, FALSE, TRUE);
1012 SG_BIGNUM_SET_SIGN(z, -1);
1013 bignum_2scmpl(z);
1014 return Sg_NormalizeBignum(z);
1015 }
1016 } else {
1017 if (ysign > 0) {
1018 /* (-,+) */
1019 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
1020 bignum_copy(xx, x);
1021 bignum_2scmpl(xx);
1022 z = bignum_xor(make_bignum(zsize), xx, y, TRUE, FALSE);
1023 SG_BIGNUM_SET_SIGN(z, -1);
1024 bignum_2scmpl(z);
1025 return Sg_NormalizeBignum(z);
1026 } else {
1027 /* (-,-) */
1028 ALLOC_TEMP_BIGNUM_REC(xx, xsize);
1029 bignum_copy(xx, x);
1030 bignum_2scmpl(xx);
1031 ALLOC_TEMP_BIGNUM_REC(yy, ysize);
1032 bignum_copy(yy, y);
1033 bignum_2scmpl(yy);
1034
1035 z = bignum_xor(make_bignum(zsize), xx, yy, TRUE, TRUE);
1036 return Sg_NormalizeBignum(z);
1037 }
1038 }
1039 }
1040
1041 #define DEF_SI_LOGOP(op) \
1042 SgObject SG_CPP_CAT(op, SI)(SgBignum *x, long y) \
1043 { \
1044 SgBignum *by; \
1045 ALLOC_TEMP_BIGNUM(by, 1); \
1046 if (y == 0) { \
1047 SG_BIGNUM_SET_SIGN(by, 0); \
1048 } else if (y == LONG_MIN) { \
1049 by->elements[0] = (unsigned long)LONG_MAX + 1; \
1050 SG_BIGNUM_SET_SIGN(by, -1); \
1051 } else if (y < 0) { \
1052 by->elements[0] = -y; \
1053 SG_BIGNUM_SET_SIGN(by, -1); \
1054 } else { \
1055 by->elements[0] = y; \
1056 SG_BIGNUM_SET_SIGN(by, 1); \
1057 } \
1058 return op(x, by); \
1059 }
1060
1061 DEF_SI_LOGOP(Sg_BignumLogAnd)
DEF_SI_LOGOP(Sg_BignumLogIor)1062 DEF_SI_LOGOP(Sg_BignumLogIor)
1063 DEF_SI_LOGOP(Sg_BignumLogXor)
1064
1065 static long bignum_safe_size_for_add(SgBignum *x, SgBignum *y)
1066 {
1067 long xsize = SG_BIGNUM_GET_COUNT(x);
1068 long ysize = SG_BIGNUM_GET_COUNT(y);
1069 return mp_safe_size_for_add(x->elements, xsize, y->elements, ysize);
1070 }
1071
bignum_add_int(SgBignum * br,SgBignum * bx,SgBignum * by)1072 static SgBignum* bignum_add_int(SgBignum *br, SgBignum *bx, SgBignum *by)
1073 {
1074 long rsize = SG_BIGNUM_GET_COUNT(br);
1075 long xsize = SG_BIGNUM_GET_COUNT(bx);
1076 long ysize = SG_BIGNUM_GET_COUNT(by);
1077 long i;
1078
1079 /* handle 0 first */
1080 if (xsize == 0) {
1081 if (ysize == 0) {
1082 SG_BIGNUM_SET_ZERO(br);
1083 return br;
1084 }
1085 for (i = 0; i < ysize; i++) {
1086 br->elements[i] = by->elements[i];
1087 }
1088 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(by));
1089 } else if (ysize == 0) {
1090 for (i = 0; i < xsize; i++) {
1091 br->elements[i] = bx->elements[i];
1092 }
1093 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1094 } else {
1095 /* if x is shorter, swap it */
1096 if (xsize < ysize) {
1097 SgBignum *t = bx;
1098 ulong tsize = xsize;
1099 bx = by; by = t;
1100 xsize = ysize; ysize = tsize;
1101 }
1102 mp_add(br->elements, rsize, bx->elements, xsize, by->elements, ysize);
1103 }
1104 return br;
1105 }
1106
bignum_sub_int(SgBignum * br,SgBignum * bx,SgBignum * by)1107 static SgBignum* bignum_sub_int(SgBignum *br, SgBignum *bx, SgBignum *by)
1108 {
1109 long rsize = SG_BIGNUM_GET_COUNT(br);
1110 long xsize = SG_BIGNUM_GET_COUNT(bx);
1111 long ysize = SG_BIGNUM_GET_COUNT(by);
1112 long i;
1113
1114 /* handle 0 first */
1115 if (xsize == 0) {
1116 if (ysize == 0) {
1117 SG_BIGNUM_SET_ZERO(br);
1118 return br;
1119 }
1120 for (i = 0; i < ysize; i++) {
1121 br->elements[i] = by->elements[i];
1122 }
1123 /* 0 - n must be - */
1124 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(by));
1125 } else if (ysize == 0) {
1126 for (i = 0; i < xsize; i++) {
1127 br->elements[i] = bx->elements[i];
1128 }
1129 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1130 } else {
1131 int flip = FALSE, c;
1132 if (xsize < ysize) {
1133 SgBignum *t = bx;
1134 bx = by; by = t;
1135 xsize = SG_BIGNUM_GET_COUNT(bx);
1136 ysize = SG_BIGNUM_GET_COUNT(by);
1137 flip = TRUE;
1138 }
1139 c = mp_sub(br->elements, rsize, bx->elements, xsize, by->elements, ysize);
1140 if (flip || c) {
1141 if (c) bignum_2scmpl(br);
1142 SG_BIGNUM_SET_SIGN(br, -SG_BIGNUM_GET_SIGN(br));
1143 }
1144 }
1145 return br;
1146 }
1147
bignum_add(SgBignum * bx,SgBignum * by)1148 static SgBignum* bignum_add(SgBignum *bx, SgBignum *by)
1149 {
1150 long rsize = bignum_safe_size_for_add(bx, by);
1151 SgBignum *br = make_bignum(rsize);
1152 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1153 if (SG_BIGNUM_GET_SIGN(bx) == SG_BIGNUM_GET_SIGN(by)) {
1154 bignum_add_int(br, bx, by);
1155 } else {
1156 bignum_sub_int(br, bx, by);
1157 }
1158 return br;
1159 }
bignum_sub(SgBignum * bx,SgBignum * by)1160 static SgBignum* bignum_sub(SgBignum *bx, SgBignum *by)
1161 {
1162 long rsize = bignum_safe_size_for_add(bx, by);
1163 SgBignum *br = make_bignum(rsize);
1164 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1165 if (SG_BIGNUM_GET_SIGN(bx) == SG_BIGNUM_GET_SIGN(by)) {
1166 bignum_sub_int(br, bx, by);
1167 } else {
1168 bignum_add_int(br, bx, by);
1169 if (SG_BIGNUM_GET_SIGN(bx) == 0) {
1170 /*flip sign*/
1171 SG_BIGNUM_SET_SIGN(br, -SG_BIGNUM_GET_SIGN(by));
1172 }
1173 }
1174 return br;
1175 }
1176
bignum_add_si(SgBignum * bx,long y)1177 static SgBignum* bignum_add_si(SgBignum *bx, long y)
1178 {
1179 ulong yabs;
1180 int ysign;
1181 SgBignum *br;
1182
1183 if (y == 0) return bx; /* short cut */
1184
1185 yabs = ((y < 0) ? -y : y);
1186 ysign = ((y < 0) ? -1: 1);
1187 br = make_bignum(bx->size + 1);
1188 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1189 /* dispatch here */
1190 if (SG_BIGNUM_GET_SIGN(bx) == ysign) {
1191 mp_add_ul(br->elements, br->size, bx->elements, bx->size, yabs);
1192 } else {
1193 mp_sub_ul(br->elements, br->size, bx->elements, bx->size, yabs);
1194 }
1195 return br;
1196 }
1197
Sg_BignumAdd(SgBignum * a,SgBignum * b)1198 SgObject Sg_BignumAdd(SgBignum *a, SgBignum *b)
1199 {
1200 return Sg_NormalizeBignum(bignum_add(a, b));
1201 }
1202
1203
Sg_BignumAddSI(SgBignum * a,long b)1204 SgObject Sg_BignumAddSI(SgBignum *a, long b)
1205 {
1206 return Sg_NormalizeBignum(bignum_add_si(a, b));
1207 }
1208
Sg_BignumSub(SgBignum * a,SgBignum * b)1209 SgObject Sg_BignumSub(SgBignum *a, SgBignum *b)
1210 {
1211 return Sg_NormalizeBignum(bignum_sub(a, b));
1212 }
1213
Sg_BignumSubSI(SgBignum * a,long b)1214 SgObject Sg_BignumSubSI(SgBignum *a, long b)
1215 {
1216 return Sg_NormalizeBignum(bignum_add_si(a, -b));
1217 }
1218
1219 /* whatever i did, it's slow. */
1220 #if 0
1221 #include <emmintrin.h>
1222
1223 typedef struct
1224 {
1225 ulong r0lo; /* r0 low */
1226 ulong r0hi; /* r0 high */
1227 ulong r1lo; /* r1 low */
1228 ulong r1hi; /* r1 high */
1229 } r4;
1230
1231 static inline SgBignum* bignum_mul_word(SgBignum *br, SgBignum *bx,
1232 unsigned long y)
1233 {
1234 int size = SG_BIGNUM_GET_COUNT(bx);
1235 /* do trivial case first */
1236 if (size == 1) {
1237 /* only one element in bx */
1238 udlong p;
1239 p = (udlong)bx->elements[0] * y;
1240 br->elements[0] = (ulong)p;
1241 br->elements[1] = (ulong)(p >> WORD_BITS);
1242 return br;
1243 } else if (size == 2) {
1244 /* only 2 elements in bx */
1245 udlong p;
1246 p = (udlong)bx->elements[0] * y;
1247 br->elements[0] = (ulong)p;
1248 p = (udlong)bx->elements[1] * y + (ulong)(p >> WORD_BITS);
1249 br->elements[1] = (ulong)p;
1250 br->elements[2] = (ulong)(p >> WORD_BITS);
1251 return br;
1252 } else {
1253 /* more than 3 elements means at least one loop */
1254 #if 1
1255 /* try use SSE as mere (64 bit) register.
1256 why is there no 64x64 multiplication? */
1257 int i;
1258 __m128i p, x, yy, t;
1259 yy = _mm_cvtsi32_si128(y);
1260 for (i = 0; i < size; i++) {
1261 x = _mm_cvtsi32_si128(bx->elements[i]);
1262 t = _mm_cvtsi32_si128(br->elements[i]); /* carry is here */
1263 p = _mm_add_epi64(_mm_mul_epu32(x, yy), t);
1264 _mm_storel_epi64(&br->elements[i], p);
1265 }
1266 /* the last carry is already stored */
1267 return br;
1268 #else
1269 /* following was actually slow... */
1270 /* debug */
1271 /* #define DEBUG_SSE2 */
1272 #ifdef DEBUG_SSE2
1273 #define dump_m128(m) \
1274 do { \
1275 r4 r_ __attribute__((aligned(16))); \
1276 _mm_store_si128((__m128i *)&r_, (m)); \
1277 fprintf(stderr, "%08lx:%08lx:%08lx:%08lx:%s\n", \
1278 r_.r0hi, r_.r0lo, r_.r1hi, r_.r1lo, #m); \
1279 } while (0)
1280 #else
1281 #define dump_m128(m) /* dummy */
1282 #endif
1283 /*
1284 p0(h,l) = x0 * y
1285 p1(h,l) = x1 * y + p0(h)
1286 p2(h,l) = x2 * y + p1(h)
1287
1288 i.e.)
1289 p0 = 0x0000002e257a5000
1290 p0(h) = 0x0000002e
1291 x1y = 0x0000002478ce03be
1292 p1 = 0x0000002478ce03ec
1293 p1(h) = 0x00000024 = high(low(x1y) + p0(h)) + high(x1y)
1294
1295 pn(h) = high(low(xn * y) + pn-1(h)) + high(xn * y)
1296
1297 c0 register holds pn-1(h)
1298 c1 register holds pn(h)
1299 */
1300 #define compute_sse(x0, x1) \
1301 do { \
1302 /* compute carry from previous (r1 register) */ \
1303 /* carry pi-1 and pi-2 */ \
1304 c1 = _mm_shuffle_epi32(p, _MM_SHUFFLE(0,1,2,3)); \
1305 c1 = _mm_and_si128(c1, ml); /* mask */ \
1306 dump_m128(c1); \
1307 xx = _mm_set_epi64x((x1), (x0)); \
1308 p = _mm_mul_epu32(xx, yy); /* (x0 * y0) & (x1 * y1) */ \
1309 dump_m128(p); \
1310 c0 = _mm_and_si128(p, ml); /* low(x1y) */ \
1311 c0 = _mm_add_epi64(c0, c1); /* low(x1y) + p0(h) */ \
1312 dump_m128(c0); \
1313 xx = _mm_srli_epi64(c0, WORD_BITS); /* high(low(x1y) + p0(h)) */ \
1314 c0 = _mm_srli_epi64(p, WORD_BITS); /* high(x1y) */ \
1315 dump_m128(xx); \
1316 xx = _mm_add_epi64(xx, c0); /* high(low(x1y) +p0(h)) + high(x1y) */ \
1317 dump_m128(xx); \
1318 /* now xx contains p1(h) in r1 so shuffle it */ \
1319 /* p1 = x*y + (ulong)(p0>>32) */ \
1320 xx = _mm_and_si128(xx, mx); /* clear r0,r1,r2 */ \
1321 xx = _mm_shuffle_epi32(xx, _MM_SHUFFLE(1,0,3,2)); \
1322 c1 = _mm_and_si128(c1, mx); \
1323 dump_m128(xx); \
1324 c0 = _mm_or_si128(xx, c1); \
1325 dump_m128(c0); \
1326 p = _mm_add_epi64(p, c0); \
1327 dump_m128(p); \
1328 } while (0)
1329
1330 int i;
1331 __m128i p, xx, yy, c0, c1, ml, mx;
1332 __attribute__((aligned(16))) r4 pr = {0, 0, 0, 0};
1333
1334 yy = _mm_set_epi64x(y, y);
1335 ml = _mm_set_epi32(0, 0xffffffffUL, 0, 0xffffffffUL);
1336 mx = _mm_set_epi32(0, 0, 0, 0xffffffffUL);
1337 /* set 0 */
1338 p = _mm_setzero_si128();
1339 for (i = 0; i < size - 1; i += 2) {
1340 compute_sse(bx->elements[i], bx->elements[i+1]);
1341 _mm_store_si128((__m128i *)&pr, p);
1342 br->elements[i] = pr.r0lo;
1343 br->elements[i+1] = pr.r1lo;
1344 }
1345 /* do the rest */
1346 if (size & 1) {
1347 compute_sse(bx->elements[i], 0);
1348 _mm_store_si128((__m128i *)&pr, p);
1349 br->elements[i] = pr.r0lo;
1350 br->elements[i+1] = pr.r0hi;
1351 } else {
1352 br->elements[i] = pr.r1hi;
1353 }
1354 return br;
1355 #undef compute_sse
1356 #endif
1357 }
1358 }
1359 #else
bignum_mul_word(SgBignum * br,SgBignum * bx,unsigned long y)1360 static inline SgBignum* bignum_mul_word(SgBignum *br, SgBignum *bx,
1361 unsigned long y)
1362 {
1363 mp_mul_ul(br->elements, br->size, bx->elements, bx->size, y);
1364 return br;
1365 }
1366 #endif
1367
bignum_mul_int(SgBignum * br,SgBignum * bx,SgBignum * by)1368 static SgBignum* bignum_mul_int(SgBignum *br, SgBignum *bx, SgBignum *by)
1369 {
1370 long xlen = SG_BIGNUM_GET_COUNT(bx);
1371 long ylen = SG_BIGNUM_GET_COUNT(by);
1372 /* early check */
1373 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx) * SG_BIGNUM_GET_SIGN(by));
1374 mp_mul(br->elements, br->size, bx->elements, xlen, by->elements, ylen);
1375 return br;
1376 }
1377
bignum_mul(SgBignum * bx,SgBignum * by)1378 static SgBignum* bignum_mul(SgBignum *bx, SgBignum *by)
1379 {
1380 long xlen = SG_BIGNUM_GET_COUNT(bx);
1381 long ylen = SG_BIGNUM_GET_COUNT(by);
1382 SgBignum *br;
1383
1384 if (xlen == 0) return bx;
1385 if (ylen == 0) return by;
1386
1387 br = make_bignum(xlen + ylen);
1388 return bignum_mul_int(br, bx, by);
1389 }
1390
bignum_mul_si(SgBignum * bx,long y)1391 static SgBignum* bignum_mul_si(SgBignum *bx, long y)
1392 {
1393 SgBignum *br;
1394 unsigned long yabs;
1395
1396 if (y == 1L) return bx;
1397 else if (y == 0) {
1398 br = make_bignum(1);
1399 SG_BIGNUM_SET_SIGN(br, 0);
1400 br->elements[0] = 0;
1401 return br;
1402 } else if (y == -1L) {
1403 br = SG_BIGNUM(Sg_BignumCopy(bx));
1404 SG_BIGNUM_SET_SIGN(br, -SG_BIGNUM_GET_SIGN(bx));
1405 return br;
1406 } else {
1407 br = make_bignum(SG_BIGNUM_GET_COUNT(bx) + 1);
1408 yabs = (y < 0) ? -y : y;
1409 SG_BIGNUM_SET_SIGN(br, SG_BIGNUM_GET_SIGN(bx));
1410 bignum_mul_word(br, bx, yabs);
1411 if (y < 0) SG_BIGNUM_SET_SIGN(br, -SG_BIGNUM_GET_SIGN(br));
1412 return br;
1413 }
1414 }
1415
1416
Sg_BignumMul(SgBignum * a,SgBignum * b)1417 SgObject Sg_BignumMul(SgBignum *a, SgBignum *b)
1418 {
1419 return Sg_NormalizeBignum(bignum_mul(a, b));
1420 }
1421
Sg_BignumMulSI(SgBignum * a,long b)1422 SgObject Sg_BignumMulSI(SgBignum *a, long b)
1423 {
1424 return Sg_NormalizeBignum(bignum_mul_si(a, b));
1425 }
1426
Sg_BignumSquare(SgBignum * bx)1427 SgObject Sg_BignumSquare(SgBignum *bx)
1428 {
1429 long len = SG_BIGNUM_GET_COUNT(bx);
1430 SgBignum *br = make_bignum(len<<1);
1431 /* square_to_len(bx->elements, len, br->elements); */
1432 mp_square(br->elements, br->size, bx->elements, len);
1433 SG_BIGNUM_SET_SIGN(br, 1);
1434 return Sg_NormalizeBignum(br);
1435 }
1436
bignum_gdiv_rem(SgBignum * dividend,SgBignum * divisor,SgBignum * quotient,SgBignum * remainder)1437 static void bignum_gdiv_rem(SgBignum *dividend, SgBignum *divisor,
1438 SgBignum *quotient, SgBignum *remainder)
1439 {
1440 ulong rsize;
1441 ulong *q = NULL, *r = NULL;
1442 ulong qs = 0, rs = 0;
1443 if (remainder) {
1444 rs = remainder->size;
1445 r = remainder->elements;
1446 }
1447 if (quotient) {
1448 q = quotient->elements;
1449 qs = quotient->size;
1450 }
1451 rsize = mp_div_rem(q, qs,
1452 r, rs,
1453 dividend->elements, dividend->size,
1454 divisor->elements, divisor->size);
1455 if (remainder) {
1456 remainder->size = rsize;
1457 remainder->sign = dividend->sign;
1458 }
1459 }
1460
bignum_gdiv_rec(SgBignum * dividend,SgBignum * divisor,SgBignum * quotient,int remainderp)1461 static SgBignum* bignum_gdiv_rec(SgBignum *dividend, SgBignum *divisor,
1462 SgBignum *quotient, int remainderp)
1463 {
1464 SgBignum *u = NULL;
1465 if (remainderp) {
1466 u = make_bignum(dividend->size + 1);
1467 }
1468 bignum_gdiv_rem(dividend, divisor, quotient, u);
1469 return u;
1470 }
1471
1472 #define bignum_gdiv(dend, dvis, quo) bignum_gdiv_rec(dend, dvis, quo, TRUE)
1473
1474
bignum_sdiv_rec(SgBignum * quotient,SgBignum * dividend,ulong divisor)1475 static ulong bignum_sdiv_rec(SgBignum *quotient,
1476 SgBignum *dividend, ulong divisor)
1477 {
1478 #ifdef USE_DLONG
1479 if (dividend->size == 1) {
1480 udlong de = dividend->elements[0];
1481 ulong q = (ulong) (de / divisor);
1482 ulong r = (ulong) (de - q * divisor);
1483 quotient->elements[0] = q;
1484 return r;
1485 } else {
1486 long n = dividend->size - 1;
1487 ulong *pu = dividend->elements;
1488 ulong *qu = quotient->elements;
1489 int shift = nlz(divisor);
1490 udlong rem = pu[n];
1491
1492 if (rem < divisor) {
1493 qu[n] = 0;
1494 } else {
1495 qu[n] = (ulong)(rem / divisor);
1496 rem = (ulong)(rem - ((udlong)qu[n] * divisor));
1497 }
1498 n--;
1499 for (; n >= 0; n--) {
1500 udlong e = (rem << WORD_BITS) | pu[n];
1501 ulong q;
1502 q = (ulong)(e/divisor);
1503 rem = (ulong)(e - ((udlong)q * divisor));
1504 qu[n] = q;
1505 }
1506 return (shift>0) ? rem % divisor: rem;
1507 }
1508 #else
1509 /* only HALF_WORD */
1510 if (divisor < HALF_WORD) {
1511 int n = dividend->size - 1;
1512 unsigned long *pu = dividend->elements;
1513 ulong q0 = 0, r0 = 0, q1, r1;
1514
1515 for (; n > 0; n--) {
1516 q1 = pu[n] / divisor + q0;
1517 r1 = ((pu[n] % divisor) << HALF_BITS) + HI(pu[n - 1]);
1518 q0 = ((r1 / divisor) << HALF_BITS);
1519 r0 = r1 % divisor;
1520 pu[n] = q1;
1521 pu[n - 1] = (r0 << HALF_BITS) + LO(pu[n - 1]);
1522 }
1523 q1 = pu[0] / divisor + q0;
1524 r1 = pu[0] % divisor;
1525 pu[0] = q1;
1526 return r1;
1527 } else {
1528 /* FIXME this doesn't work */
1529 SgBignum *bv = SG_BIGNUM(Sg_MakeBignumFromSI(divisor));
1530 SgBignum *br, *q;
1531 int i;
1532 q = make_bignum(dividend->size + 1);
1533 br = bignum_gdiv(dividend, bv, q);
1534 for (i = 0; i < q->size; i++) {
1535 dividend->elements[i] = q->elements[i];
1536 }
1537 return br->elements[0];
1538 }
1539 #endif
1540 }
1541
1542 #define bignum_sdiv(a, si) bignum_sdiv_rec(a, a, si)
1543
1544 static SgBignum *ZERO = NULL;
1545
bignum_div_rem(SgBignum * a,SgBignum * b,SgBignum ** rr)1546 static SgBignum ** bignum_div_rem(SgBignum *a, SgBignum *b, SgBignum **rr)
1547 {
1548 SgBignum *q, *r;
1549 ulong qsize = a->size - b->size + 1;
1550 if (Sg_BignumAbsCmp(a, b) < 0) {
1551 rr[0] = ZERO;
1552 rr[1] = a;
1553 return rr;
1554 }
1555 if (rr[0] && rr[0]->size >= qsize ) {
1556 q = rr[0];
1557 } else {
1558 q = make_bignum(qsize);
1559 }
1560 if (rr[1] && rr[1]->size >= a->size+1) {
1561 r = rr[1];
1562 } else {
1563 r = make_bignum(a->size+1);
1564 }
1565 /* this would help alot on sizeof(long) == 8 environment */
1566 if (b->size == 1) {
1567 ulong ur;
1568 ur = bignum_sdiv_rec(q, a, b->elements[0]);
1569 r->elements[0] = ur;
1570 r->size = 1;
1571 r->sign = a->sign;
1572 } else {
1573 bignum_gdiv_rem(a, b, q, r);
1574 }
1575 SG_BIGNUM_SET_SIGN(q, SG_BIGNUM_GET_SIGN(a) * SG_BIGNUM_GET_SIGN(b));
1576 SG_BIGNUM_SET_SIGN(r, SG_BIGNUM_GET_SIGN(a));
1577 rr[0] = bignum_normalize_rec(q, FALSE);
1578 rr[1] = bignum_normalize_rec(r, FALSE);
1579 return rr;
1580 }
1581
Sg_BignumDivRem(SgBignum * a,SgBignum * b)1582 SgObject Sg_BignumDivRem(SgBignum *a, SgBignum *b)
1583 {
1584 SgBignum *rr[2] = {NULL, };
1585 bignum_div_rem(a, b, rr);
1586 return Sg_Cons(Sg_NormalizeBignum(rr[0]), Sg_NormalizeBignum(rr[1]));
1587 }
1588
Sg_BignumModulo(SgBignum * a,SgBignum * b,int remp)1589 SgObject Sg_BignumModulo(SgBignum *a, SgBignum *b, int remp)
1590 {
1591 SgBignum *r;
1592 r = bignum_gdiv(a, b, NULL);
1593 r = Sg_NormalizeBignum(r);
1594 if (!remp
1595 && (r != SG_MAKE_INT(0))
1596 && (SG_BIGNUM_GET_SIGN(a) * SG_BIGNUM_GET_SIGN(b) < 0)) {
1597 if (SG_BIGNUMP(r)) {
1598 return Sg_BignumAdd(SG_BIGNUM(b), SG_BIGNUM(r));
1599 } else {
1600 return Sg_BignumAddSI(SG_BIGNUM(b), SG_INT_VALUE(r));
1601 }
1602 }
1603 return r;
1604 }
1605
Sg_BignumModuloSI(SgBignum * a,long b,int remp)1606 SgObject Sg_BignumModuloSI(SgBignum *a, long b, int remp)
1607 {
1608 long r = 0;
1609 long d = (b < 0) ? -b : b;
1610 long i;
1611 int sign = SG_BIGNUM_GET_SIGN(a);
1612
1613 /* ASSERT(b != 0); */
1614 for (i = SG_BIGNUM_GET_COUNT(a)-1 ; i >= 0; i--) {
1615 r = (((udlong)r<<WORD_BITS) | a->elements[i]) % d;
1616 }
1617 r = r * sign; /* got remainder */
1618
1619 if (!remp && r && sign * b < 0){
1620 return Sg_MakeIntegerFromS64((int64_t)(b + r));
1621 }
1622 return Sg_MakeInteger(r);
1623
1624 }
1625
Sg_BignumDivSI(SgBignum * a,long b,long * rem)1626 SgObject Sg_BignumDivSI(SgBignum *a, long b, long *rem)
1627 {
1628 unsigned long dd= (b < 0) ? -b : b;
1629 unsigned long rr;
1630 int d_sign = (b < 0) ? -1 : 1;
1631 SgBignum *q;
1632
1633 q = SG_BIGNUM(Sg_BignumCopy(a));
1634 rr = bignum_sdiv(q, dd);
1635
1636 if (rem) {
1637 *rem = ((SG_BIGNUM_GET_SIGN(a) < 0) ? -(signed long)rr : (signed long)rr);
1638 }
1639 SG_BIGNUM_SET_SIGN(q, SG_BIGNUM_GET_SIGN(a) * d_sign);
1640 return Sg_NormalizeBignum(q);
1641 }
1642
bn_sqrt(SgBignum * bn)1643 static SgBignum * bn_sqrt(SgBignum *bn)
1644 {
1645 long count, workpad_count, bitsize;
1646 SgBignum *s, *workpad;
1647
1648 if (SG_BIGNUM_GET_SIGN(bn) == 0) return bn;
1649 count = SG_BIGNUM_GET_COUNT(bn);
1650 ALLOC_TEMP_BIGNUM_REC(s, count);
1651 bignum_copy(s, bn);
1652 SG_BIGNUM_SET_SIGN(s, 1);
1653
1654 workpad_count = count + 1;
1655 ALLOC_TEMP_BIGNUM(workpad, workpad_count);
1656
1657 bitsize = Sg_BignumBitSize(bn);
1658 bignum_rshift(s, s, (bitsize - 1) / 2);
1659 bignum_normalize(s);
1660
1661 while (TRUE) {
1662 memset(workpad->elements, 0, sizeof(long) * workpad_count);
1663 SG_BIGNUM_SET_COUNT(workpad, workpad_count);
1664 /* only need quotient */
1665 bignum_gdiv_rec(bn, s, workpad, FALSE);
1666 if (SG_BIGNUM_GET_SIGN(workpad) == SG_BIGNUM_GET_SIGN(s)) {
1667 bignum_add_int(workpad, workpad, s);
1668 } else {
1669 bignum_sub_int(workpad, workpad, s);
1670 }
1671 bignum_rshift(workpad, workpad, 1);
1672 bignum_normalize(workpad);
1673
1674 if (Sg_BignumCmp(workpad, s) >= 0) {
1675 bignum_copy(bn, s);
1676 SG_BIGNUM_SET_SIGN(bn, 1);
1677 return bignum_normalize(bn);;
1678 }
1679 bignum_copy(s, workpad);
1680 /* count = SG_BIGNUM_GET_COUNT(workpad); */
1681 /* SG_BIGNUM_SET_COUNT(s, count); */
1682 /* memcpy(s->elements, workpad->elements, sizeof(long) * count); */
1683 }
1684 /* never reach */
1685 return bn;
1686 }
1687
Sg_BignumSqrt(SgBignum * bn)1688 SgObject Sg_BignumSqrt(SgBignum *bn)
1689 {
1690 long count;
1691 SgBignum *pad;
1692 double s;
1693
1694 count = SG_BIGNUM_GET_COUNT(bn);
1695 ALLOC_TEMP_BIGNUM_REC(pad, count);
1696 bignum_copy(pad, bn);
1697 bn_sqrt(pad);
1698 if (bn->elements[0] == pad->elements[0] * pad->elements[0]) {
1699 SgBignum *s2;
1700 ALLOC_TEMP_BIGNUM_REC(s2, SG_BIGNUM_GET_COUNT(pad)<<1);
1701 /* square_to_len(pad->elements, pad->size, s2->elements); */
1702 mp_square(s2->elements, pad->size<<1, pad->elements, pad->size);
1703 bignum_normalize(s2);
1704 /* set the same sign for temp otherwise cmp returns non 0 even the
1705 value is the same. */
1706 SG_BIGNUM_SET_SIGN(s2, SG_BIGNUM_GET_SIGN(bn));
1707 if (Sg_BignumCmp(bn, s2) == 0) {
1708 if (SG_BIGNUM_GET_SIGN(bn) > 0) return Sg_BignumToInteger(pad);
1709 return Sg_MakeComplex(SG_MAKE_INT(0), Sg_BignumToInteger(pad));
1710 }
1711 }
1712 s = Sg_BignumToDouble(bn);
1713 s = sqrt(s < 0.0 ? -s : s);
1714 if (SG_BIGNUM_GET_SIGN(bn) > 0) return Sg_MakeFlonum(s);
1715 return Sg_MakeComplex(Sg_MakeFlonum(0.0), Sg_MakeFlonum(s));
1716 }
1717
Sg_BignumSqrtApprox(SgBignum * bn)1718 SgObject Sg_BignumSqrtApprox(SgBignum *bn)
1719 {
1720 long count;
1721 SgBignum *workpad;
1722
1723 count = SG_BIGNUM_GET_COUNT(bn);
1724 ALLOC_TEMP_BIGNUM_REC(workpad, count);
1725 bignum_copy(workpad, bn);
1726 bn_sqrt(workpad);
1727 if (SG_BIGNUM_GET_SIGN(bn) > 0) return Sg_BignumToInteger(workpad);
1728 else return Sg_MakeComplex(SG_MAKE_INT(0), Sg_BignumToInteger(workpad));
1729 }
1730
Sg_MakeBignumWithSize(long size,unsigned long init)1731 SgObject Sg_MakeBignumWithSize(long size, unsigned long init)
1732 {
1733 SgBignum *b = make_bignum(size);
1734 b->elements[0] = init;
1735 return b;
1736 }
1737
Sg_BignumAccMultAddUI(SgBignum * acc,unsigned long coef,unsigned long c)1738 SgObject Sg_BignumAccMultAddUI(SgBignum *acc, unsigned long coef,
1739 unsigned long c)
1740 {
1741 SgBignum *r;
1742 unsigned long rsize = SG_BIGNUM_GET_COUNT(acc) + 1, i;
1743 unsigned long carry;
1744 ALLOC_TEMP_BIGNUM(r, rsize);
1745 r->elements[0] = c;
1746 carry = mp_mul_add(r->elements, acc->elements,
1747 SG_BIGNUM_GET_COUNT(acc), coef);
1748 r->elements[SG_BIGNUM_GET_COUNT(acc)] = carry;
1749
1750 if (r->elements[rsize - 1] == 0) {
1751 for (i = 0; i < SG_BIGNUM_GET_COUNT(acc); i++) {
1752 acc->elements[i] = r->elements[i];
1753 }
1754 return acc;
1755 } else {
1756 SgBignum *rr;
1757 rr = make_bignum(rsize + 3);
1758 SG_BIGNUM_SET_SIGN(rr, SG_BIGNUM_GET_SIGN(acc));
1759 for (i = 0; i < rsize; i++) {
1760 rr->elements[i] = r->elements[i];
1761 }
1762 return rr;
1763 }
1764 }
1765
calc_string_size(SgBignum * q,int radix)1766 static long calc_string_size(SgBignum *q, int radix)
1767 {
1768 long count = 0;
1769 long size = SG_BIGNUM_GET_COUNT(q);
1770 for (; size > 0;) {
1771 bignum_sdiv(q, radix);
1772 count++;
1773 for (; q->elements[size - 1] == 0 && size > 0; size--);
1774 }
1775 return count;
1776 }
1777
radix2_string_helper(SgObject r,ulong e,long size,long off)1778 static long radix2_string_helper(SgObject r, ulong e, long size, long off)
1779 {
1780 long i;
1781 for (i = size-1; i >= 0; i--) {
1782 int b = (e>>i) & 1;
1783 SG_STRING_VALUE_AT(r, off++) = (b) ? '1' : '0';
1784 }
1785 return off;
1786 }
radix2_string(SgBignum * b)1787 static SgObject radix2_string(SgBignum *b)
1788 {
1789 /* if the radix is 2, we can convert bit by bit */
1790 SgObject r;
1791 long count, n, off = 0, i;
1792 /* we only need to care the first element's bit size */
1793 count = n = WORD_BITS - nlz(b->elements[b->size-1]);
1794
1795 if (b->sign < 0) count++;
1796 /* compute the rest of elements bit */
1797 count += (b->size-1)*WORD_BITS;
1798
1799 /* allocate string */
1800 r = Sg_ReserveString(count, 0);
1801 if (b->sign < 0) SG_STRING_VALUE_AT(r, off++) = '-';
1802
1803 off = radix2_string_helper(r, b->elements[b->size-1], n, off);
1804 for (i = b->size-2; i >= 0; i--) {
1805 off = radix2_string_helper(r, b->elements[i], WORD_BITS, off);
1806 }
1807 return r;
1808 }
1809
1810 /* radix 8 is a bit tricky
1811 we need to split each elements 3 bit each however the word
1812 length is either 32 or 64, none of them are multiple of 3.
1813 so first we need to calculate offset of the first word then
1814 proceed.
1815 */
radix8_string(SgBignum * b)1816 static SgObject radix8_string(SgBignum *b)
1817 {
1818 static const char tab[] = "012345678";
1819 long bits = Sg_BignumBitSize(b);
1820 int ob = bits % 3; /* offset bit */
1821 long off = 0, i, fb, times, c;
1822 ulong e;
1823 SgObject r = Sg_ReserveString((bits/3)+((ob)?1:0)+((b->sign<0)?1:0), 0);
1824
1825 if (b->sign < 0) SG_STRING_VALUE_AT(r, off++) = '-';
1826 /*
1827 e.g
1828 word = 4, bits
1829 first word (fb = 7)
1830 1000100 0..0 0..0 (2 words + fb = 71)
1831 ob = 2
1832 we need to use offset -1 from first word to make total bit
1833 multiple of 3.
1834
1835 so not the first word is
1836 01000100 (fb 8)
1837 */
1838 /* handle first element */
1839 ob = (ob)? (3-ob) : 0;
1840 fb = WORD_BITS - nlz(b->elements[b->size-1]) + ob;
1841 times = fb / 3;
1842 c = fb % 3; /* for next iteration */
1843 e = b->elements[b->size-1];
1844
1845 for (i = 0; i < times; i++) {
1846 int t = (e >> ((fb-3)-(i*3))) & 7;
1847 SG_STRING_VALUE_AT(r, off++) = tab[t];
1848 }
1849 e = (e <<(3-c))&7; /* gets carry bits */
1850
1851 /* do the rest */
1852 for (i = b->size-2; i >= 0; i--) {
1853 ulong re = b->elements[i];
1854 int j;
1855
1856 times = (WORD_BITS+c)/3;
1857 for (j = 0; j < times; j++) {
1858 int t = ((re >> ((WORD_BITS-(3-c))-(j*3)))|e)&7;
1859 SG_STRING_VALUE_AT(r, off++) = tab[t];
1860 if (e) e = 0; /* clear carry */
1861 }
1862 c = (WORD_BITS+c)%3; /* next carry */
1863 /* only if the c is not 0, otherwise we'll get garbage */
1864 if (c) e = (re <<(3-c))&7; /* gets carry bits */
1865 }
1866 return r;
1867 }
1868
radix16_string(SgBignum * b,int use_upper)1869 static SgObject radix16_string(SgBignum *b, int use_upper)
1870 {
1871 /* if the radix is 16 then we can simply dump the elements */
1872 char buf[(SIZEOF_LONG<<1)+1];
1873 SgObject r;
1874 #if SIZEOF_LONG == 8
1875 char *fmt = (use_upper)? "%016lX": "%016lx";
1876 #else
1877 char *fmt = (use_upper)? "%08lX": "%08lx";
1878 #endif
1879 char *first_fmt = (use_upper)? "%lX": "%lx";
1880 long count, n, off = 0, i, j;
1881 /* this is the very first time I'm using the return value of
1882 printf related procedure...
1883 */
1884 count = n = snprintf(buf, sizeof(buf), first_fmt, b->elements[b->size-1]);
1885 if (b->sign < 0) count++;
1886 /* calculate the rest of words */
1887 count += (b->size-1) * (SIZEOF_LONG<<1);
1888
1889 r = Sg_ReserveString(count, 0);
1890 /* set the first word */
1891 if (b->sign < 0) SG_STRING_VALUE_AT(r, off++) = '-';
1892 for (i = 0; i < n; i++) {
1893 SG_STRING_VALUE_AT(r, off++) = buf[i];
1894 }
1895
1896 for (i = b->size-2; i >= 0; i--) {
1897 snprintf(buf, sizeof(buf), fmt, b->elements[i]);
1898 for (j = 0; j < sizeof(buf)-1; j++) {
1899 SG_STRING_VALUE_AT(r, off++) = buf[j];
1900 }
1901 }
1902 return r;
1903 }
1904
1905 #define SCHONEHAGE_BASE_CONVERSION_THRESHOLD 20
1906
small_bignum_to_string(SgBignum * b,int radix,int use_upper)1907 static SgObject small_bignum_to_string(SgBignum *b, int radix, int use_upper)
1908 {
1909 static const char ltab[] = "0123456789abcdefghijklmnopqrstuvwxyz";
1910 static const char utab[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
1911 const char *tab = use_upper ? utab : ltab;
1912 /* SgObject h = SG_NIL, t = SG_NIL; */
1913 SgObject rs;
1914 SgBignum *q;
1915 long rem, size, count;
1916 long i;
1917 q = SG_BIGNUM(Sg_BignumCopy(b));
1918 size = SG_BIGNUM_GET_COUNT(q);
1919
1920 count = calc_string_size(q, radix);
1921 if (SG_BIGNUM_GET_SIGN(q) < 0) count++;
1922 for (i = 0; i < size; i++) q->elements[i] = b->elements[i];
1923 rs = Sg_ReserveString(count, 0);
1924
1925 for (i = count-1; size > 0; i--) {
1926 rem = bignum_sdiv(q, radix);
1927 /* SG_APPEND1(h, t, SG_MAKE_CHAR(tab[rem])); */
1928 SG_STRING_VALUE_AT(rs, i) = tab[rem];
1929 for (; q->elements[size - 1] == 0 && size > 0; size--);
1930 }
1931 if (SG_BIGNUM_GET_SIGN(q) < 0) {
1932 SG_STRING_VALUE_AT(rs, 0) = '-';
1933 /* SG_APPEND1(h, t, SG_MAKE_CHAR('-')); */
1934 }
1935 /* return Sg_ListToString(Sg_ReverseX(h), 0, -1); */
1936 return rs;
1937 }
1938
1939 /* FIXME this is also in number.c */
roundeven(double v)1940 static inline double roundeven(double v)
1941 {
1942 double r;
1943 double frac = modf(v, &r);
1944 if (v > 0.0) {
1945 if (frac > 0.5) r += 1.0;
1946 else if (frac == 0.5) {
1947 if (fmod(r, 2.0) != 0.0) r += 1.0;
1948 }
1949 } else {
1950 if (frac < -0.5) r -= 1.0;
1951 else if (frac == -0.5) {
1952 if (fmod(r, 2.0) != 0.0) r -= 1.0;
1953 }
1954 }
1955 return r;
1956 }
1957 /* this is used here */
1958 static SgBignum * bignum_expt(SgBignum *b, long exponent);
1959
1960 #define MAX_RADIX 36
1961 static SgBignum *RADIXES[MAX_RADIX+1] = {NULL,};
1962 /* returns radix^2^exponent
1963 the exponent is usually not so big (unless the number is extremely huge,
1964 but such numbers can't be create on Sagittarius due to the memory
1965 limitation). So caluculating it each time is more expensive than
1966 allocating small space for cache.
1967 */
radix_conversion(int radix,int exponent,SgBignum *** two_exps)1968 static SgBignum* radix_conversion(int radix, int exponent, SgBignum ***two_exps)
1969 {
1970 SgObject c;
1971 SgBignum *r;
1972
1973 if (!*two_exps) {
1974 /* the first time, then create the array with exponent. don't worry,
1975 it's not so big */
1976 *two_exps = SG_NEW_ARRAY(SgBignum *, exponent+1);
1977 }
1978 r = (*two_exps)[exponent];
1979 if (r) return r;
1980
1981 c = Sg_Expt(SG_MAKE_INT(2), SG_MAKE_INT(exponent));
1982 if (!SG_INTP(c)) Sg_Error(UC("big num is too big to show"));
1983 /* keep it as bignum */
1984 r = bignum_normalize_rec(bignum_expt(RADIXES[radix], SG_INT_VALUE(c)),
1985 FALSE);
1986 (*two_exps)[exponent] = r;
1987 return r;
1988 }
1989
1990 /*
1991 it does improve the performance as long as the number is small
1992 however once it consumed all stack then the GC would get affected
1993 by the large grown stack area. I believe, this impacted performance.
1994 c.f)
1995 with stack
1996 ;; (string->number (number->string (expt 2 (vector-ref v 0))))
1997 ;; -6.915366 real 2.154000 user 0.218000 sys
1998 11.53s user 0.31s system 103% cpu 11.410 total
1999
2000 without stack
2001 ;; (string->number (number->string (expt 2 (vector-ref v 0))))
2002 ;; 0.337534 real 9.266000 user 0.437000 sys
2003 9.67s user 0.50s system 105% cpu 9.657 total
2004
2005 2 seconds are not negligible difference.
2006 */
2007 #undef USE_STACK_FOR_SCHONEHAGE
2008
2009 /* cache constant log values */
2010 /* static double DBL_LOG2 = 0.0; */
2011 #define DBL_LOG2 0.6931471805599453
2012 static double DBL_RADIXES[MAX_RADIX+1];
2013 /*
2014 Using Schonehage base conversion
2015
2016 This is much faster than before (it's now named small_bignum_to_string),
2017 however the implementation uses considerable memory if the bignum is
2018 huge. In that case it's more fight against GC than small performance
2019 tuning such as caching.
2020
2021 For future note.
2022 There's couple of things I can try to squeeze performance.
2023 - rope: string port isn't that bad but calling overhead would be redueced
2024 - non recursive conversion: then we can use stack instead of heap
2025 The second one would be a challenge.
2026 */
schonehage_to_string(SgBignum * b,SgObject out,int radix,int count,int use_upper,SgBignum *** two_exps)2027 static void schonehage_to_string(SgBignum *b, SgObject out, int radix,
2028 int count, int use_upper,
2029 /* r^2^e cache */
2030 SgBignum ***two_exps
2031 /*, int use_stack */)
2032 {
2033 long bits;
2034 int n, e;
2035 SgBignum *v;
2036 SgBignum *result[2] = {NULL, };
2037 double lr, l2;
2038 #ifdef USE_STACK_FOR_SCHONEHAGE
2039 intptr_t stack_size;
2040 ulong qsize, rsize;
2041 #endif
2042 if (b->size < SCHONEHAGE_BASE_CONVERSION_THRESHOLD) {
2043 SgObject r = small_bignum_to_string(b, radix, use_upper);
2044 long i;
2045 if (SG_STRING_SIZE(r) < count && Sg_PortPosition(out) > 0) {
2046 for (i = SG_STRING_SIZE(r); i < count; i++) Sg_PutcUnsafe(out, '0');
2047 }
2048 Sg_PutsUnsafe(out, r);
2049 return;
2050 }
2051 bits = Sg_BitSize(b);
2052 l2 = DBL_LOG2;
2053 lr = DBL_RADIXES[radix];
2054
2055 n = (int)roundeven(log(bits * l2 / lr) / lr - 1.0);
2056 v = radix_conversion(radix, n, two_exps);
2057
2058 #ifdef USE_STACK_FOR_SCHONEHAGE
2059 if (use_stack) {
2060 stack_size = Sg_AvailableStackSize((volatile uintptr_t)&bits);
2061 qsize = b->size - v->size + 1;
2062 rsize = b->size+1;
2063 if (stack_size > 0 && stack_size > BIGNUM_SIZE(qsize)) {
2064 stack_size -= BIGNUM_SIZE(qsize);
2065 ALLOC_TEMP_BIGNUM(result[0], qsize);
2066 } else {
2067 use_stack = FALSE;
2068 }
2069 if (stack_size > 0 && stack_size > BIGNUM_SIZE(rsize)) {
2070 ALLOC_TEMP_BIGNUM(result[1], rsize);
2071 } else {
2072 use_stack = FALSE;
2073 }
2074 }
2075 #endif
2076
2077 bignum_div_rem(b, v, result);
2078 e = 1<<n;
2079 schonehage_to_string(result[0], out, radix, count-e, use_upper, two_exps);
2080 schonehage_to_string(result[1], out, radix, e, use_upper, two_exps);
2081 }
2082
Sg_BignumToString(SgBignum * b,int radix,int use_upper)2083 SgObject Sg_BignumToString(SgBignum *b, int radix, int use_upper)
2084 {
2085 if (radix < 2 || radix > 36) {
2086 Sg_Error(UC("radix out of range: %d"), radix);
2087 }
2088 /* special case 0 */
2089 if (b->sign == 0 || b->size == 0) return SG_MAKE_STRING("0");
2090
2091 /* handle easily converted case */
2092 if (radix == 2) return radix2_string(b);
2093 if (radix == 8) return radix8_string(b);
2094 if (radix == 16) return radix16_string(b, use_upper);
2095
2096 /* The Art of Computer Programming Vol2 4.4 (Q 14) answer*/
2097 if (b->size < SCHONEHAGE_BASE_CONVERSION_THRESHOLD) {
2098 return small_bignum_to_string(b, radix, use_upper);
2099 } else {
2100 SgPort *out;
2101 SgStringPort sp;
2102 SgBignum **two_exponents = NULL;
2103 out = Sg_InitStringOutputPort(&sp, 1024);
2104 if (b->sign < 0) {
2105 b = SG_BIGNUM(Sg_Negate(SG_OBJ(b)));
2106 Sg_PutcUnsafe(out, '-');
2107 }
2108 schonehage_to_string(b, out, radix, 0, use_upper, &two_exponents);
2109 return Sg_GetStringFromStringPort(&sp);
2110 }
2111 }
2112
2113 /* we do this destructively */
bignum_difference(SgBignum * a,SgBignum * b)2114 static int bignum_difference(SgBignum *a, SgBignum *b)
2115 {
2116 int sign = Sg_BignumCmp(a, b);
2117
2118 if (sign == 0) return 0;
2119
2120 if (sign < 0) {
2121 SgBignum *tmp = a;
2122 a = b;
2123 b = tmp;
2124 }
2125 ASSERT(SG_BIGNUM_GET_COUNT(a) >= SG_BIGNUM_GET_COUNT(b));
2126 bignum_sub_int(a, a, b);
2127 bignum_normalize(a);
2128 return sign;
2129 }
2130
2131
gcd_fixfix(unsigned long x,unsigned long y)2132 static inline unsigned long gcd_fixfix(unsigned long x, unsigned long y)
2133 {
2134 while (y > 0) {
2135 unsigned long r = x % y;
2136 x = y;
2137 y = r;
2138 }
2139 return x;
2140 }
2141
2142 /* avoid memory allocation as much as possible */
binary_gcd(SgBignum * bx,SgBignum * by)2143 static SgObject binary_gcd(SgBignum *bx, SgBignum *by)
2144 {
2145 /* Algorithm B from Knuth section 4.5.2 */
2146 SgBignum *u = SG_BIGNUM(Sg_Abs(Sg_BignumCopy(bx)));
2147 SgBignum *v = SG_BIGNUM(Sg_Abs(Sg_BignumCopy(by)));
2148 SgObject ret;
2149 /* step B1 */
2150 long s1 = Sg_BignumFirstBitSet(u);
2151 long s2 = Sg_BignumFirstBitSet(v);
2152 long k = min(s1, s2);
2153 /* these are for step B2 */
2154 int uOdd, tsign;
2155 long lb;
2156 SgBignum *t;
2157 /* Sg_Printf(Sg_StandardErrorPort(), UC(";; x=%A, y=%A, k=%d\n"), u, v, k); */
2158 if (k != 0) {
2159 bignum_rshift(u, u, k);
2160 bignum_rshift(v, v, k);
2161 bignum_normalize(u);
2162 bignum_normalize(v);
2163 }
2164 /* Step B2 */
2165 uOdd = (k==s1);
2166 t = uOdd ? v: u;
2167 tsign = uOdd ? -1 : 1;
2168
2169 while ((lb = Sg_BignumFirstBitSet(t)) >= 0) {
2170 /* Step B3 and B4 */
2171 /* Sg_Printf(Sg_StandardErrorPort(), UC("(print (= (>> %A %d) "), t, lb); */
2172 bignum_rshift(t, t, lb);
2173 bignum_normalize(t);
2174 /* Sg_Printf(Sg_StandardErrorPort(), UC("%A))\n"), t); */
2175 /* Step B5 */
2176 if (tsign > 0) u = t;
2177 else v = t;
2178 /* special case one word numbers */
2179 if (SG_BIGNUM_GET_COUNT(u) < 2 && SG_BIGNUM_GET_COUNT(v) < 2) {
2180 unsigned long x = u->elements[0], y = v->elements[0];
2181 x = (x >= y) ? gcd_fixfix(x, y) : gcd_fixfix(y, x);
2182 ret = Sg_MakeIntegerU(x);
2183 if (k > 0) {
2184 /* Sg_Printf(Sg_StandardErrorPort(), UC(";; r=%A, k=%d\n"), ret, k); */
2185 ret = Sg_Ash(ret, k);
2186 }
2187 goto end;
2188 }
2189 /* Step B6 */
2190 /*
2191 Sg_Printf(Sg_StandardErrorPort(), UC("(print (= (abs (- %A %A))"), u, v);
2192 */
2193 if ((tsign = bignum_difference(u, v)) == 0) break;
2194 t = (tsign >= 0) ? u : v;
2195 /* Sg_Printf(Sg_StandardErrorPort(), UC("%A))\n"), t); */
2196 }
2197 if (k > 0) {
2198 /* Sg_Printf(Sg_StandardErrorPort(), UC(";; u=%A, k=%d\n"), u, k); */
2199 ret = Sg_BignumShiftLeft(u, k);
2200 } else {
2201 ret = Sg_NormalizeBignum(u);
2202 }
2203 end:
2204 /* Sg_Printf(Sg_StandardErrorPort(), UC(";; return=%A\n"), ret); */
2205 return ret;
2206 }
2207
Sg_BignumGcd(SgBignum * bx,SgBignum * by)2208 SgObject Sg_BignumGcd(SgBignum *bx, SgBignum *by)
2209 {
2210 /* return binary_gcd(bx, by); */
2211 /* hybrid was faster. */
2212 while (SG_BIGNUM_GET_COUNT(by) != 0) {
2213 SgBignum *r;
2214 SgObject t;
2215 if (labs((long)SG_BIGNUM_GET_COUNT(bx)-(long)SG_BIGNUM_GET_COUNT(by)) < 2)
2216 return binary_gcd(bx, by);
2217 r = bignum_gdiv(bx, by, NULL);
2218 t = Sg_NormalizeBignum(r);
2219 if (SG_INTP(t)) return Sg_Gcd(by, t);
2220 bx = by;
2221 by = r;
2222 }
2223 return bx;
2224 }
2225
inverse_mod_long(ulong val)2226 static long inverse_mod_long(ulong val)
2227 {
2228 long t = val;
2229 t *= 2 - val*t;
2230 t *= 2 - val*t;
2231 t *= 2 - val*t;
2232 t *= 2 - val*t;
2233 #if SIZEOF_LONG == 8
2234 t *= 2 - val*t;
2235 t *= 2 - val*t;
2236 t *= 2 - val*t;
2237 t *= 2 - val*t;
2238 #endif
2239 return t;
2240 }
2241
2242 #define BIGNUM_ZEROP(bn) (SG_BIGNUM_GET_SIGN(bn) == 0)
2243 #define BIGNUM_ONEP(bn) (SG_BIGNUM_GET_COUNT(bn) == 1 && (bn)->elements[0] == 1)
2244
bignum_mod(SgBignum * a,SgBignum * b,SgBignum * q)2245 static SgBignum * bignum_mod(SgBignum *a, SgBignum *b, SgBignum *q)
2246 {
2247 SgBignum *r = bignum_gdiv(a, b, q);
2248 bignum_normalize(r);
2249 if (BIGNUM_ZEROP(r)
2250 && (SG_BIGNUM_GET_SIGN(a) * SG_BIGNUM_GET_SIGN(b) < 0)) {
2251 return bignum_normalize(bignum_add(b, r));
2252 }
2253 return r;
2254 }
2255
2256 /*
2257 Using extended Euclidean algorithm.
2258
2259 The complexity of the compuation is O(log(m)^2).
2260
2261 For now this is actually good enough since both x and m are not
2262 really huge. If we got a huge number (e.g. 2048 bits) we can think
2263 about other algorithm.
2264 Alternatives:
2265 - The Montgomery inverse and its applications, Burton S. Kaliski:
2266
2267 - Constant Time Modular Inversion, Joppe W Bos:
2268 http://www.joppebos.com/files/CTInversion.pdf
2269 */
bignum_mod_inverse(SgBignum * x,SgBignum * m)2270 static SgBignum * bignum_mod_inverse(SgBignum *x, SgBignum *m)
2271 {
2272 SgBignum *u1, *u3, *v1, *v3, *q;
2273 int sign = 1;
2274 long size;
2275 if (SG_BIGNUM_GET_COUNT(x) > SG_BIGNUM_GET_COUNT(m)) {
2276 size = SG_BIGNUM_GET_COUNT(x) + 1;
2277 } else {
2278 size = SG_BIGNUM_GET_COUNT(m) + 1;
2279 }
2280
2281 if (SG_BIGNUM_GET_SIGN(m) < 0) {
2282 Sg_Error(UC("modulus not positive %S"), m);
2283 }
2284 u1 = Sg_MakeBignumFromSI(1);
2285 u3 = x;
2286 v1 = Sg_MakeBignumFromSI(0);
2287 v3 = m;
2288 ALLOC_TEMP_BIGNUM(q, size);
2289 while (!BIGNUM_ZEROP(v3)) {
2290 SgBignum *t3, *w, *t1;
2291 t3 = bignum_mod(u3, v3, q);
2292 bignum_normalize(q);
2293 w = bignum_normalize(bignum_mul(q, v1));
2294 t1 = bignum_normalize(bignum_add(u1, w));
2295 u1 = v1; v1 = t1; u3 = v3; v3 = t3;
2296 sign = -sign;
2297 /* reset buffer */
2298 /* clear the content as well since bignum_gdiv may not clear the
2299 array. (got the bug due to the non cleared stack area...)
2300
2301 NB: element size number is sufficient since it's initialised to 0.
2302 */
2303 memset(q->elements, 0, SG_BIGNUM_GET_COUNT(q) * sizeof(long));
2304 SG_BIGNUM_SET_COUNT(q, size);
2305 SG_BIGNUM_SET_SIGN(q, 1);
2306 }
2307 /* Sg_Printf(Sg_StandardErrorPort(), UC("x : %A\n"), x); */
2308 /* Sg_Printf(Sg_StandardErrorPort(), UC("m : %A\n"), m); */
2309 /* Sg_Printf(Sg_StandardErrorPort(), UC("u1: %A\n"), u1); */
2310 /* Sg_Printf(Sg_StandardErrorPort(), UC("v1: %A\n"), v1); */
2311 if (sign < 0) {
2312 return bignum_normalize(bignum_sub(m, u1));
2313 } else {
2314 return bignum_normalize(u1);
2315 }
2316 }
2317
2318 #define EXPMOD_MAX_WINDOWS 7
2319 static ulong exp_mod_threadh_table[EXPMOD_MAX_WINDOWS] = {
2320 7, 25, 81, 241, 673, 1793, (ulong)-1L
2321 };
2322
odd_mod_expt_rec(SgBignum * x,SgBignum * exp,SgBignum * mod,ulong inv,ulong * a,ulong * b)2323 static ulong * odd_mod_expt_rec(SgBignum *x, SgBignum *exp, SgBignum *mod,
2324 ulong inv, ulong *a, ulong *b)
2325 {
2326 SgBignum *tr, *tb;
2327 ulong *table[1 << EXPMOD_MAX_WINDOWS], *prod, *t;
2328 long i, j, tblmask, wbits, ebits;
2329 int isone = TRUE;
2330 long modlen = SG_BIGNUM_GET_COUNT(mod);
2331 long elen = SG_BIGNUM_GET_COUNT(exp);
2332 unsigned int buf = 0;
2333 long multpos;
2334 ulong *mult, *e = exp->elements + elen -1;
2335 ulong bitpos;
2336 long modlen2 = modlen<<1;
2337
2338 wbits = 0;
2339 ebits = Sg_BignumBitSize(exp);
2340 if ((ebits != 17) ||
2341 (SG_BIGNUM_GET_COUNT(exp) != 1 && exp->elements[0] != 65537)) {
2342 while ((unsigned long)ebits > exp_mod_threadh_table[wbits]) {
2343 wbits++;
2344 }
2345 }
2346
2347 tblmask = 1u << wbits;
2348 /* convert x to montgomery form */
2349 ALLOC_TEMP_BIGNUM(tb, SG_BIGNUM_GET_COUNT(x) + modlen);
2350 /* BIGNUM_CLEAR_LEFT(tb, modlen); */
2351 for (i = 0; i < SG_BIGNUM_GET_COUNT(x); i++) {
2352 tb->elements[i+modlen] = x->elements[i];
2353 }
2354 tr = bignum_normalize(bignum_gdiv(tb, mod, NULL));
2355 if (SG_BIGNUM_GET_COUNT(tr) < modlen) {
2356 /* table[0] won't be initialised, we need to do it here */
2357 ALLOC_TEMP_BUFFER(table[0], ulong, modlen);
2358 for (i = 0; i < (int)SG_BIGNUM_GET_COUNT(tr); i++) {
2359 table[0][i] = tr->elements[i];
2360 }
2361 } else {
2362 /* when tr has the same length as modlen, we can reuse it. */
2363 table[0] = tr->elements;
2364 }
2365 /* set b to the square of the base */
2366 /* blen = modlen<<1 */
2367 /* square_to_len(table[0], modlen, b); */
2368 mp_square(b, modlen2, table[0], modlen);
2369 mp_mont_reduce(b, modlen2, mod->elements, modlen, inv);
2370 /* Use high hald of b to initialise the table */
2371 t = b+modlen;
2372
2373 /* Fill int the table with odd powers of the base */
2374 for (i = 1; i < tblmask; i++) {
2375 ALLOC_TEMP_BUFFER(prod, ulong, modlen);
2376 /* alen = modlen<<1 */
2377 mp_mul(a, modlen2, t, modlen, table[i-1], modlen);
2378 mp_mont_reduce(a, modlen2, mod->elements, modlen, inv);
2379
2380 for (j = 0; j < modlen; j++) {
2381 prod[j] = a[j+modlen];
2382 }
2383 table[i] = prod;
2384 }
2385 /* hmm, let me make a new scope... */
2386 bitpos = 1UL << ((ebits-1) & (WORD_BITS-1));
2387
2388 for (i = 0; i <= wbits; i++) {
2389 buf = (buf << 1) | ((*e & bitpos) != 0);
2390 bitpos >>= 1;
2391 if (!bitpos) {
2392 e--;
2393 bitpos = 1UL << (WORD_BITS-1);
2394 elen--;
2395 }
2396 }
2397 /* The first iteration, which is hoisted out of the main loop */
2398 ebits--;
2399 multpos = ebits - wbits;
2400 while ((buf & 1) == 0) {
2401 buf >>= 1;
2402 multpos++;
2403 }
2404 mult = table[buf >> 1];
2405 buf = 0;
2406 if (multpos == ebits) {
2407 isone = FALSE;
2408 }
2409 /* the main loop */
2410 for (;;) {
2411 ebits--;
2412 buf <<= 1;
2413 if (elen) {
2414 buf |= ((*e & bitpos) != 0);
2415 bitpos >>= 1;
2416 if (!bitpos) {
2417 e--;
2418 bitpos = 1UL << (WORD_BITS-1);
2419 elen--;
2420 }
2421 }
2422 /* examine the window for pending multiplies */
2423 if (buf & tblmask) {
2424 multpos = ebits - wbits;
2425 while ((buf & 1) == 0) {
2426 buf >>= 1;
2427 multpos++;
2428 }
2429 mult = table[buf >> 1];
2430 buf = 0;
2431 }
2432 /* perform multiply */
2433 if (ebits == multpos) {
2434 /* t will be reused anyway, so make it shared */
2435 t = b+modlen;
2436 if (isone) {
2437 for (i = 0; i < modlen; i++) {
2438 t[i] = mult[i];
2439 }
2440 isone = FALSE;
2441 } else {
2442 /* alen = modlen<<1 */
2443 mp_mul(a, modlen2, t, modlen, mult, modlen);
2444 mp_mont_reduce(a, modlen2, mod->elements, modlen, inv);
2445 t = a; a = b; b = t;
2446 }
2447 }
2448 /* check if done */
2449 if (!ebits) return b;
2450 /* square the input */
2451 if (!isone) {
2452 t = b+modlen;
2453 /* square_to_len(t, modlen, a); */
2454 mp_square(a, modlen2, t, modlen);
2455 mp_mont_reduce(a, modlen2, mod->elements, modlen, inv);
2456 t = a; a = b; b = t;
2457 }
2458 }
2459 }
2460
odd_mod_expt(SgBignum * x,SgBignum * exp,SgBignum * mod)2461 static SgObject odd_mod_expt(SgBignum *x, SgBignum *exp, SgBignum *mod)
2462 {
2463 long modlen, i, buflen;
2464 SgBignum *r;
2465 ulong *a, *b, *t;
2466 ulong inv;
2467 if (BIGNUM_ONEP(exp)) return x;
2468 if (BIGNUM_ZEROP(x)) return x;
2469
2470 modlen = SG_BIGNUM_GET_COUNT(mod);
2471 inv = -inverse_mod_long(mod->elements[0]);
2472
2473 /* pre allocate */
2474 r = make_bignum(modlen);
2475
2476 buflen = modlen * 2;
2477 /* allocate buffers */
2478 ALLOC_TEMP_BUFFER(a, ulong, buflen);
2479 ALLOC_TEMP_BUFFER(b, ulong, buflen);
2480
2481 /* do main loop */
2482 b = odd_mod_expt_rec(x, exp, mod, inv, a, b);
2483 /* convert result out of montgomery form and return */
2484 t = b+modlen;
2485 for (i = 0; i < modlen; i++) {
2486 b[i] = t[i];
2487 t[i] = 0;
2488 }
2489 mp_mont_reduce(b, buflen, mod->elements, modlen, inv);
2490
2491 for (i = 0; i < modlen; i++) {
2492 r->elements[i] = t[i];
2493 }
2494 return bignum_normalize(r);
2495 }
2496
2497 /* mod(2^p */
bignum_mod2(SgBignum * x,long p)2498 static SgBignum * bignum_mod2(SgBignum *x, long p)
2499 {
2500 long numInts, excessBits, i;
2501 SgBignum *r;
2502 if (Sg_BignumBitSize(x) <= p) return x;
2503 numInts = (p + (WORD_BITS-1)) >> SHIFT_MAGIC;
2504 r = make_bignum(numInts);
2505 for (i = 0; i < numInts; i++) {
2506 r->elements[i] = x->elements[i];
2507 }
2508 excessBits = (numInts << SHIFT_MAGIC) - p;
2509 r->elements[numInts - 1] &= (1UL << (WORD_BITS - excessBits)) - 1;
2510 return bignum_normalize(r);
2511 }
2512
bignum_mod_expt2(SgBignum * x,SgBignum * e,long p)2513 static SgBignum * bignum_mod_expt2(SgBignum *x, SgBignum *e, long p)
2514 {
2515 SgBignum *base = bignum_mod2(x, p), *result = Sg_MakeBignumFromSI(1);
2516 long exp_offset = 0, limit = Sg_BignumBitSize(e);
2517 if (x->elements[0] & 1) {
2518 limit = ((p-1) < limit) ? (p-1) : limit;
2519 }
2520 while (exp_offset < limit) {
2521 if (bignum_test_bit(e, exp_offset)) {
2522 result = bignum_mul(result, base);
2523 bignum_normalize(result);
2524 result = bignum_mod2(result, p);
2525 }
2526 exp_offset++;
2527 if (exp_offset < limit) {
2528 base = bignum_mul(base, base);
2529 bignum_normalize(base);
2530 base = bignum_mod2(base, p);
2531 }
2532 }
2533 return result;
2534 }
2535
bignum_mod_expt(SgBignum * bx,SgBignum * be,SgBignum * bm)2536 static SgObject bignum_mod_expt(SgBignum *bx, SgBignum *be, SgBignum *bm)
2537 {
2538 int invertp = SG_BIGNUM_GET_SIGN(be) < 0;
2539 SgBignum *base, *result;
2540 if (invertp) {
2541 /* keep it bignum */
2542 SgBignum *b = Sg_BignumCopy(be);
2543 SG_BIGNUM_SET_SIGN(b, -SG_BIGNUM_GET_SIGN(be));
2544 be = b;
2545 }
2546 base = (SG_BIGNUM_GET_SIGN(bx) < 0 || Sg_BignumCmp(bx, bm) >= 0)
2547 ? bignum_mod(bx, bm, NULL) : bx;
2548
2549 if (bm->elements[0] & 1) { /* odd modulus */
2550 result = odd_mod_expt(base, be, bm);
2551 } else {
2552 /* even modulus. Tear it into an "odd part" (m1) and power of two (m2),
2553 exponentiate mod m1, manually exponentiate mod m2, and use Chinese
2554 Remainder Theorem to combine results.
2555 */
2556 long p = Sg_BignumFirstBitSet(bm);
2557 long lsize = 1 + (p + WORD_BITS - 1) / WORD_BITS;
2558 long rsize = SG_BIGNUM_GET_COUNT(bm) + (-p)/WORD_BITS;
2559 SgBignum *m1, *m2, *base2, *a1, *a2, *y1, *y2;
2560 SgBignum *ONE;
2561 ALLOC_TEMP_BIGNUM(ONE, 1);
2562 ONE->elements[0] = 1;
2563
2564 if (rsize < 1) {
2565 if (SG_BIGNUM_GET_SIGN(bm) < 0) {
2566 m1 = Sg_MakeBignumFromSI(-1);
2567 } else {
2568 m1 = Sg_MakeBignumFromSI(0);
2569 }
2570 } else {
2571 /* m must be positive */
2572 /* m1 = make_bignum(rsize); */
2573 ALLOC_TEMP_BIGNUM(m1, rsize);
2574 bignum_rshift(m1, bm, p);
2575 bignum_normalize(m1);
2576 }
2577 /* m2 = make_bignum(lsize); */
2578 ALLOC_TEMP_BIGNUM(m2, lsize);
2579 bignum_lshift(m2, ONE, p);
2580 bignum_normalize(m2);
2581
2582 base2 = (SG_BIGNUM_GET_SIGN(bx) || Sg_BignumCmp(bx, m1) >= 0)
2583 ? bignum_mod(bx, m1, NULL) : bx;
2584 a1 = (BIGNUM_ONEP(m1))
2585 ? Sg_MakeBignumFromSI(0)
2586 : odd_mod_expt(base2, be, m1);
2587 a2 = bignum_mod_expt2(base, be, p);
2588 y1 = bignum_mod_inverse(m2, m1);
2589 y2 = bignum_mod_inverse(m1, m2);
2590 {
2591 SgBignum *t1, *t2;
2592 t1 = bignum_normalize(bignum_mul(a1, m2));
2593 t1 = bignum_normalize(bignum_mul(t1, y1));
2594 t2 = bignum_normalize(bignum_mul(a2, m1));
2595 t2 = bignum_normalize(bignum_mul(t2, y2));
2596 result = bignum_mod(bignum_normalize(bignum_add(t1, t2)),
2597 bm, NULL);
2598 }
2599 }
2600 return (invertp) ? bignum_mod_inverse(result, bm) : result;
2601 }
2602
Sg_BignumModExpt(SgBignum * bx,SgBignum * be,SgBignum * bm)2603 SgObject Sg_BignumModExpt(SgBignum *bx, SgBignum *be, SgBignum *bm)
2604 {
2605 if (BIGNUM_ZEROP(be)) {
2606 return (BIGNUM_ONEP(bm)) ? SG_MAKE_INT(0) : SG_MAKE_INT(1);
2607 }
2608 if (BIGNUM_ONEP(bx)) {
2609 return (BIGNUM_ONEP(bm)) ? SG_MAKE_INT(0) : SG_MAKE_INT(1);
2610 }
2611 if (BIGNUM_ZEROP(bx)) {
2612 return (BIGNUM_ONEP(bm)) ? SG_MAKE_INT(0) : SG_MAKE_INT(1);
2613 }
2614 return Sg_NormalizeBignum(bignum_mod_expt(bx, be, bm));
2615 }
2616
Sg_BignumModInverse(SgBignum * bx,SgBignum * bm)2617 SgObject Sg_BignumModInverse(SgBignum *bx, SgBignum *bm)
2618 {
2619 return Sg_NormalizeBignum(bignum_mod_inverse(bx, bm));
2620 }
2621
2622
2623 /* using slide window algorithm. */
bfffo(ulong x)2624 static int bfffo(ulong x)
2625 {
2626 static int tabshi[16]={ 4, 3, 2, 2, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0 };
2627 int value = WORD_BITS - 4;
2628 ulong arg1 = x;
2629 #if SIZEOF_LONG == 8
2630 if (arg1 & ~0xffffffffUL) {value -= 32; arg1 >>= 32;}
2631 #endif
2632 if (arg1 & ~0xffffUL) {value -= 16; arg1 >>= 16;}
2633 if (arg1 & ~0x00ffUL) {value -= 8; arg1 >>= 8;}
2634 if (arg1 & ~0x000fUL) {value -= 4; arg1 >>= 4;}
2635 return value + tabshi[arg1];
2636 }
2637
2638 /* helper to avoid memory allocation. */
compute_lr_binary_size(long j,long m,long len)2639 static long compute_lr_binary_size(long j, long m, long len)
2640 {
2641 long r = len;
2642 for (; j; m <<= 1, j--) {
2643 r <<= 1;
2644 if (m < 0) {
2645 r += len;
2646 }
2647 }
2648 return r;
2649 }
2650
2651 #define copy_element(src, dst, size) \
2652 do { \
2653 long i; \
2654 for (i = 0; i < size; i++) (src)[i] = (dst)[i]; \
2655 } while (0)
2656
leftright_binray_expt(SgBignum * b,long e)2657 static SgBignum * leftright_binray_expt(SgBignum *b, long e)
2658 {
2659 SgBignum *y, *ye, *prod, *t;
2660 long m = e, j = 1 + bfffo(m);
2661 long size;
2662 long ylen = SG_BIGNUM_GET_COUNT(b);
2663 m <<= j;
2664 j = WORD_BITS-j;
2665
2666 /* setup buffer */
2667 size = compute_lr_binary_size(j, m, SG_BIGNUM_GET_COUNT(b));
2668 y = make_bignum(size);
2669 bignum_copy(y, b);
2670 ye = y;
2671 ALLOC_TEMP_BIGNUM_REC(prod, size);
2672 for (; j; m <<= 1, j--) {
2673 /* y = bignum_mul(y, y); */
2674 /* TODO make this for bignum */
2675 /* square_to_len(ye->elements, ylen, prod->elements); */
2676 mp_square(prod->elements, ylen<<1, ye->elements, ylen);
2677 ylen <<= 1;
2678 t = ye;
2679 ye = prod;
2680 prod = t;
2681 /* set ylen to ye */
2682 SG_BIGNUM_SET_COUNT(ye, ylen);
2683 if (m < 0) {
2684 bignum_mul_int(prod, ye, b);
2685 ylen += SG_BIGNUM_GET_COUNT(b);
2686 t = ye;
2687 ye = prod;
2688 prod = t;
2689 }
2690 }
2691 if (y != ye) {
2692 copy_element(y->elements, ye->elements, ylen);
2693 }
2694 SG_BIGNUM_SET_COUNT(y, ylen);
2695 return y;
2696 }
2697
vals(ulong z)2698 static long vals(ulong z)
2699 {
2700 static char tab[64] =
2701 {-1, 0, 1, 12, 2, 6, -1, 13, 3, -1, 7, -1, -1, -1, -1, 14,
2702 10, 4, -1, -1, 8, -1, -1, 25, -1, -1, -1, -1, -1, 21, 27, 15,
2703 31, 11, 5, -1, -1, -1, -1, -1, 9, -1, -1, 24, -1, -1, 20, 26,
2704 30, -1, -1, -1, -1, 23, -1, 19, 29, -1, 22, 18, 28, 17, 16, -1};
2705
2706 #if SIZEOF_INT == 8
2707 long s;
2708 #endif
2709
2710 if (!z) return -1;
2711 #if SIZEOF_INT == 8
2712 if (! (z&0xffffffff)) { s = 32; z >>=32; } else s = 0;
2713 #endif
2714 z |= ~z + 1;
2715 z += z << 4;
2716 z += z << 6;
2717 z ^= z << 16; /* or z -= z<<16 */
2718 #if SIZEOF_INT == 8
2719 return s + tab[(z&0xffffffff)>>26];
2720 #else
2721 return tab[z>>26];
2722 #endif
2723 }
2724
2725 #define BIGNUM_CLEAR_LEFT(bn, size) \
2726 do { \
2727 long i, len = (size); \
2728 long off = SG_BIGNUM_GET_COUNT(bn)-1; \
2729 for (i = 0; i < len; i++) (bn)->elements[off-i] = 0; \
2730 } while (0)
2731 /* I'm not sure if I can write this as a static function since
2732 it's using alloca to setup table. */
2733 #define setup_table(table, u, b) \
2734 do { \
2735 SgBignum *x2; \
2736 long i_, x2len =SG_BIGNUM_GET_COUNT(b)<<1; \
2737 table[0] = b; \
2738 ALLOC_TEMP_BIGNUM_REC(x2, x2len); \
2739 BIGNUM_CLEAR_LEFT(x2, 2); \
2740 mp_square(x2->elements, b->size<<1, b->elements, b->size); \
2741 for (i_ = 1; i_ < (u); i_++) { \
2742 long size = SG_BIGNUM_GET_COUNT(table[i_-1]) + x2len; \
2743 ALLOC_TEMP_BIGNUM_REC(table[i_], size); \
2744 BIGNUM_CLEAR_LEFT(table[i_], 2); \
2745 bignum_mul_int(table[i_], table[i_-1], x2); \
2746 } \
2747 } while (0)
2748
2749 /* compute result size of sliding window.
2750 this is actually the same algorithm without computing */
compute_sw_size(long l,long e,long n,SgBignum ** table)2751 static long compute_sw_size(long l, long e, long n, SgBignum **table)
2752 {
2753 long z = 0;
2754 long tw;
2755 long w, v, i;
2756 while (l >= 0) {
2757 if (e > l+1) e = l + 1;
2758 w = (n >>(l+1-e)) & ((1UL<<2)-1);
2759 v = vals(w);
2760 l -= e;
2761 tw = SG_BIGNUM_GET_COUNT(table[w>>(v+1)]);
2762 if (z) {
2763 for (i = 1; i <= e-v; i++) z <<= 1; /* square z */
2764 z += tw; /* mul */
2765 } else {
2766 z = tw;
2767 }
2768 for (i = 1; i <= v; i++) z <<= 1; /* square z */
2769 while (l >= 0) {
2770 if (n & (1UL<<l)) break;
2771 z <<= 1; /* square z */
2772 l--;
2773 }
2774 }
2775 return z;
2776 }
2777
2778 /* unfortunately, this wasn't rgith from the begining
2779 but didn't come up so that never reached here... */
sliding_window_expt(SgBignum * b,long n,long e)2780 static SgBignum * sliding_window_expt(SgBignum *b, long n, long e)
2781 {
2782 /* max must be 4 */
2783 long volatile zsize, zlen = 0;
2784 long i, l = (WORD_BITS-1) - (long)bfffo(n), u = 1UL<<(e-1);
2785 long w, v;
2786 SgBignum *r, *table[1UL<<(3-1)], *tw, *z = NULL, *prod1, *prod2;
2787
2788 setup_table(table, u, b);
2789 zsize = compute_sw_size(l, e, n, table);
2790 /* we don't need to clear the buffer here. */
2791 r = make_bignum_rec(zsize, FALSE);
2792 SG_BIGNUM_SET_SIGN(r, 1);
2793
2794 prod1 = r;
2795 if (Sg_AvailableStackSize((volatile uintptr_t)&zsize) >= zsize) {
2796 ALLOC_TEMP_BIGNUM_REC(prod2, zsize);
2797 } else {
2798 prod2 = make_bignum_rec(zsize, FALSE);
2799 }
2800
2801 while (l >= 0) {
2802 int index;
2803 if (e > l+1) e = l + 1;
2804 w = (n >>(l+1-e)) & ((1UL<<2)-1);
2805 v = vals(w);
2806 l -= e;
2807 index = (int)(w>>(v+1));
2808 tw = table[index];
2809 if (z) {
2810 SgBignum *t;
2811 ulong len = e-v;
2812 for (i = 1; i <= len; i++) {
2813 /* z = bignum_mul(z, z); */
2814 /* square_to_len(z->elements, zlen, prod1->elements); */
2815 mp_square(prod1->elements, zlen<<1, z->elements, zlen);
2816 zlen <<= 1;
2817 /* swap buffer */
2818 t = z;
2819 z = prod1;
2820 prod1 = t;
2821 SG_BIGNUM_SET_COUNT(z, zlen);
2822 }
2823 /* z = bignum_mul(z, tw); */
2824 bignum_mul_int(prod1, z, tw);
2825 zlen += SG_BIGNUM_GET_COUNT(tw);
2826 /* swap buffer */
2827 t = z;
2828 z = prod1;
2829 prod1 = t;
2830 SG_BIGNUM_SET_COUNT(z, zlen);
2831 } else {
2832 /* z = tw; */
2833 /* setup buffer */
2834 copy_element(prod1->elements, tw->elements, SG_BIGNUM_GET_COUNT(tw));
2835 z = prod1;
2836 prod1 = prod2;
2837 zlen = SG_BIGNUM_GET_COUNT(tw);
2838 SG_BIGNUM_SET_COUNT(z, zlen);
2839 }
2840 for (i = 1; i <= v; i++) {
2841 SgBignum *t;
2842 /* square_to_len(z->elements, zlen, prod1->elements); */
2843 mp_square(prod1->elements, zlen<<1, z->elements, zlen);
2844 zlen <<= 1;
2845 t = z;
2846 z = prod1;
2847 prod1 = t;
2848 SG_BIGNUM_SET_COUNT(z, zlen);
2849 }
2850 while (l >= 0) {
2851 SgBignum *t;
2852 if (n & (1UL<<l)) break;
2853 /* z = bignum_mul(z, z); */
2854 /* square_to_len(z->elements, zlen, prod1->elements); */
2855 mp_square(prod1->elements, zlen<<1, z->elements, zlen);
2856 zlen <<= 1;
2857 /* swap buffer */
2858 t = z;
2859 z = prod1;
2860 prod1 = t;
2861 SG_BIGNUM_SET_COUNT(z, zlen);
2862 l--;
2863 }
2864 }
2865 /* i'm not sure which would the proper one so check */
2866 if (r != z)
2867 copy_element(r->elements, z->elements, zsize);
2868 SG_BIGNUM_SET_COUNT(r, zsize);
2869 return r;
2870 }
2871
bignum_expt(SgBignum * b,long exponent)2872 static SgBignum * bignum_expt(SgBignum *b, long exponent)
2873 {
2874 /* exponent > 0 */
2875 long l = (WORD_BITS-1) - (long)bfffo((ulong)exponent);
2876 if (l <= 8) return leftright_binray_expt(b, exponent);
2877 else if (l <= 24) return sliding_window_expt(b, exponent, 2);
2878 else return sliding_window_expt(b, exponent, 3);
2879 }
2880
2881
Sg_BignumExpt(SgBignum * b,long exponent)2882 SgObject Sg_BignumExpt(SgBignum *b, long exponent)
2883 {
2884 /* let's try the easiest one */
2885 ASSERT(exponent >= 0);
2886 /* let's handle the rare case, this must be handled by Sg_Expt */
2887 if (SG_BIGNUM_GET_SIGN(b) == 0) {
2888 return SG_MAKE_INT(0);
2889 } else if (!exponent) {
2890 return SG_MAKE_INT(1);
2891 } else if (exponent == 1) {
2892 return SG_OBJ(b);
2893 } else {
2894 return Sg_NormalizeBignum(bignum_expt(b, exponent));
2895 }
2896 }
2897
Sg__InitBignum()2898 void Sg__InitBignum()
2899 {
2900 int i;
2901 ZERO = Sg_MakeBignumFromSI(0);
2902 RADIXES[0] = ZERO;
2903 DBL_RADIXES[0] = log(0);
2904 for (i = 1; i <= MAX_RADIX; i++) {
2905 RADIXES[i] = Sg_MakeBignumFromSI(i);
2906 DBL_RADIXES[i] = log(i);
2907 }
2908 /* DBL_LOG2 = DBL_RADIXES[2]; */
2909 }
2910
2911
2912 /*
2913 end of file
2914 Local Variables:
2915 coding: utf-8-unix
2916 End:
2917 */
2918