3 * $Id: mpx.c,v 1.4 1999/11/17 18:04:09 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.4 1999/11/17 18:04:09 mdw
34 * Add two's-complement functionality. Improve mpx_udiv a little by
35 * performing the multiplication of the divisor by q with the subtraction
38 * Revision 1.3 1999/11/13 01:57:31 mdw
39 * Remove stray debugging code.
41 * Revision 1.2 1999/11/13 01:50:59 mdw
42 * Multiprecision routines finished and tested.
44 * Revision 1.1 1999/09/03 08:41:12 mdw
49 /*----- Header files ------------------------------------------------------*/
56 #include <mLib/bits.h>
61 /*----- Loading and storing -----------------------------------------------*/
63 /* --- @mpx_storel@ --- *
65 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
66 * @void *pp@ = pointer to octet array
67 * @size_t sz@ = size of octet array
71 * Use: Stores an MP in an octet array, least significant octet
72 * first. High-end octets are silently discarded if there
73 * isn't enough space for them.
76 void mpx_storel(const mpw
*v
, const mpw
*vl
, void *pp
, size_t sz
)
79 octet
*p
= pp
, *q
= p
+ sz
;
89 *p
++ = U8(w
| n
<< bits
);
101 /* --- @mpx_loadl@ --- *
103 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
104 * @const void *pp@ = pointer to octet array
105 * @size_t sz@ = size of octet array
109 * Use: Loads an MP in an octet array, least significant octet
110 * first. High-end octets are ignored if there isn't enough
114 void mpx_loadl(mpw
*v
, mpw
*vl
, const void *pp
, size_t sz
)
118 const octet
*p
= pp
, *q
= p
+ sz
;
127 if (bits
>= MPW_BITS
) {
129 w
= n
>> (MPW_BITS
- bits
+ 8);
139 /* --- @mpx_storeb@ --- *
141 * Arguments: @const mpw *v, *vl@ = base and limit of source vector
142 * @void *pp@ = pointer to octet array
143 * @size_t sz@ = size of octet array
147 * Use: Stores an MP in an octet array, most significant octet
148 * first. High-end octets are silently discarded if there
149 * isn't enough space for them.
152 void mpx_storeb(const mpw
*v
, const mpw
*vl
, void *pp
, size_t sz
)
155 octet
*p
= pp
, *q
= p
+ sz
;
165 *--q
= U8(w
| n
<< bits
);
167 bits
+= MPW_BITS
- 8;
177 /* --- @mpx_loadb@ --- *
179 * Arguments: @mpw *v, *vl@ = base and limit of destination vector
180 * @const void *pp@ = pointer to octet array
181 * @size_t sz@ = size of octet array
185 * Use: Loads an MP in an octet array, most significant octet
186 * first. High-end octets are ignored if there isn't enough
190 void mpx_loadb(mpw
*v
, mpw
*vl
, const void *pp
, size_t sz
)
194 const octet
*p
= pp
, *q
= p
+ sz
;
203 if (bits
>= MPW_BITS
) {
205 w
= n
>> (MPW_BITS
- bits
+ 8);
215 /*----- Logical shifting --------------------------------------------------*/
217 /* --- @mpx_lsl@ --- *
219 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
220 * @const mpw *av, *avl@ = source vector base and limit
221 * @size_t n@ = number of bit positions to shift by
225 * Use: Performs a logical shift left operation on an integer.
228 void mpx_lsl(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, size_t n
)
233 /* --- Trivial special case --- */
236 MPX_COPY(dv
, dvl
, av
, avl
);
238 /* --- Single bit shifting --- */
247 *dv
++ = MPW((t
<< 1) | w
);
248 w
= t
>> (MPW_BITS
- 1);
257 /* --- Break out word and bit shifts for more sophisticated work --- */
262 /* --- Handle a shift by a multiple of the word size --- */
265 MPX_COPY(dv
+ nw
, dvl
, av
, avl
);
266 memset(dv
, 0, MPWS(nw
));
269 /* --- And finally the difficult case --- *
271 * This is a little convoluted, because I have to start from the end and
272 * work backwards to avoid overwriting the source, if they're both the same
278 size_t nr
= MPW_BITS
- nb
;
279 size_t dvn
= dvl
- dv
;
280 size_t avn
= avl
- av
;
287 if (dvn
> avn
+ nw
) {
288 size_t off
= avn
+ nw
+ 1;
289 MPX_ZERO(dv
+ off
, dvl
);
299 *--dvl
= (t
>> nr
) | w
;
310 /* --- @mpx_lsr@ --- *
312 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
313 * @const mpw *av, *avl@ = source vector base and limit
314 * @size_t n@ = number of bit positions to shift by
318 * Use: Performs a logical shift right operation on an integer.
321 void mpx_lsr(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
, size_t n
)
326 /* --- Trivial special case --- */
329 MPX_COPY(dv
, dvl
, av
, avl
);
331 /* --- Single bit shifting --- */
340 *dv
++ = MPW((t
<< (MPW_BITS
- 1)) | w
);
350 /* --- Break out word and bit shifts for more sophisticated work --- */
355 /* --- Handle a shift by a multiple of the word size --- */
358 MPX_COPY(dv
, dvl
, av
+ nw
, avl
);
360 /* --- And finally the difficult case --- */
364 size_t nr
= MPW_BITS
- nb
;
373 *dv
++ = MPW((w
>> nb
) | (t
<< nr
));
377 *dv
++ = MPW(w
>> nb
);
385 /*----- Unsigned arithmetic -----------------------------------------------*/
387 /* --- @mpx_2c@ --- *
389 * Arguments: @mpw *dv, *dvl@ = destination vector
390 * @const mpw *v, *vl@ = source vector
394 * Use: Calculates the two's complement of @v@.
397 void mpx_2c(mpw
*dv
, mpw
*dvl
, const mpw
*v
, const mpw
*vl
)
400 while (dv
< dvl
&& v
< vl
)
401 *dv
++ = c
= MPW(~*v
++);
408 MPX_UADDN(dv
, dvl
, 1);
411 /* --- @mpx_ucmp@ --- *
413 * Arguments: @const mpw *av, *avl@ = first argument vector base and limit
414 * @const mpw *bv, *bvl@ = second argument vector base and limit
416 * Returns: Less than, equal to, or greater than zero depending on
417 * whether @a@ is less than, equal to or greater than @b@,
420 * Use: Performs an unsigned integer comparison.
423 int mpx_ucmp(const mpw
*av
, const mpw
*avl
, const mpw
*bv
, const mpw
*bvl
)
428 if (avl
- av
> bvl
- bv
)
430 else if (avl
- av
< bvl
- bv
)
432 else while (avl
> av
) {
433 mpw a
= *--avl
, b
= *--bvl
;
442 /* --- @mpx_uadd@ --- *
444 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
445 * @const mpw *av, *avl@ = first addend vector base and limit
446 * @const mpw *bv, *bvl@ = second addend vector base and limit
450 * Use: Performs unsigned integer addition. If the result overflows
451 * the destination vector, high-order bits are discarded. This
452 * means that two's complement addition happens more or less for
453 * free, although that's more a side-effect than anything else.
454 * The result vector may be equal to either or both source
455 * vectors, but may not otherwise overlap them.
458 void mpx_uadd(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
459 const mpw
*bv
, const mpw
*bvl
)
463 while (av
< avl
|| bv
< bvl
) {
468 a
= (av
< avl
) ?
*av
++ : 0;
469 b
= (bv
< bvl
) ?
*bv
++ : 0;
470 x
= (mpd
)a
+ (mpd
)b
+ c
;
480 /* --- @mpx_usub@ --- *
482 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
483 * @const mpw *av, *avl@ = first argument vector base and limit
484 * @const mpw *bv, *bvl@ = second argument vector base and limit
488 * Use: Performs unsigned integer subtraction. If the result
489 * overflows the destination vector, high-order bits are
490 * discarded. This means that two's complement subtraction
491 * happens more or less for free, althuogh that's more a side-
492 * effect than anything else. The result vector may be equal to
493 * either or both source vectors, but may not otherwise overlap
497 void mpx_usub(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
498 const mpw
*bv
, const mpw
*bvl
)
502 while (av
< avl
|| bv
< bvl
) {
507 a
= (av
< avl
) ?
*av
++ : 0;
508 b
= (bv
< bvl
) ?
*bv
++ : 0;
509 x
= (mpd
)a
- (mpd
)b
- c
;
522 /* --- @mpx_umul@ --- *
524 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
525 * @const mpw *av, *avl@ = multiplicand vector base and limit
526 * @const mpw *bv, *bvl@ = multiplier vector base and limit
530 * Use: Performs unsigned integer multiplication. If the result
531 * overflows the desination vector, high-order bits are
532 * discarded. The result vector may not overlap the argument
533 * vectors in any way.
536 void mpx_umul(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
,
537 const mpw
*bv
, const mpw
*bvl
)
539 /* --- This is probably worthwhile on a multiply --- */
544 /* --- Deal with a multiply by zero --- */
551 /* --- Do the initial multiply and initialize the accumulator --- */
553 MPX_UMULN(dv
, dvl
, av
, avl
, *bv
++);
555 /* --- Do the remaining multiply/accumulates --- */
557 while (dv
< dvl
&& bv
< bvl
) {
567 x
= (mpd
)*dvv
+ (mpd
)m
* (mpd
)*avv
++ + c
;
571 MPX_UADDN(dvv
, dvl
, c
);
576 /* --- @mpx_usqr@ --- *
578 * Arguments: @mpw *dv, *dvl@ = destination vector base and limit
579 * @const mpw *av, *av@ = source vector base and limit
583 * Use: Performs unsigned integer squaring. The result vector must
584 * not overlap the source vector in any way.
587 void mpx_usqr(mpw
*dv
, mpw
*dvl
, const mpw
*av
, const mpw
*avl
)
591 /* --- Main loop --- */
599 /* --- Stop if I've run out of destination --- */
604 /* --- Work out the square at this point in the proceedings --- */
607 mpd x
= (mpd
)a
* (mpd
)a
+ *dvv
;
609 c
= MPW(x
>> MPW_BITS
);
612 /* --- Now fix up the rest of the vector upwards --- */
615 while (dvv
< dvl
&& avv
< avl
) {
616 mpd x
= (mpd
)a
* (mpd
)*avv
++;
617 mpd y
= ((x
<< 1) & MPW_MAX
) + c
+ *dvv
;
618 c
= (x
>> (MPW_BITS
- 1)) + (y
>> MPW_BITS
);
621 while (dvv
< dvl
&& c
) {
627 /* --- Get ready for the next round --- */
634 /* --- @mpx_udiv@ --- *
636 * Arguments: @mpw *qv, *qvl@ = quotient vector base and limit
637 * @mpw *rv, *rvl@ = dividend/remainder vector base and limit
638 * @const mpw *dv, *dvl@ = divisor vector base and limit
639 * @mpw *sv, *svl@ = scratch workspace
643 * Use: Performs unsigned integer division. If the result overflows
644 * the quotient vector, high-order bits are discarded. (Clearly
645 * the remainder vector can't overflow.) The various vectors
646 * may not overlap in any way. Yes, I know it's a bit odd
647 * requiring the dividend to be in the result position but it
648 * does make some sense really. The remainder must have
649 * headroom for at least two extra words. The scratch space
650 * must be at least one word larger than the divisor.
653 void mpx_udiv(mpw
*qv
, mpw
*qvl
, mpw
*rv
, mpw
*rvl
,
654 const mpw
*dv
, const mpw
*dvl
,
661 /* --- Initialize the quotient --- */
665 /* --- Perform some sanity checks --- */
668 assert(((void)"division by zero in mpx_udiv", dv
< dvl
));
670 /* --- Normalize the divisor --- *
672 * The algorithm requires that the divisor be at least two digits long.
673 * This is easy to fix.
680 for (b
= MPW_BITS
/ 2; b
; b
>>= 1) {
681 if (d
< (MPW_MAX
>> b
)) {
690 /* --- Normalize the dividend/remainder to match --- */
693 mpx_lsl(rv
, rvl
, rv
, rvl
, norm
);
694 mpx_lsl(sv
, svl
, dv
, dvl
, norm
);
704 /* --- Work out the relative scales --- */
707 size_t rvn
= rvl
- rv
;
708 size_t dvn
= dvl
- dv
;
710 /* --- If the divisor is clearly larger, notice this --- */
713 mpx_lsr(rv
, rvl
, rv
, rvl
, norm
);
720 /* --- Calculate the most significant quotient digit --- *
722 * Because the divisor has its top bit set, this can only happen once. The
723 * pointer arithmetic is a little contorted, to make sure that the
724 * behaviour is defined.
727 if (MPX_UCMP(rv
+ scale
, rvl
, >=, dv
, dvl
)) {
728 mpx_usub(rv
+ scale
, rvl
, rv
+ scale
, rvl
, dv
, dvl
);
729 if (qvl
- qv
> scale
)
733 /* --- Now for the main loop --- */
742 /* --- Get an estimate for the next quotient digit --- */
749 rh
= ((mpd
)r
<< MPW_BITS
) | rr
;
755 /* --- Refine the estimate --- */
759 mpd yl
= (mpd
)dd
* q
;
762 yh
+= yl
>> MPW_BITS
;
766 while (yh
> rh
|| (yh
== rh
&& yl
> rrr
)) {
777 /* --- Remove a chunk from the dividend --- */
784 /* --- Calculate the size of the chunk --- *
786 * This does the whole job of calculating @r >> scale - qd@.
789 for (svv
= rv
+ scale
, dvv
= dv
;
790 dvv
< dvl
&& svv
< rvl
;
792 mpd x
= (mpd
)*dvv
* (mpd
)q
+ mc
;
794 x
= (mpd
)*svv
- MPW(x
) - sc
;
803 mpd x
= (mpd
)*svv
- mc
- sc
;
813 /* --- Fix if the quotient was too large --- *
815 * This doesn't seem to happen very often.
818 if (rvl
[-1] > MPW_MAX
/ 2) {
819 mpx_uadd(rv
+ scale
, rvl
, rv
+ scale
, rvl
, dv
, dvl
);
824 /* --- Done for another iteration --- */
826 if (qvl
- qv
> scale
)
833 /* --- Now fiddle with unnormalizing and things --- */
835 mpx_lsr(rv
, rvl
, rv
, rvl
, norm
);
838 /*----- That's all, folks -------------------------------------------------*/