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