divexact.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:8k
源码类别:

数学计算

开发平台:

Unix_Linux

  1. /* mpn_divexact(qp,np,nn,dp,dn,tp) -- Divide N = {np,nn} by D = {dp,dn} storing
  2.    the result in Q = {qp,nn-dn+1} expecting no remainder.  Overlap allowed
  3.    between Q and N; all other overlap disallowed.
  4.    Contributed to the GNU project by Torbjorn Granlund.
  5.    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
  6.    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  7.    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
  8. Copyright 2006, 2007, 2009 Free Software Foundation, Inc.
  9. This file is part of the GNU MP Library.
  10. The GNU MP Library is free software; you can redistribute it and/or modify
  11. it under the terms of the GNU Lesser General Public License as published by
  12. the Free Software Foundation; either version 3 of the License, or (at your
  13. option) any later version.
  14. The GNU MP Library is distributed in the hope that it will be useful, but
  15. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17. License for more details.
  18. You should have received a copy of the GNU Lesser General Public License
  19. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20. #include "gmp.h"
  21. #include "gmp-impl.h"
  22. #include "longlong.h"
  23. #if 1
  24. void
  25. mpn_divexact (mp_ptr qp,
  26.       mp_srcptr np, mp_size_t nn,
  27.       mp_srcptr dp, mp_size_t dn)
  28. {
  29.   unsigned shift;
  30.   mp_size_t qn;
  31.   mp_ptr tp, wp;
  32.   TMP_DECL;
  33.   ASSERT (dn > 0);
  34.   ASSERT (nn >= dn);
  35.   ASSERT (dp[dn-1] > 0);
  36.   while (dp[0] == 0)
  37.     {
  38.       ASSERT (np[0] == 0);
  39.       dp++;
  40.       np++;
  41.       dn--;
  42.       nn--;
  43.     }
  44.   if (dn == 1)
  45.     {
  46.       MPN_DIVREM_OR_DIVEXACT_1 (qp, np, nn, dp[0]);
  47.       return;
  48.     }
  49.   TMP_MARK;
  50.   qn = nn + 1 - dn;
  51.   count_trailing_zeros (shift, dp[0]);
  52.   if (shift > 0)
  53.     {
  54.       mp_size_t ss = (dn > qn) ? qn + 1 : dn;
  55.       tp = TMP_ALLOC_LIMBS (ss);
  56.       mpn_rshift (tp, dp, ss, shift);
  57.       dp = tp;
  58.       /* Since we have excluded dn == 1, we have nn > qn, and we need
  59.  to shift one limb beyond qn. */
  60.       wp = TMP_ALLOC_LIMBS (qn + 1);
  61.       mpn_rshift (wp, np, qn + 1, shift);
  62.     }
  63.   else
  64.     {
  65.       wp = TMP_ALLOC_LIMBS (qn);
  66.       MPN_COPY (wp, np, qn);
  67.     }
  68.   if (dn > qn)
  69.     dn = qn;
  70.   tp = TMP_ALLOC_LIMBS (mpn_bdiv_q_itch (qn, dn));
  71.   mpn_bdiv_q (qp, wp, qn, dp, dn, tp);
  72.   TMP_FREE;
  73. }
  74. #else
  75. /* We use the Jebelean's bidirectional exact division algorithm.  This is
  76.    somewhat naively implemented, with equal quotient parts done by 2-adic
  77.    division and truncating division.  Since 2-adic division is faster, it
  78.    should be used for a larger chunk.
  79.    This code is horrendously ugly, in all sorts of ways.
  80.    * It was hacked without much care or thought, but with a testing program.
  81.    * It handles scratch space frivolously, and furthermore the itch function
  82.      is broken.
  83.    * Doesn't provide any measures to deal with mu_divappr_q's +3 error.  We
  84.      have yet to provoke an error due to this, though.
  85.    * Algorithm selection leaves a lot to be desired.  In particular, the choice
  86.      between DC and MU isn't a point, but we treat it like one.
  87.    * It makes the msb part 1 or 2 limbs larger than the lsb part, in spite of
  88.      that the latter is faster.  We should at least reverse this, but perhaps
  89.      we should make the lsb part considerably larger.  (How do we tune this?)
  90. */
  91. mp_size_t
  92. mpn_divexact_itch (mp_size_t nn, mp_size_t dn)
  93. {
  94.   return nn + dn; /* FIXME this is not right */
  95. }
  96. void
  97. mpn_divexact (mp_ptr qp,
  98.       mp_srcptr np, mp_size_t nn,
  99.       mp_srcptr dp, mp_size_t dn,
  100.       mp_ptr scratch)
  101. {
  102.   mp_size_t qn;
  103.   mp_size_t nn0, qn0;
  104.   mp_size_t nn1, qn1;
  105.   mp_ptr tp;
  106.   mp_limb_t qml;
  107.   mp_limb_t qh;
  108.   int cnt;
  109.   mp_ptr xdp;
  110.   mp_limb_t di;
  111.   mp_limb_t cy;
  112.   gmp_pi1_t dinv;
  113.   TMP_DECL;
  114.   TMP_MARK;
  115.   qn = nn - dn + 1;
  116.   /* For small divisors, and small quotients, don't use Jebelean's algorithm. */
  117.   if (dn < DIVEXACT_JEB_THRESHOLD || qn < DIVEXACT_JEB_THRESHOLD)
  118.     {
  119.       tp = scratch;
  120.       MPN_COPY (tp, np, qn);
  121.       binvert_limb (di, dp[0]);  di = -di;
  122.       dn = MIN (dn, qn);
  123.       mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di);
  124.       TMP_FREE;
  125.       return;
  126.     }
  127.   qn0 = ((nn - dn) >> 1) + 1; /* low quotient size */
  128.   /* If quotient is much larger than the divisor, the bidirectional algorithm
  129.      does not work as currently implemented.  Fall back to plain bdiv.  */
  130.   if (qn0 > dn)
  131.     {
  132.       if (BELOW_THRESHOLD (dn, DC_BDIV_Q_THRESHOLD))
  133. {
  134.   tp = scratch;
  135.   MPN_COPY (tp, np, qn);
  136.   binvert_limb (di, dp[0]);  di = -di;
  137.   dn = MIN (dn, qn);
  138.   mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di);
  139. }
  140.       else if (BELOW_THRESHOLD (dn, MU_BDIV_Q_THRESHOLD))
  141. {
  142.   tp = scratch;
  143.   MPN_COPY (tp, np, qn);
  144.   binvert_limb (di, dp[0]);  di = -di;
  145.   mpn_dcpi1_bdiv_q (qp, tp, qn, dp, dn, di);
  146. }
  147.       else
  148. {
  149.   mpn_mu_bdiv_q (qp, np, qn, dp, dn, scratch);
  150. }
  151.       TMP_FREE;
  152.       return;
  153.     }
  154.   nn0 = qn0 + qn0;
  155.   nn1 = nn0 - 1 + ((nn-dn) & 1);
  156.   qn1 = qn0;
  157.   if (LIKELY (qn0 != dn))
  158.     {
  159.       nn1 = nn1 + 1;
  160.       qn1 = qn1 + 1;
  161.       if (UNLIKELY (dp[dn - 1] == 1 && qn1 != dn))
  162. {
  163.   /* If the leading divisor limb == 1, i.e. has just one bit, we have
  164.      to include an extra limb in order to get the needed overlap.  */
  165.   /* FIXME: Now with the mu_divappr_q function, we should really need
  166.      more overlap. That indicates one of two things: (1) The test code
  167.      is not good. (2) We actually overlap too much by default.  */
  168.   nn1 = nn1 + 1;
  169.   qn1 = qn1 + 1;
  170. }
  171.     }
  172.   tp = TMP_ALLOC_LIMBS (nn1 + 1);
  173.   count_leading_zeros (cnt, dp[dn - 1]);
  174.   /* Normalize divisor, store into tmp area.  */
  175.   if (cnt != 0)
  176.     {
  177.       xdp = TMP_ALLOC_LIMBS (qn1);
  178.       mpn_lshift (xdp, dp + dn - qn1, qn1, cnt);
  179.     }
  180.   else
  181.     {
  182.       xdp = (mp_ptr) dp + dn - qn1;
  183.     }
  184.   /* Shift dividend according to the divisor normalization.  */
  185.   /* FIXME: We compute too much here for XX_divappr_q, but these functions'
  186.      interfaces want a pointer to the imaginative least significant limb, not
  187.      to the least significant *used* limb.  Of course, we could leave nn1-qn1
  188.      rubbish limbs in the low part, to save some time.  */
  189.   if (cnt != 0)
  190.     {
  191.       cy = mpn_lshift (tp, np + nn - nn1, nn1, cnt);
  192.       if (cy != 0)
  193. {
  194.   tp[nn1] = cy;
  195.   nn1++;
  196. }
  197.     }
  198.   else
  199.     {
  200.       /* FIXME: This copy is not needed for mpn_mu_divappr_q, except when the
  201.  mpn_sub_n right before is executed.  */
  202.       MPN_COPY (tp, np + nn - nn1, nn1);
  203.     }
  204.   invert_pi1 (dinv, xdp[qn1 - 1], xdp[qn1 - 2]);
  205.   if (BELOW_THRESHOLD (qn1, DC_DIVAPPR_Q_THRESHOLD))
  206.     {
  207.       qp[qn0 - 1 + nn1 - qn1] = mpn_sbpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, dinv.inv32);
  208.     }
  209.   else if (BELOW_THRESHOLD (qn1, MU_DIVAPPR_Q_THRESHOLD))
  210.     {
  211.       qp[qn0 - 1 + nn1 - qn1] = mpn_dcpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, &dinv);
  212.     }
  213.   else
  214.     {
  215.       /* FIXME: mpn_mu_divappr_q doesn't handle qh != 0.  Work around it with a
  216.  conditional subtraction here.  */
  217.       qh = mpn_cmp (tp + nn1 - qn1, xdp, qn1) >= 0;
  218.       if (qh)
  219. mpn_sub_n (tp + nn1 - qn1, tp + nn1 - qn1, xdp, qn1);
  220.       mpn_mu_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, scratch);
  221.       qp[qn0 - 1 + nn1 - qn1] = qh;
  222.     }
  223.   qml = qp[qn0 - 1];
  224.   binvert_limb (di, dp[0]);  di = -di;
  225.   if (BELOW_THRESHOLD (qn0, DC_BDIV_Q_THRESHOLD))
  226.     {
  227.       MPN_COPY (tp, np, qn0);
  228.       mpn_sbpi1_bdiv_q (qp, tp, qn0, dp, qn0, di);
  229.     }
  230.   else if (BELOW_THRESHOLD (qn0, MU_BDIV_Q_THRESHOLD))
  231.     {
  232.       MPN_COPY (tp, np, qn0);
  233.       mpn_dcpi1_bdiv_q (qp, tp, qn0, dp, qn0, di);
  234.     }
  235.   else
  236.     {
  237.       mpn_mu_bdiv_q (qp, np, qn0, dp, qn0, scratch);
  238.     }
  239.   if (qml < qp[qn0 - 1])
  240.     mpn_decr_u (qp + qn0, 1);
  241.   TMP_FREE;
  242. }
  243. #endif