libhebimath.h - libzahl - big integer library | |
git clone git://git.suckless.org/libzahl | |
Log | |
Files | |
Refs | |
README | |
LICENSE | |
--- | |
libhebimath.h (11909B) | |
--- | |
1 #include <hebimath.h> | |
2 | |
3 #include <setjmp.h> | |
4 #include <stddef.h> | |
5 #include <stdint.h> | |
6 #include <stdio.h> | |
7 #include <stdlib.h> | |
8 | |
9 #define BIGINT_LIBRARY "hebimath" | |
10 #define HEBIMATH | |
11 | |
12 typedef hebi_int z_t[1]; | |
13 | |
14 static z_t _0, _1, _2, _add_sub_a, _add_sub_b, _cmp_b, _pow_a, _pow_m; | |
15 static z_t _trunc_btest_a, _bits_lsb_a, _str_a, _str_b, _str_m, _bset_a; | |
16 static z_t _logic_a, _logic_b, _logic_x, _not_a, _not_b, _gcd_a, _gcd_b; | |
17 static z_t _shift_b, _div_a, _div_b, _divmod; | |
18 static int error; | |
19 static jmp_buf jbuf; | |
20 | |
21 #define FAST_RANDOM 0 | |
22 #define SECURE_RANDOM 0 | |
23 #define DEFAULT_RANDOM 0 | |
24 #define FASTEST_RANDOM 0 | |
25 #define LIBC_RAND_RANDOM 0 | |
26 #define LIBC_RANDOM_RANDOM 0 | |
27 #define LIBC_RAND48_RANDOM 0 | |
28 #define QUASIUNIFORM 0 | |
29 #define UNIFORM 1 | |
30 #define MODUNIFORM 2 | |
31 | |
32 #ifdef ZAHL_UNSAFE | |
33 # define try(expr) (expr) | |
34 #else | |
35 # define try(expr) do if ((error = (expr))) longjmp(jbuf, 1); while (0) | |
36 #endif | |
37 | |
38 #define zsave(a, s) zstr(a, s, sizeof(s) - 1) | |
39 #define zload(a, s) zsets(a, s) | |
40 | |
41 static inline void | |
42 zsetup(jmp_buf env) | |
43 { | |
44 *jbuf = *env; | |
45 hebi_init(_0); | |
46 hebi_init(_1); | |
47 hebi_init(_2); | |
48 hebi_init(_add_sub_a); | |
49 hebi_init(_add_sub_b); | |
50 hebi_init(_cmp_b); | |
51 hebi_init(_pow_a); | |
52 hebi_init(_pow_m); | |
53 hebi_init(_trunc_btest_a); | |
54 hebi_init(_bits_lsb_a); | |
55 hebi_init(_str_a); | |
56 hebi_init(_str_b); | |
57 hebi_init(_str_m); | |
58 hebi_init(_bset_a); | |
59 hebi_init(_logic_a); | |
60 hebi_init(_logic_b); | |
61 hebi_init(_logic_x); | |
62 hebi_init(_not_a); | |
63 hebi_init(_not_b); | |
64 hebi_init(_gcd_a); | |
65 hebi_init(_gcd_b); | |
66 hebi_init(_shift_b); | |
67 hebi_init(_div_a); | |
68 hebi_init(_div_b); | |
69 hebi_init(_divmod); | |
70 try(hebi_set_u(_0, 0)); | |
71 try(hebi_set_u(_1, 1)); | |
72 try(hebi_set_u(_2, 2)); | |
73 } | |
74 | |
75 static inline void | |
76 zunsetup(void) | |
77 { | |
78 hebi_destroy(_0); | |
79 hebi_destroy(_1); | |
80 hebi_destroy(_2); | |
81 hebi_destroy(_add_sub_a); | |
82 hebi_destroy(_add_sub_b); | |
83 hebi_destroy(_cmp_b); | |
84 hebi_destroy(_pow_a); | |
85 hebi_destroy(_pow_m); | |
86 hebi_destroy(_trunc_btest_a); | |
87 hebi_destroy(_bits_lsb_a); | |
88 hebi_destroy(_str_a); | |
89 hebi_destroy(_str_b); | |
90 hebi_destroy(_str_m); | |
91 hebi_destroy(_bset_a); | |
92 hebi_destroy(_logic_a); | |
93 hebi_destroy(_logic_b); | |
94 hebi_destroy(_logic_x); | |
95 hebi_destroy(_not_a); | |
96 hebi_destroy(_not_b); | |
97 hebi_destroy(_gcd_a); | |
98 hebi_destroy(_gcd_b); | |
99 hebi_destroy(_shift_b); | |
100 hebi_destroy(_div_a); | |
101 hebi_destroy(_div_b); | |
102 hebi_destroy(_divmod); | |
103 } | |
104 | |
105 static inline void | |
106 zperror(const char *str) | |
107 { | |
108 const char *serr; | |
109 switch (error) { | |
110 case hebi_badrange: serr = "hebi_badrange"; break; | |
111 case hebi_badvalue: serr = "hebi_badvalue"; break; | |
112 case hebi_nomem: serr = "hebi_nomem"; break; | |
113 default: serr = "unknown error"; break; | |
114 } | |
115 if (str && *str) | |
116 fprintf(stderr, "%s: %s\n", str, serr); | |
117 else | |
118 fprintf(stderr, "%s\n", serr); | |
119 } | |
120 | |
121 static inline void | |
122 zinit(z_t a) | |
123 { | |
124 hebi_init(a); | |
125 } | |
126 | |
127 static inline void | |
128 zfree(z_t a) | |
129 { | |
130 hebi_destroy(a); | |
131 } | |
132 | |
133 static inline void | |
134 zset(z_t r, const z_t a) | |
135 { | |
136 try(hebi_set(r, a)); | |
137 } | |
138 | |
139 static inline void | |
140 zsetu(z_t r, unsigned long long int a) | |
141 { | |
142 try(hebi_set_u(r, a)); | |
143 } | |
144 | |
145 static inline void | |
146 zseti(z_t r, long long int a) | |
147 { | |
148 try(hebi_set_i(r, a)); | |
149 } | |
150 | |
151 static inline void | |
152 zneg(z_t r, const z_t a) | |
153 { | |
154 try(hebi_neg(r, a)); | |
155 } | |
156 | |
157 static inline void | |
158 zabs(z_t r, const z_t a) | |
159 { | |
160 if (hebi_sign(a) < 0) | |
161 zneg(r, a); | |
162 else | |
163 zset(r, a); | |
164 } | |
165 | |
166 static inline void | |
167 zadd(z_t r, const z_t a, const z_t b) | |
168 { | |
169 try(hebi_add(r, a, b)); | |
170 } | |
171 | |
172 static inline void | |
173 zsub(z_t r, const z_t a, const z_t b) | |
174 { | |
175 try(hebi_sub(r, a, b)); | |
176 } | |
177 | |
178 static inline void | |
179 zadd_unsigned(z_t r, const z_t a, const z_t b) | |
180 { | |
181 zabs(_add_sub_a, a); | |
182 zabs(_add_sub_b, b); | |
183 zadd(r, _add_sub_a, _add_sub_b); | |
184 } | |
185 | |
186 static inline void | |
187 zsub_unsigned(z_t r, const z_t a, const z_t b) | |
188 { | |
189 zabs(_add_sub_a, a); | |
190 zabs(_add_sub_b, b); | |
191 zsub(r, _add_sub_a, _add_sub_b); | |
192 } | |
193 | |
194 static inline int | |
195 zeven(const z_t a) | |
196 { | |
197 return ~hebi_get_u(a) & 1; | |
198 } | |
199 | |
200 static inline int | |
201 zodd(const z_t a) | |
202 { | |
203 return hebi_get_u(a) & 1; | |
204 } | |
205 | |
206 static inline int | |
207 zeven_nonzero(const z_t a) | |
208 { | |
209 return zeven(a); | |
210 } | |
211 | |
212 static inline int | |
213 zodd_nonzero(const z_t a) | |
214 { | |
215 return zodd(a); | |
216 } | |
217 | |
218 static inline void | |
219 zswap(z_t a, z_t b) | |
220 { | |
221 hebi_swap(a, b); | |
222 } | |
223 | |
224 static inline int | |
225 zcmpmag(const z_t a, const z_t b) | |
226 { | |
227 return hebi_cmp_mag(a, b); | |
228 } | |
229 | |
230 static inline int | |
231 zcmp(const z_t a, const z_t b) | |
232 { | |
233 return hebi_cmp(a, b); | |
234 } | |
235 | |
236 static inline int | |
237 zcmpi(const z_t a, long long int b) | |
238 { | |
239 zseti(_cmp_b, b); | |
240 return zcmp(a, _cmp_b); | |
241 } | |
242 | |
243 static inline int | |
244 zcmpu(const z_t a, long long int b) | |
245 { | |
246 zsetu(_cmp_b, b); | |
247 return zcmp(a, _cmp_b); | |
248 } | |
249 | |
250 static inline int | |
251 zsignum(const z_t a) | |
252 { | |
253 return zcmp(a, _0); | |
254 } | |
255 | |
256 static inline int | |
257 zzero(const z_t a) | |
258 { | |
259 return !zsignum(a); | |
260 } | |
261 | |
262 static inline void | |
263 zmul(z_t r, const z_t a, const z_t b) | |
264 { | |
265 try(hebi_mul(r, a, b)); | |
266 } | |
267 | |
268 static inline void | |
269 zsqr(z_t r, const z_t a) | |
270 { | |
271 zmul(r, a, a); | |
272 } | |
273 | |
274 static inline void | |
275 zsets(z_t a, const char *s) | |
276 { | |
277 try(hebi_set_str(a, s, 0, 10)); | |
278 } | |
279 | |
280 static inline void | |
281 zdivmod(z_t q, z_t r, const z_t a, const z_t b) | |
282 { | |
283 try(hebi_div(q, r, a, b)); | |
284 } | |
285 | |
286 static inline void | |
287 zdiv(z_t q, const z_t a, const z_t b) | |
288 { | |
289 zdivmod(q, 0, a, b); | |
290 } | |
291 | |
292 static inline void | |
293 zmod(z_t r, const z_t a, const z_t b) | |
294 { | |
295 zdivmod(_divmod, r, a, b); | |
296 } | |
297 | |
298 static inline void | |
299 zmodmul(z_t r, const z_t a, const z_t b, const z_t m) | |
300 { | |
301 zmul(r, a, b); | |
302 zmod(r, r, m); | |
303 } | |
304 | |
305 static inline void | |
306 zmodsqr(z_t r, const z_t a, const z_t m) | |
307 { | |
308 zsqr(r, a); | |
309 zmod(r, r, m); | |
310 } | |
311 | |
312 static inline void | |
313 zpowu(z_t r, const z_t a, unsigned long long int b) | |
314 { | |
315 if (!b) { | |
316 if (zzero(a)) { | |
317 error = hebi_badvalue; | |
318 longjmp(jbuf, 1); | |
319 } | |
320 zsetu(r, 1); | |
321 return; | |
322 } | |
323 | |
324 zset(_pow_a, a); | |
325 zsetu(r, 1); | |
326 for (; b; b >>= 1) { | |
327 if (b & 1) | |
328 zmul(r, r, _pow_a); | |
329 zsqr(_pow_a, _pow_a); | |
330 } | |
331 } | |
332 | |
333 static inline void | |
334 zpow(z_t r, const z_t a, const z_t b) | |
335 { | |
336 zpowu(r, a, hebi_get_u(b)); | |
337 } | |
338 | |
339 static inline void | |
340 zmodpowu(z_t r, const z_t a, unsigned long long int b, const z_t m) | |
341 { | |
342 if (!b) { | |
343 if (zzero(a) || zzero(m)) { | |
344 error = hebi_badvalue; | |
345 longjmp(jbuf, 1); | |
346 } | |
347 zsetu(r, 1); | |
348 return; | |
349 } else if (zzero(m)) { | |
350 error = hebi_badvalue; | |
351 longjmp(jbuf, 1); | |
352 } | |
353 | |
354 zmod(_pow_a, a, m); | |
355 zset(_pow_m, m); | |
356 zsetu(r, 1); | |
357 for (; b; b >>= 1) { | |
358 if (b & 1) | |
359 zmodmul(r, r, _pow_a, _pow_m); | |
360 zmodsqr(_pow_a, _pow_a, _pow_m); | |
361 } | |
362 } | |
363 | |
364 static inline void | |
365 zmodpow(z_t r, const z_t a, const z_t b, const z_t m) | |
366 { | |
367 zmodpowu(r, a, hebi_get_u(b), m); | |
368 } | |
369 | |
370 static inline void | |
371 zlsh(z_t r, const z_t a, size_t b) | |
372 { | |
373 try(hebi_shl(r, a, b)); | |
374 } | |
375 | |
376 static inline void | |
377 zrsh(z_t r, const z_t a, size_t b) | |
378 { | |
379 try(hebi_shr(r, a, b)); | |
380 } | |
381 | |
382 static inline void | |
383 ztrunc(z_t r, const z_t a, size_t b) | |
384 { | |
385 zrsh(_trunc_btest_a, a, b); | |
386 zlsh(_trunc_btest_a, _trunc_btest_a, b); | |
387 zsub(r, a, _trunc_btest_a); | |
388 } | |
389 | |
390 static inline int | |
391 zbtest(z_t a, size_t b) | |
392 { | |
393 zlsh(_trunc_btest_a, a, b); | |
394 return zodd(a); | |
395 } | |
396 | |
397 static inline size_t | |
398 zbits(const z_t a) | |
399 { | |
400 hebi_uword x = x; | |
401 size_t rc = 0; | |
402 zset(_bits_lsb_a, a); | |
403 while (zsignum(_bits_lsb_a)) { | |
404 x = hebi_get_u(_bits_lsb_a); | |
405 zrsh(_bits_lsb_a, _bits_lsb_a, 8 * sizeof(x)); | |
406 rc += 8 * sizeof(x); | |
407 } | |
408 if (rc) | |
409 rc -= 8 * sizeof(x); | |
410 while (x) { | |
411 x >>= 1; | |
412 rc += 1; | |
413 } | |
414 return rc; | |
415 } | |
416 | |
417 static inline size_t | |
418 zlsb(const z_t a) | |
419 { | |
420 hebi_uword x; | |
421 size_t rc = 0; | |
422 if (zzero(a)) | |
423 return SIZE_MAX; | |
424 zset(_bits_lsb_a, a); | |
425 while (!(x = hebi_get_u(_bits_lsb_a))) { | |
426 zrsh(_bits_lsb_a, _bits_lsb_a, 8 * sizeof(x)); | |
427 rc += 8 * sizeof(x); | |
428 } | |
429 while (~x & 1) { | |
430 x >>= 1; | |
431 rc += 1; | |
432 } | |
433 return rc; | |
434 } | |
435 | |
436 static inline void | |
437 zptest(z_t w, const z_t a, int t) | |
438 { | |
439 static int gave_up = 0; | |
440 if (!gave_up) { | |
441 gave_up = 1; | |
442 printf("I'm sorry, primality test has not been implement… | |
443 } | |
444 (void) w; | |
445 (void) a; | |
446 (void) t; | |
447 } | |
448 | |
449 static inline size_t | |
450 zstr_length(const z_t a, unsigned long long int radix) | |
451 { | |
452 size_t size_total = 1, size_temp; | |
453 zset(_str_a, a); | |
454 while (!zzero(_str_a)) { | |
455 zsetu(_str_m, radix); | |
456 zset(_str_b, _str_m); | |
457 size_temp = 1; | |
458 while (zcmpmag(_str_m, _str_a) <= 0) { | |
459 zset(_str_b, _str_m); | |
460 zsqr(_str_m, _str_m); | |
461 size_temp <<= 1; | |
462 } | |
463 size_temp >>= 1; | |
464 size_total += size_temp; | |
465 zdiv(_str_a, _str_a, _str_b); | |
466 } | |
467 return size_total + (zsignum(a) < 0); | |
468 } | |
469 | |
470 static inline char * | |
471 zstr(const z_t a, char *s, size_t n) | |
472 { | |
473 if (!n) | |
474 n = zstr_length(a, 10) * 2 + 1; | |
475 try(hebi_get_str(s, n, a, 10)); | |
476 return s; | |
477 } | |
478 | |
479 static inline void | |
480 zsplit(z_t high, z_t low, const z_t a, size_t brk) | |
481 { | |
482 if (low == a) { | |
483 zrsh(high, a, brk); | |
484 ztrunc(low, a, brk); | |
485 } else { | |
486 ztrunc(low, a, brk); | |
487 zrsh(high, a, brk); | |
488 } | |
489 } | |
490 | |
491 static inline void | |
492 zbset(z_t r, const z_t a, size_t bit, int mode) | |
493 { | |
494 zrsh(_bset_a, a, bit); | |
495 if (mode && zeven(_bset_a)) { | |
496 zlsh(_bset_a, _1, bit); | |
497 zadd(r, a, _bset_a); | |
498 } else if (mode <= 0 && zodd(_bset_a)) { | |
499 zlsh(_bset_a, _1, bit); | |
500 zsub(r, a, _bset_a); | |
501 } else { | |
502 zset(r, a); | |
503 } | |
504 } | |
505 | |
506 static inline void | |
507 zrand(z_t r, int dev, int dist, const z_t n) | |
508 { | |
509 static int gave_up[] = {0, 0, 0}; | |
510 if (!gave_up[dist]) { | |
511 gave_up[dist] = 1; | |
512 printf("I'm sorry, prng has not been implemented.\n"); | |
513 } | |
514 (void) r; | |
515 (void) dev; | |
516 (void) n; | |
517 } | |
518 | |
519 static inline void | |
520 zand(z_t r, const z_t a, const z_t b) | |
521 { | |
522 int neg = hebi_sign(a) < 0 && hebi_sign(b) < 0; | |
523 hebi_uword x; | |
524 size_t i = 0; | |
525 zset(_logic_a, a); | |
526 zset(_logic_b, b); | |
527 zsetu(r, 0); | |
528 while (zsignum(_logic_a) && zsignum(_logic_b)) { | |
529 x = hebi_get_u(_logic_a) & hebi_get_u(_logic_b); | |
530 zsetu(_logic_x, x); | |
531 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); | |
532 zadd(r, r, _logic_x); | |
533 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); | |
534 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); | |
535 } | |
536 if (neg) | |
537 zneg(r, r); | |
538 } | |
539 | |
540 static inline void | |
541 zor(z_t r, const z_t a, const z_t b) | |
542 { | |
543 int neg = hebi_sign(a) < 0 || hebi_sign(b) < 0; | |
544 hebi_uword x; | |
545 size_t i = 0; | |
546 zset(_logic_a, a); | |
547 zset(_logic_b, b); | |
548 zsetu(r, 0); | |
549 while (zsignum(_logic_a) || zsignum(_logic_b)) { | |
550 x = hebi_get_u(_logic_a) | hebi_get_u(_logic_b); | |
551 zsetu(_logic_x, x); | |
552 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); | |
553 zadd(r, r, _logic_x); | |
554 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); | |
555 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); | |
556 } | |
557 if (neg) | |
558 zneg(r, r); | |
559 } | |
560 | |
561 static inline void | |
562 zxor(z_t r, const z_t a, const z_t b) | |
563 { | |
564 int neg = (hebi_sign(a) < 0) ^ (hebi_sign(b) < 0); | |
565 hebi_uword x; | |
566 size_t i = 0; | |
567 zset(_logic_a, a); | |
568 zset(_logic_b, b); | |
569 zsetu(r, 0); | |
570 while (zsignum(_logic_a) || zsignum(_logic_b)) { | |
571 x = hebi_get_u(_logic_a) ^ hebi_get_u(_logic_b); | |
572 zsetu(_logic_x, x); | |
573 zlsh(_logic_x, _logic_x, i++ * 8 * sizeof(x)); | |
574 zadd(r, r, _logic_x); | |
575 zrsh(_logic_a, _logic_a, 8 * sizeof(x)); | |
576 zrsh(_logic_b, _logic_b, 8 * sizeof(x)); | |
577 } | |
578 if (neg) | |
579 zneg(r, r); | |
580 } | |
581 | |
582 static inline void | |
583 zgcd(z_t r, const z_t a, const z_t b) | |
584 { | |
585 size_t shifts, a_lsb, b_lsb; | |
586 int neg, cmpmag; | |
587 | |
588 if (zzero(a)) { | |
589 if (r != b) | |
590 zset(r, b); | |
591 return; | |
592 } | |
593 if (zzero(b)) { | |
594 if (r != a) | |
595 zset(r, a); | |
596 return; | |
597 } | |
598 | |
599 neg = hebi_sign(a) < 0 && hebi_sign(b) < 0; | |
600 | |
601 a_lsb = zlsb(a); | |
602 b_lsb = zlsb(b); | |
603 shifts = a_lsb < b_lsb ? a_lsb : b_lsb; | |
604 zrsh(_gcd_a, a, a_lsb); | |
605 zrsh(_gcd_b, b, b_lsb); | |
606 | |
607 for (;;) { | |
608 if ((cmpmag = zcmpmag(_gcd_a, _gcd_b)) >= 0) { | |
609 if (cmpmag == 0) | |
610 break; | |
611 zswap(_gcd_a, _gcd_b); | |
612 } | |
613 zsub(_gcd_b, _gcd_b, _gcd_a); | |
614 zrsh(_gcd_b, _gcd_b, zlsb(_gcd_b)); | |
615 } | |
616 | |
617 zlsh(r, _gcd_a, shifts); | |
618 if (neg) | |
619 zneg(r, r); | |
620 } | |
621 | |
622 static inline void | |
623 znot(z_t r, const z_t a) | |
624 { | |
625 size_t bits = zbits(a); | |
626 zsetu(_not_b, 0); | |
627 zsetu(_not_a, 1); | |
628 zlsh(_not_a, _not_a, bits); | |
629 zadd(_not_b, _not_b, _logic_a); | |
630 zsub(_not_b, _not_b, _1); | |
631 zxor(r, a, _not_b); | |
632 zneg(r, r); | |
633 } | |
634 | |
635 /* Prototype declared, but implementation missing, in hebimath */ | |
636 | |
637 int | |
638 hebi_shl(hebi_int *r, const hebi_int *a, unsigned int b) | |
639 { | |
640 zpowu(_shift_b, _2, b); | |
641 zmul(r, a, _shift_b); | |
642 return hebi_success; | |
643 } | |
644 | |
645 int | |
646 hebi_shr(hebi_int *r, const hebi_int *a, unsigned int b) | |
647 { | |
648 zpowu(_shift_b, _2, b); | |
649 zdiv(r, a, _shift_b); | |
650 return hebi_success; | |
651 } | |
652 | |
653 int | |
654 hebi_div(hebi_int *q, hebi_int *r, const hebi_int *a, const hebi_int *b) | |
655 { | |
656 int neg = zsignum(a) < 0; | |
657 zset(q, _0); | |
658 zabs(_div_a, a); | |
659 zabs(_div_b, b); | |
660 if (zzero(b)) { | |
661 error = hebi_badvalue; | |
662 longjmp(jbuf, 1); | |
663 } | |
664 while (zcmpmag(_div_a, _div_b) >= 0) { | |
665 zadd(q, q, _1); | |
666 zsub(_div_a, _div_a, _div_b); | |
667 } | |
668 if (neg) | |
669 zneg(q, q); | |
670 if (r) | |
671 zset(r, _div_a); | |
672 return hebi_success; | |
673 } |