title: Newton-Raphson Division With Big Integers
body: <p>I'm making a BigInt class as a programming exercise. It uses a vector of 2's complement signed ints in base-65536 (so that 32-bit multiplications don't overflow. I will increase the base once I get it fully working).</p>
  
  <p>All of the basic math operations are coded, with one problem: division is <em>painfully</em> slow with the basic algorithm I was able to create. (It kind of works like binary division for each digit of the quotient... I'm not going to post it unless someone wants to see it....)</p>
  
  <p>Instead of my slow algorithm, I want to use Newton-Raphson to find the (shifted) reciprocal and then multiply (and shift). I think I have my head around the basics: you give the formula <strong>(x1 = x0(2 - x0 * divisor))</strong> a good initial guess, and then after some amount of iterations, x converges to the reciprocal. This part seems easy enough... but I am running into some problems when trying to apply this formula to big integers:</p>
  
  <p><strong>Problem 1:</strong></p>
  
  <p>Because I am working with integers... well... I can't use fractions. This seems to cause x to always diverge (x0 * divisor must be &lt;2 it seems?). My intuition tells me there should be some modification to the equation that would allow it to work integers (to some accuracy) but I am really struggling to find out what it is. (My lack of math skills are beating me up here....) I think I need to find some equivalent equation where instead of <strong>d</strong> there is <strong>d*[base^somePower]</strong>? Can there be some equation like <strong>(x1 = x0(2 - x0 * d))</strong> that works with whole numbers?</p>
  
  <p><strong>Problem 2:</strong></p>
  
  <p>When I use Newton's formula to find the reciprocal of some numbers, the result ends up being just a small faction below what the answer should be... ex. when trying to find reciprocal of 4 (in decimal):</p>
  
  <pre><code>x0 = 0.3
  x1 = 0.24
  x2 = 0.2496
  x3 = 0.24999936
  x4 = 0.2499999999983616
  x5 = 0.24999999999999999999998926258176
  </code></pre>
  
  <p>If I were representing numbers in base-10, I would want a result of 25 (and to remember to right shift product by 2). With some reciprocals such as 1/3, you can simply truncate the result after you know you have enough accuracy. But how can I pull out the correct reciprocal from the above result?</p>
  
  <p>Sorry if this is all too vague or if I'm asking for too much. I looked through Wikipedia and all of the research papers I could find on Google, but I feel like I'm banging my head against a wall. I appreciate any help anyone can give me!</p>
  
  <p>...</p>
  
  <p>Edit: Got the algorithm working, although it is much slower than I expected. I actually lost a lot of speed compared to my old algorithm, even on numbers with thousands of digits... I'm still missing something. It's not a problem with multiplication, which is very fast. (I am indeed using Karatsuba's algorithm).</p>
  
  <p>For anyone interested, here is my current iteration of the Newton-Raphson algorithm:</p>
  
  <pre><code>bigint operator/(const bigint&amp; lhs, const bigint&amp; rhs) {
      if (rhs == 0) throw overflow_error("Divide by zero exception");
      bigint dividend = lhs;
      bigint divisor = rhs;
  
      bool negative = 0;
      if (dividend &lt; 0) {
          negative = !negative;
          dividend.invert();
      }
      if (divisor &lt; 0) {
          negative = !negative;
          divisor.invert();
      }
  
      int k = dividend.numBits() + divisor.numBits();
      bigint pow2 = 1;
      pow2 &lt;&lt;= k + 1;
  
      bigint x = dividend - divisor;
      bigint lastx = 0;
      bigint lastlastx = 0;
      while (1) {
          x = (x * (pow2 - x * divisor)) &gt;&gt; k;
          if (x == lastx || x == lastlastx) break;
          lastlastx = lastx;
          lastx = x;
      }
      bigint quotient = dividend * x &gt;&gt; k;
      if (dividend - (quotient * divisor) &gt;= divisor) quotient++;
      if (negative)quotient.invert();
      return quotient;
  }
  </code></pre>
  
  <p>And here is my (really ugly) old algorithm that is faster:</p>
  
  <pre><code>bigint operator/(const bigint&amp; lhs, const bigint &amp; rhs) {
      if (rhs == 0) throw overflow_error("Divide by zero exception");
      bigint dividend = lhs;
      bigint divisor = rhs;
  
      bool negative = 0;
      if (dividend &lt; 0) {
          negative = !negative;
          dividend.invert();
      }
      if (divisor &lt; 0) {
          negative = !negative;
          divisor.invert();
      }
  
      bigint remainder = 0;
      bigint quotient = 0;
      while (dividend.value.size() &gt; 0) {
          remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1));
          remainder.value.push_back(0);
          remainder.unPad();
          dividend.value.pop_back();
  
          if (divisor &gt; remainder) {
              quotient.value.push_back(0);
          } else {
              int count = 0;
              int i = MSB;
              bigint value = 0;
              while (i &gt; 0) {
                  bigint increase = divisor * i;
                  bigint next = value + increase;
                  if (next &lt;= remainder) {
                      value = next;
                      count += i;
                  }
                  i &gt;&gt;= 1;
              }
              quotient.value.push_back(count);
              remainder -= value;
          }
      }
  
      for (int i = 0; i &lt; quotient.value.size() / 2; i++) {
          int swap = quotient.value.at(i);
          quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i);
          quotient.value.at(quotient.value.size() - 1 - i) = swap;
      }
  
      if (negative)quotient.invert();
      quotient.unPad();
      return quotient;
  }
  </code></pre>