Commit 244c698ecd69be18db21d09b58093ec580351855
2019-11-24T05:37:43
corrected startvalue for sigma and cutoff in mp_div
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
diff --git a/mp_div.c b/mp_div.c
index cbc52e8..b092d7b 100644
--- a/mp_div.c
+++ b/mp_div.c
@@ -26,7 +26,7 @@ mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *c, mp_int *d)
}
if (MP_HAS(S_MP_DIV_RECURSIVE)
- && (b->used > MP_MUL_KARATSUBA_CUTOFF)
+ && (b->used > (2 * MP_MUL_KARATSUBA_CUTOFF))
&& (b->used <= ((a->used/3)*2))) {
err = s_mp_div_recursive(a, b, c, d);
} else if (MP_HAS(S_MP_DIV_SCHOOL)) {
diff --git a/s_mp_div_recursive.c b/s_mp_div_recursive.c
index df9a297..d719c4e 100644
--- a/s_mp_div_recursive.c
+++ b/s_mp_div_recursive.c
@@ -20,7 +20,7 @@ static mp_err s_recursion(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
mp_int A1, A2, B1, B0, Q1, Q0, R1, R0, t;
int m = a->used - b->used, k = m/2;
- if (m < MP_MUL_KARATSUBA_CUTOFF) {
+ if (m < (MP_MUL_KARATSUBA_CUTOFF)) {
return s_mp_div_school(a, b, q, r);
}
@@ -43,15 +43,16 @@ static mp_err s_recursion(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
if ((err = mp_sub(&A1, &t, &A1)) != MP_OKAY) goto LBL_ERR;
/* while A1 < 0 do Q1 = Q1 - 1, A1 = A1 + (beta^k * B) */
- if ((err = mp_mul_2d(b, k * MP_DIGIT_BIT, &t)) != MP_OKAY) goto LBL_ERR;
- while (mp_cmp_d(&A1, 0uL) == MP_LT) {
- if ((err = mp_decr(&Q1)) != MP_OKAY) goto LBL_ERR;
- if ((err = mp_add(&A1, &t, &A1)) != MP_OKAY) goto LBL_ERR;
+ if (mp_cmp_d(&A1, 0uL) == MP_LT) {
+ if ((err = mp_mul_2d(b, k * MP_DIGIT_BIT, &t)) != MP_OKAY) goto LBL_ERR;
+ do {
+ if ((err = mp_decr(&Q1)) != MP_OKAY) goto LBL_ERR;
+ if ((err = mp_add(&A1, &t, &A1)) != MP_OKAY) goto LBL_ERR;
+ } while (mp_cmp_d(&A1, 0uL) == MP_LT);
}
-
/* (Q0, R0) = RecursiveDivRem(A1 / beta^(k), B1) */
if ((err = mp_div_2d(&A1, k * MP_DIGIT_BIT, &A1, &t)) != MP_OKAY) goto LBL_ERR;
- if ((err = s_recursion(&A1, &B1, &Q0, &R0)) != MP_OKAY) goto LBL_ERR;
+ if ((err = s_recursion(&A1, &B1, &Q0, &R0)) != MP_OKAY) goto LBL_ERR;
/* A2 = (R0*beta^k) + (A1 % beta^k) - (Q0*B0) */
if ((err = mp_lshd(&R0, k)) != MP_OKAY) goto LBL_ERR;
@@ -65,13 +66,10 @@ static mp_err s_recursion(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
if ((err = mp_add(&A2, b, &A2)) != MP_OKAY) goto LBL_ERR;
}
/* return q = (Q1*beta^k) + Q0, r = A2 */
- if (q != NULL) {
- if ((err = mp_lshd(&Q1, k)) != MP_OKAY) goto LBL_ERR;
- if ((err = mp_add(&Q1, &Q0, q)) != MP_OKAY) goto LBL_ERR;
- }
- if (r != NULL) {
- if ((err = mp_copy(&A2, r)) != MP_OKAY) goto LBL_ERR;
- }
+ if ((err = mp_lshd(&Q1, k)) != MP_OKAY) goto LBL_ERR;
+ if ((err = mp_add(&Q1, &Q0, q)) != MP_OKAY) goto LBL_ERR;
+
+ if ((err = mp_copy(&A2, r)) != MP_OKAY) goto LBL_ERR;
LBL_ERR:
mp_clear_multi(&A1, &A2, &B1, &B0, &Q1, &Q0, &R1, &R0, &t, NULL);
@@ -94,24 +92,7 @@ mp_err s_mp_div_recursive(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
/* most significant bit of a limb */
/* assumes MP_DIGIT_MAX < (sizeof(mp_digit) * CHAR_BIT) */
msb = (MP_DIGIT_MAX + (mp_digit)(1)) >> 1;
-
- /*
- Method to compute sigma shamelessly stolen from
-
- J. Burnikel and J. Ziegler, "Fast recursive division", Research Report
- MPI-I-98-1-022, Max-Planck-Institut fuer Informatik, Saarbruecken, Germany,
- October 1998. (available online)
-
- Vid. section 2.3.
- */
- m = MP_MUL_KARATSUBA_CUTOFF;
- while (m <= b->used) {
- m <<= 1;
- }
- j = (b->used + m - 1) / m;
- n = j * m;
-
- sigma = MP_DIGIT_BIT * (n - b->used);
+ sigma = 0;
msb_b = b->dp[b->used - 1];
while (msb_b < msb) {
sigma++;
@@ -139,7 +120,7 @@ mp_err s_mp_div_recursive(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
/* Q = 0 */
mp_zero(&Q);
while (m > n) {
- /* (q, r) = RecursveDivRem(A / (beta^(m-n)), B) */
+ /* (q, r) = RecursiveDivRem(A / (beta^(m-n)), B) */
j = (m - n) * MP_DIGIT_BIT;
if ((err = mp_div_2d(&A, j, &A_div, &A_mod)) != MP_OKAY) goto LBL_ERR;
if ((err = s_recursion(&A_div, &B, &Q1, &R)) != MP_OKAY) goto LBL_ERR;
@@ -152,7 +133,7 @@ mp_err s_mp_div_recursive(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r
/* m = m - n */
m = m - n;
}
- /* (q, r) = RecursveDivRem(A, B) */
+ /* (q, r) = RecursiveDivRem(A, B) */
if ((err = s_recursion(&A, &B, &Q1, &R)) != MP_OKAY) goto LBL_ERR;
/* Q = (Q * beta^m) + q, R = r */
if ((err = mp_mul_2d(&Q, m * MP_DIGIT_BIT, &Q)) != MP_OKAY) goto LBL_ERR;