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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_tdiv_qr -- Divide the numerator (np,nn) by the denominator (dp,dn) and
  2.    write the nn-dn+1 quotient limbs at qp and the dn remainder limbs at rp.  If
  3.    qxn is non-zero, generate that many fraction limbs and append them after the
  4.    other quotient limbs, and update the remainder accordingly.  The input
  5.    operands are unaffected.
  6.    Preconditions:
  7.    1. The most significant limb of of the divisor must be non-zero.
  8.    2. nn >= dn, even if qxn is non-zero.  (??? relax this ???)
  9.    The time complexity of this is O(qn*qn+M(dn,qn)), where M(m,n) is the time
  10.    complexity of multiplication.
  11. Copyright 1997, 2000, 2001, 2002, 2005, 2009 Free Software Foundation, Inc.
  12. This file is part of the GNU MP Library.
  13. The GNU MP Library is free software; you can redistribute it and/or modify
  14. it under the terms of the GNU Lesser General Public License as published by
  15. the Free Software Foundation; either version 3 of the License, or (at your
  16. option) any later version.
  17. The GNU MP Library is distributed in the hope that it will be useful, but
  18. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  19. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  20. License for more details.
  21. You should have received a copy of the GNU Lesser General Public License
  22. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  23. #include "gmp.h"
  24. #include "gmp-impl.h"
  25. #include "longlong.h"
  26. void
  27. mpn_tdiv_qr (mp_ptr qp, mp_ptr rp, mp_size_t qxn,
  28.      mp_srcptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
  29. {
  30.   ASSERT_ALWAYS (qxn == 0);
  31.   ASSERT (nn >= 0);
  32.   ASSERT (dn >= 0);
  33.   ASSERT (dn == 0 || dp[dn - 1] != 0);
  34.   ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, np, nn));
  35.   ASSERT (! MPN_OVERLAP_P (qp, nn - dn + 1 + qxn, dp, dn));
  36.   switch (dn)
  37.     {
  38.     case 0:
  39.       DIVIDE_BY_ZERO;
  40.     case 1:
  41.       {
  42. rp[0] = mpn_divrem_1 (qp, (mp_size_t) 0, np, nn, dp[0]);
  43. return;
  44.       }
  45.     case 2:
  46.       {
  47. mp_ptr n2p, d2p;
  48. mp_limb_t qhl, cy;
  49. TMP_DECL;
  50. TMP_MARK;
  51. if ((dp[1] & GMP_NUMB_HIGHBIT) == 0)
  52.   {
  53.     int cnt;
  54.     mp_limb_t dtmp[2];
  55.     count_leading_zeros (cnt, dp[1]);
  56.     cnt -= GMP_NAIL_BITS;
  57.     d2p = dtmp;
  58.     d2p[1] = (dp[1] << cnt) | (dp[0] >> (GMP_NUMB_BITS - cnt));
  59.     d2p[0] = (dp[0] << cnt) & GMP_NUMB_MASK;
  60.     n2p = TMP_ALLOC_LIMBS (nn + 1);
  61.     cy = mpn_lshift (n2p, np, nn, cnt);
  62.     n2p[nn] = cy;
  63.     qhl = mpn_divrem_2 (qp, 0L, n2p, nn + (cy != 0), d2p);
  64.     if (cy == 0)
  65.       qp[nn - 2] = qhl; /* always store nn-2+1 quotient limbs */
  66.     rp[0] = (n2p[0] >> cnt)
  67.       | ((n2p[1] << (GMP_NUMB_BITS - cnt)) & GMP_NUMB_MASK);
  68.     rp[1] = (n2p[1] >> cnt);
  69.   }
  70. else
  71.   {
  72.     d2p = (mp_ptr) dp;
  73.     n2p = TMP_ALLOC_LIMBS (nn);
  74.     MPN_COPY (n2p, np, nn);
  75.     qhl = mpn_divrem_2 (qp, 0L, n2p, nn, d2p);
  76.     qp[nn - 2] = qhl; /* always store nn-2+1 quotient limbs */
  77.     rp[0] = n2p[0];
  78.     rp[1] = n2p[1];
  79.   }
  80. TMP_FREE;
  81. return;
  82.       }
  83.     default:
  84.       {
  85. int adjust;
  86. gmp_pi1_t dinv;
  87. TMP_DECL;
  88. TMP_MARK;
  89. adjust = np[nn - 1] >= dp[dn - 1]; /* conservative tests for quotient size */
  90. if (nn + adjust >= 2 * dn)
  91.   {
  92.     mp_ptr n2p, d2p;
  93.     mp_limb_t cy;
  94.     int cnt;
  95.     qp[nn - dn] = 0;   /* zero high quotient limb */
  96.     if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0) /* normalize divisor */
  97.       {
  98. count_leading_zeros (cnt, dp[dn - 1]);
  99. cnt -= GMP_NAIL_BITS;
  100. d2p = TMP_ALLOC_LIMBS (dn);
  101. mpn_lshift (d2p, dp, dn, cnt);
  102. n2p = TMP_ALLOC_LIMBS (nn + 1);
  103. cy = mpn_lshift (n2p, np, nn, cnt);
  104. n2p[nn] = cy;
  105. nn += adjust;
  106.       }
  107.     else
  108.       {
  109. cnt = 0;
  110. d2p = (mp_ptr) dp;
  111. n2p = TMP_ALLOC_LIMBS (nn + 1);
  112. MPN_COPY (n2p, np, nn);
  113. n2p[nn] = 0;
  114. nn += adjust;
  115.       }
  116.     invert_pi1 (dinv, d2p[dn - 1], d2p[dn - 2]);
  117.     if (BELOW_THRESHOLD (dn, DC_DIV_QR_THRESHOLD))
  118.       mpn_sbpi1_div_qr (qp, n2p, nn, d2p, dn, dinv.inv32);
  119.     else if (BELOW_THRESHOLD (dn, MUPI_DIV_QR_THRESHOLD) ||   /* fast condition */
  120.      BELOW_THRESHOLD (nn, 2 * MU_DIV_QR_THRESHOLD) || /* fast condition */
  121.      (double) (2 * (MU_DIV_QR_THRESHOLD - MUPI_DIV_QR_THRESHOLD)) * dn /* slow... */
  122.      + (double) MUPI_DIV_QR_THRESHOLD * nn > (double) dn * nn)    /* ...condition */
  123.       mpn_dcpi1_div_qr (qp, n2p, nn, d2p, dn, &dinv);
  124.     else
  125.       {
  126. mp_size_t itch = mpn_mu_div_qr_itch (nn, dn, 0);
  127. mp_ptr scratch = TMP_ALLOC_LIMBS (itch);
  128. mpn_mu_div_qr (qp, rp, n2p, nn, d2p, dn, scratch);
  129. n2p = rp;
  130.       }
  131.     if (cnt != 0)
  132.       mpn_rshift (rp, n2p, dn, cnt);
  133.     else
  134.       MPN_COPY (rp, n2p, dn);
  135.     TMP_FREE;
  136.     return;
  137.   }
  138. /* When we come here, the numerator/partial remainder is less
  139.    than twice the size of the denominator.  */
  140.   {
  141.     /* Problem:
  142.        Divide a numerator N with nn limbs by a denominator D with dn
  143.        limbs forming a quotient of qn=nn-dn+1 limbs.  When qn is small
  144.        compared to dn, conventional division algorithms perform poorly.
  145.        We want an algorithm that has an expected running time that is
  146.        dependent only on qn.
  147.        Algorithm (very informally stated):
  148.        1) Divide the 2 x qn most significant limbs from the numerator
  149.   by the qn most significant limbs from the denominator.  Call
  150.   the result qest.  This is either the correct quotient, but
  151.   might be 1 or 2 too large.  Compute the remainder from the
  152.   division.  (This step is implemented by a mpn_divrem call.)
  153.        2) Is the most significant limb from the remainder < p, where p
  154.   is the product of the most significant limb from the quotient
  155.   and the next(d)?  (Next(d) denotes the next ignored limb from
  156.   the denominator.)  If it is, decrement qest, and adjust the
  157.   remainder accordingly.
  158.        3) Is the remainder >= qest?  If it is, qest is the desired
  159.   quotient.  The algorithm terminates.
  160.        4) Subtract qest x next(d) from the remainder.  If there is
  161.   borrow out, decrement qest, and adjust the remainder
  162.   accordingly.
  163.        5) Skip one word from the denominator (i.e., let next(d) denote
  164.   the next less significant limb.  */
  165.     mp_size_t qn;
  166.     mp_ptr n2p, d2p;
  167.     mp_ptr tp;
  168.     mp_limb_t cy;
  169.     mp_size_t in, rn;
  170.     mp_limb_t quotient_too_large;
  171.     unsigned int cnt;
  172.     qn = nn - dn;
  173.     qp[qn] = 0; /* zero high quotient limb */
  174.     qn += adjust; /* qn cannot become bigger */
  175.     if (qn == 0)
  176.       {
  177. MPN_COPY (rp, np, dn);
  178. TMP_FREE;
  179. return;
  180.       }
  181.     in = dn - qn; /* (at least partially) ignored # of limbs in ops */
  182.     /* Normalize denominator by shifting it to the left such that its
  183.        most significant bit is set.  Then shift the numerator the same
  184.        amount, to mathematically preserve quotient.  */
  185.     if ((dp[dn - 1] & GMP_NUMB_HIGHBIT) == 0)
  186.       {
  187. count_leading_zeros (cnt, dp[dn - 1]);
  188. cnt -= GMP_NAIL_BITS;
  189. d2p = TMP_ALLOC_LIMBS (qn);
  190. mpn_lshift (d2p, dp + in, qn, cnt);
  191. d2p[0] |= dp[in - 1] >> (GMP_NUMB_BITS - cnt);
  192. n2p = TMP_ALLOC_LIMBS (2 * qn + 1);
  193. cy = mpn_lshift (n2p, np + nn - 2 * qn, 2 * qn, cnt);
  194. if (adjust)
  195.   {
  196.     n2p[2 * qn] = cy;
  197.     n2p++;
  198.   }
  199. else
  200.   {
  201.     n2p[0] |= np[nn - 2 * qn - 1] >> (GMP_NUMB_BITS - cnt);
  202.   }
  203.       }
  204.     else
  205.       {
  206. cnt = 0;
  207. d2p = (mp_ptr) dp + in;
  208. n2p = TMP_ALLOC_LIMBS (2 * qn + 1);
  209. MPN_COPY (n2p, np + nn - 2 * qn, 2 * qn);
  210. if (adjust)
  211.   {
  212.     n2p[2 * qn] = 0;
  213.     n2p++;
  214.   }
  215.       }
  216.     /* Get an approximate quotient using the extracted operands.  */
  217.     if (qn == 1)
  218.       {
  219. mp_limb_t q0, r0;
  220. udiv_qrnnd (q0, r0, n2p[1], n2p[0] << GMP_NAIL_BITS, d2p[0] << GMP_NAIL_BITS);
  221. n2p[0] = r0 >> GMP_NAIL_BITS;
  222. qp[0] = q0;
  223.       }
  224.     else if (qn == 2)
  225.       mpn_divrem_2 (qp, 0L, n2p, 4L, d2p); /* FIXME: obsolete function */
  226.     else
  227.       {
  228. invert_pi1 (dinv, d2p[qn - 1], d2p[qn - 2]);
  229. if (BELOW_THRESHOLD (qn, DC_DIV_QR_THRESHOLD))
  230.   mpn_sbpi1_div_qr (qp, n2p, 2 * qn, d2p, qn, dinv.inv32);
  231. else if (BELOW_THRESHOLD (qn, MU_DIV_QR_THRESHOLD))
  232.   mpn_dcpi1_div_qr (qp, n2p, 2 * qn, d2p, qn, &dinv);
  233. else
  234.   {
  235.     mp_size_t itch = mpn_mu_div_qr_itch (2 * qn, qn, 0);
  236.     mp_ptr scratch = TMP_ALLOC_LIMBS (itch);
  237.     mp_ptr r2p = rp;
  238.     if (np == r2p) /* If N and R share space, put ... */
  239.       r2p += nn - qn; /* intermediate remainder at N's upper end. */
  240.     mpn_mu_div_qr (qp, r2p, n2p, 2 * qn, d2p, qn, scratch);
  241.     MPN_COPY (n2p, r2p, qn);
  242.   }
  243.       }
  244.     rn = qn;
  245.     /* Multiply the first ignored divisor limb by the most significant
  246.        quotient limb.  If that product is > the partial remainder's
  247.        most significant limb, we know the quotient is too large.  This
  248.        test quickly catches most cases where the quotient is too large;
  249.        it catches all cases where the quotient is 2 too large.  */
  250.     {
  251.       mp_limb_t dl, x;
  252.       mp_limb_t h, dummy;
  253.       if (in - 2 < 0)
  254. dl = 0;
  255.       else
  256. dl = dp[in - 2];
  257. #if GMP_NAIL_BITS == 0
  258.       x = (dp[in - 1] << cnt) | ((dl >> 1) >> ((~cnt) % GMP_LIMB_BITS));
  259. #else
  260.       x = (dp[in - 1] << cnt) & GMP_NUMB_MASK;
  261.       if (cnt != 0)
  262. x |= dl >> (GMP_NUMB_BITS - cnt);
  263. #endif
  264.       umul_ppmm (h, dummy, x, qp[qn - 1] << GMP_NAIL_BITS);
  265.       if (n2p[qn - 1] < h)
  266. {
  267.   mp_limb_t cy;
  268.   mpn_decr_u (qp, (mp_limb_t) 1);
  269.   cy = mpn_add_n (n2p, n2p, d2p, qn);
  270.   if (cy)
  271.     {
  272.       /* The partial remainder is safely large.  */
  273.       n2p[qn] = cy;
  274.       ++rn;
  275.     }
  276. }
  277.     }
  278.     quotient_too_large = 0;
  279.     if (cnt != 0)
  280.       {
  281. mp_limb_t cy1, cy2;
  282. /* Append partially used numerator limb to partial remainder.  */
  283. cy1 = mpn_lshift (n2p, n2p, rn, GMP_NUMB_BITS - cnt);
  284. n2p[0] |= np[in - 1] & (GMP_NUMB_MASK >> cnt);
  285. /* Update partial remainder with partially used divisor limb.  */
  286. cy2 = mpn_submul_1 (n2p, qp, qn, dp[in - 1] & (GMP_NUMB_MASK >> cnt));
  287. if (qn != rn)
  288.   {
  289.     ASSERT_ALWAYS (n2p[qn] >= cy2);
  290.     n2p[qn] -= cy2;
  291.   }
  292. else
  293.   {
  294.     n2p[qn] = cy1 - cy2; /* & GMP_NUMB_MASK; */
  295.     quotient_too_large = (cy1 < cy2);
  296.     ++rn;
  297.   }
  298. --in;
  299.       }
  300.     /* True: partial remainder now is neutral, i.e., it is not shifted up.  */
  301.     tp = TMP_ALLOC_LIMBS (dn);
  302.     if (in < qn)
  303.       {
  304. if (in == 0)
  305.   {
  306.     MPN_COPY (rp, n2p, rn);
  307.     ASSERT_ALWAYS (rn == dn);
  308.     goto foo;
  309.   }
  310. mpn_mul (tp, qp, qn, dp, in);
  311.       }
  312.     else
  313.       mpn_mul (tp, dp, in, qp, qn);
  314.     cy = mpn_sub (n2p, n2p, rn, tp + in, qn);
  315.     MPN_COPY (rp + in, n2p, dn - in);
  316.     quotient_too_large |= cy;
  317.     cy = mpn_sub_n (rp, np, tp, in);
  318.     cy = mpn_sub_1 (rp + in, rp + in, rn, cy);
  319.     quotient_too_large |= cy;
  320.   foo:
  321.     if (quotient_too_large)
  322.       {
  323. mpn_decr_u (qp, (mp_limb_t) 1);
  324. mpn_add_n (rp, rp, dp, dn);
  325.       }
  326.   }
  327. TMP_FREE;
  328. return;
  329.       }
  330.     }
  331. }