Home | History | Annotate | Line # | Download | only in libtommath 
      1 /*	$NetBSD: bn_mp_n_root.c,v 1.2 2017/01/28 21:31:47 christos Exp $	*/
      2 
      3 #include <tommath.h>
      4 #ifdef BN_MP_N_ROOT_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 (at) gmail.com, http://libtom.org
     18  */
     19 
     20 /* find the n'th root of an integer
     21  *
     22  * Result found such that (c)**b <= a and (c+1)**b > a
     23  *
     24  * This algorithm uses Newton's approximation
     25  * x[i+1] = x[i] - f(x[i])/f'(x[i])
     26  * which will find the root in log(N) time where
     27  * each step involves a fair bit.  This is not meant to
     28  * find huge roots [square and cube, etc].
     29  */
     30 int mp_n_root (mp_int * a, mp_digit b, mp_int * c)
     31 {
     32   mp_int  t1, t2, t3;
     33   int     res, neg;
     34 
     35   /* input must be positive if b is even */
     36   if ((b & 1) == 0 && a->sign == MP_NEG) {
     37     return MP_VAL;
     38   }
     39 
     40   if ((res = mp_init (&t1)) != MP_OKAY) {
     41     return res;
     42   }
     43 
     44   if ((res = mp_init (&t2)) != MP_OKAY) {
     45     goto LBL_T1;
     46   }
     47 
     48   if ((res = mp_init (&t3)) != MP_OKAY) {
     49     goto LBL_T2;
     50   }
     51 
     52   /* if a is negative fudge the sign but keep track */
     53   neg     = a->sign;
     54   a->sign = MP_ZPOS;
     55 
     56   /* t2 = 2 */
     57   mp_set (&t2, 2);
     58 
     59   do {
     60     /* t1 = t2 */
     61     if ((res = mp_copy (&t2, &t1)) != MP_OKAY) {
     62       goto LBL_T3;
     63     }
     64 
     65     /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */
     66 
     67     /* t3 = t1**(b-1) */
     68     if ((res = mp_expt_d (&t1, b - 1, &t3)) != MP_OKAY) {
     69       goto LBL_T3;
     70     }
     71 
     72     /* numerator */
     73     /* t2 = t1**b */
     74     if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
     75       goto LBL_T3;
     76     }
     77 
     78     /* t2 = t1**b - a */
     79     if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {
     80       goto LBL_T3;
     81     }
     82 
     83     /* denominator */
     84     /* t3 = t1**(b-1) * b  */
     85     if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
     86       goto LBL_T3;
     87     }
     88 
     89     /* t3 = (t1**b - a)/(b * t1**(b-1)) */
     90     if ((res = mp_div (&t2, &t3, &t3, NULL)) != MP_OKAY) {
     91       goto LBL_T3;
     92     }
     93 
     94     if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
     95       goto LBL_T3;
     96     }
     97   }  while (mp_cmp (&t1, &t2) != MP_EQ);
     98 
     99   /* result can be off by a few so check */
    100   for (;;) {
    101     if ((res = mp_expt_d (&t1, b, &t2)) != MP_OKAY) {
    102       goto LBL_T3;
    103     }
    104 
    105     if (mp_cmp (&t2, a) == MP_GT) {
    106       if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
    107          goto LBL_T3;
    108       }
    109     } else {
    110       break;
    111     }
    112   }
    113 
    114   /* reset the sign of a first */
    115   a->sign = neg;
    116 
    117   /* set the result */
    118   mp_exch (&t1, c);
    119 
    120   /* set the sign of the result */
    121   c->sign = neg;
    122 
    123   res = MP_OKAY;
    124 
    125 LBL_T3:mp_clear (&t3);
    126 LBL_T2:mp_clear (&t2);
    127 LBL_T1:mp_clear (&t1);
    128   return res;
    129 }
    130 #endif
    131 
    132 /* Source: /cvs/libtom/libtommath/bn_mp_n_root.c,v  */
    133 /* Revision: 1.4  */
    134 /* Date: 2006/12/28 01:25:13  */
    135