1: /* 2: * Copyright (c) 1980 Regents of the University of California. 3: * All rights reserved. The Berkeley software License Agreement 4: * specifies the terms and conditions for redistribution. 5: */ 6: 7: #ifndef lint 8: static char sccsid[] = "@(#)mdiv.c 5.1 (Berkeley) 4/30/85"; 9: #endif not lint 10: 11: #include <mp.h> 12: mdiv(a,b,q,r) MINT *a,*b,*q,*r; 13: { MINT x,y; 14: int sign; 15: sign=1; 16: x.val=a->val; 17: y.val=b->val; 18: x.len=a->len; 19: if(x.len<0) {sign= -1; x.len= -x.len;} 20: y.len=b->len; 21: if(y.len<0) {sign= -sign; y.len= -y.len;} 22: xfree(q); 23: xfree(r); 24: m_div(&x,&y,q,r); 25: if(sign==-1) 26: { q->len= -q->len; 27: r->len = - r->len; 28: } 29: return; 30: } 31: m_dsb(q,n,a,b) short *a,*b; 32: { long int x,qx; 33: int borrow,j; 34: short u; 35: qx=q; 36: borrow=0; 37: for(j=0;j<n;j++) 38: { x=borrow-a[j]*qx+b[j]; 39: b[j]=x&077777; 40: borrow=x>>15; 41: } 42: x=borrow+b[j]; 43: b[j]=x&077777; 44: if(x>>15 ==0) { return(0);} 45: borrow=0; 46: for(j=0;j<n;j++) 47: { u=a[j]+b[j]+borrow; 48: if(u<0) borrow=1; 49: else borrow=0; 50: b[j]=u&077777; 51: } 52: { return(1);} 53: } 54: m_trq(v1,v2,u1,u2,u3) 55: { long int d; 56: long int x1; 57: if(u1==v1) d=077777; 58: else d=(u1*0100000L+u2)/v1; 59: while(1) 60: { x1=u1*0100000L+u2-v1*d; 61: x1=x1*0100000L+u3-v2*d; 62: if(x1<0) d=d-1; 63: else {return(d);} 64: } 65: } 66: m_div(a,b,q,r) MINT *a,*b,*q,*r; 67: { MINT u,v,x,w; 68: short d,*qval; 69: int qq,n,v1,v2,j; 70: u.len=v.len=x.len=w.len=0; 71: if(b->len==0) { fatal("mdiv divide by zero"); return;} 72: if(b->len==1) 73: { r->val=xalloc(1,"m_div1"); 74: sdiv(a,b->val[0],q,r->val); 75: if(r->val[0]==0) 76: { shfree(r->val); 77: r->len=0; 78: } 79: else r->len=1; 80: return; 81: } 82: if(a->len < b->len) 83: { q->len=0; 84: r->len=a->len; 85: r->val=xalloc(r->len,"m_div2"); 86: for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq]; 87: return; 88: } 89: x.len=1; 90: x.val = &d; 91: n=b->len; 92: d=0100000L/(b->val[n-1]+1L); 93: mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */ 94: mult(b,&x,&v); 95: v1=v.val[n-1]; 96: v2=v.val[n-2]; 97: qval=xalloc(a->len-n+1,"m_div3"); 98: for(j=a->len-n;j>=0;j--) 99: { qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]); 100: if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1; 101: qval[j]=qq; 102: } 103: x.len=n; 104: x.val=u.val; 105: mcan(&x); 106: sdiv(&x,d,&w,(short *)&qq); 107: r->len=w.len; 108: r->val=w.val; 109: q->val=qval; 110: qq=a->len-n+1; 111: if(qq>0 && qval[qq-1]==0) qq -= 1; 112: q->len=qq; 113: if(qq==0) shfree(qval); 114: if(x.len!=0) xfree(&u); 115: xfree(&v); 116: return; 117: }