c_cpp/etc/calc/qmod.c

   1 /*
   2  * qmod - modular arithmetic routines for normal numbers and REDC numbers
   3  *
   4  * Copyright (C) 1999-2007  David I. Bell and Ernest Bowen
   5  *
   6  * Primary author:  David I. Bell
   7  *
   8  * Calc is open software; you can redistribute it and/or modify it under
   9  * the terms of the version 2.1 of the GNU Lesser General Public License
  10  * as published by the Free Software Foundation.
  11  *
  12  * Calc is distributed in the hope that it will be useful, but WITHOUT
  13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  15  * Public License for more details.
  16  *
  17  * A copy of version 2.1 of the GNU Lesser General Public License is
  18  * distributed with calc under the filename COPYING-LGPL.  You should have
  19  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  20  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  21  *
  22  * @(#) $Revision: 30.1 $
  23  * @(#) $Id: qmod.c,v 30.1 2007/03/16 11:09:46 chongo Exp $
  24  * @(#) $Source: /usr/local/src/bin/calc/RCS/qmod.c,v $
  25  *
  26  * Under source code control:   1991/05/22 23:15:07
  27  * File existed as early as:    1991
  28  *
  29  * Share and enjoy!  :-)        http://www.isthe.com/chongo/tech/comp/calc/
  30  */
  31
  32
  33 #include <stdio.h>
  34 #include "qmath.h"
  35 #include "config.h"
  36
  37
  38 /*
  39  * Structure used for caching REDC information.
  40  */
  41 typedef struct  {
  42         NUMBER  *rnum;          /* modulus being cached */
  43         REDC    *redc;          /* REDC information for modulus */
  44         long    age;            /* age counter for reallocation */
  45 } REDC_CACHE;
  46
  47
  48 STATIC long redc_age;                   /* current age counter */
  49 STATIC REDC_CACHE redc_cache[MAXREDC];  /* cached REDC info */
  50
  51
  52 S_FUNC REDC *qfindredc(NUMBER *q);
  53
  54
  55 /*
  56  * qmod(q1, q2, rnd) returns zero if q1 is a multiple of q2; it
  57  * q1 if q2 is zero.  For other q1 and q2, it returns one of
  58  * the two remainders with absolute value less than abs(q2)
  59  * when q1 is divided by q2;  which remainder is returned is
  60  * determined by rnd and the signs and relative sizes of q1 and q2.
  61  */
  62 NUMBER *
  63 qmod(NUMBER *q1, NUMBER *q2, long rnd)
  64 {
  65         ZVALUE tmp, tmp1, tmp2;
  66         NUMBER *q;
  67
  68         if (qiszero(q2)) return qlink(q1);
  69         if (qiszero(q1)) return qlink(&_qzero_);
  70         if (qisint(q1) && qisint(q2)) {         /* easy case */
  71                 zmod(q1->num, q2->num, &tmp, rnd);
  72                 if (ziszero(tmp)) {
  73                         zfree(tmp);
  74                         return qlink(&_qzero_);
  75                 }
  76                 if(zisone(tmp)) {
  77                         zfree(tmp);
  78                         return qlink(&_qone_);
  79                 }
  80                 q = qalloc();
  81                 q->num = tmp;
  82                 return q;
  83         }
  84         zmul(q1->num, q2->den, &tmp1);
  85         zmul(q2->num, q1->den, &tmp2);
  86         zmod(tmp1, tmp2, &tmp, rnd);
  87         zfree(tmp1);
  88         zfree(tmp2);
  89         if (ziszero(tmp)) {
  90                 zfree(tmp);
  91                 return qlink(&_qzero_);
  92         }
  93         zmul(q1->den, q2->den, &tmp1);
  94         q = qalloc();
  95         zreduce(tmp, tmp1, &q->num, &q->den);
  96         zfree(tmp1);
  97         zfree(tmp);
  98         return q;
  99 }
 100
 101
 102 /*
 103  * Given two numbers q1, q2, qquomod(q1, q2, quo, mod, rnd) assigns
 104  * quo(q1, q2, rnd) to quo and mod(q1, q2, rnd) to mod. The quotient
 105  * quo is always an integer, q1 = q2 * quo + mod, and if q2 is non-zero,
 106  * abs(mod) < abs(q2). If q2 is zero, quo is assigned the value zero and
 107  * mod = q1.
 108  * Which of the two possible quotient-remainder pairs is assigned
 109  * to quo and mod is determined by the fifth argument rnd.
 110  * If rnd is zero, the remainder has the sign of q2
 111  * and the quotient is rounded towards zero.
 112  *
 113  * The function returns TRUE or FALSE according as the remainder is
 114  * nonzero or zero
 115  *
 116  * Given:
 117  *      q1              number to be divided
 118  *      q2              divisor
 119  *      quo             quotient
 120  *      mod             remainder
 121  *      rnd             rounding-type specifier
 122  */
 123 BOOL
 124 qquomod(NUMBER *q1, NUMBER *q2, NUMBER **quo, NUMBER **mod, long rnd)
 125 {
 126         NUMBER *qq, *qm;
 127         ZVALUE tmp1, tmp2, tmp3, tmp4;
 128
 129         if (qiszero(q2)) {                      /* zero divisor case */
 130                 qq = qlink(&_qzero_);
 131                 qm = qlink(q1);
 132         } else if (qisint(q1) && qisint(q2)) {  /* integer args case */
 133                 zdiv(q1->num, q2->num, &tmp1, &tmp2, rnd);
 134                 if (ziszero(tmp1)) {
 135                         zfree(tmp1);
 136                         zfree(tmp2);
 137                         qq = qlink(&_qzero_);
 138                         qm = qlink(q1);
 139                 } else {
 140                         qq = qalloc();
 141                         qq->num = tmp1;
 142                         if (ziszero(tmp2)) {
 143                                 zfree(tmp2);
 144                                 qm = qlink(&_qzero_);
 145                         } else {
 146                                 qm = qalloc();
 147                                 qm->num = tmp2;
 148                         }
 149                 }
 150         } else {                                        /* fractional case */
 151                 zmul(q1->num, q2->den, &tmp1);
 152                 zmul(q2->num, q1->den, &tmp2);
 153                 zdiv(tmp1, tmp2, &tmp3, &tmp4, rnd);
 154                 zfree(tmp1);
 155                 zfree(tmp2);
 156                 if (ziszero(tmp3)) {
 157                         zfree(tmp3);
 158                         zfree(tmp4);
 159                         qq = qlink(&_qzero_);
 160                         qm = qlink(q1);
 161                 } else {
 162                         qq = qalloc();
 163                         qq->num = tmp3;
 164                         if (ziszero(tmp4)) {
 165                                 zfree(tmp4);
 166                                 qm = qlink(&_qzero_);
 167                         } else {
 168                                 qm = qalloc();
 169                                 zmul(q1->den, q2->den, &tmp1);
 170                                 zreduce(tmp4, tmp1, &qm->num, &qm->den);
 171                                 zfree(tmp1);
 172                                 zfree(tmp4);
 173                         }
 174                 }
 175         }
 176         *quo = qq;
 177         *mod = qm;
 178         return !qiszero(qm);
 179 }
 180
 181
 182 /*
 183  * Return whether or not two integers are congruent modulo a third integer.
 184  * Returns TRUE if the numbers are not congruent, and FALSE if they are.
 185  */
 186 BOOL
 187 qcmpmod(NUMBER *q1, NUMBER *q2, NUMBER *q3)
 188 {
 189         if (qisneg(q3) || qiszero(q3))
 190                 math_error("Non-positive modulus");
 191         if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
 192                 math_error("Non-integers for qcmpmod");
 193         if (q1 == q2)
 194                 return FALSE;
 195         return zcmpmod(q1->num, q2->num, q3->num);
 196 }
 197
 198
 199 /*
 200  * Convert an integer into REDC format for use in faster modular arithmetic.
 201  * The number can be negative or out of modulus range.
 202  *
 203  * given:
 204  *      q1              number to convert into REDC format
 205  *      q2              modulus
 206  */
 207 NUMBER *
 208 qredcin(NUMBER *q1, NUMBER *q2)
 209 {
 210         REDC *rp;               /* REDC information */
 211         NUMBER *r;              /* result */
 212
 213         if (qisfrac(q1))
 214                 math_error("Non-integer for qredcin");
 215         rp = qfindredc(q2);
 216         r = qalloc();
 217         zredcencode(rp, q1->num, &r->num);
 218         if (qiszero(r)) {
 219                 qfree(r);
 220                 return qlink(&_qzero_);
 221         }
 222         return r;
 223 }
 224
 225
 226 /*
 227  * Convert a REDC format number back into a normal integer.
 228  * The resulting number is in the range 0 to the modulus - 1.
 229  *
 230  * given:
 231  *      q1              number to convert into REDC format
 232  *      q2              modulus
 233  */
 234 NUMBER *
 235 qredcout(NUMBER *q1, NUMBER *q2)
 236 {
 237         REDC *rp;               /* REDC information */
 238         NUMBER *r;              /* result */
 239
 240         if (qisfrac(q1))
 241                 math_error("Non-integer argument for rcout");
 242         rp = qfindredc(q2);
 243         if (qiszero(q1) || qisunit(q2))
 244                 return qlink(&_qzero_);
 245         r = qalloc();
 246         zredcdecode(rp, q1->num, &r->num);
 247         if (zisunit(r->num)) {
 248                 qfree(r);
 249                 r = qlink(&_qone_);
 250         }
 251         return r;
 252 }
 253
 254
 255 /*
 256  * Multiply two REDC format numbers together producing a REDC format result.
 257  * This multiplication is done modulo the specified modulus.
 258  *
 259  * given:
 260  *      q1              REDC numbers to be multiplied
 261  *      q2              REDC numbers to be multiplied
 262  *      q3              modulus
 263  */
 264 NUMBER *
 265 qredcmul(NUMBER *q1, NUMBER *q2, NUMBER *q3)
 266 {
 267         REDC *rp;               /* REDC information */
 268         NUMBER *r;              /* result */
 269
 270         if (qisfrac(q1) || qisfrac(q2))
 271                 math_error("Non-integer argument for rcmul");
 272         rp = qfindredc(q3);
 273         if (qiszero(q1) || qiszero(q2) || qisunit(q3))
 274                 return qlink(&_qzero_);
 275         r = qalloc();
 276         zredcmul(rp, q1->num, q2->num, &r->num);
 277         return r;
 278 }
 279
 280
 281 /*
 282  * Square a REDC format number to produce a REDC format result.
 283  * This squaring is done modulo the specified modulus.
 284  *
 285  * given:
 286  *      q1              REDC numbers to be squared
 287  *      q2              modulus
 288  */
 289 NUMBER *
 290 qredcsquare(NUMBER *q1, NUMBER *q2)
 291 {
 292         REDC *rp;               /* REDC information */
 293         NUMBER *r;              /* result */
 294
 295         if (qisfrac(q1))
 296                 math_error("Non-integer argument for rcsq");
 297         rp = qfindredc(q2);
 298         if (qiszero(q1) || qisunit(q2))
 299                 return qlink(&_qzero_);
 300         r = qalloc();
 301         zredcsquare(rp, q1->num, &r->num);
 302         return r;
 303 }
 304
 305
 306 /*
 307  * Raise a REDC format number to the indicated power producing a REDC
 308  * format result.  This is done modulo the specified modulus.  The
 309  * power to be raised to is a normal number.
 310  *
 311  * given:
 312  *      q1              REDC number to be raised
 313  *      q2              power to be raised to
 314  *      q3              modulus
 315  */
 316 NUMBER *
 317 qredcpower(NUMBER *q1, NUMBER *q2, NUMBER *q3)
 318 {
 319         REDC *rp;               /* REDC information */
 320         NUMBER *r;              /* result */
 321
 322         if (qisfrac(q1) || qisfrac(q2) || qisfrac(q2))
 323                 math_error("Non-integer argument for rcpow");
 324         if (qisneg(q2))
 325                 math_error("Negative exponent argument for rcpow");
 326         rp = qfindredc(q3);
 327         r = qalloc();
 328         zredcpower(rp, q1->num, q2->num, &r->num);
 329         return r;
 330 }
 331
 332
 333 /*
 334  * Search for and return the REDC information for the specified number.
 335  * The information is cached into a local table so that future calls
 336  * for this information will be quick.  If the table fills up, then
 337  * the oldest cached entry is reused.
 338  *
 339  * given:
 340  *      q               modulus to find REDC information of
 341  */
 342 S_FUNC REDC *
 343 qfindredc(NUMBER *q)
 344 {
 345         register REDC_CACHE *rcp;
 346         REDC_CACHE *bestrcp;
 347
 348         /*
 349          * First try for an exact pointer match in the table.
 350          */
 351         for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) {
 352                 if (q == rcp->rnum) {
 353                         rcp->age = ++redc_age;
 354                         return rcp->redc;
 355                 }
 356         }
 357
 358         /*
 359          * Search the table again looking for a value which matches.
 360          */
 361         for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) {
 362                 if (rcp->age && (qcmp(q, rcp->rnum) == 0)) {
 363                         rcp->age = ++redc_age;
 364                         return rcp->redc;
 365                 }
 366         }
 367
 368         /*
 369          * Must invalidate an existing entry in the table.
 370          * Find the oldest (or first unused) entry.
 371          * But first make sure the modulus will be reasonable.
 372          */
 373         if (qisfrac(q) || qisneg(q)) {
 374                 math_error("REDC modulus must be positive odd integer");
 375                 /*NOTREACHED*/
 376         }
 377
 378         bestrcp = NULL;
 379         for (rcp = redc_cache; rcp <= &redc_cache[MAXREDC-1]; rcp++) {
 380                 if ((bestrcp == NULL) || (rcp->age < bestrcp->age))
 381                         bestrcp = rcp;
 382         }
 383
 384         /*
 385          * Found the best entry.
 386          * Free the old information for the entry if necessary,
 387          * then initialize it.
 388          */
 389         rcp = bestrcp;
 390         if (rcp->age) {
 391                 rcp->age = 0;
 392                 qfree(rcp->rnum);
 393                 zredcfree(rcp->redc);
 394         }
 395
 396         rcp->redc = zredcalloc(q->num);
 397         rcp->rnum = qlink(q);
 398         rcp->age = ++redc_age;
 399         return rcp->redc;
 400 }
 401
 402 void
 403 showredcdata(void)
 404 {
 405         REDC_CACHE *rcp;
 406         long i;
 407
 408         for (i = 0, rcp = redc_cache; i < MAXREDC; i++, rcp++) {
 409                 if (rcp->age > 0) {
 410                         printf("%-8ld%-8ld", i, rcp->age);
 411                         qprintnum(rcp->rnum, 0);
 412                         printf("\n");
 413                 }
 414         }
 415 }
 416
 417 void
 418 freeredcdata(void)
 419 {
 420         REDC_CACHE *rcp;
 421         long i;
 422
 423         for (i = 0, rcp = redc_cache; i < MAXREDC; i++, rcp++) {
 424                 if (rcp->age > 0) {
 425                         rcp->age = 0;
 426                         qfree(rcp->rnum);
 427                         zredcfree(rcp->redc);
 428                 }
 429         }
 430 }