chiark / gitweb /
math/fgoldi.c: Add support for Hamburg's `Goldilocks' field.
[catacomb] / math / mpx.h
CommitLineData
d03ab969 1/* -*-c-*-
d03ab969 2 *
3 * Low level multiprecision arithmetic
4 *
5 * (c) 1999 Straylight/Edgeware
6 */
7
45c0fd36 8/*----- Licensing notice --------------------------------------------------*
d03ab969 9 *
10 * This file is part of Catacomb.
11 *
12 * Catacomb is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU Library General Public License as
14 * published by the Free Software Foundation; either version 2 of the
15 * License, or (at your option) any later version.
45c0fd36 16 *
d03ab969 17 * Catacomb is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Library General Public License for more details.
45c0fd36 21 *
d03ab969 22 * You should have received a copy of the GNU Library General Public
23 * License along with Catacomb; if not, write to the Free
24 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
25 * MA 02111-1307, USA.
26 */
27
a86e33af 28#ifndef CATACOMB_MPX_H
29#define CATACOMB_MPX_H
d03ab969 30
31#ifdef __cplusplus
32 extern "C" {
33#endif
34
35/*----- The idea ----------------------------------------------------------*
36 *
37 * This file provides functions and macros which work on vectors of words as
38 * unsigned multiprecision integers. The interface works in terms of base
39 * and limit pointers (i.e., a pointer to the start of a vector, and a
40 * pointer just past its end) rather than base pointer and length, because
41 * that requires more arithmetic and state to work on.
42 *
43 * The interfaces are slightly bizarre in other ways. Try to use the
44 * higher-level functions where you can: they're rather better designed to
45 * actually be friendly and useful.
46 */
47
48/*----- Header files ------------------------------------------------------*/
49
50#include <string.h>
51
a86e33af 52#ifndef CATACOMB_MPW_H
3c9ede17 53# include "mpw.h"
d03ab969 54#endif
55
56/*----- General manipulation ----------------------------------------------*/
57
58/* --- @MPX_SHRINK@ --- *
59 *
60 * Arguments: @const mpw *v@ = pointer to vector of words
61 * @const mpw *vl@ = (updated) current limit of vector
62 *
63 * Use: Shrinks down the limit of a multiprecision integer vector.
64 */
65
66#define MPX_SHRINK(v, vl) do { \
3c9ede17 67 const mpw *_vv = (v), *_vvl = (vl); \
68 while (_vvl > _vv && !_vvl[-1]) \
69 _vvl--; \
70 (vl) = (mpw *)_vvl; \
d03ab969 71} while (0)
72
73/* --- @MPX_BITS@ --- *
74 *
75 * Arguments: @unsigned long b@ = result variable
76 * @const mpw *v@ = pointer to array of words
77 * @const mpw *vl@ = limit of vector (from @MPX_SHRINK@)
78 *
79 * Use: Calculates the number of bits in a multiprecision value.
80 */
81
82#define MPX_BITS(b, v, vl) do { \
83 const mpw *_v = (v), *_vl = (vl); \
3c9ede17 84 MPX_SHRINK(_v, _vl); \
d03ab969 85 if (_v == _vl) \
86 (b) = 0; \
87 else { \
45c0fd36 88 unsigned long _b = MPW_BITS * (_vl - _v - 1) + 1; \
d03ab969 89 mpw _w = _vl[-1]; \
c29970a7 90 unsigned _k = MPW_P2; \
d03ab969 91 while (_k) { \
92 if (_w >> _k) { \
93 _w >>= _k; \
94 _b += _k; \
95 } \
96 _k >>= 1; \
97 } \
98 (b) = _b; \
99 } \
100} while (0)
101
102/* --- @MPX_OCTETS@ --- *
103 *
104 * Arguments: @size_t o@ = result variable
3c9ede17 105 * @const mpw *v, *vl@ = pointer to array of words
d03ab969 106 *
107 * Use: Calculates the number of octets in a multiprecision value.
108 */
109
3c9ede17 110#define MPX_OCTETS(o, v, vl) do { \
f09e814a 111 unsigned long _bb; \
112 MPX_BITS(_bb, (v), (vl)); \
113 (o) = (_bb + 7) >> 3; \
114} while (0)
115
116/* --- @MPX_OCTETS2C@ --- *
117 *
118 * Arguments: @size_t o@ = result variable
119 * @const mpw *v, *vl@ = pointer to array of words
120 *
121 * Use: Calculates the number of octets in a multiprecision value, if
122 * you represent it as two's complement.
123 */
124
125#define MPX_OCTETS2C(o, v, vl) do { \
126 unsigned long _bb; \
127 MPX_BITS(_bb, (v), (vl)); \
128 (o) = (_bb >> 3) + 1; \
d03ab969 129} while (0)
130
131/* --- @MPX_COPY@ --- *
132 *
133 * Arguments: @dv, dvl@ = destination vector base and limit
134 * @av, avl@ = source vector base and limit
135 *
136 * Use: Copies a multiprecision integer.
137 */
138
3c9ede17 139#define MPX_COPY(dv, dvl, av, avl) do { \
140 mpw *_dv = (dv), *_dvl = (dvl); \
141 size_t _dn = _dvl - _dv; \
142 const mpw *_av = (av), *_avl = (avl); \
143 size_t _an = _avl - _av; \
d03ab969 144 if (_av == _dv) { \
145 if (_dvl > _avl) \
3c9ede17 146 memset(_dv, 0, MPWS(_dn - _an)); \
d03ab969 147 } else if (_an >= _dn) \
148 memmove(_dv, _av, MPWS(_dn)); \
149 else { \
150 memmove(_dv, _av, MPWS(_an)); \
151 memset(_dv + _an, 0, MPWS(_dn - _an)); \
152 } \
153} while (0)
154
155/* --- @MPX_ZERO@ --- *
156 *
157 * Arguments: @v, vl@ = base and limit of vector to clear
158 *
159 * Use: Zeroes the area between the two vector pointers.
160 */
161
3c9ede17 162#define MPX_ZERO(v, vl) do { \
d03ab969 163 mpw *_v = (v), *_vl = (vl); \
3c9ede17 164 if (_v < _vl) \
165 memset(_v, 0, MPWS(_vl - _v)); \
d03ab969 166} while (0)
167
81578196 168/* --- @MPX_ONE@ --- *
169 *
170 * Arguments: @v, vl@ = base and limit of vector to clear
171 *
172 * Use: Fills the area between the two vector pointers with ones.
173 */
174
175#define MPX_ONE(v, vl) do { \
176 mpw * _v = (v); \
177 const mpw *_vl = (vl); \
178 while (_v < _vl) \
179 *_v++ = MPW_MAX; \
180} while (0)
181
d03ab969 182/*----- Loading and storing -----------------------------------------------*/
183
184/* --- @mpx_storel@ --- *
185 *
186 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
3c9ede17 187 * @void *p@ = pointer to octet array
d03ab969 188 * @size_t sz@ = size of octet array
189 *
190 * Returns: ---
191 *
192 * Use: Stores an MP in an octet array, least significant octet
193 * first. High-end octets are silently discarded if there
194 * isn't enough space for them.
195 */
196
197extern void mpx_storel(const mpw */*v*/, const mpw */*vl*/,
3c9ede17 198 void */*p*/, size_t /*sz*/);
d03ab969 199
200/* --- @mpx_loadl@ --- *
201 *
202 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
3c9ede17 203 * @const void *p@ = pointer to octet array
d03ab969 204 * @size_t sz@ = size of octet array
205 *
206 * Returns: ---
207 *
208 * Use: Loads an MP in an octet array, least significant octet
209 * first. High-end octets are ignored if there isn't enough
210 * space for them.
211 */
212
213extern void mpx_loadl(mpw */*v*/, mpw */*vl*/,
3c9ede17 214 const void */*p*/, size_t /*sz*/);
d03ab969 215
216/* --- @mpx_storeb@ --- *
217 *
218 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
3c9ede17 219 * @void *p@ = pointer to octet array
d03ab969 220 * @size_t sz@ = size of octet array
221 *
222 * Returns: ---
223 *
224 * Use: Stores an MP in an octet array, most significant octet
225 * first. High-end octets are silently discarded if there
226 * isn't enough space for them.
227 */
228
229extern void mpx_storeb(const mpw */*v*/, const mpw */*vl*/,
3c9ede17 230 void */*p*/, size_t /*sz*/);
d03ab969 231
232/* --- @mpx_loadb@ --- *
233 *
234 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
3c9ede17 235 * @const void *p@ = pointer to octet array
d03ab969 236 * @size_t sz@ = size of octet array
237 *
238 * Returns: ---
239 *
240 * Use: Loads an MP in an octet array, most significant octet
241 * first. High-end octets are ignored if there isn't enough
242 * space for them.
243 */
244
245extern void mpx_loadb(mpw */*v*/, mpw */*vl*/,
3c9ede17 246 const void */*p*/, size_t /*sz*/);
d03ab969 247
f09e814a 248/* --- @mpx_storel2cn@ --- *
249 *
250 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
251 * @void *pp@ = pointer to octet array
252 * @size_t sz@ = size of octet array
253 *
254 * Returns: ---
255 *
256 * Use: Stores a negative MP in an octet array, least significant
257 * octet first, as two's complement. High-end octets are
258 * silently discarded if there isn't enough space for them.
259 * This obviously makes the output bad.
260 */
261
262extern void mpx_storel2cn(const mpw */*v*/, const mpw */*vl*/,
263 void */*p*/, size_t /*sz*/);
264
265/* --- @mpx_loadl2cn@ --- *
266 *
267 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
268 * @const void *pp@ = pointer to octet array
269 * @size_t sz@ = size of octet array
270 *
271 * Returns: ---
272 *
273 * Use: Loads a negative MP in an octet array, least significant
274 * octet first, as two's complement. High-end octets are
275 * ignored if there isn't enough space for them. This probably
276 * means you made the wrong choice coming here.
277 */
278
279extern void mpx_loadl2cn(mpw */*v*/, mpw */*vl*/,
280 const void */*p*/, size_t /*sz*/);
281
282/* --- @mpx_storeb2cn@ --- *
283 *
284 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
285 * @void *pp@ = pointer to octet array
286 * @size_t sz@ = size of octet array
287 *
288 * Returns: ---
289 *
290 * Use: Stores a negative MP in an octet array, most significant
291 * octet first, as two's complement. High-end octets are
292 * silently discarded if there isn't enough space for them,
293 * which probably isn't what you meant.
294 */
295
296extern void mpx_storeb2cn(const mpw */*v*/, const mpw */*vl*/,
297 void */*p*/, size_t /*sz*/);
298
299/* --- @mpx_loadb2cn@ --- *
300 *
301 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
302 * @const void *pp@ = pointer to octet array
303 * @size_t sz@ = size of octet array
304 *
305 * Returns: ---
306 *
307 * Use: Loads a negative MP in an octet array, most significant octet
308 * first as two's complement. High-end octets are ignored if
309 * there isn't enough space for them. This probably means you
310 * chose this function wrongly.
311 */
312
313extern void mpx_loadb2cn(mpw */*v*/, mpw */*vl*/,
314 const void */*p*/, size_t /*sz*/);
315
316
d03ab969 317/*----- Logical shifting --------------------------------------------------*/
318
319/* --- @mpx_lsl@ --- *
320 *
321 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
322 * @const mpw *av, *avl@ = source vector base and limit
323 * @size_t n@ = number of bit positions to shift by
324 *
325 * Returns: ---
326 *
327 * Use: Performs a logical shift left operation on an integer.
328 */
329
330extern void mpx_lsl(mpw */*dv*/, mpw */*dvl*/,
331 const mpw */*av*/, const mpw */*avl*/,
332 size_t /*n*/);
333
81578196 334/* --- @mpx_lslc@ --- *
335 *
336 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
337 * @const mpw *av, *avl@ = source vector base and limit
338 * @size_t n@ = number of bit positions to shift by
339 *
340 * Returns: ---
341 *
342 * Use: Performs a logical shift left operation on an integer, only
343 * it fills in the bits with ones instead of zeroes.
344 */
345
346extern void mpx_lslc(mpw */*dv*/, mpw */*dvl*/,
347 const mpw */*av*/, const mpw */*avl*/,
348 size_t /*n*/);
349
d03ab969 350/* --- @mpx_lsr@ --- *
351 *
352 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
353 * @const mpw *av, *avl@ = source vector base and limit
354 * @size_t n@ = number of bit positions to shift by
355 *
356 * Returns: ---
357 *
358 * Use: Performs a logical shift right operation on an integer.
359 */
360
361extern void mpx_lsr(mpw */*dv*/, mpw */*dvl*/,
362 const mpw */*av*/, const mpw */*avl*/,
363 size_t /*n*/);
364
0f32e0f8 365/*----- Bitwise operations ------------------------------------------------*/
366
f09e814a 367/* --- @mpx_bitop@ --- *
0f32e0f8 368 *
369 * Arguments: @mpw *dv, *dvl@ = destination vector
370 * @const mpw *av, *avl@ = first source vector
371 * @const mpw *bv, *bvl@ = second source vector
372 *
373 * Returns: ---
374 *
f09e814a 375 * Use: Provide the dyadic boolean functions. The functions are
376 * named after the truth table they generate:
377 *
378 * a: 0011
379 * b: 0101
380 * @mpx_bitXXXX@
0f32e0f8 381 */
382
f09e814a 383#define MPX_DOBIN(what) \
384 what(0000) what(0001) what(0010) what(0011) \
385 what(0100) what(0101) what(0110) what(0111) \
386 what(1000) what(1001) what(1010) what(1011) \
387 what(1100) what(1101) what(1110) what(1111)
0f32e0f8 388
f09e814a 389#define MPX_BITDECL(string) \
390 extern void mpx_bit##string(mpw */*dv*/, mpw */*dvl*/, \
391 const mpw */*av*/, const mpw */*avl*/, \
392 const mpw */*bv*/, const mpw */*bvl*/);
393MPX_DOBIN(MPX_BITDECL)
0f32e0f8 394
f09e814a 395/* --- @mpx_[n]and@, @mpx_[n]or@, @mpx_xor@ --- *
396 *
397 * Synonyms for the commonly-used functions above.
398 */
399
45c0fd36
MW
400#define mpx_and mpx_bit0001
401#define mpx_or mpx_bit0111
f09e814a 402#define mpx_nand mpx_bit1110
45c0fd36
MW
403#define mpx_nor mpx_bit1000
404#define mpx_xor mpx_bit0110
f09e814a 405
406/* --- @mpx_not@ --- *
407 *
408 * Arguments: @mpw *dv, *dvl@ = destination vector
409 * @const mpw *av, *avl@ = first source vector
410 *
411 * Returns: ---
412 *
413 * Use; Bitwise NOT.
414 */
0f32e0f8 415
416extern void mpx_not(mpw */*dv*/, mpw */*dvl*/,
417 const mpw */*av*/, const mpw */*avl*/);
418
d03ab969 419/*----- Unsigned arithmetic -----------------------------------------------*/
420
7c13f461 421/* --- @mpx_2c@ --- *
422 *
423 * Arguments: @mpw *dv, *dvl@ = destination vector
424 * @const mpw *v, *vl@ = source vector
425 *
426 * Returns: ---
427 *
428 * Use: Calculates the two's complement of @v@.
429 */
430
431extern void mpx_2c(mpw */*dv*/, mpw */*dvl*/,
432 const mpw */*v*/, const mpw */*vl*/);
433
1a05a8ef 434/* --- @mpx_ueq@ --- *
435 *
436 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
437 * @const mpw *bv, *bvl@ = second argument vector base and limit
438 *
439 * Returns: Nonzero if the two vectors are equal.
440 *
441 * Use: Performs an unsigned integer test for equality.
442 */
443
444extern int mpx_ueq(const mpw */*av*/, const mpw */*avl*/,
445 const mpw */*bv*/, const mpw */*bvl*/);
446
d03ab969 447/* --- @mpx_ucmp@ --- *
448 *
449 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
450 * @const mpw *bv, *bvl@ = second argument vector base and limit
451 *
452 * Returns: Less than, equal to, or greater than zero depending on
453 * whether @a@ is less than, equal to or greater than @b@,
454 * respectively.
455 *
456 * Use: Performs an unsigned integer comparison.
457 */
458
459#define MPX_UCMP(av, avl, op, dv, dvl) \
460 (mpx_ucmp((av), (avl), (dv), (dvl)) op 0)
461
462extern int mpx_ucmp(const mpw */*av*/, const mpw */*avl*/,
463 const mpw */*bv*/, const mpw */*bvl*/);
464
465/* --- @mpx_uadd@ --- *
466 *
467 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
468 * @const mpw *av, *avl@ = first addend vector base and limit
469 * @const mpw *bv, *bvl@ = second addend vector base and limit
470 *
471 * Returns: ---
472 *
473 * Use: Performs unsigned integer addition. If the result overflows
474 * the destination vector, high-order bits are discarded. This
475 * means that two's complement addition happens more or less for
476 * free, although that's more a side-effect than anything else.
477 * The result vector may be equal to either or both source
478 * vectors, but may not otherwise overlap them.
479 */
480
481extern void mpx_uadd(mpw */*dv*/, mpw */*dvl*/,
482 const mpw */*av*/, const mpw */*avl*/,
483 const mpw */*bv*/, const mpw */*bvl*/);
484
dd517851 485/* --- @mpx_uaddn@ --- *
486 *
487 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
488 * @mpw n@ = other addend
3c9ede17 489 *
dd517851 490 * Returns: ---
3c9ede17 491 *
492 * Use: Adds a small integer to a multiprecision number.
493 */
494
495#define MPX_UADDN(dv, dvl, n) do { \
496 mpw *_ddv = (dv), *_ddvl = (dvl); \
497 mpw _c = (n); \
498 \
499 while (_c && _ddv < _ddvl) { \
500 mpd _x = (mpd)*_ddv + (mpd)_c; \
501 *_ddv++ = MPW(_x); \
502 _c = _x >> MPW_BITS; \
503 } \
504} while (0)
505
dd517851 506extern void mpx_uaddn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
507
f46efa79 508/* --- @mpx_uaddnlsl@ --- *
509 *
510 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
511 * @mpw a@ = second argument
512 * @unsigned o@ = offset in bits
513 *
514 * Returns: ---
515 *
516 * Use: Computes %$d + 2^o a$%. If the result overflows then
517 * high-order bits are discarded, as usual. We must have
518 * @0 < o < MPW_BITS@.
519 */
520
521extern void mpx_uaddnlsl(mpw */*dv*/, mpw */*dvl*/,
522 mpw /*a*/, unsigned /*o*/);
523
d03ab969 524/* --- @mpx_usub@ --- *
525 *
526 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
527 * @const mpw *av, *avl@ = first argument vector base and limit
528 * @const mpw *bv, *bvl@ = second argument vector base and limit
529 *
530 * Returns: ---
531 *
532 * Use: Performs unsigned integer subtraction. If the result
533 * overflows the destination vector, high-order bits are
534 * discarded. This means that two's complement subtraction
3c9ede17 535 * happens more or less for free, although that's more a side-
d03ab969 536 * effect than anything else. The result vector may be equal to
537 * either or both source vectors, but may not otherwise overlap
538 * them.
539 */
540
541extern void mpx_usub(mpw */*dv*/, mpw */*dvl*/,
542 const mpw */*av*/, const mpw */*avl*/,
543 const mpw */*bv*/, const mpw */*bvl*/);
544
dd517851 545/* --- @mpx_usubn@ --- *
3c9ede17 546 *
dd517851 547 * Arguments: @mpw *dv, *dvl@ = source and destination base and limit
548 * @n@ = subtrahend
549 *
550 * Returns: ---
3c9ede17 551 *
552 * Use: Subtracts a small integer from a multiprecision number.
553 */
554
555#define MPX_USUBN(dv, dvl, n) do { \
556 mpw *_ddv = (dv), *_ddvl = (dvl); \
557 mpw _c = (n); \
558 \
559 while (_ddv < _ddvl) { \
560 mpd _x = (mpd)*_ddv - (mpd)_c; \
561 *_ddv++ = MPW(_x); \
562 if (_x >> MPW_BITS) \
563 _c = 1; \
564 else \
565 break; \
566 } \
567} while (0)
568
dd517851 569extern void mpx_usubn(mpw */*dv*/, mpw */*dvl*/, mpw /*n*/);
570
f46efa79 571/* --- @mpx_usubnlsl@ --- *
572 *
573 * Arguments: @mpw *dv, *dvl@ = destination and first argument vector
574 * @mpw a@ = second argument
575 * @unsigned o@ = offset in bits
576 *
577 * Returns: ---
578 *
579 * Use: Computes %$d - 2^o a$%. If the result overflows then
580 * high-order bits are discarded, as usual, so you get two's
581 * complement. Which might be what you wanted... We must have
582 * @0 < o < MPW_BITS@.
583 */
584
585extern void mpx_usubnlsl(mpw */*dv*/, mpw */*dvl*/,
586 mpw /*a*/, unsigned /*o*/);
587
3c9ede17 588/* --- @mpx_umul@ --- *
589 *
590 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
591 * @const mpw *av, *avl@ = multiplicand vector base and limit
592 * @const mpw *bv, *bvl@ = multiplier vector base and limit
593 *
594 * Returns: ---
595 *
596 * Use: Performs unsigned integer multiplication. If the result
597 * overflows the desination vector, high-order bits are
598 * discarded. The result vector may not overlap the argument
599 * vectors in any way.
600 */
601
602extern void mpx_umul(mpw */*dv*/, mpw */*dvl*/,
603 const mpw */*av*/, const mpw */*avl*/,
604 const mpw */*bv*/, const mpw */*bvl*/);
605
dd517851 606/* --- @mpx_umuln@ --- *
d03ab969 607 *
dd517851 608 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
609 * @const mpw *av, *avl@ = multiplicand vector base and limit
610 * @mpw m@ = multiplier
611 *
612 * Returns: ---
d03ab969 613 *
614 * Use: Multiplies a multiprecision integer by a single-word value.
615 * The destination and source may be equal. The destination
616 * is completely cleared after use.
617 */
618
619#define MPX_UMULN(dv, dvl, av, avl, m) do { \
620 mpw *_dv = (dv), *_dvl = (dvl); \
621 const mpw *_av = (av), *_avl = (avl); \
622 mpw _c = 0; \
623 mpd _m = (m); \
624 \
625 while (_av < _avl) { \
626 mpd _x; \
627 if (_dv >= _dvl) \
628 break; \
3c9ede17 629 _x = (mpd)_m * (mpd)*_av++ + _c; \
d03ab969 630 *_dv++ = MPW(_x); \
631 _c = _x >> MPW_BITS; \
632 } \
633 if (_dv < _dvl) { \
634 *_dv++ = MPW(_c); \
635 MPX_ZERO(_dv, _dvl); \
636 } \
637} while (0)
638
dd517851 639extern void mpx_umuln(mpw */*dv*/, mpw */*dvl*/,
640 const mpw */*av*/, const mpw */*avl*/, mpw m);
641
642/* --- @mpx_umlan@ --- *
d03ab969 643 *
dd517851 644 * Arguments: @mpw *dv, *dvl@ = destination/accumulator base and limit
645 * @const mpw *av, *avl@ = multiplicand vector base and limit
646 * @mpw m@ = multiplier
647 *
648 * Returns: ---
d03ab969 649 *
650 * Use: Multiplies a multiprecision integer by a single-word value
651 * and adds the result to an accumulator.
652 */
653
654#define MPX_UMLAN(dv, dvl, av, avl, m) do { \
655 mpw *_dv = (dv), *_dvl = (dvl); \
656 const mpw *_av = (av), *_avl = (avl); \
7c13f461 657 mpw _cc = 0; \
d03ab969 658 mpd _m = (m); \
659 \
5bf74dea 660 while (_dv < _dvl && _av < _avl) { \
d03ab969 661 mpd _x; \
7c13f461 662 _x = (mpd)*_dv + (mpd)_m * (mpd)*_av++ + _cc; \
d03ab969 663 *_dv++ = MPW(_x); \
7c13f461 664 _cc = _x >> MPW_BITS; \
d03ab969 665 } \
7c13f461 666 MPX_UADDN(_dv, _dvl, _cc); \
d03ab969 667} while (0)
668
dd517851 669extern void mpx_umlan(mpw */*dv*/, mpw */*dvl*/,
670 const mpw */*av*/, const mpw */*avl*/, mpw m);
671
3c9ede17 672/* --- @mpx_usqr@ --- *
d03ab969 673 *
674 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
3c9ede17 675 * @const mpw *av, *av@ = source vector base and limit
d03ab969 676 *
677 * Returns: ---
678 *
3c9ede17 679 * Use: Performs unsigned integer squaring. The result vector must
680 * not overlap the source vector in any way.
d03ab969 681 */
682
3c9ede17 683extern void mpx_usqr(mpw */*dv*/, mpw */*dvl*/,
684 const mpw */*av*/, const mpw */*avl*/);
d03ab969 685
5bf74dea 686/* --- @mpx_udiv@ --- *
687 *
688 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
689 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
690 * @const mpw *dv, *dvl@ = divisor vector base and limit
45c0fd36 691 * @mpw *sv, *svl@ = scratch workspace
5bf74dea 692 *
693 * Returns: ---
694 *
695 * Use: Performs unsigned integer division. If the result overflows
696 * the quotient vector, high-order bits are discarded. (Clearly
697 * the remainder vector can't overflow.) The various vectors
698 * may not overlap in any way. Yes, I know it's a bit odd
699 * requiring the dividend to be in the result position but it
700 * does make some sense really. The remainder must have
701 * headroom for at least two extra words. The scratch space
702 * must be at least one word larger than the divisor.
703 */
704
705extern void mpx_udiv(mpw */*qv*/, mpw */*qvl*/, mpw */*rv*/, mpw */*rvl*/,
706 const mpw */*dv*/, const mpw */*dvl*/,
707 mpw */*sv*/, mpw */*svl*/);
708
698bd937 709/* --- @mpx_udivn@ --- *
710 *
711 * Arguments: @mpw *qv, *qvl@ = storage for the quotient (may overlap
712 * dividend)
713 * @const mpw *rv, *rvl@ = dividend
714 * @mpw d@ = single-precision divisor
715 *
716 * Returns: Remainder after divison.
717 *
718 * Use: Performs a single-precision division operation.
719 */
720
721extern mpw mpx_udivn(mpw */*qv*/, mpw */*qvl*/,
722 const mpw */*rv*/, const mpw */*rvl*/, mpw /*d*/);
723
5bf74dea 724/*----- Karatsuba multiplication algorithms -------------------------------*/
725
52cdaca9 726/* --- @MPK_THRESH@ --- *
5bf74dea 727 *
728 * This is the limiting length for using Karatsuba algorithms. It's best to
729 * use the simpler classical multiplication method on numbers smaller than
dd22938e 730 * this. It is unsafe to make this constant less than four (i.e., the
731 * algorithms will fail).
5bf74dea 732 */
733
d6b9dc04 734#define MPK_THRESH 32
5bf74dea 735
a86e33af 736/* --- @mpx_kmul@ --- *
737 *
738 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
739 * @const mpw *av, *avl@ = pointer to first argument
740 * @const mpw *bv, *bvl@ = pointer to second argument
741 * @mpw *sv, *svl@ = pointer to scratch workspace
742 *
743 * Returns: ---
744 *
745 * Use: Multiplies two multiprecision integers using Karatsuba's
746 * algorithm. This is rather faster than traditional long
747 * multiplication (e.g., @mpx_umul@) on large numbers, although
748 * more expensive on small ones.
749 *
dd22938e 750 * The destination must be three times as large as the larger
751 * argument. The scratch space must be five times as large as
752 * the larger argument.
a86e33af 753 */
754
a86e33af 755extern void mpx_kmul(mpw */*dv*/, mpw */*dvl*/,
756 const mpw */*av*/, const mpw */*avl*/,
757 const mpw */*bv*/, const mpw */*bvl*/,
758 mpw */*sv*/, mpw */*svl*/);
759
5bf74dea 760/* --- @mpx_ksqr@ --- *
d03ab969 761 *
5bf74dea 762 * Arguments: @mpw *dv, *dvl@ = pointer to destination buffer
763 * @const mpw *av, *avl@ = pointer to first argument
764 * @mpw *sv, *svl@ = pointer to scratch workspace
d03ab969 765 *
766 * Returns: ---
767 *
5bf74dea 768 * Use: Squares a multiprecision integers using something similar to
769 * Karatsuba's multiplication algorithm. This is rather faster
770 * than traditional long multiplication (e.g., @mpx_umul@) on
771 * large numbers, although more expensive on small ones, and
772 * rather simpler than full-blown Karatsuba multiplication.
773 *
dd22938e 774 * The destination must be three times as large as the larger
775 * argument. The scratch space must be five times as large as
776 * the larger argument.
d03ab969 777 */
778
5bf74dea 779extern void mpx_ksqr(mpw */*dv*/, mpw */*dvl*/,
780 const mpw */*av*/, const mpw */*avl*/,
3c9ede17 781 mpw */*sv*/, mpw */*svl*/);
d03ab969 782
783/*----- That's all, folks -------------------------------------------------*/
784
785#ifdef __cplusplus
786 }
787#endif
788
789#endif