From e10a1c48270678096d6891ed049ce0e2d1656405 Mon Sep 17 00:00:00 2001 From: ssinai1 Date: Tue, 23 Aug 2022 13:52:34 +0200 Subject: [PATCH 1/4] Improve mp_root_n code --- mp_root_n.c | 31 +++++++------------------------ 1 file changed, 7 insertions(+), 24 deletions(-) diff --git a/mp_root_n.c b/mp_root_n.c index d904df883..47bb2f535 100644 --- a/mp_root_n.c +++ b/mp_root_n.c @@ -35,7 +35,7 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) a_ = *a; a_.sign = MP_ZPOS; - /* Compute seed: 2^(log_2(n)/b + 2)*/ + /* Compute seed: 2^(log_2(n)/b + 1)*/ ilog2 = mp_count_bits(a); /* @@ -65,8 +65,13 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) goto LBL_ERR; } /* Start value must be larger than root */ - ilog2 += 2; + ilog2 += 1; if ((err = mp_2expt(&t2,ilog2)) != MP_OKAY) goto LBL_ERR; + /* + Due to the convexity of the region and the upward rounding of the x[i] values, + the series will be monotonically decreasing and x[i] remain larger than + the exact value of the root (if the sarting value was larger). + */ do { /* t1 = t2 */ if ((err = mp_copy(&t2, &t1)) != MP_OKAY) goto LBL_ERR; @@ -92,31 +97,9 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) if ((err = mp_sub(&t1, &t3, &t2)) != MP_OKAY) goto LBL_ERR; - /* - Number of rounds is at most log_2(root). If it is more it - got stuck, so break out of the loop and do the rest manually. - */ - if (ilog2-- == 0) { - break; - } } while (mp_cmp(&t1, &t2) != MP_EQ); /* result can be off by a few so check */ - /* Loop beneath can overshoot by one if found root is smaller than actual root */ - for (;;) { - mp_ord cmp; - if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; - cmp = mp_cmp(&t2, &a_); - if (cmp == MP_EQ) { - err = MP_OKAY; - goto LBL_ERR; - } - if (cmp == MP_LT) { - if ((err = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY) goto LBL_ERR; - } else { - break; - } - } /* correct overshoot from above or from recurrence */ for (;;) { if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; From 36a055e3e41530701aec9323646f682c32f479bb Mon Sep 17 00:00:00 2001 From: ssinai1 Date: Wed, 24 Aug 2022 13:57:06 +0200 Subject: [PATCH 2/4] for loop put back (for developability) and fixes --- mp_root_n.c | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/mp_root_n.c b/mp_root_n.c index 47bb2f535..b0cd2a290 100644 --- a/mp_root_n.c +++ b/mp_root_n.c @@ -57,7 +57,7 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) err = MP_OKAY; goto LBL_ERR; } - ilog2 = ilog2 / b; + ilog2 = (ilog2 - 1) / b; if (ilog2 == 0) { mp_set(c, 1uL); c->sign = a->sign; @@ -100,6 +100,23 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) } while (mp_cmp(&t1, &t2) != MP_EQ); /* result can be off by a few so check */ + /* Loop beneath can overshoot by one if found root is smaller than actual root */ + for (;;) { + mp_ord cmp; + if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; + cmp = mp_cmp(&t2, &a_); + if (cmp == MP_EQ) { + err = MP_OKAY; + mp_exch(&t1, c); + c->sign = a->sign; + goto LBL_ERR; + } + if (cmp == MP_LT) { + if ((err = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY) goto LBL_ERR; + } else { + break; + } + } /* correct overshoot from above or from recurrence */ for (;;) { if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; From b36deb0cd8cd9979f496fc1ebe0a481e2b0a3527 Mon Sep 17 00:00:00 2001 From: ssinai1 Date: Thu, 25 Aug 2022 17:30:55 +0200 Subject: [PATCH 3/4] uses less exponentiation --- mp_root_n.c | 61 +++++++++++++++++++++++++---------------------------- 1 file changed, 29 insertions(+), 32 deletions(-) diff --git a/mp_root_n.c b/mp_root_n.c index b0cd2a290..aeba58863 100644 --- a/mp_root_n.c +++ b/mp_root_n.c @@ -14,7 +14,7 @@ */ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) { - mp_int t1, t2, t3, a_; + mp_int t1, t2, t3, a_, d; int ilog2; mp_err err; @@ -27,7 +27,7 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) return MP_VAL; } - if ((err = mp_init_multi(&t1, &t2, &t3, NULL)) != MP_OKAY) { + if ((err = mp_init_multi(&t1, &t2, &t3, &d, NULL)) != MP_OKAY) { return err; } @@ -66,17 +66,12 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) } /* Start value must be larger than root */ ilog2 += 1; - if ((err = mp_2expt(&t2,ilog2)) != MP_OKAY) goto LBL_ERR; - /* - Due to the convexity of the region and the upward rounding of the x[i] values, - the series will be monotonically decreasing and x[i] remain larger than - the exact value of the root (if the sarting value was larger). - */ + if ((err = mp_2expt(&t1, ilog2)) != MP_OKAY) goto LBL_ERR; + do { - /* t1 = t2 */ - if ((err = mp_copy(&t2, &t1)) != MP_OKAY) goto LBL_ERR; + mp_ord cmp; - /* t2 = t1 - ((t1**b - a) / (b * t1**(b-1))) */ + /* t2 = t1 - ceiling(((t1**b - a) / (b * t1**(b-1)))) */ /* t3 = t1**(b-1) */ if ((err = mp_expt_n(&t1, b - 1, &t3)) != MP_OKAY) goto LBL_ERR; @@ -85,6 +80,14 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) /* t2 = t1**b */ if ((err = mp_mul(&t3, &t1, &t2)) != MP_OKAY) goto LBL_ERR; + cmp = mp_cmp(&t2, &a_); + if (cmp == MP_EQ || cmp == MP_LT) { + err = MP_OKAY; + mp_exch(&t1, c); + c->sign = a->sign; + goto LBL_ERR; + } + /* t2 = t1**b - a */ if ((err = mp_sub(&t2, &a_, &t2)) != MP_OKAY) goto LBL_ERR; @@ -93,30 +96,24 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) if ((err = mp_mul_d(&t3, (mp_digit)b, &t3)) != MP_OKAY) goto LBL_ERR; /* t3 = (t1**b - a)/(b * t1**(b-1)) */ - if ((err = mp_div(&t2, &t3, &t3, NULL)) != MP_OKAY) goto LBL_ERR; + if ((err = mp_div(&t2, &t3, &t3, &d)) != MP_OKAY) goto LBL_ERR; + /* round up t3 - so t1 will be rounded down */ + if(mp_cmp_d(&d, 0) != MP_EQ) { + if ((err = mp_add_d(&t3, 1uL, &t3)) != MP_OKAY) goto LBL_ERR; + } + cmp = mp_cmp_d(&t3, 0); + if (cmp == MP_EQ || cmp == MP_LT) { + /* this should never happen */ + err = MP_ERR; + goto LBL_ERR; + } - if ((err = mp_sub(&t1, &t3, &t2)) != MP_OKAY) goto LBL_ERR; + /* t1 = t1 - t3 */ + if ((err = mp_sub(&t1, &t3, &t1)) != MP_OKAY) goto LBL_ERR; - } while (mp_cmp(&t1, &t2) != MP_EQ); + } while (mp_cmp_d(&t3, 1uL) == MP_GT); /* result can be off by a few so check */ - /* Loop beneath can overshoot by one if found root is smaller than actual root */ - for (;;) { - mp_ord cmp; - if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; - cmp = mp_cmp(&t2, &a_); - if (cmp == MP_EQ) { - err = MP_OKAY; - mp_exch(&t1, c); - c->sign = a->sign; - goto LBL_ERR; - } - if (cmp == MP_LT) { - if ((err = mp_add_d(&t1, 1uL, &t1)) != MP_OKAY) goto LBL_ERR; - } else { - break; - } - } /* correct overshoot from above or from recurrence */ for (;;) { if ((err = mp_expt_n(&t1, b, &t2)) != MP_OKAY) goto LBL_ERR; @@ -134,7 +131,7 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) c->sign = a->sign; LBL_ERR: - mp_clear_multi(&t1, &t2, &t3, NULL); + mp_clear_multi(&t1, &t2, &t3, &d, NULL); return err; } From 1c287f879da023e572fc008c0d5435c55905b818 Mon Sep 17 00:00:00 2001 From: ssinai1 Date: Fri, 26 Aug 2022 19:34:49 +0200 Subject: [PATCH 4/4] mp_cmp_d replaced --- mp_root_n.c | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/mp_root_n.c b/mp_root_n.c index aeba58863..e3c862b31 100644 --- a/mp_root_n.c +++ b/mp_root_n.c @@ -98,20 +98,15 @@ mp_err mp_root_n(const mp_int *a, int b, mp_int *c) /* t3 = (t1**b - a)/(b * t1**(b-1)) */ if ((err = mp_div(&t2, &t3, &t3, &d)) != MP_OKAY) goto LBL_ERR; /* round up t3 - so t1 will be rounded down */ - if(mp_cmp_d(&d, 0) != MP_EQ) { + if(!mp_iszero(&d)) { if ((err = mp_add_d(&t3, 1uL, &t3)) != MP_OKAY) goto LBL_ERR; } - cmp = mp_cmp_d(&t3, 0); - if (cmp == MP_EQ || cmp == MP_LT) { - /* this should never happen */ - err = MP_ERR; - goto LBL_ERR; - } /* t1 = t1 - t3 */ if ((err = mp_sub(&t1, &t3, &t1)) != MP_OKAY) goto LBL_ERR; - } while (mp_cmp_d(&t3, 1uL) == MP_GT); + /* while t3 != 1 */ + } while (!((t3.used == 1u) && (t3.dp[0] == 1u))); /* result can be off by a few so check */ /* correct overshoot from above or from recurrence */