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