Return to dgcd.c CVS log

Up to [Research Unix] / researchv10no / cmd / lfactor

Annotation of researchv10no/cmd/lfactor/dgcd.c, revision 1.1.1.1

1.1       root        1: /*
                      2: double
                      3: dgcd(a,b)
                      4: double a,b;
                      5: {
                      6:        double n,d,q,r;
                      7: 
                      8:        n = a;
                      9:        d = b;
                     10:        while(d != 0){
                     11:                modf(n/d, &q);
                     12:                r = n - q*d;
                     13:                n = d;
                     14:                d = r;
                     15:        }
                     16:        if(n<0) n = -n;
                     17:        return(n);
                     18: }
                     19: */
                     20: 
                     21: /*
                     22:  * dinv computes and returns the g.c.d g = (a,b)
                     23:  * and then provides numbers c,d such that
                     24:  * a*c + b*d = g
                     25:  *
                     26:  * cf. Gauss D.A. Sect. 27.
                     27:  */
                     28: double
                     29: dinv(a, b, c, d)
                     30: double a, b;
                     31: double *c, *d;
                     32: {
                     33:        double num, den, q, r;
                     34:        double temp;
                     35:        double oldx = 0;
                     36:        double newx = 1;
                     37:        double oldy = 1;
                     38:        double newy = 0;
                     39: 
                     40:        num = b;
                     41:        den = a;
                     42:        while(den != 0){
                     43:                modf(num/den, &q);
                     44:                r = num - q*den;
                     45:                temp = newx;
                     46:                newx = oldx - q*newx;
                     47:                oldx = temp;
                     48:                temp = newy;
                     49:                newy = oldy - q*newy;
                     50:                oldy = temp;
                     51:                num = den;
                     52:                den = r;
                     53:        }
                     54:        if(num <0){
                     55:                num = -num;
                     56:                oldx = -oldx;
                     57:                oldy = -oldy;
                     58:        }
                     59:        *c = oldx;
                     60:        *d = oldy;
                     61:        return(num);
                     62: }

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