Introduction
Introduction Statistics Contact Development Disclaimer Help
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 }
You are viewing proxied material from suckless.org. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.