--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/tomsfastmath/src/divide/fp_div_d.c	Wed Nov 23 18:10:20 2011 +0700
@@ -0,0 +1,98 @@
+/* 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>
+
+static int s_is_power_of_two(fp_digit b, int *p)
+{
+   int x;
+
+   /* fast return if no power of two */
+   if ((b==0) || (b & (b-1))) {
+      return 0;
+   }
+
+   for (x = 0; x < DIGIT_BIT; x++) {
+      if (b == (((fp_digit)1)<<x)) {
+         *p = x;
+         return 1;
+      }
+   }
+   return 0;
+}
+
+/* a/b => cb + d == a */
+int fp_div_d(fp_int *a, fp_digit b, fp_int *c, fp_digit *d)
+{
+  fp_int   q;
+  fp_word  w;
+  fp_digit t;
+  int      ix;
+
+  /* cannot divide by zero */
+  if (b == 0) {
+     return FP_VAL;
+  }
+
+  /* quick outs */
+  if (b == 1 || fp_iszero(a) == 1) {
+     if (d != NULL) {
+        *d = 0;
+     }
+     if (c != NULL) {
+        fp_copy(a, c);
+     }
+     return FP_OKAY;
+  }
+
+  /* power of two ? */
+  if (s_is_power_of_two(b, &ix) == 1) {
+     if (d != NULL) {
+        *d = a->dp[0] & ((((fp_digit)1)<<ix) - 1);
+     }
+     if (c != NULL) {
+        fp_div_2d(a, ix, c, NULL);
+     }
+     return FP_OKAY;
+  }
+
+  /* no easy answer [c'est la vie].  Just division */
+  fp_init(&q);
+  
+  q.used = a->used;
+  q.sign = a->sign;
+  w = 0;
+  for (ix = a->used - 1; ix >= 0; ix--) {
+     w = (w << ((fp_word)DIGIT_BIT)) | ((fp_word)a->dp[ix]);
+     
+     if (w >= b) {
+        t = (fp_digit)(w / b);
+        w -= ((fp_word)t) * ((fp_word)b);
+      } else {
+        t = 0;
+      }
+      q.dp[ix] = (fp_digit)t;
+  }
+  
+  if (d != NULL) {
+     *d = (fp_digit)w;
+  }
+  
+  if (c != NULL) {
+     fp_clamp(&q);
+     fp_copy(&q, c);
+  }
+ 
+  return FP_OKAY;
+}
+
+
+/* $Source$ */
+/* $Revision$ */
+/* $Date$ */