2.9BSD - /usr/src/lib/ape/msqrt.c

   1: #ifndef lint
   2: static char sccs_id[] = "@(#)msqrt.c	2.1	7/6/82";
   3: #endif  lint
   4: 
   5: #include <ape.h>
   6: msqrt(a,b,r) /* b = square root of a - r, where r is the smallest
   7: 		possible positive integral remainder.
   8: 		The algorithm used is Newton's method.
   9: 		The value returned is the length of the
  10: 		remainder. */
  11: MINT *a,*b,*r;
  12: {   MINT x,junk,y;
  13:     int j;
  14: 
  15:     x.len=junk.len=y.len=0;
  16:     if (a->len<0) aperror("msqrt: neg arg");
  17:     if (a->len==0)
  18:         {
  19:         b->len=0;
  20:         r->len=0;
  21:         return(0);
  22:         }
  23: /* Generate the initial approximation: */
  24:     if (a->len%2==1) x.len=(1+a->len)/2;
  25:     else x.len=1+a->len/2;
  26:     x.val=xalloc(x.len,"msqrt");
  27:     for (j = 0; j < x.len; x.val[j++] = 0)
  28:         ;
  29:     if (a->len%2==1) x.val[x.len-1]=0400;
  30:     else x.val[x.len-1]=1;
  31: loop:
  32:     mdiv(a,&x,&y,&junk);
  33:     xfree(&junk);
  34:     madd(&x,&y,&y);
  35:     sdiv(&y,2,&y,(short *)&j);
  36:     if (mcmp(&x,&y)>0)
  37:         {
  38:         move(&y,&x);
  39:         xfree(&y);
  40:         goto loop;
  41:         }
  42:     xfree(&y);
  43:     move(&x,b);
  44:     mult(&x,&x,&x);
  45:     msub(a,&x,r);
  46:     xfree(&x);
  47:     return(r->len);
  48: }

Defined functions

Defined variables