libs/libtommath/bn_mp_n_root.c

   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@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 $ */