1 | /* |
2 | * Copyright (c) 2008-2020 Stefan Krah. All rights reserved. |
3 | * |
4 | * Redistribution and use in source and binary forms, with or without |
5 | * modification, are permitted provided that the following conditions |
6 | * are met: |
7 | * |
8 | * 1. Redistributions of source code must retain the above copyright |
9 | * notice, this list of conditions and the following disclaimer. |
10 | * |
11 | * 2. Redistributions in binary form must reproduce the above copyright |
12 | * notice, this list of conditions and the following disclaimer in the |
13 | * documentation and/or other materials provided with the distribution. |
14 | * |
15 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND |
16 | * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
17 | * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
18 | * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE |
19 | * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
20 | * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS |
21 | * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
22 | * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT |
23 | * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY |
24 | * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF |
25 | * SUCH DAMAGE. |
26 | */ |
27 | |
28 | |
29 | #ifndef LIBMPDEC_UMODARITH_H_ |
30 | #define LIBMPDEC_UMODARITH_H_ |
31 | |
32 | |
33 | #include "mpdecimal.h" |
34 | |
35 | #include "constants.h" |
36 | #include "typearith.h" |
37 | |
38 | |
39 | /* Bignum: Low level routines for unsigned modular arithmetic. These are |
40 | used in the fast convolution functions for very large coefficients. */ |
41 | |
42 | |
43 | /**************************************************************************/ |
44 | /* ANSI modular arithmetic */ |
45 | /**************************************************************************/ |
46 | |
47 | |
48 | /* |
49 | * Restrictions: a < m and b < m |
50 | * ACL2 proof: umodarith.lisp: addmod-correct |
51 | */ |
52 | static inline mpd_uint_t |
53 | addmod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
54 | { |
55 | mpd_uint_t s; |
56 | |
57 | s = a + b; |
58 | s = (s < a) ? s - m : s; |
59 | s = (s >= m) ? s - m : s; |
60 | |
61 | return s; |
62 | } |
63 | |
64 | /* |
65 | * Restrictions: a < m and b < m |
66 | * ACL2 proof: umodarith.lisp: submod-2-correct |
67 | */ |
68 | static inline mpd_uint_t |
69 | submod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
70 | { |
71 | mpd_uint_t d; |
72 | |
73 | d = a - b; |
74 | d = (a < b) ? d + m : d; |
75 | |
76 | return d; |
77 | } |
78 | |
79 | /* |
80 | * Restrictions: a < 2m and b < 2m |
81 | * ACL2 proof: umodarith.lisp: section ext-submod |
82 | */ |
83 | static inline mpd_uint_t |
84 | ext_submod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
85 | { |
86 | mpd_uint_t d; |
87 | |
88 | a = (a >= m) ? a - m : a; |
89 | b = (b >= m) ? b - m : b; |
90 | |
91 | d = a - b; |
92 | d = (a < b) ? d + m : d; |
93 | |
94 | return d; |
95 | } |
96 | |
97 | /* |
98 | * Reduce double word modulo m. |
99 | * Restrictions: m != 0 |
100 | * ACL2 proof: umodarith.lisp: section dw-reduce |
101 | */ |
102 | static inline mpd_uint_t |
103 | dw_reduce(mpd_uint_t hi, mpd_uint_t lo, mpd_uint_t m) |
104 | { |
105 | mpd_uint_t r1, r2, w; |
106 | |
107 | _mpd_div_word(&w, &r1, hi, m); |
108 | _mpd_div_words(&w, &r2, r1, lo, m); |
109 | |
110 | return r2; |
111 | } |
112 | |
113 | /* |
114 | * Subtract double word from a. |
115 | * Restrictions: a < m |
116 | * ACL2 proof: umodarith.lisp: section dw-submod |
117 | */ |
118 | static inline mpd_uint_t |
119 | dw_submod(mpd_uint_t a, mpd_uint_t hi, mpd_uint_t lo, mpd_uint_t m) |
120 | { |
121 | mpd_uint_t d, r; |
122 | |
123 | r = dw_reduce(hi, lo, m); |
124 | d = a - r; |
125 | d = (a < r) ? d + m : d; |
126 | |
127 | return d; |
128 | } |
129 | |
130 | #ifdef CONFIG_64 |
131 | |
132 | /**************************************************************************/ |
133 | /* 64-bit modular arithmetic */ |
134 | /**************************************************************************/ |
135 | |
136 | /* |
137 | * A proof of the algorithm is in literature/mulmod-64.txt. An ACL2 |
138 | * proof is in umodarith.lisp: section "Fast modular reduction". |
139 | * |
140 | * Algorithm: calculate (a * b) % p: |
141 | * |
142 | * a) hi, lo <- a * b # Calculate a * b. |
143 | * |
144 | * b) hi, lo <- R(hi, lo) # Reduce modulo p. |
145 | * |
146 | * c) Repeat step b) until 0 <= hi * 2**64 + lo < 2*p. |
147 | * |
148 | * d) If the result is less than p, return lo. Otherwise return lo - p. |
149 | */ |
150 | |
151 | static inline mpd_uint_t |
152 | x64_mulmod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
153 | { |
154 | mpd_uint_t hi, lo, x, y; |
155 | |
156 | |
157 | _mpd_mul_words(&hi, &lo, a, b); |
158 | |
159 | if (m & (1ULL<<32)) { /* P1 */ |
160 | |
161 | /* first reduction */ |
162 | x = y = hi; |
163 | hi >>= 32; |
164 | |
165 | x = lo - x; |
166 | if (x > lo) hi--; |
167 | |
168 | y <<= 32; |
169 | lo = y + x; |
170 | if (lo < y) hi++; |
171 | |
172 | /* second reduction */ |
173 | x = y = hi; |
174 | hi >>= 32; |
175 | |
176 | x = lo - x; |
177 | if (x > lo) hi--; |
178 | |
179 | y <<= 32; |
180 | lo = y + x; |
181 | if (lo < y) hi++; |
182 | |
183 | return (hi || lo >= m ? lo - m : lo); |
184 | } |
185 | else if (m & (1ULL<<34)) { /* P2 */ |
186 | |
187 | /* first reduction */ |
188 | x = y = hi; |
189 | hi >>= 30; |
190 | |
191 | x = lo - x; |
192 | if (x > lo) hi--; |
193 | |
194 | y <<= 34; |
195 | lo = y + x; |
196 | if (lo < y) hi++; |
197 | |
198 | /* second reduction */ |
199 | x = y = hi; |
200 | hi >>= 30; |
201 | |
202 | x = lo - x; |
203 | if (x > lo) hi--; |
204 | |
205 | y <<= 34; |
206 | lo = y + x; |
207 | if (lo < y) hi++; |
208 | |
209 | /* third reduction */ |
210 | x = y = hi; |
211 | hi >>= 30; |
212 | |
213 | x = lo - x; |
214 | if (x > lo) hi--; |
215 | |
216 | y <<= 34; |
217 | lo = y + x; |
218 | if (lo < y) hi++; |
219 | |
220 | return (hi || lo >= m ? lo - m : lo); |
221 | } |
222 | else { /* P3 */ |
223 | |
224 | /* first reduction */ |
225 | x = y = hi; |
226 | hi >>= 24; |
227 | |
228 | x = lo - x; |
229 | if (x > lo) hi--; |
230 | |
231 | y <<= 40; |
232 | lo = y + x; |
233 | if (lo < y) hi++; |
234 | |
235 | /* second reduction */ |
236 | x = y = hi; |
237 | hi >>= 24; |
238 | |
239 | x = lo - x; |
240 | if (x > lo) hi--; |
241 | |
242 | y <<= 40; |
243 | lo = y + x; |
244 | if (lo < y) hi++; |
245 | |
246 | /* third reduction */ |
247 | x = y = hi; |
248 | hi >>= 24; |
249 | |
250 | x = lo - x; |
251 | if (x > lo) hi--; |
252 | |
253 | y <<= 40; |
254 | lo = y + x; |
255 | if (lo < y) hi++; |
256 | |
257 | return (hi || lo >= m ? lo - m : lo); |
258 | } |
259 | } |
260 | |
261 | static inline void |
262 | x64_mulmod2c(mpd_uint_t *a, mpd_uint_t *b, mpd_uint_t w, mpd_uint_t m) |
263 | { |
264 | *a = x64_mulmod(*a, w, m); |
265 | *b = x64_mulmod(*b, w, m); |
266 | } |
267 | |
268 | static inline void |
269 | x64_mulmod2(mpd_uint_t *a0, mpd_uint_t b0, mpd_uint_t *a1, mpd_uint_t b1, |
270 | mpd_uint_t m) |
271 | { |
272 | *a0 = x64_mulmod(*a0, b0, m); |
273 | *a1 = x64_mulmod(*a1, b1, m); |
274 | } |
275 | |
276 | static inline mpd_uint_t |
277 | x64_powmod(mpd_uint_t base, mpd_uint_t exp, mpd_uint_t umod) |
278 | { |
279 | mpd_uint_t r = 1; |
280 | |
281 | while (exp > 0) { |
282 | if (exp & 1) |
283 | r = x64_mulmod(r, base, umod); |
284 | base = x64_mulmod(base, base, umod); |
285 | exp >>= 1; |
286 | } |
287 | |
288 | return r; |
289 | } |
290 | |
291 | /* END CONFIG_64 */ |
292 | #else /* CONFIG_32 */ |
293 | |
294 | |
295 | /**************************************************************************/ |
296 | /* 32-bit modular arithmetic */ |
297 | /**************************************************************************/ |
298 | |
299 | #if defined(ANSI) |
300 | #if !defined(LEGACY_COMPILER) |
301 | /* HAVE_UINT64_T */ |
302 | static inline mpd_uint_t |
303 | std_mulmod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
304 | { |
305 | return ((mpd_uuint_t) a * b) % m; |
306 | } |
307 | |
308 | static inline void |
309 | std_mulmod2c(mpd_uint_t *a, mpd_uint_t *b, mpd_uint_t w, mpd_uint_t m) |
310 | { |
311 | *a = ((mpd_uuint_t) *a * w) % m; |
312 | *b = ((mpd_uuint_t) *b * w) % m; |
313 | } |
314 | |
315 | static inline void |
316 | std_mulmod2(mpd_uint_t *a0, mpd_uint_t b0, mpd_uint_t *a1, mpd_uint_t b1, |
317 | mpd_uint_t m) |
318 | { |
319 | *a0 = ((mpd_uuint_t) *a0 * b0) % m; |
320 | *a1 = ((mpd_uuint_t) *a1 * b1) % m; |
321 | } |
322 | /* END HAVE_UINT64_T */ |
323 | #else |
324 | /* LEGACY_COMPILER */ |
325 | static inline mpd_uint_t |
326 | std_mulmod(mpd_uint_t a, mpd_uint_t b, mpd_uint_t m) |
327 | { |
328 | mpd_uint_t hi, lo, q, r; |
329 | _mpd_mul_words(&hi, &lo, a, b); |
330 | _mpd_div_words(&q, &r, hi, lo, m); |
331 | return r; |
332 | } |
333 | |
334 | static inline void |
335 | std_mulmod2c(mpd_uint_t *a, mpd_uint_t *b, mpd_uint_t w, mpd_uint_t m) |
336 | { |
337 | *a = std_mulmod(*a, w, m); |
338 | *b = std_mulmod(*b, w, m); |
339 | } |
340 | |
341 | static inline void |
342 | std_mulmod2(mpd_uint_t *a0, mpd_uint_t b0, mpd_uint_t *a1, mpd_uint_t b1, |
343 | mpd_uint_t m) |
344 | { |
345 | *a0 = std_mulmod(*a0, b0, m); |
346 | *a1 = std_mulmod(*a1, b1, m); |
347 | } |
348 | /* END LEGACY_COMPILER */ |
349 | #endif |
350 | |
351 | static inline mpd_uint_t |
352 | std_powmod(mpd_uint_t base, mpd_uint_t exp, mpd_uint_t umod) |
353 | { |
354 | mpd_uint_t r = 1; |
355 | |
356 | while (exp > 0) { |
357 | if (exp & 1) |
358 | r = std_mulmod(r, base, umod); |
359 | base = std_mulmod(base, base, umod); |
360 | exp >>= 1; |
361 | } |
362 | |
363 | return r; |
364 | } |
365 | #endif /* ANSI CONFIG_32 */ |
366 | |
367 | |
368 | /**************************************************************************/ |
369 | /* Pentium Pro modular arithmetic */ |
370 | /**************************************************************************/ |
371 | |
372 | /* |
373 | * A proof of the algorithm is in literature/mulmod-ppro.txt. The FPU |
374 | * control word must be set to 64-bit precision and truncation mode |
375 | * prior to using these functions. |
376 | * |
377 | * Algorithm: calculate (a * b) % p: |
378 | * |
379 | * p := prime < 2**31 |
380 | * pinv := (long double)1.0 / p (precalculated) |
381 | * |
382 | * a) n = a * b # Calculate exact product. |
383 | * b) qest = n * pinv # Calculate estimate for q = n / p. |
384 | * c) q = (qest+2**63)-2**63 # Truncate qest to the exact quotient. |
385 | * d) r = n - q * p # Calculate remainder. |
386 | * |
387 | * Remarks: |
388 | * |
389 | * - p = dmod and pinv = dinvmod. |
390 | * - dinvmod points to an array of three uint32_t, which is interpreted |
391 | * as an 80 bit long double by fldt. |
392 | * - Intel compilers prior to version 11 do not seem to handle the |
393 | * __GNUC__ inline assembly correctly. |
394 | * - random tests are provided in tests/extended/ppro_mulmod.c |
395 | */ |
396 | |
397 | #if defined(PPRO) |
398 | #if defined(ASM) |
399 | |
400 | /* Return (a * b) % dmod */ |
401 | static inline mpd_uint_t |
402 | ppro_mulmod(mpd_uint_t a, mpd_uint_t b, double *dmod, uint32_t *dinvmod) |
403 | { |
404 | mpd_uint_t retval; |
405 | |
406 | __asm__ ( |
407 | "fildl %2\n\t" |
408 | "fildl %1\n\t" |
409 | "fmulp %%st, %%st(1)\n\t" |
410 | "fldt (%4)\n\t" |
411 | "fmul %%st(1), %%st\n\t" |
412 | "flds %5\n\t" |
413 | "fadd %%st, %%st(1)\n\t" |
414 | "fsubrp %%st, %%st(1)\n\t" |
415 | "fldl (%3)\n\t" |
416 | "fmulp %%st, %%st(1)\n\t" |
417 | "fsubrp %%st, %%st(1)\n\t" |
418 | "fistpl %0\n\t" |
419 | : "=m" (retval) |
420 | : "m" (a), "m" (b), "r" (dmod), "r" (dinvmod), "m" (MPD_TWO63) |
421 | : "st" , "memory" |
422 | ); |
423 | |
424 | return retval; |
425 | } |
426 | |
427 | /* |
428 | * Two modular multiplications in parallel: |
429 | * *a0 = (*a0 * w) % dmod |
430 | * *a1 = (*a1 * w) % dmod |
431 | */ |
432 | static inline void |
433 | ppro_mulmod2c(mpd_uint_t *a0, mpd_uint_t *a1, mpd_uint_t w, |
434 | double *dmod, uint32_t *dinvmod) |
435 | { |
436 | __asm__ ( |
437 | "fildl %2\n\t" |
438 | "fildl (%1)\n\t" |
439 | "fmul %%st(1), %%st\n\t" |
440 | "fxch %%st(1)\n\t" |
441 | "fildl (%0)\n\t" |
442 | "fmulp %%st, %%st(1) \n\t" |
443 | "fldt (%4)\n\t" |
444 | "flds %5\n\t" |
445 | "fld %%st(2)\n\t" |
446 | "fmul %%st(2)\n\t" |
447 | "fadd %%st(1)\n\t" |
448 | "fsub %%st(1)\n\t" |
449 | "fmull (%3)\n\t" |
450 | "fsubrp %%st, %%st(3)\n\t" |
451 | "fxch %%st(2)\n\t" |
452 | "fistpl (%0)\n\t" |
453 | "fmul %%st(2)\n\t" |
454 | "fadd %%st(1)\n\t" |
455 | "fsubp %%st, %%st(1)\n\t" |
456 | "fmull (%3)\n\t" |
457 | "fsubrp %%st, %%st(1)\n\t" |
458 | "fistpl (%1)\n\t" |
459 | : : "r" (a0), "r" (a1), "m" (w), |
460 | "r" (dmod), "r" (dinvmod), |
461 | "m" (MPD_TWO63) |
462 | : "st" , "memory" |
463 | ); |
464 | } |
465 | |
466 | /* |
467 | * Two modular multiplications in parallel: |
468 | * *a0 = (*a0 * b0) % dmod |
469 | * *a1 = (*a1 * b1) % dmod |
470 | */ |
471 | static inline void |
472 | ppro_mulmod2(mpd_uint_t *a0, mpd_uint_t b0, mpd_uint_t *a1, mpd_uint_t b1, |
473 | double *dmod, uint32_t *dinvmod) |
474 | { |
475 | __asm__ ( |
476 | "fildl %3\n\t" |
477 | "fildl (%2)\n\t" |
478 | "fmulp %%st, %%st(1)\n\t" |
479 | "fildl %1\n\t" |
480 | "fildl (%0)\n\t" |
481 | "fmulp %%st, %%st(1)\n\t" |
482 | "fldt (%5)\n\t" |
483 | "fld %%st(2)\n\t" |
484 | "fmul %%st(1), %%st\n\t" |
485 | "fxch %%st(1)\n\t" |
486 | "fmul %%st(2), %%st\n\t" |
487 | "flds %6\n\t" |
488 | "fldl (%4)\n\t" |
489 | "fxch %%st(3)\n\t" |
490 | "fadd %%st(1), %%st\n\t" |
491 | "fxch %%st(2)\n\t" |
492 | "fadd %%st(1), %%st\n\t" |
493 | "fxch %%st(2)\n\t" |
494 | "fsub %%st(1), %%st\n\t" |
495 | "fxch %%st(2)\n\t" |
496 | "fsubp %%st, %%st(1)\n\t" |
497 | "fxch %%st(1)\n\t" |
498 | "fmul %%st(2), %%st\n\t" |
499 | "fxch %%st(1)\n\t" |
500 | "fmulp %%st, %%st(2)\n\t" |
501 | "fsubrp %%st, %%st(3)\n\t" |
502 | "fsubrp %%st, %%st(1)\n\t" |
503 | "fxch %%st(1)\n\t" |
504 | "fistpl (%2)\n\t" |
505 | "fistpl (%0)\n\t" |
506 | : : "r" (a0), "m" (b0), "r" (a1), "m" (b1), |
507 | "r" (dmod), "r" (dinvmod), |
508 | "m" (MPD_TWO63) |
509 | : "st" , "memory" |
510 | ); |
511 | } |
512 | /* END PPRO GCC ASM */ |
513 | #elif defined(MASM) |
514 | |
515 | /* Return (a * b) % dmod */ |
516 | static inline mpd_uint_t __cdecl |
517 | ppro_mulmod(mpd_uint_t a, mpd_uint_t b, double *dmod, uint32_t *dinvmod) |
518 | { |
519 | mpd_uint_t retval; |
520 | |
521 | __asm { |
522 | mov eax, dinvmod |
523 | mov edx, dmod |
524 | fild b |
525 | fild a |
526 | fmulp st(1), st |
527 | fld TBYTE PTR [eax] |
528 | fmul st, st(1) |
529 | fld MPD_TWO63 |
530 | fadd st(1), st |
531 | fsubp st(1), st |
532 | fld QWORD PTR [edx] |
533 | fmulp st(1), st |
534 | fsubp st(1), st |
535 | fistp retval |
536 | } |
537 | |
538 | return retval; |
539 | } |
540 | |
541 | /* |
542 | * Two modular multiplications in parallel: |
543 | * *a0 = (*a0 * w) % dmod |
544 | * *a1 = (*a1 * w) % dmod |
545 | */ |
546 | static inline mpd_uint_t __cdecl |
547 | ppro_mulmod2c(mpd_uint_t *a0, mpd_uint_t *a1, mpd_uint_t w, |
548 | double *dmod, uint32_t *dinvmod) |
549 | { |
550 | __asm { |
551 | mov ecx, dmod |
552 | mov edx, a1 |
553 | mov ebx, dinvmod |
554 | mov eax, a0 |
555 | fild w |
556 | fild DWORD PTR [edx] |
557 | fmul st, st(1) |
558 | fxch st(1) |
559 | fild DWORD PTR [eax] |
560 | fmulp st(1), st |
561 | fld TBYTE PTR [ebx] |
562 | fld MPD_TWO63 |
563 | fld st(2) |
564 | fmul st, st(2) |
565 | fadd st, st(1) |
566 | fsub st, st(1) |
567 | fmul QWORD PTR [ecx] |
568 | fsubp st(3), st |
569 | fxch st(2) |
570 | fistp DWORD PTR [eax] |
571 | fmul st, st(2) |
572 | fadd st, st(1) |
573 | fsubrp st(1), st |
574 | fmul QWORD PTR [ecx] |
575 | fsubp st(1), st |
576 | fistp DWORD PTR [edx] |
577 | } |
578 | } |
579 | |
580 | /* |
581 | * Two modular multiplications in parallel: |
582 | * *a0 = (*a0 * b0) % dmod |
583 | * *a1 = (*a1 * b1) % dmod |
584 | */ |
585 | static inline void __cdecl |
586 | ppro_mulmod2(mpd_uint_t *a0, mpd_uint_t b0, mpd_uint_t *a1, mpd_uint_t b1, |
587 | double *dmod, uint32_t *dinvmod) |
588 | { |
589 | __asm { |
590 | mov ecx, dmod |
591 | mov edx, a1 |
592 | mov ebx, dinvmod |
593 | mov eax, a0 |
594 | fild b1 |
595 | fild DWORD PTR [edx] |
596 | fmulp st(1), st |
597 | fild b0 |
598 | fild DWORD PTR [eax] |
599 | fmulp st(1), st |
600 | fld TBYTE PTR [ebx] |
601 | fld st(2) |
602 | fmul st, st(1) |
603 | fxch st(1) |
604 | fmul st, st(2) |
605 | fld DWORD PTR MPD_TWO63 |
606 | fld QWORD PTR [ecx] |
607 | fxch st(3) |
608 | fadd st, st(1) |
609 | fxch st(2) |
610 | fadd st, st(1) |
611 | fxch st(2) |
612 | fsub st, st(1) |
613 | fxch st(2) |
614 | fsubrp st(1), st |
615 | fxch st(1) |
616 | fmul st, st(2) |
617 | fxch st(1) |
618 | fmulp st(2), st |
619 | fsubp st(3), st |
620 | fsubp st(1), st |
621 | fxch st(1) |
622 | fistp DWORD PTR [edx] |
623 | fistp DWORD PTR [eax] |
624 | } |
625 | } |
626 | #endif /* PPRO MASM (_MSC_VER) */ |
627 | |
628 | |
629 | /* Return (base ** exp) % dmod */ |
630 | static inline mpd_uint_t |
631 | ppro_powmod(mpd_uint_t base, mpd_uint_t exp, double *dmod, uint32_t *dinvmod) |
632 | { |
633 | mpd_uint_t r = 1; |
634 | |
635 | while (exp > 0) { |
636 | if (exp & 1) |
637 | r = ppro_mulmod(r, base, dmod, dinvmod); |
638 | base = ppro_mulmod(base, base, dmod, dinvmod); |
639 | exp >>= 1; |
640 | } |
641 | |
642 | return r; |
643 | } |
644 | #endif /* PPRO */ |
645 | #endif /* CONFIG_32 */ |
646 | |
647 | |
648 | #endif /* LIBMPDEC_UMODARITH_H_ */ |
649 | |