Return to util.c CVS log

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

Annotation of researchv10no/cmd/lfactor/util.c, revision 1.1

1.1     ! root        1: #define testing 0
        !             2: /*
        !             3:  * compos performs the composition of two reduced_______ forms.
        !             4:  * Thus we are entitled to assume that
        !             5:  * |b| <= a <= c
        !             6:  * and
        !             7:  * a,b < 1e16
        !             8:  */
        !             9: #include "mp.h"
        !            10: #include "form.h"
        !            11: form
        !            12: compos(f1,f2)
        !            13: form f1, f2;
        !            14: {
        !            15:        form f3;
        !            16:        mint m, r, y;
        !            17:        mint yq, yr;
        !            18:        mint v1, v2;
        !            19:        mint d, d1, s, n;
        !            20:        mint x1, x2, y1, y2;
        !            21:        mint temp, temp1, temp2;
        !            22:        mint zero;
        !            23:        double dinv();
        !            24:        double modf();
        !            25: 
        !            26:        zero = itom(0);
        !            27:        n.high = s.high = 0;
        !            28:        s.low = (f1.b.low + f2.b.low)*0.5;
        !            29:        n.low = s.low - f2.b.low;
        !            30: #if testing
        !            31:        temp1 = msub(mult(f1.c,f1.a),mult(f2.c,f2.a));
        !            32:        if(mcmp(temp1,mult(n,s)) != 0){
        !            33:                printf("compos: unequal determinants!\n");
        !            34:                tobort(5);
        !            35:        }
        !            36: #endif
        !            37:        d.high = 0;
        !            38:        x1.high = y1.high = 0;
        !            39:        d.low = dinv(f1.a.low,f2.a.low,&x1.low,&y1.low);
        !            40:        d1.high = 0;
        !            41:        x2.high = 0;
        !            42:        y2.high = 0;
        !            43:        v1.high = v2.high = 0;
        !            44:        if(d.low == 1){
        !            45:                d1.low = 1;
        !            46:                v1.low = f1.a.low;
        !            47:                v2.low = f2.a.low;
        !            48:                m = mult(y1,n);
        !            49:        }else{
        !            50:                d1.low = dinv(s.low,d.low,&x2.low,&y2.low);
        !            51:                v1.low = f1.a.low/d1.low;
        !            52:                v2.low = f2.a.low/d1.low;
        !            53:                mdiv(mult(y1,y2), v1, &temp);
        !            54:                temp1 = mult(temp, n);
        !            55:                temp2 = mult(x2, f2.c);
        !            56:                m = msub(temp1, temp2);
        !            57:        }
        !            58:        mdiv(m, v1, &r);
        !            59:        if((r.low + r.low) > v1.low)
        !            60:                r.low -= v1.low;
        !            61:        else if((r.low+r.low) < -v1.low)
        !            62:                r.low += v1.low;
        !            63:        yq = mdiv(mult(r,r), v1, &yr);
        !            64:        temp = mult(f2.c, d1);
        !            65:        temp = madd(temp, mult(r,f2.b));
        !            66:        y = madd(temp,mult(v2,yr));
        !            67:        temp = mdiv(y, v1, &temp2);
        !            68:        if(mcmp(temp2,zero) != 0) tobort(10);
        !            69:        f3.a = mult(v1, v2);
        !            70:        f3.b = madd(f2.b, smult(mult(v2,r),2));
        !            71:        f3.c = madd(temp, mult(yq,v2));
        !            72:        reduce(&f3);
        !            73:        return(f3);;
        !            74: }
        !            75: void
        !            76: reduce(f)
        !            77: form *f;
        !            78: {
        !            79:        mint babs;
        !            80:        mint zero, one;
        !            81:        mint mtemp, k;
        !            82:        mint ka;
        !            83:        double dt1, dt2, dtemp;
        !            84: 
        !            85:        zero = itom(0);
        !            86:        one = itom(1);
        !            87: back:
        !            88:        if(mcmp(f->a,f->c) > 0){
        !            89:                mtemp = f->a;
        !            90:                f->a = f->c;
        !            91:                f->c = mtemp;
        !            92:                f->b = mchs(f->b);
        !            93:        }
        !            94:        babs = f->b;
        !            95:        if(mcmp(f->b,zero) < 0)
        !            96:                babs = mchs(babs);
        !            97:        if(mcmp(f->a,babs) >= 0){
        !            98:                if((mcmp(f->a,f->c)==0) || (mcmp(f->a,babs)==0)){
        !            99:                        f->b = babs;
        !           100:                }
        !           101:                return;
        !           102:        }
        !           103:        dt1 = f->b.high*1e16 + f->b.low;
        !           104:        dt2 = f->a.high*1e16 + f->a.low;
        !           105:        modf(dt1/dt2, &dtemp);
        !           106:        if(dtemp>0) dtemp += 1;
        !           107:        else if(dtemp<0) dtemp -= 1;
        !           108:        modf(dtemp/2, &dtemp);
        !           109:        modf(dtemp/1e16, &k.high);
        !           110:        k.low = dtemp - 1e16*k.high;
        !           111:        ka = mult(k,f->a);
        !           112:        f->b = msub(f->b,madd(ka,ka));
        !           113:        f->c = msub(f->c,mult(k,madd(ka,f->b)));
        !           114:        goto back;
        !           115: }
        !           116: tobort(n)
        !           117: int n;
        !           118: {
        !           119:        printf("Error %d\n", n);
        !           120:        abort();
        !           121: }
        !           122: 
        !           123: form
        !           124: formpow(f,h)
        !           125: form f;
        !           126: double h;
        !           127: {
        !           128:        form ftemp;
        !           129:        int i,n;
        !           130:        int bit[120];
        !           131:        double junk;
        !           132:        double temph;
        !           133:        double modf();
        !           134: 
        !           135:        if(h == 0.){
        !           136:                ftemp = f;
        !           137:                ftemp.b = mchs(ftemp.b);
        !           138:                ftemp = compos(f, ftemp);
        !           139:                return(ftemp);
        !           140:        }
        !           141: 
        !           142:        i = 0;
        !           143:        temph = h;
        !           144:        while(temph != 0){
        !           145:                modf(temph/2, &junk);
        !           146:                if(2.*junk != temph)
        !           147:                        bit[i++] = 1;
        !           148:                else
        !           149:                        bit[i++] = 0;
        !           150:                temph = junk;
        !           151:        }
        !           152:        n = i-1;
        !           153: 
        !           154:        ftemp = f;
        !           155:        for(i=n; i>0; i--){
        !           156:                ftemp = compos(ftemp, ftemp);
        !           157:                if(bit[i-1] != 0){
        !           158:                        ftemp = compos(f, ftemp);
        !           159:                }
        !           160:        }
        !           161:        return(ftemp);
        !           162: }
        !           163: void
        !           164: formout(f)
        !           165: form f;
        !           166: {
        !           167:        printf("(");
        !           168:        mout1(f.a);
        !           169:        printf(",");
        !           170:        mout1(f.b);
        !           171:        printf(",");
        !           172:        mout1(f.c);
        !           173:        printf(")");
        !           174:        return;
        !           175: }

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