Addition of fast division (recursive divrem only)

demo/test.c

@@ -2327,6 +2327,139 @@ static int test_mp_radix_size(void)
 }
 
 
+/* Some larger values to test the fast division algorithm */
+static int test_s_mp_div_recursive(void)
+{
+   mp_int a, b, c_q, c_r, d_q, d_r;
+   int size, err;
+
+   if ((err = mp_init_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+
+   for (size = MP_KARATSUBA_MUL_CUTOFF; size < 3 * MP_KARATSUBA_MUL_CUTOFF; size += 10) {
+      fprintf(stderr,"sizes = %d / %d\n", 10 * size, size);
+      /* Relation 10:1 */
+      if ((err = mp_rand(&a, 10 * size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_rand(&b, size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_recursive(&a, &b, &c_q, &c_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_school(&a, &b, &d_q, &d_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_q, &d_q) != MP_EQ) {
+         fprintf(stderr, "1. Recursive division failed at sizes %d / %d, wrong quotient\n",
+                 10 * size, size);
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_r, &d_r) != MP_EQ) {
+         fprintf(stderr, "1. Recursive division failed at sizes %d / %d, wrong remainder\n",
+                 10 * size, size);
+         goto LBL_ERR;
+      }
+      fprintf(stderr,"sizes = %d / %d\n", 2 * size, size);
+      /* Relation 2:1 */
+      if ((err = mp_rand(&a, 2 * size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_rand(&b, size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_recursive(&a, &b, &c_q, &c_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_school(&a, &b, &d_q, &d_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_q, &d_q) != MP_EQ) {
+         fprintf(stderr, "2. Recursive division failed at sizes %d / %d, wrong quotient\n",
+                 2 * size, size);
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_r, &d_r) != MP_EQ) {
+         fprintf(stderr, "2. Recursive division failed at sizes %d / %d, wrong remainder\n",
+                 2 * size, size);
+         goto LBL_ERR;
+      }
+      fprintf(stderr,"sizes = %d / %d\n", 3 * size, 2 * size);
+      /* Upper limit 3:2 */
+      if ((err = mp_rand(&a, 3 * size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_rand(&b, 2 * size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_recursive(&a, &b, &c_q, &c_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_school(&a, &b, &d_q, &d_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_q, &d_q) != MP_EQ) {
+         fprintf(stderr, "3. Recursive division failed at sizes %d / %d, wrong quotient\n",
+                 3 * size, 2 * size);
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_r, &d_r) != MP_EQ) {
+         fprintf(stderr, "3. Recursive division failed at sizes %d / %d, wrong remainder\n",
+                 3 * size, 2 * size);
+         goto LBL_ERR;
+      }
+   }
+
+   mp_clear_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL);
+   return EXIT_SUCCESS;
+LBL_ERR:
+   mp_clear_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL);
+   return EXIT_FAILURE;
+}
+
+static int test_s_mp_div_small(void)
+{
+   mp_int a, b, c_q, c_r, d_q, d_r;
+   int size, err;
+
+   if ((err = mp_init_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL)) != MP_OKAY) {
+      goto LBL_ERR;
+   }
+   for (size = 1; size < MP_KARATSUBA_MUL_CUTOFF; size += 10) {
+      fprintf(stderr,"sizes = %d / %d\n", 2 * size, size);
+      /* Relation 10:1 */
+      if ((err = mp_rand(&a, 2 * size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = mp_rand(&b, size)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_small(&a, &b, &c_q, &c_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if ((err = s_mp_div_school(&a, &b, &d_q, &d_r)) != MP_OKAY) {
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_q, &d_q) != MP_EQ) {
+         fprintf(stderr, "1. Small division failed at sizes %d / %d, wrong quotient\n",
+                 2 * size, size);
+         goto LBL_ERR;
+      }
+      if (mp_cmp(&c_r, &d_r) != MP_EQ) {
+         fprintf(stderr, "1. Small division failed at sizes %d / %d, wrong remainder\n",
+                 2 * size, size);
+         goto LBL_ERR;
+      }
+   }
+   mp_clear_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL);
+   return EXIT_SUCCESS;
+LBL_ERR:
+   mp_clear_multi(&a, &b, &c_q, &c_r, &d_q, &d_r, NULL);
+   return EXIT_FAILURE;
+}
+
 
 static int test_mp_read_write_ubin(void)
 {
@@ -2500,6 +2633,8 @@ static int unit_tests(int argc, char **argv)
       T1(mp_sqrt, MP_SQRT),
       T1(mp_sqrtmod_prime, MP_SQRTMOD_PRIME),
       T1(mp_xor, MP_XOR),
+      T2(s_mp_div_recursive, S_MP_DIV_RECURSIVE, S_MP_DIV_SCHOOL),
+      T2(s_mp_div_small, S_MP_DIV_SMALL, S_MP_DIV_SCHOOL),
       T1(s_mp_balance_mul, S_MP_BALANCE_MUL),
       T1(s_mp_karatsuba_mul, S_MP_KARATSUBA_MUL),
       T1(s_mp_karatsuba_sqr, S_MP_KARATSUBA_SQR),

libtommath_VS2008.vcproj

@@ -852,6 +852,18 @@
 			RelativePath="s_mp_balance_mul.c"
 			>
 		</File>
+		<File
+			RelativePath="s_mp_div_recursive.c"
+			>
+		</File>
+		<File
+			RelativePath="s_mp_div_school.c"
+			>
+		</File>
+		<File
+			RelativePath="s_mp_div_small.c"
+			>
+		</File>
 		<File
 			RelativePath="s_mp_exptmod.c"
 			>

makefile

@@ -44,12 +44,12 @@ mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.
 mp_root_u32.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o mp_set_l.o \
 mp_set_ll.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_set_ull.o mp_shrink.o mp_signed_rsh.o mp_sqr.o \
 mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o \
-mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o s_mp_exptmod.o \
-s_mp_exptmod_fast.o s_mp_get_bit.o s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o \
-s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o \
-s_mp_mul_digs_fast.o s_mp_mul_high_digs.o s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o \
-s_mp_rand_jenkins.o s_mp_rand_platform.o s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o \
-s_mp_toom_mul.o s_mp_toom_sqr.o
+mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_get_bit.o \
+s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o \
+s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o s_mp_mul_digs_fast.o s_mp_mul_high_digs.o \
+s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o s_mp_rand_jenkins.o s_mp_rand_platform.o \
+s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o s_mp_toom_mul.o s_mp_toom_sqr.o
 
 #END_INS
 

makefile.mingw

@@ -47,12 +47,12 @@ mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.
 mp_root_u32.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o mp_set_l.o \
 mp_set_ll.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_set_ull.o mp_shrink.o mp_signed_rsh.o mp_sqr.o \
 mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o \
-mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o s_mp_exptmod.o \
-s_mp_exptmod_fast.o s_mp_get_bit.o s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o \
-s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o \
-s_mp_mul_digs_fast.o s_mp_mul_high_digs.o s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o \
-s_mp_rand_jenkins.o s_mp_rand_platform.o s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o \
-s_mp_toom_mul.o s_mp_toom_sqr.o
+mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_get_bit.o \
+s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o \
+s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o s_mp_mul_digs_fast.o s_mp_mul_high_digs.o \
+s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o s_mp_rand_jenkins.o s_mp_rand_platform.o \
+s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o s_mp_toom_mul.o s_mp_toom_sqr.o
 
 HEADERS_PUB=tommath.h
 HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)

makefile.msvc

@@ -39,12 +39,12 @@ mp_reduce_2k_setup.obj mp_reduce_2k_setup_l.obj mp_reduce_is_2k.obj mp_reduce_is
 mp_root_u32.obj mp_rshd.obj mp_sbin_size.obj mp_set.obj mp_set_double.obj mp_set_i32.obj mp_set_i64.obj mp_set_l.obj \
 mp_set_ll.obj mp_set_u32.obj mp_set_u64.obj mp_set_ul.obj mp_set_ull.obj mp_shrink.obj mp_signed_rsh.obj mp_sqr.obj \
 mp_sqrmod.obj mp_sqrt.obj mp_sqrtmod_prime.obj mp_sub.obj mp_sub_d.obj mp_submod.obj mp_to_radix.obj mp_to_sbin.obj \
-mp_to_ubin.obj mp_ubin_size.obj mp_unpack.obj mp_xor.obj mp_zero.obj s_mp_add.obj s_mp_balance_mul.obj s_mp_exptmod.obj \
-s_mp_exptmod_fast.obj s_mp_get_bit.obj s_mp_invmod_fast.obj s_mp_invmod_slow.obj s_mp_karatsuba_mul.obj \
-s_mp_karatsuba_sqr.obj s_mp_log.obj s_mp_log_d.obj s_mp_montgomery_reduce_fast.obj s_mp_mul_digs.obj \
-s_mp_mul_digs_fast.obj s_mp_mul_high_digs.obj s_mp_mul_high_digs_fast.obj s_mp_prime_is_divisible.obj \
-s_mp_rand_jenkins.obj s_mp_rand_platform.obj s_mp_reverse.obj s_mp_sqr.obj s_mp_sqr_fast.obj s_mp_sub.obj \
-s_mp_toom_mul.obj s_mp_toom_sqr.obj
+mp_to_ubin.obj mp_ubin_size.obj mp_unpack.obj mp_xor.obj mp_zero.obj s_mp_add.obj s_mp_balance_mul.obj \
+s_mp_div_recursive.obj s_mp_div_school.obj s_mp_div_small.obj s_mp_exptmod.obj s_mp_exptmod_fast.obj s_mp_get_bit.obj \
+s_mp_invmod_fast.obj s_mp_invmod_slow.obj s_mp_karatsuba_mul.obj s_mp_karatsuba_sqr.obj s_mp_log.obj s_mp_log_d.obj \
+s_mp_montgomery_reduce_fast.obj s_mp_mul_digs.obj s_mp_mul_digs_fast.obj s_mp_mul_high_digs.obj \
+s_mp_mul_high_digs_fast.obj s_mp_prime_is_divisible.obj s_mp_rand_jenkins.obj s_mp_rand_platform.obj \
+s_mp_reverse.obj s_mp_sqr.obj s_mp_sqr_fast.obj s_mp_sub.obj s_mp_toom_mul.obj s_mp_toom_sqr.obj
 
 HEADERS_PUB=tommath.h
 HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)

makefile.shared

@@ -41,12 +41,12 @@ mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.
 mp_root_u32.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o mp_set_l.o \
 mp_set_ll.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_set_ull.o mp_shrink.o mp_signed_rsh.o mp_sqr.o \
 mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o \
-mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o s_mp_exptmod.o \
-s_mp_exptmod_fast.o s_mp_get_bit.o s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o \
-s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o \
-s_mp_mul_digs_fast.o s_mp_mul_high_digs.o s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o \
-s_mp_rand_jenkins.o s_mp_rand_platform.o s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o \
-s_mp_toom_mul.o s_mp_toom_sqr.o
+mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_get_bit.o \
+s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o \
+s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o s_mp_mul_digs_fast.o s_mp_mul_high_digs.o \
+s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o s_mp_rand_jenkins.o s_mp_rand_platform.o \
+s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o s_mp_toom_mul.o s_mp_toom_sqr.o
 
 #END_INS
 

makefile.unix

@@ -48,12 +48,12 @@ mp_reduce_2k_setup.o mp_reduce_2k_setup_l.o mp_reduce_is_2k.o mp_reduce_is_2k_l.
 mp_root_u32.o mp_rshd.o mp_sbin_size.o mp_set.o mp_set_double.o mp_set_i32.o mp_set_i64.o mp_set_l.o \
 mp_set_ll.o mp_set_u32.o mp_set_u64.o mp_set_ul.o mp_set_ull.o mp_shrink.o mp_signed_rsh.o mp_sqr.o \
 mp_sqrmod.o mp_sqrt.o mp_sqrtmod_prime.o mp_sub.o mp_sub_d.o mp_submod.o mp_to_radix.o mp_to_sbin.o \
-mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o s_mp_exptmod.o \
-s_mp_exptmod_fast.o s_mp_get_bit.o s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o \
-s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o \
-s_mp_mul_digs_fast.o s_mp_mul_high_digs.o s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o \
-s_mp_rand_jenkins.o s_mp_rand_platform.o s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o \
-s_mp_toom_mul.o s_mp_toom_sqr.o
+mp_to_ubin.o mp_ubin_size.o mp_unpack.o mp_xor.o mp_zero.o s_mp_add.o s_mp_balance_mul.o \
+s_mp_div_recursive.o s_mp_div_school.o s_mp_div_small.o s_mp_exptmod.o s_mp_exptmod_fast.o s_mp_get_bit.o \
+s_mp_invmod_fast.o s_mp_invmod_slow.o s_mp_karatsuba_mul.o s_mp_karatsuba_sqr.o s_mp_log.o s_mp_log_d.o \
+s_mp_montgomery_reduce_fast.o s_mp_mul_digs.o s_mp_mul_digs_fast.o s_mp_mul_high_digs.o \
+s_mp_mul_high_digs_fast.o s_mp_prime_is_divisible.o s_mp_rand_jenkins.o s_mp_rand_platform.o \
+s_mp_reverse.o s_mp_sqr.o s_mp_sqr_fast.o s_mp_sub.o s_mp_toom_mul.o s_mp_toom_sqr.o
 
 HEADERS_PUB=tommath.h
 HEADERS=tommath_private.h tommath_class.h tommath_superclass.h tommath_cutoffs.h $(HEADERS_PUB)