mroot.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Bignum root
   4 //
   5 //      Returns null pointer if not perfect root.
   6 //
   7 //      The sign of the radicand is ignored.
   8 //
   9 //-----------------------------------------------------------------------------
  10
  11 #include "stdafx.h"
  12 #include "defs.h"
  13
  14 unsigned int *
  15 mroot(unsigned int *n, unsigned int index)
  16 {
  17         int i, j, k;
  18         unsigned int m, *x, *y;
  19
  20         if (index == 0)
  21                 stop("root index is zero");
  22
  23         // count number of bits
  24
  25         k = 32 * (MLENGTH(n) - 1);
  26
  27         m = n[MLENGTH(n) - 1];
  28
  29         while (m) {
  30                 m >>= 1;
  31                 k++;
  32         }
  33
  34         if (k == 0)
  35                 return mint(0);
  36
  37         // initial guess
  38
  39         k = (k - 1) / index;
  40
  41         j = k / 32 + 1;
  42         x = mnew(j);
  43         MSIGN(x) = 1;
  44         MLENGTH(x) = j;
  45         for (i = 0; i < j; i++)
  46                 x[i] = 0;
  47
  48         while (k >= 0) {
  49                 mp_set_bit(x, k);
  50                 y = mpow(x, index);
  51                 switch (mcmp(y, n)) {
  52                 case -1:
  53                         break;
  54                 case 0:
  55                         mfree(y);
  56                         return x;
  57                 case 1:
  58                         mp_clr_bit(x, k);
  59                         break;
  60                 }
  61                 mfree(y);
  62                 k--;
  63         }
  64
  65         mfree(x);
  66
  67         return 0;
  68 }