/* * Copyright (c) 1993 David I. Bell * Permission is granted to use, distribute, or modify this source, * provided that this copyright notice remains intact. * * Modular arithmetic routines for normal numbers, and also using * the faster REDC algorithm. */ #include <stdio.h> #include "math.h" /* * Structure used for caching REDC information. */ typedef struct { NUMBER *num; /* modulus being cached */ REDC *redc; /* REDC information for modulus */ long age; /* age counter for reallocation */ } REDC_CACHE; static long redc_age; /* current age counter */ static REDC_CACHE redc_cache[MAXREDC]; /* cached REDC info */ static REDC *qfindredc(); /* * Return the remainder when one number is divided by another. * The second argument cannot be negative. The result is normalized * to lie in the range 0 to q2, and works for fractional values as: * qmod(q1, q2) = q1 - (int(q1 / q2) * q2). */ NUMBER * qmod(q1, q2) register NUMBER *q1, *q2; { ZVALUE res; NUMBER *q, *tmp; if (qisneg(q2) || qiszero(q2)) error("Non-positive modulus"); if (qisint(q1) && qisint(q2)) { /* easy case */ zmod(q1->num, q2->num, &res); if (iszero(res)) { freeh(res.v); return qlink(&_qzero_); } if (isone(res)) { freeh(res.v); return qlink(&_qone_); } q = qalloc(); q->num = res; return q; } q = qquo(q1, q2); tmp = qmul(q, q2); qfree(q); q = qsub(q1, tmp); qfree(tmp); if (qisneg(q)) { tmp = qadd(q2, q); qfree(q); q = tmp; } return q; } /* * Given two numbers (a and b), calculate the quotient (c) and remainder (d) * when a is divided by b. This is defined so 0 <= d < b, and c is integral, * and a = b * c + d. The results are returned indirectly through pointers. * This works for integers or fractions. Returns whether or not there is a * remainder. Examples: * qquomod(11, 4, &x, &y) sets x to 2, y to 3, and returns TRUE. * qquomod(-7, 3, &x, &y) sets x to -3, y to 2, and returns TRUE. */ BOOL qquomod(q1, q2, retqdiv, retqmod) NUMBER *q1, *q2; /* numbers to do quotient with */ NUMBER **retqdiv; /* returned quotient */ NUMBER **retqmod; /* returned modulo */ { NUMBER *qq, *qm, *tmp; if (qisneg(q2) || qiszero(q2)) error("Non-positive modulus"); if (qisint(q1) && qisint(q2)) { /* integer case */ qq = qalloc(); qm = qalloc(); zdiv(q1->num, q2->num, &qq->num, &qm->num); if (!qisneg(q1) || qiszero(qm)) { *retqdiv = qq; *retqmod = qm; return !qiszero(qm); } /* * Need to fix up negative results. */ tmp = qdec(qq); qfree(qq); qq = tmp; tmp = qsub(q2, qm); qfree(qm); qm = tmp; *retqdiv = qq; *retqmod = qm; return TRUE; } /* * Fractional case. */ qq = qquo(q1, q2); tmp = qmul(qq, q2); qm = qsub(q1, tmp); qfree(tmp); if (qisneg(qm)) { tmp = qadd(qm, q2); qfree(qm); qm = tmp; tmp = qdec(qq); qfree(qq); qq = tmp; } *retqdiv = qq; *retqmod = qm; return !qiszero(qm); } #if 0 /* * Return the product of two integers modulo a third integer. * The result is in the range 0 to q3 - 1 inclusive. * q4 = (q1 * q2) mod q3. */ NUMBER * qmulmod(q1, q2, q3) NUMBER *q1, *q2, *q3; { NUMBER *q; if (qisneg(q3) || qiszero(q3)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) error("Non-integers for qmulmod"); if (qiszero(q1) || qiszero(q2) || qisunit(q3)) return qlink(&_qzero_); q = qalloc(); zmulmod(q1->num, q2->num, q3->num, &q->num); return q; } /* * Return the square of an integer modulo another integer. * The result is in the range 0 to q2 - 1 inclusive. * q2 = (q1^2) mod q2. */ NUMBER * qsquaremod(q1, q2) NUMBER *q1, *q2; { NUMBER *q; if (qisneg(q2) || qiszero(q2)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2)) error("Non-integers for qsquaremod"); if (qiszero(q1) || qisunit(q2)) return qlink(&_qzero_); if (qisunit(q1)) return qlink(&_qone_); q = qalloc(); zsquaremod(q1->num, q2->num, &q->num); return q; } /* * Return the sum of two integers modulo a third integer. * The result is in the range 0 to q3 - 1 inclusive. * q4 = (q1 + q2) mod q3. */ NUMBER * qaddmod(q1, q2, q3) NUMBER *q1, *q2, *q3; { NUMBER *q; if (qisneg(q3) || qiszero(q3)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) error("Non-integers for qaddmod"); q = qalloc(); zaddmod(q1->num, q2->num, q3->num, &q->num); return q; } /* * Return the difference of two integers modulo a third integer. * The result is in the range 0 to q3 - 1 inclusive. * q4 = (q1 - q2) mod q3. */ NUMBER * qsubmod(q1, q2, q3) NUMBER *q1, *q2, *q3; { NUMBER *q; if (qisneg(q3) || qiszero(q3)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) error("Non-integers for qsubmod"); if (q1 == q2) return qlink(&_qzero_); q = qalloc(); zsubmod(q1->num, q2->num, q3->num, &q->num); return q; } /* * Return the negative of an integer modulo another integer. * The result is in the range 0 to q2 - 1 inclusive. * q2 = (-q1) mod q2. */ NUMBER * qnegmod(q1, q2) NUMBER *q1, *q2; { NUMBER *q; if (qisneg(q2) || qiszero(q2)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2)) error("Non-integers for qnegmod"); if (qiszero(q1) || qisunit(q2)) return qlink(&_qzero_); q = qalloc(); znegmod(q1->num, q2->num, &q->num); return q; } #endif /* * Return the integer congruent to an integer whose absolute value is smallest. * This is a unique integer in the range int((q2-1)/2 to int(q2/2), inclusive. * For example, for a modulus of 7, returned values are [-3, 3], and for a * modulus of 8, returned values are [-3, 4]. */ NUMBER * qminmod(q1, q2) NUMBER *q1, *q2; { NUMBER *q; if (qisneg(q2) || qiszero(q2)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2)) error("Non-integers for qminmod"); if (qiszero(q1) || (q1->num.len < q2->num.len - 1)) return qlink(q1); q = qalloc(); zminmod(q1->num, q2->num, &q->num); return q; } /* * Return whether or not two integers are congruent modulo a third integer. * Returns TRUE if the numbers are not congruent, and FALSE if they are. */ BOOL qcmpmod(q1, q2, q3) NUMBER *q1, *q2, *q3; { if (qisneg(q3) || qiszero(q3)) error("Non-positive modulus"); if (qisfrac(q1) || qisfrac(q2) || qisfrac(q3)) error("Non-integers for qcmpmod"); if (q1 == q2) return FALSE; return zcmpmod(q1->num, q2->num, q3->num); } /* * Convert an integer into REDC format for use in faster modular arithmetic. * The number can be negative or out of modulus range. */ NUMBER * qredcin(q1, q2) NUMBER *q1; /* number to convert into REDC format */ NUMBER *q2; /* modulus */ { REDC *rp; /* REDC information */ NUMBER *r; /* result */ if (qisfrac(q1)) error("Non-integer for qredcin"); rp = qfindredc(q2); if (qiszero(q1)) return qlink(&_qzero_); r = qalloc(); zredcencode(rp, q1->num, &r->num); return r; } /* * Convert a REDC format number back into a normal integer. * The resulting number is in the range 0 to the modulus - 1. */ NUMBER * qredcout(q1, q2) NUMBER *q1; /* number to convert out of REDC format */ NUMBER *q2; /* modulus */ { REDC *rp; /* REDC information */ NUMBER *r; /* result */ if (qisfrac(q1) || qisneg(q1)) error("Non-positive integer required for qredcout"); rp = qfindredc(q2); if (qiszero(q1)) return qlink(&_qzero_); r = qalloc(); zredcdecode(rp, q1->num, &r->num); if (isunit(r->num)) { qfree(r); r = qlink(&_qone_); } return r; } /* * Multiply two REDC format numbers together producing a REDC format result. * This multiplication is done modulo the specified modulus. */ NUMBER * qredcmul(q1, q2, q3) NUMBER *q1, *q2; /* REDC numbers to be multiplied */ NUMBER *q3; /* modulus */ { REDC *rp; /* REDC information */ NUMBER *r; /* result */ if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2)) error("Non-positive integers required for qredcmul"); rp = qfindredc(q3); if (qiszero(q1) || qiszero(q2)) return qlink(&_qzero_); r = qalloc(); zredcmul(rp, q1->num, q2->num, &r->num); return r; } /* * Square a REDC format number to produce a REDC format result. * This squaring is done modulo the specified modulus. */ NUMBER * qredcsquare(q1, q2) NUMBER *q1; /* REDC number to be squared */ NUMBER *q2; /* modulus */ { REDC *rp; /* REDC information */ NUMBER *r; /* result */ if (qisfrac(q1) || qisneg(q1)) error("Non-positive integer required for qredcsquare"); rp = qfindredc(q2); if (qiszero(q1)) return qlink(&_qzero_); r = qalloc(); zredcsquare(rp, q1->num, &r->num); return r; } /* * Raise a REDC format number to the indicated power producing a REDC * format result. This is done modulo the specified modulus. The * power to be raised to is a normal number. */ NUMBER * qredcpower(q1, q2, q3) NUMBER *q1; /* REDC number to be raised */ NUMBER *q2; /* power to be raised to */ NUMBER *q3; /* modulus */ { REDC *rp; /* REDC information */ NUMBER *r; /* result */ if (qisfrac(q1) || qisneg(q1) || qisfrac(q2) || qisneg(q2)) error("Non-positive integers required for qredcpower"); rp = qfindredc(q3); if (qiszero(q1) || qisunit(q3)) return qlink(&_qzero_); r = qalloc(); zredcpower(rp, q1->num, q2->num, &r->num); return r; } /* * Search for and return the REDC information for the specified number. * The information is cached into a local table so that future calls * for this information will be quick. If the table fills up, then * the oldest cached entry is reused. */ static REDC * qfindredc(q) NUMBER *q; /* modulus to find REDC information of */ { register REDC_CACHE *rcp; REDC_CACHE *bestrcp; /* * First try for an exact pointer match in the table. */ for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { if (q == rcp->num) { rcp->age = ++redc_age; return rcp->redc; } } /* * Search the table again looking for a value which matches. */ for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { if (rcp->age && (qcmp(q, rcp->num) == 0)) { rcp->age = ++redc_age; return rcp->redc; } } /* * Must invalidate an existing entry in the table. * Find the oldest (or first unused) entry. * But first make sure the modulus will be reasonable. */ if (qisfrac(q) || qiseven(q) || qisneg(q)) error("REDC modulus must be positive odd integer"); bestrcp = NULL; for (rcp = redc_cache; rcp < &redc_cache[MAXREDC]; rcp++) { if ((bestrcp == NULL) || (rcp->age < bestrcp->age)) bestrcp = rcp; } /* * Found the best entry. * Free the old information for the entry if necessary, * then initialize it. */ rcp = bestrcp; if (rcp->age) { rcp->age = 0; qfree(rcp->num); zredcfree(rcp->redc); } rcp->redc = zredcalloc(q->num); rcp->num = qlink(q); rcp->age = ++redc_age; return rcp->redc; } /* END CODE */