1 files changed, 63 insertions, 0 deletions
diff --git a/libtommath/bn_s_mp_div_3.c b/libtommath/bn_s_mp_div_3.c
new file mode 100644
index 0000000..e0aeefc
--- /dev/null
+++ b/libtommath/bn_s_mp_div_3.c
@@ -0,0 +1,63 @@
+#include "tommath_private.h"
+#ifdef BN_S_MP_DIV_3_C
+/* LibTomMath, multiple-precision integer library -- Tom St Denis */
+/* SPDX-License-Identifier: Unlicense */
+
+/* divide by three (based on routine from MPI and the GMP manual) */
+mp_err s_mp_div_3(const mp_int *a, mp_int *c, mp_digit *d)
+{
+   mp_int   q;
+   mp_word  w, t;
+   mp_digit b;
+   mp_err   err;
+   int      ix;
+
+   /* b = 2**MP_DIGIT_BIT / 3 */
+   b = ((mp_word)1 << (mp_word)MP_DIGIT_BIT) / (mp_word)3;
+
+   if ((err = mp_init_size(&q, a->used)) != MP_OKAY) {
+      return err;
+   }
+
+   q.used = a->used;
+   q.sign = a->sign;
+   w = 0;
+   for (ix = a->used - 1; ix >= 0; ix--) {
+      w = (w << (mp_word)MP_DIGIT_BIT) | (mp_word)a->dp[ix];
+
+      if (w >= 3u) {
+         /* multiply w by [1/3] */
+         t = (w * (mp_word)b) >> (mp_word)MP_DIGIT_BIT;
+
+         /* now subtract 3 * [w/3] from w, to get the remainder */
+         w -= t+t+t;
+
+         /* fixup the remainder as required since
+          * the optimization is not exact.
+          */
+         while (w >= 3u) {
+            t += 1u;
+            w -= 3u;
+         }
+      } else {
+         t = 0;
+      }
+      q.dp[ix] = (mp_digit)t;
+   }
+
+   /* [optional] store the remainder */
+   if (d != NULL) {
+      *d = (mp_digit)w;
+   }
+
+   /* [optional] store the quotient */
+   if (c != NULL) {
+      mp_clamp(&q);
+      mp_exch(&q, c);
+   }
+   mp_clear(&q);
+
+   return err;
+}
+
+#endif