cheng/mpir
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Improve tail of sb_divappr_q and fix bug.

Browse Source
This commit is contained in:
William Hart 2014-02-21 13:41:37 +00:00
parent e41d78e02b
commit e9cd157e50
1 changed files with 85 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

119
mpn/generic/sb_divappr_q.c
View File

@ -32,7 +32,7 @@ along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
#include "gmp-impl.h"
#include "longlong.h"
#define SB_DIVAPPR_Q_SMALL_THRESHOLD 6
#define SB_DIVAPPR_Q_SMALL_THRESHOLD 30
void __divappr_helper(mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t qn)
{
@ -46,40 +46,19 @@ void __divappr_helper(mp_ptr qp, mp_ptr np, mp_srcptr dp, mp_size_t qn)
}
}
mp_limb_t
mpn_sb_divappr_q (mp_ptr qp,
void
mpn_sb_divappr_q_small (mp_ptr qp,
mp_ptr np, mp_size_t nn,
mp_srcptr dp, mp_size_t dn,
mp_limb_t dinv, mp_limb_t d1inv)
mp_limb_t d1inv, mp_size_t qn)
{
mp_limb_t qh;
mp_size_t qn, i;
mp_size_t i;
mp_limb_t n1, n0;
mp_limb_t d1, d0, r1, r2;
mp_limb_t cy, cy1;
mp_limb_t q;
mp_limb_t flag;
ASSERT (dn > 2);
ASSERT (nn >= dn);
ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
np += nn;
qn = nn - dn;
if (qn + 1 < dn)
{
dp += dn - (qn + 1);
dn = qn + 1;
}
qh = mpn_cmp (np - dn, dp, dn) >= 0;
if (qh != 0)
mpn_sub_n (np - dn, np - dn, dp, dn);
if (dn <= SB_DIVAPPR_Q_SMALL_THRESHOLD)
{
/* Reduce until dn - 2 >= qn */
/* Reduce until dn - 2 >= qn */
for (qn--, np--; qn > dn - 2; qn--)
{
/* fetch next word */
@ -117,7 +96,7 @@ mpn_sb_divappr_q (mp_ptr qp,
if (cy > dp[qn] || (cy == dp[qn] && mpn_cmp(np - qn + 1, dp, qn) >= 0))
{
__divappr_helper(qp, np - qn, dp, qn);
return qh;
return;
}
mpir_divapprox32_preinv2(q, cy, np[0], d1inv);
@ -138,8 +117,42 @@ mpn_sb_divappr_q (mp_ptr qp,
qp[qn - 1] = q;
dp++;
}
np[1] = cy;
}
np[1] = cy;
}
mp_limb_t
mpn_sb_divappr_q (mp_ptr qp,
mp_ptr np, mp_size_t nn,
mp_srcptr dp, mp_size_t dn,
mp_limb_t dinv, mp_limb_t d1inv)
{
mp_limb_t qh;
mp_size_t qn, i;
mp_limb_t n1, n0;
mp_limb_t d1, d0, r1, r2;
mp_limb_t cy, cy1;
mp_limb_t q;
ASSERT (dn > 2);
ASSERT (nn >= dn);
ASSERT ((dp[dn-1] & GMP_NUMB_HIGHBIT) != 0);
np += nn;
qn = nn - dn;
if (qn + 1 < dn)
{
dp += dn - (qn + 1);
dn = qn + 1;
}
qh = mpn_cmp (np - dn, dp, dn) >= 0;
if (qh != 0)
mpn_sub_n (np - dn, np - dn, dp, dn);
if (dn <= SB_DIVAPPR_Q_SMALL_THRESHOLD)
mpn_sb_divappr_q_small (qp, np, nn, dp, dn, d1inv, qn);
else
{
d1 = dp[dn - 1];
@ -179,7 +192,7 @@ mpn_sb_divappr_q (mp_ptr qp,
np--;
qn--;
for ( ; qn >= 0; qn--)
for ( ; qn > 0; qn--)
{
/* fetch next word */
cy = np[1];
@ -206,7 +219,7 @@ mpn_sb_divappr_q (mp_ptr qp,
udiv_qr_3by2(q, np[1], np[0], cy, np[1], np[0], d1, d0, dinv);
/* np -= dp*q */
cy1 = qn == 0 ? 0 : mpn_submul_1(np - qn, dp, qn, q);
cy1 = mpn_submul_1(np - qn, dp, qn, q);
sub_333(cy, np[1], np[0], 0, np[1], np[0], 0, 0, cy1);
}
@ -216,7 +229,7 @@ mpn_sb_divappr_q (mp_ptr qp,
udiv_qr_3by2(q, np[1], np[0], cy, np[1], np[0], d1, d0, dinv);
/* np -= dp*q */
cy1 = qn == 0 ? 0 : mpn_submul_1(np - qn, dp, qn, q);
cy1 = mpn_submul_1(np - qn, dp, qn, q);
sub_333(cy, np[1], np[0], 0, np[1], np[0], 0, 0, cy1);
}
@ -230,6 +243,44 @@ mpn_sb_divappr_q (mp_ptr qp,
qp[qn] = q;
dp++;
}
/* fetch next word */
cy = np[1];
np--;
/* rare case where truncation ruins normalisation */
if (UNLIKELY(cy >= d1))
{
if (cy > d1 || (cy == d1 && mpn_cmp(np - qn + 1, dp, qn + 1) >= 0))
{
__divappr_helper(qp, np - qn, dp, qn + 1);
return qh;
}
if (np[1] >= d0)
{
q = ~CNST_LIMB(0);
/* np -= dp*q */
cy -= mpn_submul_1(np - qn, dp, qn + 2, q);
} else
udiv_qr_3by2(q, np[1], np[0], cy, np[1], np[0], d1, d0, dinv);
}
else
udiv_qr_3by2(q, np[1], np[0], cy, np[1], np[0], d1, d0, dinv);
/* correct if quotient is too large */
if (UNLIKELY(cy != 0))
{
q--;
cy += mpn_add_n(np - qn, np - qn, dp, qn + 2);
}
qp[qn] = q;
dp++;
np[2] = cy;
}