1*48370Sbostic /*- 2*48370Sbostic * %sccs.include.proprietary.c% 3*48370Sbostic */ 4*48370Sbostic 519819Sdist #ifndef lint 6*48370Sbostic static char sccsid[] = "@(#)mdiv.c 5.3 (Berkeley) 04/19/91"; 7*48370Sbostic #endif /* not lint */ 819819Sdist 99946Ssam #include <mp.h> 109946Ssam mdiv(a,b,q,r) MINT *a,*b,*q,*r; 119946Ssam { MINT x,y; 129946Ssam int sign; 139946Ssam sign=1; 149946Ssam x.val=a->val; 159946Ssam y.val=b->val; 169946Ssam x.len=a->len; 179946Ssam if(x.len<0) {sign= -1; x.len= -x.len;} 189946Ssam y.len=b->len; 199946Ssam if(y.len<0) {sign= -sign; y.len= -y.len;} 209946Ssam xfree(q); 219946Ssam xfree(r); 229946Ssam m_div(&x,&y,q,r); 239946Ssam if(sign==-1) 249946Ssam { q->len= -q->len; 259946Ssam r->len = - r->len; 269946Ssam } 279946Ssam return; 289946Ssam } 299946Ssam m_dsb(q,n,a,b) short *a,*b; 309946Ssam { long int x,qx; 3114332Ssam int borrow,j; 3214332Ssam short u; 339946Ssam qx=q; 349946Ssam borrow=0; 359946Ssam for(j=0;j<n;j++) 369946Ssam { x=borrow-a[j]*qx+b[j]; 379946Ssam b[j]=x&077777; 389946Ssam borrow=x>>15; 399946Ssam } 409946Ssam x=borrow+b[j]; 419946Ssam b[j]=x&077777; 429946Ssam if(x>>15 ==0) { return(0);} 439946Ssam borrow=0; 449946Ssam for(j=0;j<n;j++) 459946Ssam { u=a[j]+b[j]+borrow; 469946Ssam if(u<0) borrow=1; 479946Ssam else borrow=0; 489946Ssam b[j]=u&077777; 499946Ssam } 509946Ssam { return(1);} 519946Ssam } 529946Ssam m_trq(v1,v2,u1,u2,u3) 539946Ssam { long int d; 549946Ssam long int x1; 559946Ssam if(u1==v1) d=077777; 569946Ssam else d=(u1*0100000L+u2)/v1; 579946Ssam while(1) 589946Ssam { x1=u1*0100000L+u2-v1*d; 599946Ssam x1=x1*0100000L+u3-v2*d; 609946Ssam if(x1<0) d=d-1; 619946Ssam else {return(d);} 629946Ssam } 639946Ssam } 649946Ssam m_div(a,b,q,r) MINT *a,*b,*q,*r; 659946Ssam { MINT u,v,x,w; 669946Ssam short d,*qval; 679946Ssam int qq,n,v1,v2,j; 689946Ssam u.len=v.len=x.len=w.len=0; 699946Ssam if(b->len==0) { fatal("mdiv divide by zero"); return;} 709946Ssam if(b->len==1) 719946Ssam { r->val=xalloc(1,"m_div1"); 729946Ssam sdiv(a,b->val[0],q,r->val); 739946Ssam if(r->val[0]==0) 749946Ssam { shfree(r->val); 759946Ssam r->len=0; 769946Ssam } 779946Ssam else r->len=1; 789946Ssam return; 799946Ssam } 809946Ssam if(a->len < b->len) 819946Ssam { q->len=0; 829946Ssam r->len=a->len; 839946Ssam r->val=xalloc(r->len,"m_div2"); 849946Ssam for(qq=0;qq<r->len;qq++) r->val[qq]=a->val[qq]; 859946Ssam return; 869946Ssam } 879946Ssam x.len=1; 889946Ssam x.val = &d; 899946Ssam n=b->len; 909946Ssam d=0100000L/(b->val[n-1]+1L); 919946Ssam mult(a,&x,&u); /*subtle: relies on fact that mult allocates extra space */ 929946Ssam mult(b,&x,&v); 939946Ssam v1=v.val[n-1]; 949946Ssam v2=v.val[n-2]; 959946Ssam qval=xalloc(a->len-n+1,"m_div3"); 969946Ssam for(j=a->len-n;j>=0;j--) 979946Ssam { qq=m_trq(v1,v2,u.val[j+n],u.val[j+n-1],u.val[j+n-2]); 989946Ssam if(m_dsb(qq,n,v.val,&(u.val[j]))) qq -= 1; 999946Ssam qval[j]=qq; 1009946Ssam } 1019946Ssam x.len=n; 1029946Ssam x.val=u.val; 1039946Ssam mcan(&x); 1049946Ssam sdiv(&x,d,&w,(short *)&qq); 1059946Ssam r->len=w.len; 1069946Ssam r->val=w.val; 1079946Ssam q->val=qval; 1089946Ssam qq=a->len-n+1; 1099946Ssam if(qq>0 && qval[qq-1]==0) qq -= 1; 1109946Ssam q->len=qq; 1119946Ssam if(qq==0) shfree(qval); 1129946Ssam if(x.len!=0) xfree(&u); 1139946Ssam xfree(&v); 1149946Ssam return; 1159946Ssam } 116