3 * $Id: mpx.c,v 1.1 1999/09/03 08:41:12 mdw Exp $
5 * Low-level multiprecision arithmetic
7 * (c) 1999 Straylight/Edgeware
10 /*----- Licensing notice --------------------------------------------------*
12 * This file is part of Catacomb.
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.
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.
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,
30 /*----- Revision history --------------------------------------------------*
33 * Revision 1.1 1999/09/03 08:41:12 mdw
38 /*----- Header files ------------------------------------------------------*/
44 #include <mLib/bits.h>
49 /*----- Loading and storing -----------------------------------------------*/
51 /* --- @mpx_storel@ --- *
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
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.
64 void mpx_storel(const mpw
*v
, const mpw
*vl
, octet
*p
, size_t sz
)
77 *p
++ = U8(w
| n
<< bits
);
89 /* --- @mpx_loadl@ --- *
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
97 * Use: Loads an MP in an octet array, least significant octet
98 * first. High-end octets are ignored if there isn't enough
102 void mpx_loadl(mpw
*v
, const mpw
*vl
, const octet
*p
, size_t sz
)
105 const octet
*q
= p
+ sz
;
114 if (bits
>= MPW_BITS
) {
116 w
= n
>> (MPW_BITS
- bits
+ 8);
126 /* --- @mpx_storeb@ --- *
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
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.
139 void mpx_storeb(const mpw
*v
, const mpw
*vl
, octet
*p
, size_t sz
);
152 *--q
= U8(w
| n
<< bits
);
154 bits
+= MPW_BITS
- 8;
164 /* --- @mpx_loadb@ --- *
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
172 * Use: Loads an MP in an octet array, most significant octet
173 * first. High-end octets are ignored if there isn't enough
177 void mpx_loadb(mpw
*v
, const mpw
*vl
, const octet
*p
, size_t sz
)
180 const octet
*q
= p
+ sz
;
189 if (bits
>= MPW_BITS
) {
191 w
= n
>> (MPW_BITS
- bits
+ 8);
201 /*----- Logical shifting --------------------------------------------------*/
203 /* --- @mpx_lsl@ --- *
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
211 * Use: Performs a logical shift left operation on an integer.
214 void mpx_lsl(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, size_t n
)
219 /* --- Trivial special case --- */
222 MPX_COPY(dv
, dvl
, av
, avl
);
224 /* --- Single bit shifting --- */
233 *dv
++ = MPW((t
<< 1) | w
);
234 w
= t
>> (MPW_BITS
- 1);
242 /* --- Break out word and bit shifts for more sophisticated work --- */
247 /* --- Handle a shift by a multiple of the word size --- */
250 MPX_COPY(dv
+ nw
, dvl
, av
, avl
);
251 memset(dv
, 0, MPWS(nw
));
254 /* --- And finally the difficult case --- */
258 size_t nr
= MPW_BITS
- nb
;
260 if (dv
+ nw
>= dvl
) {
264 memset(dv
, 0, MPWS(nw
));
273 *dv
++ = MPW((w
>> nr
) | (t
<< nb
));
278 *dv
++ = MPW(w
>> nr
);
286 /* --- @mpx_lsr@ --- *
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
294 * Use: Performs a logical shift right operation on an integer.
297 void mpx_lsr(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, size_t n
)
302 /* --- Trivial special case --- */
305 MPX_COPY(dv
, dvl
, av
, avl
);
307 /* --- Single bit shifting --- */
316 *dv
++ = MPW((t
<< (MPW_BITS
- 1)) | w
);
325 /* --- Break out word and bit shifts for more sophisticated work --- */
330 /* --- Handle a shift by a multiple of the word size --- */
333 MPX_COPY(dv
, dvl
, av
+ nw
, avl
);
335 /* --- And finally the difficult case --- */
339 size_t nr
= MPW_BITS
- nb
;
348 *dv
++ = MPW((w
>> nb
) | (t
<< nr
));
352 *dv
++ = MPW(w
>> nb
);
360 /*----- Unsigned arithmetic -----------------------------------------------*/
362 /* --- @mpx_ucmp@ --- *
364 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
365 * @const mpw *bv, *bvl@ = second argument vector base and limit
367 * Returns: Less than, equal to, or greater than zero depending on
368 * whether @a@ is less than, equal to or greater than @b@,
371 * Use: Performs an unsigned integer comparison.
374 int mpx_ucmp(const mpw
*av
, const mpw
*avl
, const mpw
*bv
, const mpw
*bvl
)
379 if (avl
- av
> bvl
- bv
)
381 else if (avl
- av
< bvl
- bv
)
383 else while (avl
> av
) {
384 mpw a
= *--avl
, b
= *--bvl
;
393 /* --- @mpx_uadd@ --- *
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
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.
409 void mpx_uadd(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
410 const mpw
*bv
, const mpw
*bvl
)
414 while (av
< avl
|| bv
< bvl
) {
419 a
= (av
< avl
) ?
*av
++ : 0;
420 b
= (bv
< bvl
) ?
*bv
++ : 0;
421 x
= (mpd
)a
+ (mpd
)b
+ c
;
431 /* --- @mpx_usub@ --- *
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
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
448 void mpx_usub(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
449 const mpw
*bv
, const mpw
*bvl
)
453 while (av
< avl
|| bv
< bvl
) {
458 a
= (av
< avl
) ?
*av
++ : 0;
459 b
= (bv
< bvl
) ?
*bv
++ : 0;
460 x
= (mpd
)a
- (mpd
)b
+ c
;
470 /* --- @mpx_umul@ --- *
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
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.
484 void mpx_umul(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
485 const mpw
*bv
, const mpw
*bvl
)
487 /* --- This is probably worthwhile on a multiply --- */
492 /* --- Deal with a multiply by zero --- */
495 MPX_COPY(dv
, dvl
, bv
, bvl
);
499 /* --- Do the initial multiply and initialize the accumulator --- */
501 MPX_UMULN(dv
, dvl
, av
, avl
, *bv
++);
503 /* --- Do the remaining multiply/accumulates --- */
515 x
= *dvv
+ m
* *av
++ + c
;
525 /* --- @mpx_udiv@ --- *
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
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.
542 void mpx_udiv(mpw
*qv
, mpw
*qvl
, mpw
*rv
, mpw
*rvl
,
543 const mpw
*dv
, const mpw
*dvl
)
550 /* --- Initialize the quotient --- */
554 /* --- Normalize the divisor --- *
556 * The algorithm requires that the divisor be at least two digits long.
557 * This is easy to fix.
562 assert(((void)"division by zero in mpx_udiv", dv
< dvl
));
573 while (d
< MPW_MAX
/ 2) {
579 /* --- Normalize the dividend/remainder to match --- */
581 mpx_lsl(rv
, rvl
, rv
, rvl
, norm
);
584 /* --- Work out the relative scales --- */
587 size_t rvn
= rvl
- rv
;
588 size_t dvn
= dvn
- dv
;
590 /* --- If the divisor is clearly larger, notice this --- */
593 mpx_lsr(rv
, rvl
, rv
, rvl
, norm
);
600 /* --- Calculate the most significant quotient digit --- *
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.
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
)
613 /* --- Now for the main loop --- */
626 /* --- Get an estimate for the next quotient digit --- */
632 mpd rx
= (r
<< MPW_BITS
) | rr
;
636 /* --- Refine the estimate --- */
640 mpd yl
= (mpd
)dd
* q
;
644 /*----- That's all, folks -------------------------------------------------*/