projects / u / mdw / catacomb / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
blame | history | raw | HEAD
Initial import.
[u/mdw/catacomb] / mpx.c
1 /* -*-c-*-
2 *
3 * $Id: mpx.c,v 1.1 1999/09/03 08:41:12 mdw Exp $
4 *
5 * Low-level multiprecision arithmetic
6 *
7 * (c) 1999 Straylight/Edgeware
8 */
9
10 /*----- Licensing notice --------------------------------------------------*
11 *
12 * This file is part of Catacomb.
13 *
14 * Catacomb is free software; you can redistribute it and/or modify
15 * it under the terms of the GNU Library General Public License as
16 * published by the Free Software Foundation; either version 2 of the
17 * License, or (at your option) any later version.
18 *
19 * Catacomb is distributed in the hope that it will be useful,
20 * but WITHOUT ANY WARRANTY; without even the implied warranty of
21 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 * GNU Library General Public License for more details.
23 *
24 * You should have received a copy of the GNU Library General Public
25 * License along with Catacomb; if not, write to the Free
26 * Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
27 * MA 02111-1307, USA.
28 */
29
30 /*----- Revision history --------------------------------------------------*
31 *
32 * $Log: mpx.c,v $
33 * Revision 1.1 1999/09/03 08:41:12 mdw
34 * Initial import.
35 *
36 */
37
38 /*----- Header files ------------------------------------------------------*/
39
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <string.h>
43
44 #include <mLib/bits.h>
45
46 #include "mptypes.h"
47 #include "mpx.h"
48
49 /*----- Loading and storing -----------------------------------------------*/
50
51 /* --- @mpx_storel@ --- *
52 *
53 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
54 * @octet *p@ = pointer to octet array
55 * @size_t sz@ = size of octet array
56 *
57 * Returns: ---
58 *
59 * Use: Stores an MP in an octet array, least significant octet
60 * first. High-end octets are silently discarded if there
61 * isn't enough space for them.
62 */
63
64 void mpx_storel(const mpw *v, const mpw *vl, octet *p, size_t sz)
65 {
66 mpw n, w = 0;
67 octet *q = p + sz;
68 unsigned bits = 0;
69
70 while (p < q) {
71 if (bits < 8) {
72 if (v >= vl) {
73 *p++ = U8(w);
74 break;
75 }
76 n = *v++;
77 *p++ = U8(w | n << bits);
78 w = n >> (8 - bits);
79 bits += MPW_BITS - 8;
80 } else {
81 *p++ = U8(w);
82 w >>= 8;
83 bits -= 8;
84 }
85 }
86 memset(p, 0, q - p);
87 }
88
89 /* --- @mpx_loadl@ --- *
90 *
91 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
92 * @const octet *p@ = pointer to octet array
93 * @size_t sz@ = size of octet array
94 *
95 * Returns: ---
96 *
97 * Use: Loads an MP in an octet array, least significant octet
98 * first. High-end octets are ignored if there isn't enough
99 * space for them.
100 */
101
102 void mpx_loadl(mpw *v, const mpw *vl, const octet *p, size_t sz)
103 {
104 unsigned n;
105 const octet *q = p + sz;
106 unsigned bits = 0;
107
108 if (v >= vl)
109 return;
110 while (p < q) {
111 n = U8(*p++);
112 w |= n << bits;
113 bits += 8;
114 if (bits >= MPW_BITS) {
115 *v++ = MPW(w);
116 w = n >> (MPW_BITS - bits + 8);
117 bits -= MPW_BITS;
118 if (v >= vl)
119 return;
120 }
121 }
122 *v++ = w;
123 MPX_ZERO(v, vl);
124 }
125
126 /* --- @mpx_storeb@ --- *
127 *
128 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
129 * @octet *p@ = pointer to octet array
130 * @size_t sz@ = size of octet array
131 *
132 * Returns: ---
133 *
134 * Use: Stores an MP in an octet array, most significant octet
135 * first. High-end octets are silently discarded if there
136 * isn't enough space for them.
137 */
138
139 void mpx_storeb(const mpw *v, const mpw *vl, octet *p, size_t sz);
140 {
141 mpw n, w = 0;
142 octet *q = p + sz;
143 unsigned bits = 0;
144
145 while (q > p) {
146 if (bits < 8) {
147 if (v >= vl) {
148 *--q = U8(w);
149 break;
150 }
151 n = *v++;
152 *--q = U8(w | n << bits);
153 w = n >> (8 - bits);
154 bits += MPW_BITS - 8;
155 } else {
156 *--q = U8(w);
157 w >>= 8;
158 bits -= 8;
159 }
160 }
161 memset(p, 0, q - p);
162 }
163
164 /* --- @mpx_loadb@ --- *
165 *
166 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
167 * @const octet *p@ = pointer to octet array
168 * @size_t sz@ = size of octet array
169 *
170 * Returns: ---
171 *
172 * Use: Loads an MP in an octet array, most significant octet
173 * first. High-end octets are ignored if there isn't enough
174 * space for them.
175 */
176
177 void mpx_loadb(mpw *v, const mpw *vl, const octet *p, size_t sz)
178 {
179 unsigned n;
180 const octet *q = p + sz;
181 unsigned bits = 0;
182
183 if (v >= vl)
184 return;
185 while (q > p) {
186 n = U8(*--q);
187 w |= n << bits;
188 bits += 8;
189 if (bits >= MPW_BITS) {
190 *v++ = MPW(w);
191 w = n >> (MPW_BITS - bits + 8);
192 bits -= MPW_BITS;
193 if (v >= vl)
194 return;
195 }
196 }
197 *v++ = w;
198 MPX_ZERO(v, vl);
199 }
200
201 /*----- Logical shifting --------------------------------------------------*/
202
203 /* --- @mpx_lsl@ --- *
204 *
205 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
206 * @const mpw *av, *avl@ = source vector base and limit
207 * @size_t n@ = number of bit positions to shift by
208 *
209 * Returns: ---
210 *
211 * Use: Performs a logical shift left operation on an integer.
212 */
213
214 void mpx_lsl(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
215 {
216 size_t nw;
217 unsigned nb;
218
219 /* --- Trivial special case --- */
220
221 if (n == 0)
222 MPX_COPY(dv, dvl, av, avl);
223
224 /* --- Single bit shifting --- */
225
226 else if (n == 1) {
227 mpw w = 0;
228 while (av < avl) {
229 mpw t;
230 if (dv >= dvl)
231 goto done;
232 t = *av++;
233 *dv++ = MPW((t << 1) | w);
234 w = t >> (MPW_BITS - 1);
235 }
236 if (dv >= dvl)
237 goto done;
238 *dv++ = MPW(w);
239 MPX_ZERO(dv, dvl);
240 }
241
242 /* --- Break out word and bit shifts for more sophisticated work --- */
243
244 nw = n / MPW_BITS;
245 nb = n % MPW_BITS;
246
247 /* --- Handle a shift by a multiple of the word size --- */
248
249 if (nb == 0) {
250 MPX_COPY(dv + nw, dvl, av, avl);
251 memset(dv, 0, MPWS(nw));
252 }
253
254 /* --- And finally the difficult case --- */
255
256 else {
257 mpw w;
258 size_t nr = MPW_BITS - nb;
259
260 if (dv + nw >= dvl) {
261 MPX_ZERO(dv, dvl);
262 goto done;
263 }
264 memset(dv, 0, MPWS(nw));
265 dv += nw;
266 w = *av++;
267
268 while (av < avl) {
269 mpw t;
270 if (dv >= dvl)
271 goto done;
272 t = *av++;
273 *dv++ = MPW((w >> nr) | (t << nb));
274 w = t;
275 }
276
277 if (dv < dvl) {
278 *dv++ = MPW(w >> nr);
279 MPX_ZERO(dv, dvl);
280 }
281 }
282
283 done:;
284 }
285
286 /* --- @mpx_lsr@ --- *
287 *
288 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
289 * @const mpw *av, *avl@ = source vector base and limit
290 * @size_t n@ = number of bit positions to shift by
291 *
292 * Returns: ---
293 *
294 * Use: Performs a logical shift right operation on an integer.
295 */
296
297 void mpx_lsr(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl, size_t n)
298 {
299 size_t nw;
300 unsigned nb;
301
302 /* --- Trivial special case --- */
303
304 if (n == 0)
305 MPX_COPY(dv, dvl, av, avl);
306
307 /* --- Single bit shifting --- */
308
309 else if (n == 1) {
310 mpw w = *av++ >> 1;
311 while (av < avl) {
312 mpw t;
313 if (dv >= dvl)
314 goto done;
315 t = *av++;
316 *dv++ = MPW((t << (MPW_BITS - 1)) | w);
317 w = t >> 1;
318 }
319 if (dv >= dvl)
320 goto done;
321 *dv++ = MPW(w);
322 MPX_ZERO(dv, dvl);
323 }
324
325 /* --- Break out word and bit shifts for more sophisticated work --- */
326
327 nw = n / MPW_BITS;
328 nb = n % MPW_BITS;
329
330 /* --- Handle a shift by a multiple of the word size --- */
331
332 if (nb == 0)
333 MPX_COPY(dv, dvl, av + nw, avl);
334
335 /* --- And finally the difficult case --- */
336
337 else {
338 mpw w;
339 size_t nr = MPW_BITS - nb;
340
341 av += nw;
342 w = *av++;
343 while (av < avl) {
344 mpw t;
345 if (dv >= dvl)
346 goto done;
347 t = *av++;
348 *dv++ = MPW((w >> nb) | (t << nr));
349 w = t;
350 }
351 if (dv < dvl) {
352 *dv++ = MPW(w >> nb);
353 MPX_ZERO(dv, dvl);
354 }
355 }
356
357 done:;
358 }
359
360 /*----- Unsigned arithmetic -----------------------------------------------*/
361
362 /* --- @mpx_ucmp@ --- *
363 *
364 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
365 * @const mpw *bv, *bvl@ = second argument vector base and limit
366 *
367 * Returns: Less than, equal to, or greater than zero depending on
368 * whether @a@ is less than, equal to or greater than @b@,
369 * respectively.
370 *
371 * Use: Performs an unsigned integer comparison.
372 */
373
374 int mpx_ucmp(const mpw *av, const mpw *avl, const mpw *bv, const mpw *bvl)
375 {
376 MPX_SHRINK(av, avl);
377 MPX_SHRINK(bv, bvl);
378
379 if (avl - av > bvl - bv)
380 return (+1);
381 else if (avl - av < bvl - bv)
382 return (-1);
383 else while (avl > av) {
384 mpw a = *--avl, b = *--bvl;
385 if (a > b)
386 return (+1);
387 else if (a < b)
388 return (-1);
389 }
390 return (0);
391 }
392
393 /* --- @mpx_uadd@ --- *
394 *
395 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
396 * @const mpw *av, *avl@ = first addend vector base and limit
397 * @const mpw *bv, *bvl@ = second addend vector base and limit
398 *
399 * Returns: ---
400 *
401 * Use: Performs unsigned integer addition. If the result overflows
402 * the destination vector, high-order bits are discarded. This
403 * means that two's complement addition happens more or less for
404 * free, although that's more a side-effect than anything else.
405 * The result vector may be equal to either or both source
406 * vectors, but may not otherwise overlap them.
407 */
408
409 void mpx_uadd(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
410 const mpw *bv, const mpw *bvl)
411 {
412 mpw c = 0;
413
414 while (av < avl || bv < bvl) {
415 mpw a, b;
416 mpd x;
417 if (dv >= dvl)
418 return;
419 a = (av < avl) ? *av++ : 0;
420 b = (bv < bvl) ? *bv++ : 0;
421 x = (mpd)a + (mpd)b + c;
422 *dv++ = MPW(x);
423 c = x >> MPW_BITS;
424 }
425 if (dv < dvl) {
426 *dv++ = c;
427 MPX_ZERO(dv, dvl);
428 }
429 }
430
431 /* --- @mpx_usub@ --- *
432 *
433 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
434 * @const mpw *av, *avl@ = first argument vector base and limit
435 * @const mpw *bv, *bvl@ = second argument vector base and limit
436 *
437 * Returns: ---
438 *
439 * Use: Performs unsigned integer subtraction. If the result
440 * overflows the destination vector, high-order bits are
441 * discarded. This means that two's complement subtraction
442 * happens more or less for free, althuogh that's more a side-
443 * effect than anything else. The result vector may be equal to
444 * either or both source vectors, but may not otherwise overlap
445 * them.
446 */
447
448 void mpx_usub(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
449 const mpw *bv, const mpw *bvl)
450 {
451 mpw c = 0;
452
453 while (av < avl || bv < bvl) {
454 mpw a, b;
455 mpd x;
456 if (dv >= dvl)
457 return;
458 a = (av < avl) ? *av++ : 0;
459 b = (bv < bvl) ? *bv++ : 0;
460 x = (mpd)a - (mpd)b + c;
461 *dv++ = MPW(x);
462 if (c >> MPW_BITS)
463 c = MPW(~0u);
464 }
465 c = c ? ~0u : 0;
466 while (dv < dvl)
467 *dv++ = c
468 }
469
470 /* --- @mpx_umul@ --- *
471 *
472 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
473 * @const mpw *av, *avl@ = multiplicand vector base and limit
474 * @const mpw *bv, *bvl@ = multiplier vector base and limit
475 *
476 * Returns: ---
477 *
478 * Use: Performs unsigned integer multiplication. If the result
479 * overflows the desination vector, high-order bits are
480 * discarded. The result vector may not overlap the argument
481 * vectors in any way.
482 */
483
484 void mpx_umul(mpw *dv, mpw *dvl, const mpw *av, const mpw *avl,
485 const mpw *bv, const mpw *bvl)
486 {
487 /* --- This is probably worthwhile on a multiply --- */
488
489 MPX_SHRINK(av, avl);
490 MPX_SHRINK(bv, bvl);
491
492 /* --- Deal with a multiply by zero --- */
493
494 if (bv == bvl) {
495 MPX_COPY(dv, dvl, bv, bvl);
496 return;
497 }
498
499 /* --- Do the initial multiply and initialize the accumulator --- */
500
501 MPX_UMULN(dv, dvl, av, avl, *bv++);
502
503 /* --- Do the remaining multiply/accumulates --- */
504
505 while (bv < bvl) {
506 mpw m = *bv++;
507 mpw c = ;
508 const mpw *avv = av;
509 mpw *dvv = ++dv;
510
511 while (avv < avl) {
512 mpd x;
513 if (dvv >= dvl)
514 goto next;
515 x = *dvv + m * *av++ + c;
516 *dv++ = MPW(x);
517 c = x >> MPW_BITS;
518 }
519 if (dvv < dvl)
520 *dvv++ = MPW(c);
521 next:;
522 }
523 }
524
525 /* --- @mpx_udiv@ --- *
526 *
527 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
528 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
529 * @const mpw *dv, *dvl@ = divisor vector base and limit
530 *
531 * Returns: ---
532 *
533 * Use: Performs unsigned integer division. If the result overflows
534 * the quotient vector, high-order bits are discarded. (Clearly
535 * the remainder vector can't overflow.) The various vectors
536 * may not overlap in any way. Yes, I know it's a bit odd
537 * requiring the dividend to be in the result position but it
538 * does make some sense really. The remainder must have
539 * headroom for at least two extra words.
540 */
541
542 void mpx_udiv(mpw *qv, mpw *qvl, mpw *rv, mpw *rvl,
543 const mpw *dv, const mpw *dvl)
544 {
545 mpw spare[2];
546 unsigned norm = 0;
547 size_t scale;
548 mpw d, dd;
549
550 /* --- Initialize the quotient --- */
551
552 MPX_ZERO(qv, qvl);
553
554 /* --- Normalize the divisor --- *
555 *
556 * The algorithm requires that the divisor be at least two digits long.
557 * This is easy to fix.
558 */
559
560 MPX_SHRINK(dv, dvl);
561
562 assert(((void)"division by zero in mpx_udiv", dv < dvl));
563
564 d = dvl[-1];
565 if (dv + 1 == dvl) {
566 spare[0] = 0;
567 spare[1] = d;
568 dv = spare;
569 dvl = spare + 2;
570 norm += MPW_BITS;
571 }
572
573 while (d < MPW_MAX / 2) {
574 d <<= 1;
575 norm += 1;
576 }
577 dd = dvl[-2];
578
579 /* --- Normalize the dividend/remainder to match --- */
580
581 mpx_lsl(rv, rvl, rv, rvl, norm);
582 MPX_SHRINK(rv, rvl);
583
584 /* --- Work out the relative scales --- */
585
586 {
587 size_t rvn = rvl - rv;
588 size_t dvn = dvn - dv;
589
590 /* --- If the divisor is clearly larger, notice this --- */
591
592 if (dvn > rvn) {
593 mpx_lsr(rv, rvl, rv, rvl, norm);
594 return;
595 }
596
597 scale = rvn - dvn;
598 }
599
600 /* --- Calculate the most significant quotient digit --- *
601 *
602 * Because the divisor has its top bit set, this can only happen once. The
603 * pointer arithmetic is a little contorted, to make sure that the
604 * behaviour is defined.
605 */
606
607 if (MPX_UCMP(rv + scale, rvl, >=, dv, dvl)) {
608 mpx_usub(rv + scale, rvl, rv + scale, rvl, dv, dvl);
609 if (qvl - qv > scale)
610 qv[scale] = 1;
611 }
612
613 /* --- Now for the main loop --- */
614
615 {
616 mpw *rvv;
617 mpw r;
618
619 scale--;
620 rvv = rvl - 2;
621 r = rvv[1];
622
623 while (scale) {
624 mpw q, rr;
625
626 /* --- Get an estimate for the next quotient digit --- */
627
628 rr = *rvv--;
629 if (r == d)
630 q = MPW_MAX;
631 else {
632 mpd rx = (r << MPW_BITS) | rr;
633 q = MPW(rx / d);
634 }
635
636 /* --- Refine the estimate --- */
637
638 {
639 mpd yh = (mpd)d * q;
640 mpd yl = (mpd)dd * q;
641
642 }
643
644 /*----- That's all, folks -------------------------------------------------*/
Catacomb cryptographic library (downstream of http://git.distorted.org.uk/~mdw/catacomb/)