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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_divrem_1 -- mpn by limb division.
  2. Copyright 1991, 1993, 1994, 1996, 1998, 1999, 2000, 2002, 2003 Free Software
  3. Foundation, Inc.
  4. This file is part of the GNU MP Library.
  5. The GNU MP Library is free software; you can redistribute it and/or modify
  6. it under the terms of the GNU Lesser General Public License as published by
  7. the Free Software Foundation; either version 3 of the License, or (at your
  8. option) any later version.
  9. The GNU MP Library is distributed in the hope that it will be useful, but
  10. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  11. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  12. License for more details.
  13. You should have received a copy of the GNU Lesser General Public License
  14. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  15. #include "gmp.h"
  16. #include "gmp-impl.h"
  17. #include "longlong.h"
  18. /* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
  19.    meaning the quotient size where that should happen, the quotient size
  20.    being how many udiv divisions will be done.
  21.    The default is to use preinv always, CPUs where this doesn't suit have
  22.    tuned thresholds.  Note in particular that preinv should certainly be
  23.    used if that's the only division available (USE_PREINV_ALWAYS).  */
  24. #ifndef DIVREM_1_NORM_THRESHOLD
  25. #define DIVREM_1_NORM_THRESHOLD  0
  26. #endif
  27. #ifndef DIVREM_1_UNNORM_THRESHOLD
  28. #define DIVREM_1_UNNORM_THRESHOLD  0
  29. #endif
  30. /* If the cpu only has multiply-by-inverse division (eg. alpha), then NORM
  31.    and UNNORM thresholds are 0 and only the inversion code is included.
  32.    If multiply-by-inverse is never viable, then NORM and UNNORM thresholds
  33.    will be MP_SIZE_T_MAX and only the plain division code is included.
  34.    Otherwise mul-by-inverse is better than plain division above some
  35.    threshold, and best results are obtained by having code for both present.
  36.    The main reason for separating the norm and unnorm cases is that not all
  37.    CPUs give zero for "n0 >> GMP_LIMB_BITS" which would arise in the unnorm
  38.    code used on an already normalized divisor.
  39.    If UDIV_NEEDS_NORMALIZATION is false then plain division uses the same
  40.    non-shifting code for both the norm and unnorm cases, though with
  41.    different criteria for skipping a division, and with different thresholds
  42.    of course.  And in fact if inversion is never viable, then that simple
  43.    non-shifting division would be all that's left.
  44.    The NORM and UNNORM thresholds might not differ much, but if there's
  45.    going to be separate code for norm and unnorm then it makes sense to have
  46.    separate thresholds.  One thing that's possible is that the
  47.    mul-by-inverse might be better only for normalized divisors, due to that
  48.    case not needing variable bit shifts.
  49.    Notice that the thresholds are tested after the decision to possibly skip
  50.    one divide step, so they're based on the actual number of divisions done.
  51.    For the unnorm case, it would be possible to call mpn_lshift to adjust
  52.    the dividend all in one go (into the quotient space say), rather than
  53.    limb-by-limb in the loop.  This might help if mpn_lshift is a lot faster
  54.    than what the compiler can generate for EXTRACT.  But this is left to CPU
  55.    specific implementations to consider, especially since EXTRACT isn't on
  56.    the dependent chain.  */
  57. mp_limb_t
  58. mpn_divrem_1 (mp_ptr qp, mp_size_t qxn,
  59.       mp_srcptr up, mp_size_t un, mp_limb_t d)
  60. {
  61.   mp_size_t  n;
  62.   mp_size_t  i;
  63.   mp_limb_t  n1, n0;
  64.   mp_limb_t  r = 0;
  65.   ASSERT (qxn >= 0);
  66.   ASSERT (un >= 0);
  67.   ASSERT (d != 0);
  68.   /* FIXME: What's the correct overlap rule when qxn!=0? */
  69.   ASSERT (MPN_SAME_OR_SEPARATE_P (qp+qxn, up, un));
  70.   n = un + qxn;
  71.   if (n == 0)
  72.     return 0;
  73.   d <<= GMP_NAIL_BITS;
  74.   qp += (n - 1);   /* Make qp point at most significant quotient limb */
  75.   if ((d & GMP_LIMB_HIGHBIT) != 0)
  76.     {
  77.       if (un != 0)
  78. {
  79.   /* High quotient limb is 0 or 1, skip a divide step. */
  80.   mp_limb_t q;
  81.   r = up[un - 1] << GMP_NAIL_BITS;
  82.   q = (r >= d);
  83.   *qp-- = q;
  84.   r -= (d & -q);
  85.   r >>= GMP_NAIL_BITS;
  86.   n--;
  87.   un--;
  88. }
  89.       if (BELOW_THRESHOLD (n, DIVREM_1_NORM_THRESHOLD))
  90. {
  91. plain:
  92.   for (i = un - 1; i >= 0; i--)
  93.     {
  94.       n0 = up[i] << GMP_NAIL_BITS;
  95.       udiv_qrnnd (*qp, r, r, n0, d);
  96.       r >>= GMP_NAIL_BITS;
  97.       qp--;
  98.     }
  99.   for (i = qxn - 1; i >= 0; i--)
  100.     {
  101.       udiv_qrnnd (*qp, r, r, CNST_LIMB(0), d);
  102.       r >>= GMP_NAIL_BITS;
  103.       qp--;
  104.     }
  105.   return r;
  106. }
  107.       else
  108. {
  109.   /* Multiply-by-inverse, divisor already normalized. */
  110.   mp_limb_t dinv;
  111.   invert_limb (dinv, d);
  112.   for (i = un - 1; i >= 0; i--)
  113.     {
  114.       n0 = up[i] << GMP_NAIL_BITS;
  115.       udiv_qrnnd_preinv (*qp, r, r, n0, d, dinv);
  116.       r >>= GMP_NAIL_BITS;
  117.       qp--;
  118.     }
  119.   for (i = qxn - 1; i >= 0; i--)
  120.     {
  121.       udiv_qrnnd_preinv (*qp, r, r, CNST_LIMB(0), d, dinv);
  122.       r >>= GMP_NAIL_BITS;
  123.       qp--;
  124.     }
  125.   return r;
  126. }
  127.     }
  128.   else
  129.     {
  130.       /* Most significant bit of divisor == 0.  */
  131.       int norm;
  132.       /* Skip a division if high < divisor (high quotient 0).  Testing here
  133.  before normalizing will still skip as often as possible.  */
  134.       if (un != 0)
  135. {
  136.   n1 = up[un - 1] << GMP_NAIL_BITS;
  137.   if (n1 < d)
  138.     {
  139.       r = n1 >> GMP_NAIL_BITS;
  140.       *qp-- = 0;
  141.       n--;
  142.       if (n == 0)
  143. return r;
  144.       un--;
  145.     }
  146. }
  147.       if (! UDIV_NEEDS_NORMALIZATION
  148.   && BELOW_THRESHOLD (n, DIVREM_1_UNNORM_THRESHOLD))
  149. goto plain;
  150.       count_leading_zeros (norm, d);
  151.       d <<= norm;
  152.       r <<= norm;
  153.       if (UDIV_NEEDS_NORMALIZATION
  154.   && BELOW_THRESHOLD (n, DIVREM_1_UNNORM_THRESHOLD))
  155. {
  156.   if (un != 0)
  157.     {
  158.       n1 = up[un - 1] << GMP_NAIL_BITS;
  159.       r |= (n1 >> (GMP_LIMB_BITS - norm));
  160.       for (i = un - 2; i >= 0; i--)
  161. {
  162.   n0 = up[i] << GMP_NAIL_BITS;
  163.   udiv_qrnnd (*qp, r, r,
  164.       (n1 << norm) | (n0 >> (GMP_NUMB_BITS - norm)),
  165.       d);
  166.   r >>= GMP_NAIL_BITS;
  167.   qp--;
  168.   n1 = n0;
  169. }
  170.       udiv_qrnnd (*qp, r, r, n1 << norm, d);
  171.       r >>= GMP_NAIL_BITS;
  172.       qp--;
  173.     }
  174.   for (i = qxn - 1; i >= 0; i--)
  175.     {
  176.       udiv_qrnnd (*qp, r, r, CNST_LIMB(0), d);
  177.       r >>= GMP_NAIL_BITS;
  178.       qp--;
  179.     }
  180.   return r >> norm;
  181. }
  182.       else
  183. {
  184.   mp_limb_t  dinv;
  185.   invert_limb (dinv, d);
  186.   if (un != 0)
  187.     {
  188.       n1 = up[un - 1] << GMP_NAIL_BITS;
  189.       r |= (n1 >> (GMP_LIMB_BITS - norm));
  190.       for (i = un - 2; i >= 0; i--)
  191. {
  192.   n0 = up[i] << GMP_NAIL_BITS;
  193.   udiv_qrnnd_preinv (*qp, r, r,
  194.      ((n1 << norm) | (n0 >> (GMP_NUMB_BITS - norm))),
  195.      d, dinv);
  196.   r >>= GMP_NAIL_BITS;
  197.   qp--;
  198.   n1 = n0;
  199. }
  200.       udiv_qrnnd_preinv (*qp, r, r, n1 << norm, d, dinv);
  201.       r >>= GMP_NAIL_BITS;
  202.       qp--;
  203.     }
  204.   for (i = qxn - 1; i >= 0; i--)
  205.     {
  206.       udiv_qrnnd_preinv (*qp, r, r, CNST_LIMB(0), d, dinv);
  207.       r >>= GMP_NAIL_BITS;
  208.       qp--;
  209.     }
  210.   return r >> norm;
  211. }
  212.     }
  213. }