1#include <tommath.h>
2#ifdef BN_MP_SQRT_C
3/* LibTomMath, multiple-precision integer library -- Tom St Denis
4 *
5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality.
7 *
8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place.
11 *
12 * The library is free for all purposes without any express
13 * guarantee it works.
14 *
15 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16 */
17
18/* this function is less generic than mp_n_root, simpler and faster */
19int mp_sqrt(mp_int *arg, mp_int *ret)
20{
21  int res;
22  mp_int t1,t2;
23
24  /* must be positive */
25  if (arg->sign == MP_NEG) {
26    return MP_VAL;
27  }
28
29  /* easy out */
30  if (mp_iszero(arg) == MP_YES) {
31    mp_zero(ret);
32    return MP_OKAY;
33  }
34
35  if ((res = mp_init_copy(&t1, arg)) != MP_OKAY) {
36    return res;
37  }
38
39  if ((res = mp_init(&t2)) != MP_OKAY) {
40    goto E2;
41  }
42
43  /* First approx. (not very bad for large arg) */
44  mp_rshd (&t1,t1.used/2);
45
46  /* t1 > 0  */
47  if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
48    goto E1;
49  }
50  if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
51    goto E1;
52  }
53  if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
54    goto E1;
55  }
56  /* And now t1 > sqrt(arg) */
57  do {
58    if ((res = mp_div(arg,&t1,&t2,NULL)) != MP_OKAY) {
59      goto E1;
60    }
61    if ((res = mp_add(&t1,&t2,&t1)) != MP_OKAY) {
62      goto E1;
63    }
64    if ((res = mp_div_2(&t1,&t1)) != MP_OKAY) {
65      goto E1;
66    }
67    /* t1 >= sqrt(arg) >= t2 at this point */
68  } while (mp_cmp_mag(&t1,&t2) == MP_GT);
69
70  mp_exch(&t1,ret);
71
72E1: mp_clear(&t2);
73E2: mp_clear(&t1);
74  return res;
75}
76
77#endif
78
79/* $Source$ */
80/* $Revision$ */
81/* $Date$ */
82