1 /* Some functions for modular arithmetic with residues and modulus in
2 unsigned long variables. Due to inlining, this file must be included
3 in the caller's source code with #include */
4
5 /* Naming convention: all function start with modul, for
6 MODulus Unsigned Long, followed by underscore, functionality of function
7 (add, mul, etc), and possibly underscore and specification of what argument
8 types the function takes (_ul, etc). */
9
10 #ifndef MOD_UL_H
11 #define MOD_UL_H
12
13 #include "cado_config.h" // just because we're a header.
14 /**********************************************************************/
15 #include <limits.h> // for LONG_BIT, ULONG_MAX
16 #include <stdint.h> // for uint64_t, int64_t
17 #include <stdlib.h> // for size_t, llabs
18 #include <stdio.h>
19 #include "macros.h" // for MAYBE_UNUSED, ASSERT_EXPENSIVE, ASSERT, ASS...
20 #include "ularith.h" // for ularith_mul_ul_ul_2ul, ularith_div_2ul_ul_ul_r
21
22
23 /*********************************************************************/
24 /* Helper macros, see also ularith.h */
25
26 /* A macro for function renaming. All functions here start with modul_ */
27 #define MODUL_RENAME(x) modul_##x
28
29 #define MODUL_SIZE 1
30 #define MODUL_MAXBITS LONG_BIT
31
32 typedef unsigned long residueul_t[MODUL_SIZE];
33 typedef unsigned long modintul_t[MODUL_SIZE];
34 typedef unsigned long modulusul_t[MODUL_SIZE];
35
36 /* ================= Functions that are part of the API ================= */
37
38 /* Some functions for integers of the same width as the modulus */
39
40 MAYBE_UNUSED
41 static inline void
modul_intinit(modintul_t r)42 modul_intinit (modintul_t r)
43 {
44 r[0] = 0;
45 }
46
47
48 MAYBE_UNUSED
49 static inline void
modul_intclear(modintul_t r MAYBE_UNUSED)50 modul_intclear (modintul_t r MAYBE_UNUSED)
51 {
52 return;
53 }
54
55
56 MAYBE_UNUSED
57 static inline void
modul_intset(modintul_t r,const modintul_t s)58 modul_intset (modintul_t r, const modintul_t s)
59 {
60 r[0] = s[0];
61 }
62
63
64 MAYBE_UNUSED
65 static inline void
modul_intset_ul(modintul_t r,const unsigned long s)66 modul_intset_ul (modintul_t r, const unsigned long s)
67 {
68 r[0] = s;
69 }
70
71
72 /* The two mod*_uls() functions import/export modint_t from/to an array of
73 unsigned longs. For modul_intset_ul, the size of the array is passed as
74 a parameter n. For mod_intget_uls(), the required array size can be
75 determined via mod_intbits(); if the modint_t is zero, mod_intget_uls()
76 writes 0 to the first output unsigned long. It returns the number of
77 unsigned longs written. */
78 MAYBE_UNUSED
79 static inline void
modul_intset_uls(modintul_t r,const unsigned long * s,const size_t n)80 modul_intset_uls (modintul_t r, const unsigned long *s, const size_t n)
81 {
82 ASSERT_ALWAYS(n <= MODUL_SIZE);
83 if (n == 0)
84 r[0] = 0;
85 else
86 r[0] = s[0];
87 }
88
89
90 MAYBE_UNUSED
91 static inline unsigned long
modul_intget_ul(const modintul_t s)92 modul_intget_ul (const modintul_t s)
93 {
94 return s[0];
95 }
96
97
98 MAYBE_UNUSED
99 static inline size_t
modul_intget_uls(unsigned long * r,const modintul_t s)100 modul_intget_uls (unsigned long *r, const modintul_t s)
101 {
102 r[0] = s[0];
103 return 1;
104 }
105
106
107 MAYBE_UNUSED
108 static inline int
modul_intequal(const modintul_t a,const modintul_t b)109 modul_intequal (const modintul_t a, const modintul_t b)
110 {
111 return (a[0] == b[0]);
112 }
113
114
115 MAYBE_UNUSED
116 static inline int
modul_intequal_ul(const modintul_t a,const unsigned long b)117 modul_intequal_ul (const modintul_t a, const unsigned long b)
118 {
119 return (a[0] == b);
120 }
121
122
123 MAYBE_UNUSED
124 static inline int
modul_intcmp(const modintul_t a,const modintul_t b)125 modul_intcmp (const modintul_t a, const modintul_t b)
126 {
127 return (a[0] < b[0]) ? -1 : (a[0] == b[0]) ? 0 : 1;
128 }
129
130
131 MAYBE_UNUSED
132 static inline int
modul_intcmp_ul(const modintul_t a,const unsigned long b)133 modul_intcmp_ul (const modintul_t a, const unsigned long b)
134 {
135 return (a[0] < b) ? -1 : (a[0] == b) ? 0 : 1;
136 }
137
138
139 MAYBE_UNUSED
140 static inline int
modul_intcmp_uint64(const modintul_t a,const uint64_t b)141 modul_intcmp_uint64 (const modintul_t a, const uint64_t b)
142 {
143 if (b > ULONG_MAX)
144 return -1;
145 return (a[0] < b) ? -1 : (a[0] == b) ? 0 : 1;
146 }
147
148
149 MAYBE_UNUSED
150 static inline int
modul_intfits_ul(const modintul_t a MAYBE_UNUSED)151 modul_intfits_ul (const modintul_t a MAYBE_UNUSED)
152 {
153 return 1;
154 }
155
156
157 MAYBE_UNUSED
158 static inline void
modul_intadd(modintul_t r,const modintul_t a,const modintul_t b)159 modul_intadd (modintul_t r, const modintul_t a, const modintul_t b)
160 {
161 r[0] = a[0] + b[0];
162 }
163
164
165 MAYBE_UNUSED
166 static inline void
modul_intsub(modintul_t r,const modintul_t a,const modintul_t b)167 modul_intsub (modintul_t r, const modintul_t a, const modintul_t b)
168 {
169 r[0] = a[0] - b[0];
170 }
171
172
173 MAYBE_UNUSED
174 static inline void
modul_intshr(modintul_t r,const modintul_t s,const int i)175 modul_intshr (modintul_t r, const modintul_t s, const int i)
176 {
177 r[0] = s[0] >> i;
178 }
179
180
181 MAYBE_UNUSED
182 static inline void
modul_intshl(modintul_t r,const modintul_t s,const int i)183 modul_intshl (modintul_t r, const modintul_t s, const int i)
184 {
185 r[0] = s[0] << i;
186 }
187
188
189 /* Returns the number of bits in a, that is, floor(log_2(n))+1.
190 For n==0 returns 0. */
191 MAYBE_UNUSED
192 static inline size_t
modul_intbits(const modintul_t a)193 modul_intbits (const modintul_t a)
194 {
195 if (a[0] == 0)
196 return 0;
197 return LONG_BIT - ularith_clz (a[0]);
198 }
199
200
201 /* r = n/d. We require d|n */
202 MAYBE_UNUSED
203 static inline void
modul_intdivexact(modintul_t r,const modintul_t n,const modintul_t d)204 modul_intdivexact (modintul_t r, const modintul_t n, const modintul_t d)
205 {
206 r[0] = n[0] / d[0];
207 }
208
209
210 /* r = n%d */
211 MAYBE_UNUSED
212 static inline void
modul_intmod(modintul_t r,const modintul_t n,const modintul_t d)213 modul_intmod (modintul_t r, const modintul_t n,
214 const modintul_t d)
215 {
216 r[0] = n[0] % d[0];
217 }
218
219
220 /* Functions for the modulus */
221
222 MAYBE_UNUSED
223 static inline void
modul_initmod_ul(modulusul_t m,const unsigned long s)224 modul_initmod_ul (modulusul_t m, const unsigned long s)
225 {
226 m[0] = s;
227 }
228
229
230 MAYBE_UNUSED
231 static inline void
modul_initmod_int(modulusul_t m,const modintul_t s)232 modul_initmod_int (modulusul_t m, const modintul_t s)
233 {
234 m[0] = s[0];
235 }
236
237
238 MAYBE_UNUSED
239 static inline unsigned long
modul_getmod_ul(const modulusul_t m)240 modul_getmod_ul (const modulusul_t m)
241 {
242 return m[0];
243 }
244
245
246 MAYBE_UNUSED
247 static inline void
modul_getmod_int(modintul_t r,const modulusul_t m)248 modul_getmod_int (modintul_t r, const modulusul_t m)
249 {
250 r[0] = m[0];
251 }
252
253
254 MAYBE_UNUSED
255 static inline void
modul_clearmod(modulusul_t m MAYBE_UNUSED)256 modul_clearmod (modulusul_t m MAYBE_UNUSED)
257 {
258 return;
259 }
260
261
262 /* Functions for residues */
263 static inline void
264 modul_neg (residueul_t, const residueul_t, const modulusul_t);
265
266 /* Initialises a residue_t type and sets it to zero */
267 MAYBE_UNUSED
268 static inline void
modul_init(residueul_t r,const modulusul_t m MAYBE_UNUSED)269 modul_init (residueul_t r, const modulusul_t m MAYBE_UNUSED)
270 {
271 r[0] = 0UL;
272 return;
273 }
274
275
276 /* Initialises a residue_t type, but does not set it to zero. For fixed length
277 residue_t types, that leaves nothing to do at all. */
278 MAYBE_UNUSED
279 static inline void
modul_init_noset0(residueul_t r MAYBE_UNUSED,const modulusul_t m MAYBE_UNUSED)280 modul_init_noset0 (residueul_t r MAYBE_UNUSED,
281 const modulusul_t m MAYBE_UNUSED)
282 {
283 return;
284 }
285
286
287 MAYBE_UNUSED
288 static inline void
modul_clear(residueul_t r MAYBE_UNUSED,const modulusul_t m MAYBE_UNUSED)289 modul_clear (residueul_t r MAYBE_UNUSED,
290 const modulusul_t m MAYBE_UNUSED)
291 {
292 return;
293 }
294
295
296 MAYBE_UNUSED
297 static inline void
modul_set(residueul_t r,const residueul_t s,const modulusul_t m MAYBE_UNUSED)298 modul_set (residueul_t r, const residueul_t s,
299 const modulusul_t m MAYBE_UNUSED)
300 {
301 ASSERT_EXPENSIVE (s[0] < m[0]);
302 r[0] = s[0];
303 }
304
305
306 MAYBE_UNUSED
307 static inline void
modul_set_ul(residueul_t r,const unsigned long s,const modulusul_t m)308 modul_set_ul (residueul_t r, const unsigned long s, const modulusul_t m)
309 {
310 r[0] = s % m[0];
311 }
312
313
314 /* Sets the residue_t to the class represented by the integer s. Assumes that
315 s is reduced (mod m), i.e. 0 <= s < m */
316
317 MAYBE_UNUSED
318 static inline void
modul_set_ul_reduced(residueul_t r,const unsigned long s,const modulusul_t m MAYBE_UNUSED)319 modul_set_ul_reduced (residueul_t r, const unsigned long s,
320 const modulusul_t m MAYBE_UNUSED)
321 {
322 ASSERT (s < m[0]);
323 r[0] = s;
324 }
325
326
327 MAYBE_UNUSED
328 static inline void
modul_set_int(residueul_t r,const modintul_t s,const modulusul_t m)329 modul_set_int (residueul_t r, const modintul_t s,
330 const modulusul_t m)
331 {
332 r[0] = s[0] % m[0];
333 }
334
335
336 MAYBE_UNUSED
337 static inline void
modul_set_int_reduced(residueul_t r,const modintul_t s,const modulusul_t m MAYBE_UNUSED)338 modul_set_int_reduced (residueul_t r, const modintul_t s,
339 const modulusul_t m MAYBE_UNUSED)
340 {
341 ASSERT (s[0] < m[0]);
342 r[0] = s[0];
343 }
344
345
346 MAYBE_UNUSED
347 static inline void
modul_set_uint64(residueul_t r,const uint64_t s,const modulusul_t m)348 modul_set_uint64 (residueul_t r, const uint64_t s, const modulusul_t m)
349 {
350 r[0] = s % m[0];
351 }
352
353
354 MAYBE_UNUSED
355 static inline void
modul_set_int64(residueul_t r,const int64_t s,const modulusul_t m)356 modul_set_int64 (residueul_t r, const int64_t s, const modulusul_t m)
357 {
358 r[0] = llabs(s) % m[0];
359 if (s < 0)
360 modul_neg(r, r, m);
361 }
362
363
364 /* These two are so trivial that we don't really require m in the
365 * interface. For 1 we might, as the internal representation might
366 * not use "1" for 1 (e.g. when using Montgomery's REDC.)
367 * For interface homogeneity we make even modul_set0 take the m parameter.
368 */
369 MAYBE_UNUSED
370 static inline void
modul_set0(residueul_t r,const modulusul_t m MAYBE_UNUSED)371 modul_set0 (residueul_t r, const modulusul_t m MAYBE_UNUSED)
372 {
373 r[0] = 0UL;
374 }
375
376
377 MAYBE_UNUSED
378 static inline void
modul_set1(residueul_t r,const modulusul_t m MAYBE_UNUSED)379 modul_set1 (residueul_t r, const modulusul_t m MAYBE_UNUSED)
380 {
381 r[0] = m[0] != 1UL;
382 }
383
384
385 /* Exchanges the values of the two arguments */
386
387 MAYBE_UNUSED
388 static inline void
modul_swap(residueul_t a,residueul_t b,const modulusul_t m MAYBE_UNUSED)389 modul_swap (residueul_t a, residueul_t b,
390 const modulusul_t m MAYBE_UNUSED)
391 {
392 unsigned long t;
393 ASSERT_EXPENSIVE (a[0] < m[0] && b[0] < m[0]);
394 t = a[0];
395 a[0] = b[0];
396 b[0] = t;
397 }
398
399
400 MAYBE_UNUSED
401 static inline unsigned long
modul_get_ul(const residueul_t s,const modulusul_t m MAYBE_UNUSED)402 modul_get_ul (const residueul_t s, const modulusul_t m MAYBE_UNUSED)
403 {
404 ASSERT_EXPENSIVE (s[0] < m[0]);
405 return s[0];
406 }
407
408
409 MAYBE_UNUSED
410 static inline void
modul_get_int(modintul_t r,const residueul_t s,const modulusul_t m MAYBE_UNUSED)411 modul_get_int (modintul_t r, const residueul_t s,
412 const modulusul_t m MAYBE_UNUSED)
413 {
414 ASSERT_EXPENSIVE (s[0] < m[0]);
415 r[0] = s[0];
416 }
417
418
419 MAYBE_UNUSED
420 static inline int
modul_equal(const residueul_t a,const residueul_t b,const modulusul_t m MAYBE_UNUSED)421 modul_equal (const residueul_t a, const residueul_t b,
422 const modulusul_t m MAYBE_UNUSED)
423 {
424 ASSERT_EXPENSIVE (a[0] < m[0] && b[0] < m[0]);
425 return (a[0] == b[0]);
426 }
427
428
429 MAYBE_UNUSED
430 static inline int
modul_is0(const residueul_t a,const modulusul_t m MAYBE_UNUSED)431 modul_is0 (const residueul_t a, const modulusul_t m MAYBE_UNUSED)
432 {
433 ASSERT_EXPENSIVE (a[0] < m[0]);
434 return (a[0] == 0UL);
435 }
436
437
438 MAYBE_UNUSED
439 static inline int
modul_is1(const residueul_t a,const modulusul_t m MAYBE_UNUSED)440 modul_is1 (const residueul_t a, const modulusul_t m MAYBE_UNUSED)
441 {
442 ASSERT_EXPENSIVE (a[0] < m[0]);
443 return (a[0] == 1UL);
444 }
445
446
447 MAYBE_UNUSED
448 static inline void
modul_add(residueul_t r,const residueul_t a,const residueul_t b,const modulusul_t m)449 modul_add (residueul_t r, const residueul_t a, const residueul_t b,
450 const modulusul_t m)
451 {
452 #ifdef MODTRACE
453 printf ("modul_add: a = %lu, b = %lu", a[0], b[0]);
454 #endif
455
456 ularith_addmod_ul_ul(r, a[0], b[0], m[0]);
457
458 #ifdef MODTRACE
459 printf (", r = %lu\n", r[0]);
460 #endif
461 }
462
463
464 MAYBE_UNUSED
465 static inline void
modul_add1(residueul_t r,const residueul_t a,const modulusul_t m)466 modul_add1 (residueul_t r, const residueul_t a, const modulusul_t m)
467 {
468 ASSERT_EXPENSIVE (a[0] < m[0]);
469 r[0] = a[0] + 1;
470 if (r[0] == m[0])
471 r[0] = 0;
472 }
473
474
475 MAYBE_UNUSED
476 static inline void
modul_add_ul(residueul_t r,const residueul_t a,const unsigned long b,const modulusul_t m)477 modul_add_ul (residueul_t r, const residueul_t a, const unsigned long b,
478 const modulusul_t m)
479 {
480 ularith_addmod_ul_ul(r, a[0], b % m[0], m[0]);
481 }
482
483 MAYBE_UNUSED
484 static inline void
modul_sub(residueul_t r,const residueul_t a,const residueul_t b,const modulusul_t m)485 modul_sub (residueul_t r, const residueul_t a, const residueul_t b,
486 const modulusul_t m)
487 {
488 #ifdef MODTRACE
489 printf ("submod_ul: a = %lu, b = %lu", a[0], b[0]);
490 #endif
491
492 ularith_submod_ul_ul(r, a[0], b[0], m[0]);
493
494 #ifdef MODTRACE
495 printf (", r = %lu\n", r[0]);
496 #endif
497 }
498
499
500 MAYBE_UNUSED
501 static inline void
modul_sub_ul(residueul_t r,const residueul_t a,const unsigned long b,const modulusul_t m)502 modul_sub_ul (residueul_t r, const residueul_t a, const unsigned long b,
503 const modulusul_t m)
504 {
505 ularith_submod_ul_ul(r, a[0], b % m[0], m[0]);
506 }
507
508
509 MAYBE_UNUSED
510 static inline void
modul_neg(residueul_t r,const residueul_t a,const modulusul_t m)511 modul_neg (residueul_t r, const residueul_t a, const modulusul_t m)
512 {
513 ASSERT_EXPENSIVE (a[0] < m[0]);
514 if (a[0] == 0UL)
515 r[0] = a[0];
516 else
517 r[0] = m[0] - a[0];
518 }
519
520
521 MAYBE_UNUSED
522 static inline void
modul_mul(residueul_t r,const residueul_t a,const residueul_t b,const modulusul_t m)523 modul_mul (residueul_t r, const residueul_t a, const residueul_t b,
524 const modulusul_t m)
525 {
526 unsigned long t1, t2;
527 #if defined(MODTRACE)
528 printf ("mulmod_ul: a = %lu, b = %lu", a[0], b[0]);
529 #endif
530
531 ASSERT_EXPENSIVE (a[0] < m[0] && b[0] < m[0]);
532 ularith_mul_ul_ul_2ul (&t1, &t2, a[0], b[0]);
533 ularith_div_2ul_ul_ul_r (r, t1, t2, m[0]);
534
535 #if defined(MODTRACE)
536 printf (", r = %lu\n", r);
537 #endif
538 }
539
540 MAYBE_UNUSED
541 static inline void
modul_sqr(residueul_t r,const residueul_t a,const modulusul_t m)542 modul_sqr (residueul_t r, const residueul_t a, const modulusul_t m)
543 {
544 unsigned long t1, t2;
545
546 ASSERT_EXPENSIVE (a[0] < m[0]);
547 ularith_mul_ul_ul_2ul (&t1, &t2, a[0], a[0]);
548 ularith_div_2ul_ul_ul_r (r, t1, t2, m[0]);
549 }
550
551 /* Computes (a * 2^wordsize) % m */
552 MAYBE_UNUSED
553 static inline void
modul_tomontgomery(residueul_t r,const residueul_t a,const modulusul_t m)554 modul_tomontgomery (residueul_t r, const residueul_t a, const modulusul_t m)
555 {
556 ASSERT_EXPENSIVE (a[0] < m[0]);
557 ularith_div_2ul_ul_ul_r (r, 0UL, a[0], m[0]);
558 }
559
560
561 /* Computes (a / 2^wordsize) % m */
562 MAYBE_UNUSED
563 static inline void
modul_frommontgomery(residueul_t r,const residueul_t a,const unsigned long invm,const modulusul_t m)564 modul_frommontgomery (residueul_t r, const residueul_t a,
565 const unsigned long invm, const modulusul_t m)
566 {
567 unsigned long t_low, t_high;
568 t_low = a[0] * invm;
569 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
570 r[0] = t_high + (a[0] != 0UL ? 1UL : 0UL);
571 }
572
573
574 /* Computes (a / 2^wordsize) % m, but result can be r = m.
575 Input a must not be equal 0 */
576 MAYBE_UNUSED
577 static inline void
modul_redcsemi_ul_not0(residueul_t r,const unsigned long a,const unsigned long invm,const modulusul_t m)578 modul_redcsemi_ul_not0 (residueul_t r, const unsigned long a,
579 const unsigned long invm, const modulusul_t m)
580 {
581 unsigned long t_low, t_high;
582
583 ASSERT (a != 0);
584
585 t_low = a * invm; /* t_low <= 2^w-1 */
586 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
587 /* t_high:t_low <= (2^w-1) * m */
588 r[0] = t_high + 1UL;
589 /* (t_high+1):t_low <= 2^w + (2^w-1) * m <= 2^w + 2^w*m - m
590 <= 2^w * (m + 1) - m */
591 /* r <= floor ((2^w * (m + 1) - m) / 2^w) <= floor((m + 1) - m/2^w)
592 <= m */
593 }
594
595
596 /* Computes ((a + b) / 2^wordsize) % m. a <= m is permissible */
597 MAYBE_UNUSED
598 static inline void
modul_addredc_ul(residueul_t r,const residueul_t a,const unsigned long b,const unsigned long invm,const modulusul_t m)599 modul_addredc_ul (residueul_t r, const residueul_t a, const unsigned long b,
600 const unsigned long invm, const modulusul_t m)
601 {
602 unsigned long s_low, s_high, t_low, t_high;
603
604 ASSERT_EXPENSIVE (a[0] <= m[0]);
605 s_low = b;
606 s_high = 0UL;
607 ularith_add_ul_2ul (&s_low, &s_high, a[0]);
608
609 t_low = s_low * invm;
610 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
611 ASSERT_EXPENSIVE (s_low + t_low == 0UL);
612 r[0] = t_high + s_high + (s_low != 0UL ? 1UL : 0UL);
613
614 /* r = ((a+b) + (((a+b)%2^w * invm) % 2^w) * m) / 2^w Use a<=m-1, b<=2^w-1
615 r <= (m + 2^w - 1 + (2^w - 1) * m) / 2^w
616 = (m - 1 + 2^w + m*2^w - m) / 2^w
617 = (- 1 + 2^w + m2^w) / 2^w
618 = m + 1 - 1/2^w
619 r <= m, since r is an integer
620 */
621 if (r[0] == m[0])
622 r[0] = 0UL;
623 }
624
625
626 /* Computes ((a + b) / 2^wordsize) % m, but result can be == m.
627 a <= m is permissible */
628 MAYBE_UNUSED
629 static inline void
modul_addredcsemi_ul(residueul_t r,const residueul_t a,const unsigned long b,const unsigned long invm,const modulusul_t m)630 modul_addredcsemi_ul (residueul_t r, const residueul_t a,
631 const unsigned long b, const unsigned long invm,
632 const modulusul_t m)
633 {
634 unsigned long s_low, s_high, t_low;
635
636 ASSERT_EXPENSIVE(a[0] <= m[0]);
637 s_low = b;
638 #ifdef HAVE_GCC_STYLE_AMD64_INLINE_ASM
639 {
640 unsigned char sb;
641 __asm__ __VOLATILE (
642 "addq %2, %0\n\t" /* cy * 2^w + s_low = a + b */
643 "setne %1\n\t" /* if (s_low != 0) sb = 1 */
644 "adcb $0, %1\n" /* sb += cy */
645 : "+&r" (s_low), "=qm" (sb)
646 : "rm" (a[0])
647 : "cc");
648 s_high = sb;
649 }
650 #elif defined(__i386__) && defined(__GNUC__)
651 {
652 unsigned char sb;
653 __asm__ __VOLATILE (
654 "addl %2, %0\n\t"
655 "setne %1\n\t"
656 "adcb $0, %1\n"
657 : "+&r" (s_low), "=qm" (sb)
658 : "rm" (a[0])
659 : "cc");
660 s_high = sb;
661 }
662 #else
663 s_high = 0UL;
664 ularith_add_ul_2ul (&s_low, &s_high, a[0]);
665 s_high += (s_low != 0UL ? 1UL : 0UL);
666 #endif
667
668 t_low = s_low * invm;
669 ularith_mul_ul_ul_2ul (&t_low, r, t_low, m[0]);
670 ASSERT_EXPENSIVE (s_low + t_low == 0UL);
671 r[0] += s_high;
672
673 /* r = ((a+b) + (((a+b)%2^w * invm) % 2^w) * m) / 2^w
674 r <= ((a+b) + (2^w - 1) * m) / 2^w
675 r <= (m + 2^w-1 + m*2^w - m) / 2^w
676 r <= (2^w -1 + p2^w) / 2^w
677 r <= p + 1 - 1/2^w
678 r <= p
679 */
680 }
681
682
683 MAYBE_UNUSED
684 static inline void
modul_mulredc(residueul_t r,const residueul_t a,const residueul_t b,const unsigned long invm,const modulusul_t m)685 modul_mulredc (residueul_t r, const residueul_t a, const residueul_t b,
686 const unsigned long invm, const modulusul_t m)
687 {
688 unsigned long p_low, p_high, t_low, t_high;
689
690 ASSERT_EXPENSIVE (m[0] % 2 != 0);
691 ASSERT_EXPENSIVE (a[0] < m[0] && b[0] < m[0]);
692 #if defined(MODTRACE)
693 printf ("(%lu * %lu / 2^%ld) %% %lu",
694 a[0], b[0], 8 * sizeof(unsigned long), m[0]);
695 #endif
696
697 ularith_mul_ul_ul_2ul (&p_low, &p_high, a[0], b[0]);
698 t_low = p_low * invm;
699 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
700 /* Let w = 2^wordsize. We know (p_high * w + p_low) + (t_high * w + t_low)
701 == 0 (mod w) so either p_low == t_low == 0, or p_low !=0 and t_low != 0.
702 In the former case we want p_high + t_high + 1, in the latter
703 p_high + t_high */
704 /* Since a <= p-1 and b <= p-1, and p <= w-1, a*b <= w^2 - 4*w + 4, so
705 adding 1 to p_high is safe */
706 #if 0
707 /* S_lower? */
708 ularith_add_ul_2ul (&p_low, &p_high, t_low);
709 #else
710 p_high += (p_low != 0UL ? 1UL : 0UL);
711 #endif
712
713 modul_add (r, &p_high, &t_high, m);
714
715 #if defined(MODTRACE)
716 printf (" == %lu /* PARI */ \n", r[0]);
717 #endif
718 }
719
720
721 /* FIXME: check for overflow if b > m */
722 MAYBE_UNUSED
723 static inline void
modul_mulredc_ul(residueul_t r,const residueul_t a,const unsigned long b,const unsigned long invm,const modulusul_t m)724 modul_mulredc_ul (residueul_t r, const residueul_t a, const unsigned long b,
725 const unsigned long invm, const modulusul_t m)
726 {
727 unsigned long p_low, p_high, t_low, t_high;
728 ASSERT_EXPENSIVE (m[0] % 2 != 0);
729 #if defined(MODTRACE)
730 printf ("(%lu * %lu / 2^%ld) %% %lu",
731 a[0], b, 8 * sizeof(unsigned long), m[0]);
732 #endif
733
734 ularith_mul_ul_ul_2ul (&p_low, &p_high, a[0], b);
735 t_low = p_low * invm;
736 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
737 p_high += (p_low != 0UL ? 1UL : 0UL);
738 r[0] = (p_high >= m[0] - t_high) ? (p_high - (m[0] - t_high)) : (p_high + t_high);
739
740 #if defined(MODTRACE)
741 printf (" == %lu /* PARI */ \n", r[0]);
742 #endif
743 }
744
745
746 /* Computes (a * b + c)/ 2^wordsize % m. Requires that
747 a * b + c < 2^wordsize * m */
748
749 MAYBE_UNUSED
750 static inline void
modul_muladdredc_ul(residueul_t r,const residueul_t a,const unsigned long b,const unsigned long c,const unsigned long invm,const modulusul_t m)751 modul_muladdredc_ul (residueul_t r, const residueul_t a, const unsigned long b,
752 const unsigned long c, const unsigned long invm,
753 const modulusul_t m)
754 {
755 unsigned long p_low, p_high, t_low, t_high;
756 ASSERT_EXPENSIVE (m[0] % 2 != 0);
757 #if defined(MODTRACE)
758 printf ("(%lu * %lu / 2^%ld) %% %lu",
759 a[0], b, 8 * sizeof(unsigned long), m[0]);
760 #endif
761
762 ularith_mul_ul_ul_2ul (&p_low, &p_high, a[0], b);
763 ularith_add_ul_2ul (&p_low, &p_high, c);
764 t_low = p_low * invm;
765 ularith_mul_ul_ul_2ul (&t_low, &t_high, t_low, m[0]);
766 p_high += (p_low != 0UL ? 1UL : 0UL);
767 r[0] = (p_high >= m[0] - t_high) ? (p_high - (m[0] - t_high)) : (p_high + t_high);
768
769 #if defined(MODTRACE)
770 printf (" == %lu /* PARI */ \n", r[0]);
771 #endif
772 }
773
774
775 MAYBE_UNUSED
776 static inline void
modul_div2(residueul_t r,const residueul_t a,const modulusul_t m)777 modul_div2 (residueul_t r, const residueul_t a, const modulusul_t m)
778 {
779 r[0] = ularith_div2mod(a[0], m[0]);
780 }
781
782
783 MAYBE_UNUSED
784 static inline int
modul_next(residueul_t r,const modulusul_t m)785 modul_next (residueul_t r, const modulusul_t m)
786 {
787 return (++r[0] == m[0]);
788 }
789
790
791 MAYBE_UNUSED
792 static inline int
modul_finished(const residueul_t r,const modulusul_t m)793 modul_finished (const residueul_t r, const modulusul_t m)
794 {
795 return (r[0] == m[0]);
796 }
797
798 #ifdef __cplusplus
799 extern "C" {
800 #endif
801
802 /* prototypes of non-inline functions */
803 int modul_div3 (residueul_t, const residueul_t, const modulusul_t);
804 int modul_div5 (residueul_t, const residueul_t, const modulusul_t);
805 int modul_div7 (residueul_t, const residueul_t, const modulusul_t);
806 int modul_div11 (residueul_t, const residueul_t, const modulusul_t);
807 int modul_div13 (residueul_t, const residueul_t, const modulusul_t);
808 void modul_gcd (modintul_t, const residueul_t, const modulusul_t);
809 void modul_pow_ul (residueul_t, const residueul_t, const unsigned long,
810 const modulusul_t);
811 void modul_2pow_ul (residueul_t, const unsigned long, const modulusul_t);
812 void modul_pow_mp (residueul_t, const residueul_t, const unsigned long *,
813 const int, const modulusul_t);
814 void modul_2pow_mp (residueul_t, const unsigned long *, const int,
815 const modulusul_t);
816 void modul_V_ul (residueul_t, const residueul_t, const unsigned long,
817 const modulusul_t);
818 void modul_V_mp (residueul_t, const residueul_t, const unsigned long *,
819 const int, const modulusul_t);
820 int modul_sprp (const residueul_t, const modulusul_t);
821 int modul_sprp2 (const modulusul_t);
822 int modul_isprime (const modulusul_t);
823 int modul_inv (residueul_t, const residueul_t, const modulusul_t);
824 int modul_inv_odd (residueul_t, const residueul_t, const modulusul_t);
825 int modul_inv_powerof2 (residueul_t, const residueul_t, const modulusul_t);
826 int modul_batchinv (residueul_t *, const residueul_t *, size_t,
827 const residueul_t, const modulusul_t);
828 int modul_jacobi (const residueul_t, const modulusul_t);
829
830 #ifdef __cplusplus
831 }
832 #endif
833
834 #endif /* MOD_UL_H */
835