--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomsfastmath/src/divide/fp_div.c	Wed Nov 23 18:10:20 2011 +0700
@@ -0,0 +1,157 @@
+/* TomsFastMath, a fast ISO C bignum library.
+ * 
+ * This project is meant to fill in where LibTomMath
+ * falls short.  That is speed ;-)
+ *
+ * This project is public domain and free for all purposes.
+ * 
+ * Tom St Denis, [email protected]
+ */
+#include <tfm.h>
+
+/* a/b => cb + d == a */
+int fp_div(fp_int *a, fp_int *b, fp_int *c, fp_int *d)
+{
+  fp_int  q, x, y, t1, t2;
+  int     n, t, i, norm, neg;
+
+  /* is divisor zero ? */
+  if (fp_iszero (b) == 1) {
+    return FP_VAL;
+  }
+
+  /* if a < b then q=0, r = a */
+  if (fp_cmp_mag (a, b) == FP_LT) {
+    if (d != NULL) {
+      fp_copy (a, d);
+    } 
+    if (c != NULL) {
+      fp_zero (c);
+    }
+    return FP_OKAY;
+  }
+
+  fp_init(&q);
+  q.used = a->used + 2;
+
+  fp_init(&t1);
+  fp_init(&t2);
+  fp_init_copy(&x, a);
+  fp_init_copy(&y, b);
+
+  /* fix the sign */
+  neg = (a->sign == b->sign) ? FP_ZPOS : FP_NEG;
+  x.sign = y.sign = FP_ZPOS;
+
+  /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
+  norm = fp_count_bits(&y) % DIGIT_BIT;
+  if (norm < (int)(DIGIT_BIT-1)) {
+     norm = (DIGIT_BIT-1) - norm;
+     fp_mul_2d (&x, norm, &x);
+     fp_mul_2d (&y, norm, &y);
+  } else {
+     norm = 0;
+  }
+
+  /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
+  n = x.used - 1;
+  t = y.used - 1;
+
+  /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
+  fp_lshd (&y, n - t);                                             /* y = y*b**{n-t} */
+
+  while (fp_cmp (&x, &y) != FP_LT) {
+    ++(q.dp[n - t]);
+    fp_sub (&x, &y, &x);
+  }
+
+  /* reset y by shifting it back down */
+  fp_rshd (&y, n - t);
+
+  /* step 3. for i from n down to (t + 1) */
+  for (i = n; i >= (t + 1); i--) {
+    if (i > x.used) {
+      continue;
+    }
+
+    /* step 3.1 if xi == yt then set q{i-t-1} to b-1, 
+     * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
+    if (x.dp[i] == y.dp[t]) {
+      q.dp[i - t - 1] = ((((fp_word)1) << DIGIT_BIT) - 1);
+    } else {
+      fp_word tmp;
+      tmp = ((fp_word) x.dp[i]) << ((fp_word) DIGIT_BIT);
+      tmp |= ((fp_word) x.dp[i - 1]);
+      tmp /= ((fp_word) y.dp[t]);
+      q.dp[i - t - 1] = (fp_digit) (tmp);
+    }
+
+    /* while (q{i-t-1} * (yt * b + y{t-1})) > 
+             xi * b**2 + xi-1 * b + xi-2 
+     
+       do q{i-t-1} -= 1; 
+    */
+    q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
+    do {
+      q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
+
+      /* find left hand */
+      fp_zero (&t1);
+      t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
+      t1.dp[1] = y.dp[t];
+      t1.used = 2;
+      fp_mul_d (&t1, q.dp[i - t - 1], &t1);
+
+      /* find right hand */
+      t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
+      t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
+      t2.dp[2] = x.dp[i];
+      t2.used = 3;
+    } while (fp_cmp_mag(&t1, &t2) == FP_GT);
+
+    /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
+    fp_mul_d (&y, q.dp[i - t - 1], &t1);
+    fp_lshd  (&t1, i - t - 1);
+    fp_sub   (&x, &t1, &x);
+
+    /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
+    if (x.sign == FP_NEG) {
+      fp_copy (&y, &t1);
+      fp_lshd (&t1, i - t - 1);
+      fp_add (&x, &t1, &x);
+      q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
+    }
+  }
+
+  /* now q is the quotient and x is the remainder 
+   * [which we have to normalize] 
+   */
+  
+  /* get sign before writing to c */
+  x.sign = x.used == 0 ? FP_ZPOS : a->sign;
+
+  if (c != NULL) {
+    fp_clamp (&q);
+    fp_copy (&q, c);
+    c->sign = neg;
+  }
+
+  if (d != NULL) {
+    fp_div_2d (&x, norm, &x, NULL);
+
+/* the following is a kludge, essentially we were seeing the right remainder but 
+   with excess digits that should have been zero
+ */
+    for (i = b->used; i < x.used; i++) {
+        x.dp[i] = 0;
+    }
+    fp_clamp(&x);
+    fp_copy (&x, d);
+  }
+
+  return FP_OKAY;
+}
+
+/* $Source$ */
+/* $Revision$ */
+/* $Date$ */