home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / unix / volume26 / calc / part05 / qmod.c < prev   
C/C++ Source or Header  |  1992-05-09  |  11KB  |  485 lines
  1. /*
  2.  * Copyright (c) 1992 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Modular arithmetic routines for normal numbers, and also using
  7.  * the faster REDC algorithm.
  8.  */
  9.  
  10. #include <stdio.h>
  11. #include "math.h"
  12.  
  13.  
  14. /*
  15.  * Structure used for caching REDC information.
  16.  */
  17. typedef struct    {
  18.     NUMBER    *num;        /* modulus being cached */
  19.     REDC    *redc;        /* REDC information for modulus */
  20.     long    age;        /* age counter for reallocation */
  21. } REDC_CACHE;
  22.  
  23.  
  24. static long redc_age;            /* current age counter */
  25. static REDC_CACHE redc_cache[MAXREDC];    /* cached REDC info */
  26.  
  27.  
  28. static REDC *qfindredc();
  29.  
  30.  
  31. /*
  32.  * Return the remainder when one number is divided by another.
  33.  * The second argument cannot be negative.  The result is normalized
  34.  * to lie in the range 0 to q2, and works for fractional values as:
  35.  *    qmod(q1, q2) = q1 - (int(q1 / q2) * q2).
  36.  */
  37. NUMBER *
  38. qmod(q1, q2)
  39.     register NUMBER *q1, *q2;
  40. {
  41.     ZVALUE res;
  42.     NUMBER *q, *tmp;
  43.  
  44.     if (qisneg(q2) || qiszero(q2))
  45.         error("Non-positive modulus");
  46.     if (qisint(q1) && qisint(q2)) {        /* easy case */
  47.         zmod(q1->num, q2->num, &res);
  48.         if (iszero(res)) {
  49.             freeh(res.v);
  50.             return qlink(&_qzero_);
  51.         }
  52.         if (isone(res)) {
  53.             freeh(res.v);
  54.             return qlink(&_qone_);
  55.         }
  56.         q = qalloc();
  57.         q->num = res;
  58.         return q;
  59.     }
  60.     q = qquo(q1, q2);
  61.     tmp = qmul(q, q2);
  62.     qfree(q);
  63.     q = qsub(q1, tmp);
  64.     qfree(tmp);
  65.     if (qisneg(q)) {
  66.         tmp = qadd(q2, q);
  67.         qfree(q);
  68.         q = tmp;
  69.     }
  70.     return q;
  71. }
  72.  
  73.  
  74. /*
  75.  * Given two numbers (a and b), calculate the quotient (c) and remainder (d)
  76.  * when a is divided by b.  This is defined so 0 <= d < b, and c is integral,
  77.  * and a = b * c + d.  The results are returned indirectly through pointers.
  78.  * This works for integers or fractions.  Returns whether or not there is a
  79.  * remainder.  Examples:
  80.  *    qquomod(11, 4, &x, &y) sets x to 2, y to 3, and returns TRUE.
  81.  *    qquomod(-7, 3, &x, &y) sets x to -3, y to 2, and returns TRUE.
  82.  */
  83. BOOL
  84. qquomod(q1, q2, retqdiv, retqmod)
  85.     NUMBER *q1, *q2;        /* numbers to do quotient with */
  86.     NUMBER **retqdiv;        /* returned quotient */
  87.     NUMBER **retqmod;        /* returned modulo */
  88. {
  89.     NUMBER *qq, *qm, *tmp;
  90.  
  91.     if (qisneg(q2) || qiszero(q2))
  92.         error("Non-positive modulus");
  93.  
  94.     if (qisint(q1) && qisint(q2)) {        /* integer case */
  95.         qq = qalloc();
  96.         qm = qalloc();
  97.         zdiv(q1->num, q2->num, &qq->num, &qm->num);
  98.         if (!qisneg(q1) || qiszero(qm)) {
  99.             *retqdiv = qq;
  100.             *retqmod = qm;
  101.             return !qiszero(qm);
  102.         }
  103.  
  104.         /*
  105.          * Need to fix up negative results.
  106.          */
  107.         tmp = qdec(qq);
  108.         qfree(qq);
  109.         qq = tmp;
  110.         tmp = qsub(q2, qm);
  111.         qfree(qm);
  112.         qm = tmp;
  113.         *retqdiv = qq;
  114.         *retqmod = qm;
  115.         return TRUE;
  116.     }
  117.  
  118.     /*
  119.      * Fractional case.
  120.      */
  121.     qq = qquo(q1, q2);
  122.     tmp = qmul(qq, q2);
  123.     qm = qsub(q1, tmp);
  124.     qfree(tmp);
  125.     if (qisneg(qm)) {
  126.         tmp = qadd(qm, q2);
  127.         qfree(qm);
  128.         qm = tmp;
  129.         tmp = qdec(qq);
  130.         qfree(qq);
  131.         qq = tmp;
  132.     }
  133.     *retqdiv = qq;
  134.     *retqmod = qm;
  135.     return !qiszero(qm);
  136. }
  137.  
  138.  
  139. #if 0
  140. /*
  141.  * Return the product of two integers modulo a third integer.
  142.  * The result is in the range 0 to q3 - 1 inclusive.
  143.  *    q4 = (q1 * q2) mod q3.
  144.  */
  145. NUMBER *
  146. qmulmod(q1, q2, q3)
  147.     NUMBER *q1, *q2, *q3;
  148. {
  149.     NUMBER *q;
  150.  
  151.     if (qisneg(q3) || qiszero(q3))
  152.         error("Non-positive modulus");
  153.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  154.         error("Non-integers for qmulmod");
  155.     if (qiszero(q1) || qiszero(q2) || qisunit(q3))
  156.         return qlink(&_qzero_);
  157.     q = qalloc();
  158.     zmulmod(q1->num, q2->num, q3->num, &q->num);
  159.     return q;
  160. }
  161.  
  162.  
  163. /*
  164.  * Return the square of an integer modulo another integer.
  165.  * The result is in the range 0 to q2 - 1 inclusive.
  166.  *    q2 = (q1^2) mod q2.
  167.  */
  168. NUMBER *
  169. qsquaremod(q1, q2)
  170.     NUMBER *q1, *q2;
  171. {
  172.     NUMBER *q;
  173.  
  174.     if (qisneg(q2) || qiszero(q2))
  175.         error("Non-positive modulus");
  176.     if (qisfrac(q1) || qisfrac(q2))
  177.         error("Non-integers for qsquaremod");
  178.     if (qiszero(q1) || qisunit(q2))
  179.         return qlink(&_qzero_);
  180.     if (qisunit(q1))
  181.         return qlink(&_qone_);
  182.     q = qalloc();
  183.     zsquaremod(q1->num, q2->num, &q->num);
  184.     return q;
  185. }
  186.  
  187.  
  188. /*
  189.  * Return the sum of two integers modulo a third integer.
  190.  * The result is in the range 0 to q3 - 1 inclusive.
  191.  *    q4 = (q1 + q2) mod q3.
  192.  */
  193. NUMBER *
  194. qaddmod(q1, q2, q3)
  195.     NUMBER *q1, *q2, *q3;
  196. {
  197.     NUMBER *q;
  198.  
  199.     if (qisneg(q3) || qiszero(q3))
  200.         error("Non-positive modulus");
  201.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  202.         error("Non-integers for qaddmod");
  203.     q = qalloc();
  204.     zaddmod(q1->num, q2->num, q3->num, &q->num);
  205.     return q;
  206. }
  207.  
  208.  
  209. /*
  210.  * Return the difference of two integers modulo a third integer.
  211.  * The result is in the range 0 to q3 - 1 inclusive.
  212.  *    q4 = (q1 - q2) mod q3.
  213.  */
  214. NUMBER *
  215. qsubmod(q1, q2, q3)
  216.     NUMBER *q1, *q2, *q3;
  217. {
  218.     NUMBER *q;
  219.  
  220.     if (qisneg(q3) || qiszero(q3))
  221.         error("Non-positive modulus");
  222.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  223.         error("Non-integers for qsubmod");
  224.     if (q1 == q2)
  225.         return qlink(&_qzero_);
  226.     q = qalloc();
  227.     zsubmod(q1->num, q2->num, q3->num, &q->num);
  228.     return q;
  229. }
  230.  
  231.  
  232. /*
  233.  * Return the negative of an integer modulo another integer.
  234.  * The result is in the range 0 to q2 - 1 inclusive.
  235.  *    q2 = (-q1) mod q2.
  236.  */
  237. NUMBER *
  238. qnegmod(q1, q2)
  239.     NUMBER *q1, *q2;
  240. {
  241.     NUMBER *q;
  242.  
  243.     if (qisneg(q2) || qiszero(q2))
  244.         error("Non-positive modulus");
  245.     if (qisfrac(q1) || qisfrac(q2))
  246.         error("Non-integers for qnegmod");
  247.     if (qiszero(q1) || qisunit(q2))
  248.         return qlink(&_qzero_);
  249.     q = qalloc();
  250.     znegmod(q1->num, q2->num, &q->num);
  251.     return q;
  252. }
  253. #endif
  254.  
  255.  
  256. /*
  257.  * Return the integer congruent to an integer whose absolute value is smallest.
  258.  * This is a unique integer in the range int((q2-1)/2 to int(q2/2), inclusive.
  259.  * For example, for a modulus of 7, returned values are [-3, 3], and for a
  260.  * modulus of 8, returned values are [-3, 4].
  261.  */
  262. NUMBER *
  263. qminmod(q1, q2)
  264.     NUMBER *q1, *q2;
  265. {
  266.     NUMBER *q;
  267.  
  268.     if (qisneg(q2) || qiszero(q2))
  269.         error("Non-positive modulus");
  270.     if (qisfrac(q1) || qisfrac(q2))
  271.         error("Non-integers for qminmod");
  272.     if (qiszero(q1) || (q1->num.len < q2->num.len - 1))
  273.         return qlink(q1);
  274.     q = qalloc();
  275.     zminmod(q1->num, q2->num, &q->num);
  276.     return q;
  277. }
  278.  
  279.  
  280. /*
  281.  * Return whether or not two integers are congruent modulo a third integer.
  282.  * Returns TRUE if the numbers are not congruent, and FALSE if they are.
  283.  */
  284. BOOL
  285. qcmpmod(q1, q2, q3)
  286.     NUMBER *q1, *q2, *q3;
  287. {
  288.     if (qisneg(q3) || qiszero(q3))
  289.         error("Non-positive modulus");
  290.     if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3))
  291.         error("Non-integers for qcmpmod");
  292.     if (q1 == q2)
  293.         return FALSE;
  294.     return zcmpmod(q1->num, q2->num, q3->num);
  295. }
  296.  
  297.  
  298. /*
  299.  * Convert an integer into REDC format for use in faster modular arithmetic.
  300.  * The number can be negative or out of modulus range.
  301.  */
  302. NUMBER *
  303. qredcin(q1, q2)
  304.     NUMBER *q1;        /* number to convert into REDC format */
  305.     NUMBER *q2;        /* modulus */
  306. {
  307.     REDC *rp;        /* REDC information */
  308.     NUMBER *r;        /* result */
  309.  
  310.     if (qisfrac(q1))
  311.         error("Non-integer for qredcin");
  312.     rp = qfindredc(q2);
  313.     if (qiszero(q1))
  314.         return qlink(&_qzero_);
  315.     r = qalloc();
  316.     zredcencode(rp, q1->num, &r->num);
  317.     return r;
  318. }
  319.  
  320.  
  321. /*
  322.  * Convert a REDC format number back into a normal integer.
  323.  * The resulting number is in the range 0 to the modulus - 1.
  324.  */
  325. NUMBER *
  326. qredcout(q1, q2)
  327.     NUMBER *q1;        /* number to convert out of REDC format */
  328.     NUMBER *q2;        /* modulus */
  329. {
  330.     REDC *rp;        /* REDC information */
  331.     NUMBER *r;        /* result */
  332.  
  333.     if (qisfrac(q1) || qisneg(q1))
  334.         error("Non-positive integer required for qredcout");
  335.     rp = qfindredc(q2);
  336.     if (qiszero(q1))
  337.         return qlink(&_qzero_);
  338.     r = qalloc();
  339.     zredcdecode(rp, q1->num, &r->num);
  340.     if (isunit(r->num)) {
  341.         qfree(r);
  342.         r = qlink(&_qone_);
  343.     }
  344.     return r;
  345. }
  346.  
  347.  
  348. /*
  349.  * Multiply two REDC format numbers together producing a REDC format result.
  350.  * This multiplication is done modulo the specified modulus.
  351.  */
  352. NUMBER *
  353. qredcmul(q1, q2, q3)
  354.     NUMBER *q1, *q2;    /* REDC numbers to be multiplied */
  355.     NUMBER *q3;        /* modulus */
  356. {
  357.     REDC *rp;        /* REDC information */
  358.     NUMBER *r;        /* result */
  359.  
  360.     if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
  361.         error("Non-positive integers required for qredcmul");
  362.     rp = qfindredc(q3);
  363.     if (qiszero(q1) || qiszero(q2))
  364.         return qlink(&_qzero_);
  365.     r = qalloc();
  366.     zredcmul(rp, q1->num, q2->num, &r->num);
  367.     return r;
  368. }
  369.  
  370.  
  371. /*
  372.  * Square a REDC format number to produce a REDC format result.
  373.  * This squaring is done modulo the specified modulus.
  374.  */
  375. NUMBER *
  376. qredcsquare(q1, q2)
  377.     NUMBER *q1;        /* REDC number to be squared */
  378.     NUMBER *q2;        /* modulus */
  379. {
  380.     REDC *rp;        /* REDC information */
  381.     NUMBER *r;        /* result */
  382.  
  383.     if (qisfrac(q1) || qisneg(q1))
  384.         error("Non-positive integer required for qredcsquare");
  385.     rp = qfindredc(q2);
  386.     if (qiszero(q1))
  387.         return qlink(&_qzero_);
  388.     r = qalloc();
  389.     zredcsquare(rp, q1->num, &r->num);
  390.     return r;
  391. }
  392.  
  393.  
  394. /*
  395.  * Raise a REDC format number to the indicated power producing a REDC
  396.  * format result.  This is done modulo the specified modulus.  The
  397.  * power to be raised to is a normal number.
  398.  */
  399. NUMBER *
  400. qredcpower(q1, q2, q3)
  401.     NUMBER *q1;        /* REDC number to be raised */
  402.     NUMBER *q2;        /* power to be raised to */
  403.     NUMBER *q3;        /* modulus */
  404. {
  405.     REDC *rp;        /* REDC information */
  406.     NUMBER *r;        /* result */
  407.  
  408.     if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2))
  409.         error("Non-positive integers required for qredcpower");
  410.     rp = qfindredc(q3);
  411.     if (qiszero(q1) || qisunit(q3))
  412.         return qlink(&_qzero_);
  413.     r = qalloc();
  414.     zredcpower(rp, q1->num, q2->num, &r->num);
  415.     return r;
  416. }
  417.  
  418.  
  419. /*
  420.  * Search for and return the REDC information for the specified number.
  421.  * The information is cached into a local table so that future calls
  422.  * for this information will be quick.  If the table fills up, then
  423.  * the oldest cached entry is reused.
  424.  */
  425. static REDC *
  426. qfindredc(q)
  427.     NUMBER *q;        /* modulus to find REDC information of */
  428. {
  429.     register REDC_CACHE *rcp;
  430.     REDC_CACHE *bestrcp;
  431.  
  432.     /*
  433.      * First try for an exact pointer match in the table.
  434.      */
  435.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  436.         if (q == rcp->num) {
  437.             rcp->age = ++redc_age;
  438.             return rcp->redc;
  439.         }
  440.     }
  441.  
  442.     /*
  443.      * Search the table again looking for a value which matches.
  444.      */
  445.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  446.         if (rcp->age && (qcmp(q, rcp->num) == 0)) {
  447.             rcp->age = ++redc_age;
  448.             return rcp->redc;
  449.         }
  450.     }
  451.  
  452.     /*
  453.      * Must invalidate an existing entry in the table.
  454.      * Find the oldest (or first unused) entry.
  455.      * But first make sure the modulus will be reasonable.
  456.      */
  457.     if (qisfrac(q) || qiseven(q) || qisneg(q))
  458.         error("REDC modulus must be positive odd integer");
  459.  
  460.     bestrcp = NULL;
  461.     for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) {
  462.         if ((bestrcp == NULL) || (rcp->age < bestrcp->age))
  463.             bestrcp = rcp;
  464.     }
  465.  
  466.     /*
  467.      * Found the best entry.
  468.      * Free the old information for the entry if necessary,
  469.      * then initialize it.
  470.      */
  471.     rcp = bestrcp;
  472.     if (rcp->age) {
  473.         rcp->age = 0;
  474.         qfree(rcp->num);
  475.         zredcfree(rcp->redc);
  476.     }
  477.  
  478.     rcp->redc = zredcalloc(q->num);
  479.     rcp->num = qlink(q);
  480.     rcp->age = ++redc_age;
  481.     return rcp->redc;
  482. }
  483.  
  484. /* END CODE */
  485.