libtfm.h - libzahl - big integer library
(HTM) git clone git://git.suckless.org/libzahl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
libtfm.h (7384B)
---
1 #include <tfm.h>
2
3 #include <setjmp.h>
4 #include <stddef.h>
5 #include <stdint.h>
6 #include <stdio.h>
7
8 #define BIGINT_LIBRARY "TomsFastMath"
9
10 #define FAST_RANDOM 0
11 #define SECURE_RANDOM 0
12 #define DEFAULT_RANDOM 0
13 #define FASTEST_RANDOM 0
14 #define LIBC_RAND_RANDOM 0
15 #define LIBC_RANDOM_RANDOM 0
16 #define LIBC_RAND48_RANDOM 0
17 #define QUASIUNIFORM 0
18 #define UNIFORM 1
19 #define MODUNIFORM 2
20
21 typedef fp_int z_t[1];
22
23 static z_t _0, _1, _a, _b;
24 static int _tmp, error;
25 static jmp_buf jbuf;
26
27 #ifdef ZAHL_UNSAFE
28 # define try(expr) (expr)
29 #else
30 # define try(expr) do if ((error = (expr))) longjmp(jbuf, 1); while (0)
31 #endif
32
33 static inline void
34 fp_set_int(z_t a, uint32_t d)
35 {
36 a->dp[0] = d;
37 a->used = !!d;
38 a->sign = 0;
39 }
40
41 static inline void
42 zsetup(jmp_buf env)
43 {
44 *jbuf = *env;
45 fp_init(_0);
46 fp_init(_1);
47 fp_init(_a);
48 fp_init(_b);
49 fp_set_int(_0, 0);
50 fp_set_int(_1, 1);
51 }
52
53 static inline void
54 zunsetup(void)
55 {
56 /* Do nothing */
57 }
58
59 static inline void
60 zperror(const char *str)
61 {
62 const char *serr;
63 switch (error) {
64 case FP_VAL: serr = "FP_VAL"; break;
65 case FP_MEM: serr = "FP_MEM"; break;
66 default: serr = "unknown error"; break;
67 }
68 if (str && *str)
69 fprintf(stderr, "%s: %s\n", str, serr);
70 else
71 fprintf(stderr, "%s\n", serr);
72 }
73
74 static inline void
75 zinit(z_t a)
76 {
77 fp_init(a);
78 }
79
80 static inline void
81 zfree(z_t a)
82 {
83 (void) a;
84 }
85
86 static inline void
87 zadd(z_t r, z_t a, z_t b)
88 {
89 fp_add(a, b, r);
90 }
91
92 static inline void
93 zsub(z_t r, z_t a, z_t b)
94 {
95 fp_sub(a, b, r);
96 }
97
98 static inline void
99 zset(z_t r, z_t a)
100 {
101 fp_copy(a, r);
102 }
103
104 static inline int
105 zeven(z_t a)
106 {
107 return fp_iseven(a);
108 }
109
110 static inline int
111 zodd(z_t a)
112 {
113 return fp_isodd(a);
114 }
115
116 static inline int
117 zeven_nonzero(z_t a)
118 {
119 return fp_iseven(a);
120 }
121
122 static inline int
123 zodd_nonzero(z_t a)
124 {
125 return fp_isodd(a);
126 }
127
128 static inline int
129 zzero(z_t a)
130 {
131 return fp_iszero(a);
132 }
133
134 static inline int
135 zsignum(z_t a)
136 {
137 return fp_cmp(a, _0);
138 }
139
140 static inline size_t
141 zbits(z_t a)
142 {
143 return fp_count_bits(a);
144 }
145
146 static inline size_t
147 zlsb(z_t a)
148 {
149 return fp_cnt_lsb(a);
150 }
151
152 static inline void
153 zlsh(z_t r, z_t a, size_t b)
154 {
155 fp_mul_2d(a, b, r);
156 }
157
158 static inline void
159 zrsh(z_t r, z_t a, size_t b)
160 {
161 fp_div_2d(a, b, r, 0);
162 }
163
164 static inline void
165 ztrunc(z_t r, z_t a, size_t b)
166 {
167 fp_mod_2d(a, b, r);
168 }
169
170 static inline void
171 zgcd(z_t r, z_t a, z_t b)
172 {
173 fp_gcd(a, b, r);
174 }
175
176 static inline void
177 zmul(z_t r, z_t a, z_t b)
178 {
179 fp_mul(a, b, r);
180 }
181
182 static inline void
183 zsqr(z_t r, z_t a)
184 {
185 fp_sqr(a, r);
186 }
187
188 static inline void
189 zmodmul(z_t r, z_t a, z_t b, z_t m)
190 {
191 try(fp_mulmod(a, b, m, r));
192 }
193
194 static inline void
195 zmodsqr(z_t r, z_t a, z_t m)
196 {
197 try(fp_sqrmod(a, m, r));
198 }
199
200 static inline void
201 zsets(z_t a, char *s)
202 {
203 try(fp_read_radix(a, s, 10));
204 }
205
206 static inline size_t
207 zstr_length(z_t a, unsigned long long int b)
208 {
209 try(fp_radix_size(a, b, &_tmp));
210 return _tmp;
211 }
212
213 static inline char *
214 zstr(z_t a, char *s, size_t n)
215 {
216 try(fp_toradix(a, s, 10));
217 return s;
218 (void) n;
219 }
220
221 static inline int
222 zptest(z_t w, z_t a, int t)
223 {
224 return fp_isprime_ex(a, t);
225 (void) w; /* Note, the witness is not returned. */
226 }
227
228 static inline void
229 zdivmod(z_t q, z_t r, z_t a, z_t b)
230 {
231 try(fp_div(a, b, q, r));
232 }
233
234 static inline void
235 zdiv(z_t q, z_t a, z_t b)
236 {
237 zdivmod(q, 0, a, b);
238 }
239
240 static inline void
241 zmod(z_t r, z_t a, z_t b)
242 {
243 try(fp_mod(a, b, r));
244 }
245
246 static inline void
247 zneg(z_t r, z_t a)
248 {
249 fp_neg(a, r);
250 }
251
252 static inline void
253 zabs(z_t r, z_t a)
254 {
255 fp_abs(a, r);
256 }
257
258 static inline void
259 zadd_unsigned(z_t r, z_t a, z_t b)
260 {
261 zabs(_a, a);
262 zabs(_b, b);
263 zadd(r, _a, _b);
264 }
265
266 static inline void
267 zsub_unsigned(z_t r, z_t a, z_t b)
268 {
269 zabs(_a, a);
270 zabs(_b, b);
271 zsub(r, _a, _b);
272 }
273
274 static inline void
275 zswap(z_t a, z_t b)
276 {
277 z_t t;
278 zset(t, a);
279 zset(a, b);
280 zset(b, t);
281 }
282
283 static inline void
284 zand(z_t r, z_t a, z_t b)
285 {
286 int i;
287 for (i = 0; i < a->used && i < b->used; i++)
288 r->dp[i] = a->dp[i] & b->dp[i];
289 r->used = i;
290 while (r->used && !r->dp[r->used])
291 r->used--;
292 r->sign = a->sign & b->sign;
293 }
294
295 static inline void
296 zor(z_t r, z_t a, z_t b)
297 {
298 int i;
299 for (i = 0; i < a->used && i < b->used; i++)
300 r->dp[i] = a->dp[i] | b->dp[i];
301 for (; i < a->used; i++)
302 r->dp[i] = a->dp[i];
303 for (; i < b->used; i++)
304 r->dp[i] = b->dp[i];
305 r->used = i;
306 r->sign = a->sign | b->sign;
307 }
308
309 static inline void
310 zxor(z_t r, z_t a, z_t b)
311 {
312 int i;
313 for (i = 0; i < a->used && i < b->used; i++)
314 r->dp[i] = a->dp[i] ^ b->dp[i];
315 for (; i < a->used; i++)
316 r->dp[i] = a->dp[i];
317 for (; i < b->used; i++)
318 r->dp[i] = b->dp[i];
319 r->used = i;
320 while (r->used && !r->dp[r->used])
321 r->used--;
322 r->sign = a->sign ^ b->sign;
323 }
324
325 static inline size_t
326 zsave(z_t a, char *s)
327 {
328 _tmp = !s ? fp_signed_bin_size(a) : (fp_to_signed_bin(a, (unsigned char *)s), 0);
329 return _tmp;
330 }
331
332 static inline size_t
333 zload(z_t a, const char *s)
334 {
335 fp_read_signed_bin(a, (unsigned char *)s, _tmp);
336 return _tmp;
337 }
338
339 static inline void
340 zsetu(z_t r, unsigned long long int val)
341 {
342 uint32_t high = (uint32_t)(val >> 32);
343 uint32_t low = (uint32_t)val;
344
345 if (high) {
346 fp_set_int(r, high);
347 fp_set_int(_a, low);
348 fp_lshd(r, 32);
349 zadd(r, r, _a);
350 } else {
351 fp_set_int(r, low);
352 }
353
354 }
355
356 static inline void
357 zseti(z_t r, long long int val)
358 {
359 if (val >= 0) {
360 zsetu(r, (unsigned long long int)val);
361 } else {
362 zsetu(r, (unsigned long long int)-val);
363 zneg(r, r);
364 }
365 }
366
367 static inline int
368 zcmpmag(z_t a, z_t b)
369 {
370 return fp_cmp_mag(a, b);
371 }
372
373 static inline int
374 zcmp(z_t a, z_t b)
375 {
376 return fp_cmp(a, b);
377 }
378
379 static inline int
380 zcmpi(z_t a, long long int b)
381 {
382 zseti(_b, b);
383 return zcmp(a, _b);
384 }
385
386 static inline int
387 zcmpu(z_t a, unsigned long long int b)
388 {
389 zsetu(_b, b);
390 return zcmp(a, _b);
391 }
392
393 static inline void
394 zsplit(z_t high, z_t low, z_t a, size_t brk)
395 {
396 if (low == a) {
397 zrsh(high, a, brk);
398 ztrunc(low, a, brk);
399 } else {
400 ztrunc(low, a, brk);
401 zrsh(high, a, brk);
402 }
403 }
404
405 static inline void
406 znot(z_t r, z_t a)
407 {
408 fp_2expt(_a, zbits(a));
409 fp_sub_d(_a, 1, _a);
410 zand(r, a, _a);
411 zneg(r, r);
412 }
413
414 static inline int
415 zbtest(z_t a, size_t bit)
416 {
417 fp_2expt(_b, bit);
418 zand(_b, a, _b);
419 return !zzero(_b);
420 }
421
422 static inline void
423 zbset(z_t r, z_t a, size_t bit, int mode)
424 {
425 if (mode > 0) {
426 fp_2expt(_b, bit);
427 zor(r, a, _b);
428 } else if (mode < 0 || zbtest(a, bit)) {
429 fp_2expt(_b, bit);
430 zxor(r, a, _b);
431 }
432 }
433
434 static inline void
435 zrand(z_t r, int dev, int dist, z_t n)
436 {
437 static int gave_up = 0;
438 int bits;
439 (void) dev;
440
441 if (zzero(n)) {
442 fp_zero(r);
443 return;
444 }
445 if (zsignum(n) < 0) {
446 return;
447 }
448
449 bits = zbits(n);
450
451 switch (dist) {
452 case QUASIUNIFORM:
453 fp_rand(r, bits);
454 zadd(r, r, _1);
455 zmul(r, r, n);
456 zrsh(r, r, bits);
457 break;
458
459 case UNIFORM:
460 if (!gave_up) {
461 gave_up = 1;
462 printf("I'm sorry, this is too difficult, I give up.\n");
463 }
464 break;
465
466 case MODUNIFORM:
467 fp_rand(r, bits);
468 if (zcmp(r, n) > 0)
469 zsub(r, r, n);
470 break;
471
472 default:
473 abort();
474 }
475 }
476
477 static inline void
478 zpowu(z_t r, z_t a, unsigned long long int b)
479 {
480 if (!b) {
481 if (zzero(a)) {
482 error = FP_VAL;
483 longjmp(jbuf, 1);
484 }
485 zsetu(r, 1);
486 return;
487 }
488 zset(_a, a);
489 zsetu(r, 1);
490 for (; b; b >>= 1) {
491 if (b & 1)
492 zmul(r, r, _a);
493 zsqr(_a, _a);
494 }
495 }
496
497 static inline void
498 zpow(z_t r, z_t a, z_t b)
499 {
500 zpowu(r, a, b->used ? b->dp[0] : 0);
501 }
502
503 static inline void
504 zmodpow(z_t r, z_t a, z_t b, z_t m)
505 {
506 try(fp_exptmod(a, b, m, r));
507 }
508
509 static inline void
510 zmodpowu(z_t r, z_t a, unsigned long long int b, z_t m)
511 {
512 fp_set_int(_b, b);
513 zmodpow(r, a, _b, m);
514 }