libhebimath.h - libzahl - big integer library
(HTM) git clone git://git.suckless.org/libzahl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) 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 implemented.\n");
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 }