Return to xgcd.c CVS log

Up to [MW Coherent from dump] / coherent / d / usr / lib / libmp

Annotation of coherent/d/usr/lib/libmp/xgcd.c, revision 1.1

1.1     ! root        1: #include "mprec.h"
        !             2: 
        !             3: 
        !             4: /*
        !             5:  *     Xgcd sets the value of the mint pointed to by "g" to the greatest
        !             6:  *     common divisor of the mints pointed to by "a" and "b".
        !             7:  *     It sets the mints "r" and "s" so that the following are true:
        !             8:  *             g = r * a + s * b
        !             9:  *     Note that the only restriction as to the distinctness
        !            10:  *     of the arguments is that "r", "s" and "g" are all different.
        !            11:  */
        !            12: 
        !            13: void
        !            14: xgcd(a, b, r, s, g)
        !            15: mint *a, *b, *r, *s, *g;
        !            16: {
        !            17:        register mint *temp;
        !            18:        mint    *r1, *s1, *t, *u;
        !            19:        mint    al, bl, quot, t1;
        !            20: 
        !            21:        /* make local copyies of arguments */
        !            22:        minit(&al);
        !            23:        mcopy(a, &al);
        !            24:        a = &al;
        !            25:        minit(&bl);
        !            26:        mcopy(b, &bl);
        !            27:        b = &bl;
        !            28: 
        !            29:        /*
        !            30:         * The following matrix equation will always hold:
        !            31:         *
        !            32:         *      r1 s1        original a     a
        !            33:         *      t  u  times  original b  =  b
        !            34:         */
        !            35: 
        !            36:        r1 = itom(1);
        !            37:        s1 = itom(0);
        !            38:        t  = itom(0);
        !            39:        u  = itom(1);
        !            40: 
        !            41:        /* calculate gcd by Euclidean algorithm */
        !            42:        minit(&quot);
        !            43:        minit(&t1);
        !            44:        while (!zerop(b)) {
        !            45:                /*
        !            46:                 * Perform following matrix assignment:
        !            47:                 *
        !            48:                 *      a  r1 s1        0      1         a  r1 s1
        !            49:                 *      b  t  u =       1  -quot  times  b  t  u
        !            50:                 *
        !            51:                 * (where quot = a/b).  This is done by the factorization:
        !            52:                 *
        !            53:                 *      0      1     0  1         1  -quot
        !            54:                 *      1  -quot  =  1  0  times  0      1
        !            55:                 */
        !            56: 
        !            57:                mdiv(a, b, &quot, a);
        !            58:                mult(&quot, t, &t1);
        !            59:                msub(r1, &t1, r1);
        !            60:                mult(&quot, u, &t1);
        !            61:                msub(s1, &t1, s1);
        !            62:                temp = r1;
        !            63:                r1 = t;
        !            64:                t = temp;
        !            65:                temp = s1;
        !            66:                s1 = u;
        !            67:                u = temp;
        !            68:                temp = a;
        !            69:                a = b;
        !            70:                b = temp;
        !            71:        }
        !            72: 
        !            73:        /* set results */
        !            74:        mpfree(g->val);
        !            75:        *g = *a;
        !            76:        mpfree(r->val);
        !            77:        *r = *r1;
        !            78:        mpfree(s->val);
        !            79:        *s = *s1;
        !            80: 
        !            81:        /* throw away garbage */
        !            82:        mpfree(quot.val);
        !            83:        mpfree(t1.val);
        !            84:        mpfree(b->val);
        !            85:        mpfree(r1);
        !            86:        mpfree(s1);
        !            87:        mintfr(t);
        !            88:        mintfr(u);
        !            89: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.