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

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.