userapps/opensource/sshd/libtommath/bn_mp_n_root.c

   1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   2  *
   3  * LibTomMath is library that provides for multiple-precision
   4  * integer arithmetic as well as number theoretic functionality.
   5  *
   6  * The library is designed directly after the MPI library by
   7  * Michael Fromberger but has been written from scratch with
   8  * additional optimizations in place.
   9  *
  10  * The library is free for all purposes without any express
  11  * guarantee it works.
  12  *
  13  * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
  14  */
  15 #include <tommath.h>
  16
  17 /* find the n'th root of an integer
  18  *
  19  * Result found such that (c)**b <= a and (c+1)**b > a
  20  *
  21  * This algorithm uses Newton's approximation
  22  * x[i+1] = x[i] - f(x[i])/f'(x[i])
  23  * which will find the root in log(N) time where
  24  * each step involves a fair bit.  This is not meant to
  25  * find huge roots [square and cube, etc].
  26  */
  27 int
  28 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 __T1;
  44   }
  45
  46   if ((res = mp_init (&t3)) != MP_OKAY) {
  47     goto __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 __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 __T3;
  68     }
  69
  70     /* numerator */
  71     /* t2 = t1**b */
  72     if ((res = mp_mul (&t3, &t1, &t2)) != MP_OKAY) {
  73       goto __T3;
  74     }
  75
  76     /* t2 = t1**b - a */
  77     if ((res = mp_sub (&t2, a, &t2)) != MP_OKAY) {
  78       goto __T3;
  79     }
  80
  81     /* denominator */
  82     /* t3 = t1**(b-1) * b  */
  83     if ((res = mp_mul_d (&t3, b, &t3)) != MP_OKAY) {
  84       goto __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 __T3;
  90     }
  91
  92     if ((res = mp_sub (&t1, &t3, &t2)) != MP_OKAY) {
  93       goto __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 __T3;
 101     }
 102
 103     if (mp_cmp (&t2, a) == MP_GT) {
 104       if ((res = mp_sub_d (&t1, 1, &t1)) != MP_OKAY) {
 105     goto __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 __T3:mp_clear (&t3);
 124 __T2:mp_clear (&t2);
 125 __T1:mp_clear (&t1);
 126   return res;
 127 }