chiark - git - mdw - catacomb/blob - mp-arith.c

   1 /* -*-c-*-
   2  *
   3  * $Id: mp-arith.c,v 1.1 1999/11/17 18:02:16 mdw Exp $
   4  *
   5  * Basic arithmetic on multiprecision integers
   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: mp-arith.c,v $
  33  * Revision 1.1  1999/11/17 18:02:16  mdw
  34  * New multiprecision integer arithmetic suite.
  35  *
  36  */
  37
  38 /*----- Header files ------------------------------------------------------*/
  39
  40 #include "mp.h"
  41
  42 /*----- Main code ---------------------------------------------------------*/
  43
  44 /* --- @mp_2c@ --- *
  45  *
  46  * Arguments:   @mp *a@ = source
  47  *
  48  * Returns:     Result, @a@ converted to two's complement notation.
  49  */
  50
  51 mp *mp_2c(mp *d, mp *a)
  52 {
  53   if (!(a->f & MP_NEG))
  54     return (MP_COPY(a));
  55
  56   MP_MODIFY(d, MP_LEN(a));
  57   mpx_2c(d->v, d->vl, a->v, a->vl);
  58   d->f = a->f & MP_BURN;
  59   MP_SHRINK(d);
  60   return (d);
  61 }
  62
  63 /* --- @mp_sm@ --- *
  64  *
  65  * Arguments:   @mp *d@ = destination
  66  *              @mp *a@ = source
  67  *
  68  * Returns:     Result, @a@ converted to the native signed-magnitude
  69  *              notation.
  70  */
  71
  72 mp *mp_sm(mp *d, mp *a)
  73 {
  74   if (!MP_LEN(a) || a->vl[-1] < MPW_MAX / 2)
  75     return (MP_COPY(a));
  76
  77   MP_MODIFY(d, MP_LEN(a));
  78   mpx_2c(d->v, d->vl, a->v, a->vl);
  79   d->f = (a->f & (MP_BURN | MP_NEG)) ^ MP_NEG;
  80   MP_SHRINK(d);
  81   return (d);
  82 }
  83
  84 /* --- @mp_lsl@ --- *
  85  *
  86  * Arguments:   @mp *d@ = destination
  87  *              @const mp *a@ = source
  88  *              @size_t n@ = number of bits to move
  89  *
  90  * Returns:     Result, @a@ shifted left by @n@.
  91  */
  92
  93 mp *mp_lsl(mp *d, const mp *a, size_t n)
  94 {
  95   MP_MODIFY(d, MP_LEN(a) + (n + MPW_BITS - 1) / MPW_BITS);
  96   mpx_lsl(d->v, d->vl, a->v, a->vl, n);
  97   d->f = a->f & (MP_NEG | MP_BURN);
  98   MP_SHRINK(d);
  99   return (d);
 100 }
 101
 102 /* --- @mp_lsr@ --- *
 103  *
 104  * Arguments:   @mp *d@ = destination
 105  *              @const mp *a@ = source
 106  *              @size_t n@ = number of bits to move
 107  *
 108  * Returns:     Result, @a@ shifted left by @n@.
 109  */
 110
 111 mp *mp_lsr(mp *d, const mp *a, size_t n)
 112 {
 113   MP_MODIFY(d, MP_LEN(a));
 114   mpx_lsr(d->v, d->vl, a->v, a->vl, n);
 115   d->f = a->f & (MP_NEG | MP_BURN);
 116   MP_SHRINK(d);
 117   return (d);
 118 }
 119
 120 /* --- @mp_cmp@ --- *
 121  *
 122  * Arguments:   @const mp *a, *b@ = two numbers
 123  *
 124  * Returns:     Less than, equal to or greater than zero, according to
 125  *              whether @a@ is less than, equal to or greater than @b@.
 126  */
 127
 128 int mp_cmp(const mp *a, const mp *b)
 129 {
 130   if (!((a->f ^ b->f) & MP_NEG))
 131     return (mpx_ucmp(a->v, a->vl, b->v, b->vl));
 132   else if (a->f & MP_NEG)
 133     return (-1);
 134   else
 135     return (+1);
 136 }
 137
 138 /* --- @mp_add@ --- *
 139  *
 140  * Arguments:   @mp *d@ = destination
 141  *              @const mp *a, *b@ = sources
 142  *
 143  * Returns:     Result, @a@ added to @b@.
 144  */
 145
 146 mp *mp_add(mp *d, const mp *a, const mp *b)
 147 {
 148   MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1);
 149   if (!((a->f ^ b->f) & MP_NEG))
 150     mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl);
 151   else {
 152     if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
 153       const mp *t = a; a = b; b = t;
 154     }
 155     mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl);
 156   }
 157   d->f = ((a->f | b->f) & MP_BURN) | (a->f & MP_NEG);
 158   MP_SHRINK(d);
 159   return (d);
 160 }
 161
 162 /* --- @mp_sub@ --- *
 163  *
 164  * Arguments:   @mp *d@ = destination
 165  *              @const mp *a, *b@ = sources
 166  *
 167  * Returns:     Result, @b@ subtracted from @a@.
 168  */
 169
 170 mp *mp_sub(mp *d, const mp *a, const mp *b)
 171 {
 172   unsigned sgn = 0;
 173   MP_MODIFY(d, (MP_LEN(a) > MP_LEN(b) ? MP_LEN(a) : MP_LEN(b)) + 1);
 174   if ((a->f ^ b->f) & MP_NEG)
 175     mpx_uadd(d->v, d->vl, a->v, a->vl, b->v, b->vl);
 176   else {
 177     if (MPX_UCMP(a->v, a->vl, <, b->v, b->vl)) {
 178       const mp *t = a; a = b; b = t;
 179       sgn = MP_NEG;
 180     }
 181     mpx_usub(d->v, d->vl, a->v, a->vl, b->v, b->vl);
 182   }
 183   d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ sgn) & MP_NEG);
 184   MP_SHRINK(d);
 185   return (d);
 186 }
 187
 188 /* --- @mp_mul@ --- *
 189  *
 190  * Arguments:   @mp *d@ = destination
 191  *              @const mp *a, *b@ = sources
 192  *
 193  * Returns:     Result, @a@ multiplied by @b@.
 194  */
 195
 196 mp *mp_mul(mp *d, const mp *a, const mp *b)
 197 {
 198   if (d == a || d == b)
 199     d = MP_NEW;
 200   MP_MODIFY(d, MP_LEN(a) + MP_LEN(b));
 201   mpx_umul(d->v, d->vl, a->v, a->vl, b->v, b->vl);
 202   d->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
 203   MP_SHRINK(d);
 204   return (d);
 205 }
 206
 207 /* --- @mp_sqr@ --- *
 208  *
 209  * Arguments:   @mp *d@ = destination
 210  *              @const mp *a@ = source
 211  *
 212  * Returns:     Result, @a@ squared.
 213  */
 214
 215 mp *mp_sqr(mp *d, const mp *a)
 216 {
 217   if (d == a)
 218     d = MP_NEW;
 219   MP_MODIFY(d, 2 * MP_LEN(a));
 220   mpx_usqr(d->v, d->vl, a->v, a->vl);
 221   d->f = a->f & MP_BURN;
 222   MP_SHRINK(d);
 223   return (d);
 224 }
 225
 226 /* --- @mp_div@ --- *
 227  *
 228  * Arguments:   @mp **qq, **rr@ = destination, quotient and remainder
 229  *              @const mp *a, *b@ = sources
 230  *
 231  * Use:         Calculates the quotient and remainder when @a@ is divided by
 232  *              @b@.  The destinations @*qq@ and @*rr@ must be distinct.
 233  *              Either of @qq@ or @rr@ may be null to indicate that the
 234  *              result is irrelevant.  (Discarding both results is silly.)
 235  *              There is a performance advantage if @a == *rr@.
 236  *
 237  *              The behaviour when @a@ and @b@ have the same sign is
 238  *              straightforward.  When the signs differ, this implementation
 239  *              chooses @r@ to have the same sign as @b@, rather than the
 240  *              more normal choice that the remainder has the same sign as
 241  *              the dividend.  This makes modular arithmetic a little more
 242  *              straightforward.
 243  */
 244
 245 void mp_div(mp **qq, mp **rr, const mp *a, const mp *b)
 246  {
 247   mp *r = rr ? *rr : MP_NEW;
 248   mp *q = qq ? *qq : MP_NEW;
 249   mpw *sv, *svl;
 250
 251   /* --- Set up some temporary workspace --- */
 252
 253   {
 254     size_t rq = MP_LEN(b) + 1;
 255     sv = MP_ALLOC(rq);
 256     svl = sv + rq;
 257   }
 258
 259   /* --- Set the remainder up right --- *
 260    *
 261    * Just in case the divisor is larger, be able to cope with this.  It's not
 262    * important in @mpx_udiv@, but it is here because of the sign correction.
 263    */
 264
 265   {
 266     size_t rq = MP_LEN(a) + 2;
 267     if (MP_LEN(b) > rq)
 268       rq = MP_LEN(b);
 269
 270     if (r == a) {
 271       MP_SPLIT(r);
 272       MP_ENSURE(r, MP_LEN(r) + 2);
 273     } else {
 274       if (r == b)
 275         r = MP_NEW;
 276       MP_MODIFY(r, MP_LEN(a) + 2);
 277       memcpy(r->v, a->v, MPWS(MP_LEN(a)));
 278       memset(r->v + MP_LEN(a), 0, MPWS(2));
 279     }
 280   }
 281
 282   /* --- Fix up the quotient too --- */
 283
 284   if (q == a || q == b)
 285     q = MP_NEW;
 286   MP_MODIFY(q, MP_LEN(a));
 287
 288   /* --- Perform the calculation --- */
 289
 290   mpx_udiv(q->v, q->vl, r->v, r->vl, b->v, b->vl, sv, svl);
 291
 292   /* --- Sort out the sign of the results --- *
 293    *
 294    * If the signs of the arguments differ, and the remainder is nonzero, I
 295    * must add one to the absolute value of the quotient and subtract the
 296    * remainder from @b@.
 297    */
 298
 299   q->f = ((a->f | b->f) & MP_BURN) | ((a->f ^ b->f) & MP_NEG);
 300   if (q->f & MP_NEG) {
 301     mpw *v = r->v;
 302     while (v < r->vl) {
 303       if (*v) {
 304         MPX_UADDN(q->v, q->vl, 1);
 305         mpx_usub(r->v, r->vl, b->v, b->vl, r->v, r->vl);
 306         break;
 307       }
 308     }
 309   }
 310
 311   r->f = ((a->f | b->f) & MP_BURN) | (b->f & MP_NEG);
 312
 313   /* --- Store the return values --- */
 314
 315   if (!qq)
 316     MP_DROP(q);
 317   else {
 318     MP_SHRINK(q);
 319     *qq = q;
 320   }
 321
 322   if (!rr)
 323     MP_DROP(r);
 324   else {
 325     MP_SHRINK(r);
 326     *rr = r;
 327   }
 328
 329   MP_FREE(sv);
 330 }
 331
 332 /*----- Test rig ----------------------------------------------------------*/
 333
 334 #ifdef TEST_RIG
 335
 336 static int verify(const char *op, mp *expect, mp *result, mp *a, mp *b)
 337 {
 338   if (MP_CMP(expect, !=, result)) {
 339     fprintf(stderr, "\n*** %s failed", op);
 340     fputs("\n*** a      = ", stderr); mp_writefile(a, stderr, 10);
 341     fputs("\n*** b      = ", stderr); mp_writefile(b, stderr, 10);
 342     fputs("\n*** result = ", stderr); mp_writefile(result, stderr, 10);
 343     fputs("\n*** expect = ", stderr); mp_writefile(expect, stderr, 10);
 344     fputc('\n', stderr);
 345     return (0);
 346   }
 347   return (1);
 348 }
 349
 350 #define RIG(name, op)                                                   \
 351   static int t ## name(dstr *v)                                         \
 352   {                                                                     \
 353     mp *a = *(mp **)v[0].buf;                                           \
 354     mpw n = *(int *)v[1].buf;                                           \
 355     mp b;                                                               \
 356     mp *r = *(mp **)v[2].buf;                                           \
 357     mp *c = op(MP_NEW, a, n);                                           \
 358     int ok;                                                             \
 359     mp_build(&b, &n, &n + 1);                                           \
 360     ok = verify(#name, r, c, a, &b);                                    \
 361     mp_drop(a); mp_drop(c); mp_drop(r);                                 \
 362     return (ok);                                                        \
 363   }
 364
 365 RIG(lsl, mp_lsl)
 366 RIG(lsr, mp_lsr)
 367
 368 #undef RIG
 369
 370 #define RIG(name, op)                                                   \
 371   static int t ## name(dstr *v)                                         \
 372   {                                                                     \
 373     mp *a = *(mp **)v[0].buf;                                           \
 374     mp *b = *(mp **)v[1].buf;                                           \
 375     mp *r = *(mp **)v[2].buf;                                           \
 376     mp *c = op(MP_NEW, a, b);                                           \
 377     int ok = verify(#name, r, c, a, b);                                 \
 378     mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(r);                     \
 379     return (ok);                                                        \
 380   }
 381
 382 RIG(add, mp_add)
 383 RIG(sub, mp_sub)
 384 RIG(mul, mp_mul)
 385
 386 #undef RIG
 387
 388 static int tdiv(dstr *v)
 389 {
 390   mp *a = *(mp **)v[0].buf;
 391   mp *b = *(mp **)v[1].buf;
 392   mp *q = *(mp **)v[2].buf;
 393   mp *r = *(mp **)v[3].buf;
 394   mp *c = MP_NEW, *d = MP_NEW;
 395   int ok = 1;
 396   mp_div(&c, &d, a, b);
 397   ok &= verify("div(quotient)", q, c, a, b);
 398   ok &= verify("div(remainder)", r, d, a, b);
 399   mp_drop(a); mp_drop(b); mp_drop(c); mp_drop(d); mp_drop(r); mp_drop(q);
 400   return (ok);
 401 }
 402
 403 static test_chunk tests[] = {
 404   { "lsl", tlsl, { &type_mp, &type_mp, &type_mp, 0 } },
 405   { "lsr", tlsr, { &type_mp, &type_mp, &type_mp, 0 } },
 406   { "add", tadd, { &type_mp, &type_mp, &type_mp, 0 } },
 407   { "sub", tsub, { &type_mp, &type_mp, &type_mp, 0 } },
 408   { "mul", tmul, { &type_mp, &type_mp, &type_mp, 0 } },
 409   { "div", tdiv, { &type_mp, &type_mp, &type_mp, &type_mp, 0 } },
 410   { 0, 0, { 0 } },
 411 };
 412
 413 int main(int argc, char *argv[])
 414 {
 415   sub_init();
 416   test_run(argc, argv, tests, SRCDIR "/tests/mp");
 417   return (0);
 418 }
 419
 420 #endif
 421
 422 /*----- That's all, folks -------------------------------------------------*/