Return to mdiv.c CVS log

Up to [Research Unix] / researchv10no / libmp

Annotation of researchv10no/libmp/mdiv.c, revision 1.1.1.1

1.1       root        1: #include "mp.h"
                      2: int mpdivdebug;
                      3: mdiv(a,b,q,r) mint *a,*b,*q,*r;
                      4: {      mint *x,*y, *z;
                      5:        int alen, blen;
                      6:        x = itom(0), y = itom(0);
                      7:        z = itom(0);
                      8:        move(a, x);
                      9:        move(b, y);
                     10:        alen = x->len;
                     11:        if(x->len<0) {x->len= -x->len;}
                     12:        blen = y->len;
                     13:        if(y->len<0) {y->len= -y->len;}
                     14:        xfree(q);
                     15:        xfree(r);
                     16:        m_div(x,y,q,r);
                     17:        if(mpdivdebug) {
                     18:                mult(y, q, z);
                     19:                madd(z, r, z);
                     20:                if(mcmp(z, x) != 0) {
                     21:                        mout(a);
                     22:                        mout(b);
                     23:                        fatal("mdiv err");
                     24:                }
                     25:        }
                     26:        if(alen < 0) {
                     27:                mint o;
                     28:                short i = 1;
                     29:                o.len = 1;
                     30:                o.val = &i;
                     31:                if(r->len == 0) {
                     32:                        if(blen > 0)
                     33:                                q->len = -q->len;
                     34:                        goto out;
                     35:                }
                     36:                msub(r, y, r);
                     37:                r->len = -r->len;
                     38:                madd(q, &o, q);
                     39:                if(blen > 0)
                     40:                        q->len = -q->len;
                     41:                goto out;
                     42:        }
                     43:        if(blen < 0)
                     44:                q->len = -q->len;
                     45: out:
                     46:        mfree(z);
                     47:        mfree(y);
                     48:        mfree(x);
                     49: }
                     50: #ifdef vax
                     51: union zz {
                     52:        long a;
                     53:        struct {
                     54:                unsigned int lo:15, hi:15;
                     55:        } b;
                     56: };
                     57: #else
                     58: union zz {
                     59:        long a;
                     60:        struct {
                     61:                unsigned int :2, hi:15,lo:15;
                     62:        } b;
                     63: };
                     64: #endif
                     65: m_dsb(q,n,a,b) short *a,*b;
                     66: {      long int qx, u;
                     67:        union zz x;
                     68:        int borrow,j;
                     69:        qx=q;
                     70:        x.a = 0;
                     71:        for(j = 0; j < n; j++) {
                     72:                x.a = qx * a[j] + x.b.hi;
                     73:                if((b[j] -= x.b.lo) < 0) {
                     74:                        b[j] += (1 << 15);
                     75:                        x.b.hi += 1;
                     76:                }
                     77:        }
                     78:        if((b[j] -= x.b.hi) >= 0)
                     79:                return(0);
                     80:        borrow=0;
                     81:        for(j=0;j<n;j++)
                     82:        {       u=a[j]+b[j]+borrow;
                     83:                if(u & 0100000)
                     84:                        borrow = 1;
                     85:                else borrow=0;
                     86:                b[j]=u&077777;
                     87:        }
                     88:        { return(1);}
                     89: }
                     90: m_trq(v1,v2,u1,u2,u3)
                     91: {      long int d;
                     92:        long int x1;
                     93:        if(u1==v1) d=077777;
                     94:        else d=(u1*0100000L+u2)/v1;
                     95:        while(1)
                     96:        {       x1=u1*0100000L+u2-v1*d;
                     97:                x1=x1*0100000L+u3-v2*d;
                     98:                if(x1<0) d=d-1;
                     99:                else {return(d);}
                    100:        }
                    101: }
                    102: m_div(a,b,q,r) mint *a,*b,*q,*r;
                    103: {      mint u,v,x,w;
                    104:        short d,*qval;
                    105:        int qq,n,v1,v2,j;
                    106:        u.len=v.len=x.len=w.len=0;
                    107:        if(b->len==0) { fatal("mdiv divide by zero"); return;}
                    108:        if(b->len==1)
                    109:        {       r->val=xalloc(1,"m_div1");
                    110:                sdiv(a,b->val[0],q,r->val);
                    111:                if(r->val[0]==0)
                    112:                {       shfree(r->val);
                    113:                        r->len=0;
                    114:                }
                    115:                else r->len=1;
                    116:                return;
                    117:        }
                    118:        if(a->len < b->len)
                    119:        {       q->len=0;
                    120:                r->len=a->len;
                    121:                r->val=xalloc(r->len,"m_div2");
                    122:                for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq];
                    123:                return;
                    124:        }
                    125:        x.len=1;
                    126:        x.val = &d;
                    127:        n=b->len;
                    128:        d=0100000L/(b->val[n-1]+1L);
                    129:        mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */
                    130:        mult(b,&x,&v);
                    131:        v1=v.val[n-1];
                    132:        v2=v.val[n-2];
                    133:        qval=xalloc(a->len-n+1,"m_div3");
                    134:        for(j=a->len-n;j>=0;j--)
                    135:        {       qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]);
                    136:                if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1;
                    137:                qval[j]=qq;
                    138:        }
                    139:        x.len=n;
                    140:        x.val=u.val;
                    141:        mcan(&x);
                    142:        sdiv(&x,d,&w,(short *)&qq);
                    143:        r->len=w.len;
                    144:        r->val=w.val;
                    145:        q->val=qval;
                    146:        qq=a->len-n+1;
                    147:        if(qq>0 && qval[qq-1]==0) qq -= 1;
                    148:        q->len=qq;
                    149:        if(qq==0) shfree(qval);
                    150:        if(x.len!=0) xfree(&u);
                    151:        xfree(&v);
                    152:        return;
                    153: }

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