usr/src/lib/libmp/common/msqrt.c

   1 /*      Copyright (c) 1984, 1986, 1987, 1988, 1989 AT&T */
   2 /*        All Rights Reserved   */
   3
   4
   5 /*
   6  * Copyright (c) 1980 Regents of the University of California.
   7  * All rights reserved.  The Berkeley software License Agreement
   8  * specifies the terms and conditions for redistribution.
   9  */
  10 /*      Portions Copyright(c) 1988, Sun Microsystems Inc.       */
  11 /*      All Rights Reserved                                     */
  12
  13 /*
  14  * Copyright (c) 1997, by Sun Microsystems, Inc.
  15  * All rights reserved.
  16  */
  17
  18 #ident  "%Z%%M% %I%     %E% SMI"        /* SVr4.0 1.1   */
  19
  20 /* LINTLIBRARY */
  21
  22 #include <mp.h>
  23 #include <sys/types.h>
  24 #include "libmp.h"
  25
  26 int
  27 mp_msqrt(MINT *a, MINT *b, MINT *r)
  28 {
  29         MINT a0, x, junk, y;
  30         int j;
  31
  32         a0.len = junk.len = y.len = 0;
  33         if (a->len < 0)
  34                 _mp_fatal("mp_msqrt: neg arg");
  35         if (a->len == 0) {
  36                 b->len = 0;
  37                 r->len = 0;
  38                 return (0);
  39         }
  40         if (a->len % 2 == 1)
  41                 x.len = (1 + a->len) / 2;
  42         else
  43                 x.len = 1 + a->len / 2;
  44         x.val = _mp_xalloc(x.len, "mp_msqrt");
  45         for (j = 0; j < x.len; x.val[j++] = 0)
  46                 ;
  47         if (a->len % 2 == 1)
  48                 x.val[x.len - 1] = 0400;
  49         else
  50                 x.val[x.len - 1] = 1;
  51         _mp_move(a, &a0);
  52         _mp_xfree(b);
  53         _mp_xfree(r);
  54 loop:
  55         mp_mdiv(&a0, &x, &y, &junk);
  56         _mp_xfree(&junk);
  57         mp_madd(&x, &y, &y);
  58         mp_sdiv(&y, 2, &y, (short *) &j);
  59         if (mp_mcmp(&x, &y) > 0) {
  60                 _mp_xfree(&x);
  61                 _mp_move(&y, &x);
  62                 _mp_xfree(&y);
  63                 goto loop;
  64         }
  65         _mp_xfree(&y);
  66         _mp_move(&x, b);
  67         mp_mult(&x, &x, &x);
  68         mp_msub(&a0, &x, r);
  69         _mp_xfree(&x);
  70         _mp_xfree(&a0);
  71         return (r->len);
  72 }