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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_mod_1(dividend_ptr, dividend_size, divisor_limb) --
  2.    Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
  3.    Return the single-limb remainder.
  4.    There are no constraints on the value of the divisor.
  5. Copyright 1991, 1993, 1994, 1999, 2000, 2002, 2007, 2008, 2009 Free
  6. Software Foundation, Inc.
  7. This file is part of the GNU MP Library.
  8. The GNU MP Library is free software; you can redistribute it and/or modify
  9. it under the terms of the GNU Lesser General Public License as published by
  10. the Free Software Foundation; either version 3 of the License, or (at your
  11. option) any later version.
  12. The GNU MP Library is distributed in the hope that it will be useful, but
  13. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15. License for more details.
  16. You should have received a copy of the GNU Lesser General Public License
  17. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  18. #include "gmp.h"
  19. #include "gmp-impl.h"
  20. #include "longlong.h"
  21. /* The size where udiv_qrnnd_preinv should be used rather than udiv_qrnnd,
  22.    meaning the quotient size where that should happen, the quotient size
  23.    being how many udiv divisions will be done.
  24.    The default is to use preinv always, CPUs where this doesn't suit have
  25.    tuned thresholds.  Note in particular that preinv should certainly be
  26.    used if that's the only division available (USE_PREINV_ALWAYS).  */
  27. #ifndef MOD_1_NORM_THRESHOLD
  28. #define MOD_1_NORM_THRESHOLD  0
  29. #endif
  30. #ifndef MOD_1_UNNORM_THRESHOLD
  31. #define MOD_1_UNNORM_THRESHOLD  0
  32. #endif
  33. #ifndef MOD_1U_TO_MOD_1_1_THRESHOLD
  34. #define MOD_1U_TO_MOD_1_1_THRESHOLD  MP_SIZE_T_MAX /* default is not to use mpn_mod_1s */
  35. #endif
  36. #ifndef MOD_1N_TO_MOD_1_1_THRESHOLD
  37. #define MOD_1N_TO_MOD_1_1_THRESHOLD  MP_SIZE_T_MAX /* default is not to use mpn_mod_1s */
  38. #endif
  39. #ifndef MOD_1_1_TO_MOD_1_2_THRESHOLD
  40. #define MOD_1_1_TO_MOD_1_2_THRESHOLD  10
  41. #endif
  42. #ifndef MOD_1_2_TO_MOD_1_4_THRESHOLD
  43. #define MOD_1_2_TO_MOD_1_4_THRESHOLD  20
  44. #endif
  45. /* The comments in mpn/generic/divrem_1.c apply here too.
  46.    As noted in the algorithms section of the manual, the shifts in the loop
  47.    for the unnorm case can be avoided by calculating r = a%(d*2^n), followed
  48.    by a final (r*2^n)%(d*2^n).  In fact if it happens that a%(d*2^n) can
  49.    skip a division where (a*2^n)%(d*2^n) can't then there's the same number
  50.    of divide steps, though how often that happens depends on the assumed
  51.    distributions of dividend and divisor.  In any case this idea is left to
  52.    CPU specific implementations to consider.  */
  53. static mp_limb_t
  54. mpn_mod_1_unnorm (mp_srcptr up, mp_size_t un, mp_limb_t d)
  55. {
  56.   mp_size_t  i;
  57.   mp_limb_t  n1, n0, r;
  58.   mp_limb_t  dummy;
  59.   int cnt;
  60.   ASSERT (un > 0);
  61.   ASSERT (d != 0);
  62.   d <<= GMP_NAIL_BITS;
  63.   /* Skip a division if high < divisor.  Having the test here before
  64.      normalizing will still skip as often as possible.  */
  65.   r = up[un - 1] << GMP_NAIL_BITS;
  66.   if (r < d)
  67.     {
  68.       r >>= GMP_NAIL_BITS;
  69.       un--;
  70.       if (un == 0)
  71. return r;
  72.     }
  73.   else
  74.     r = 0;
  75.   /* If udiv_qrnnd doesn't need a normalized divisor, can use the simple
  76.      code above. */
  77.   if (! UDIV_NEEDS_NORMALIZATION
  78.       && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
  79.     {
  80.       for (i = un - 1; i >= 0; i--)
  81. {
  82.   n0 = up[i] << GMP_NAIL_BITS;
  83.   udiv_qrnnd (dummy, r, r, n0, d);
  84.   r >>= GMP_NAIL_BITS;
  85. }
  86.       return r;
  87.     }
  88.   count_leading_zeros (cnt, d);
  89.   d <<= cnt;
  90.   n1 = up[un - 1] << GMP_NAIL_BITS;
  91.   r = (r << cnt) | (n1 >> (GMP_LIMB_BITS - cnt));
  92.   if (UDIV_NEEDS_NORMALIZATION
  93.       && BELOW_THRESHOLD (un, MOD_1_UNNORM_THRESHOLD))
  94.     {
  95.       for (i = un - 2; i >= 0; i--)
  96. {
  97.   n0 = up[i] << GMP_NAIL_BITS;
  98.   udiv_qrnnd (dummy, r, r,
  99.       (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt)),
  100.       d);
  101.   r >>= GMP_NAIL_BITS;
  102.   n1 = n0;
  103. }
  104.       udiv_qrnnd (dummy, r, r, n1 << cnt, d);
  105.       r >>= GMP_NAIL_BITS;
  106.       return r >> cnt;
  107.     }
  108.   else
  109.     {
  110.       mp_limb_t inv;
  111.       invert_limb (inv, d);
  112.       for (i = un - 2; i >= 0; i--)
  113. {
  114.   n0 = up[i] << GMP_NAIL_BITS;
  115.   udiv_qrnnd_preinv (dummy, r, r,
  116.      (n1 << cnt) | (n0 >> (GMP_NUMB_BITS - cnt)),
  117.      d, inv);
  118.   r >>= GMP_NAIL_BITS;
  119.   n1 = n0;
  120. }
  121.       udiv_qrnnd_preinv (dummy, r, r, n1 << cnt, d, inv);
  122.       r >>= GMP_NAIL_BITS;
  123.       return r >> cnt;
  124.     }
  125. }
  126. static mp_limb_t
  127. mpn_mod_1_norm (mp_srcptr up, mp_size_t un, mp_limb_t d)
  128. {
  129.   mp_size_t  i;
  130.   mp_limb_t  n0, r;
  131.   mp_limb_t  dummy;
  132.   ASSERT (un > 0);
  133.   d <<= GMP_NAIL_BITS;
  134.   ASSERT (d & GMP_LIMB_HIGHBIT);
  135.   /* High limb is initial remainder, possibly with one subtract of
  136.      d to get r<d.  */
  137.   r = up[un - 1] << GMP_NAIL_BITS;
  138.   if (r >= d)
  139.     r -= d;
  140.   r >>= GMP_NAIL_BITS;
  141.   un--;
  142.   if (un == 0)
  143.     return r;
  144.   if (BELOW_THRESHOLD (un, MOD_1_NORM_THRESHOLD))
  145.     {
  146.       for (i = un - 1; i >= 0; i--)
  147. {
  148.   n0 = up[i] << GMP_NAIL_BITS;
  149.   udiv_qrnnd (dummy, r, r, n0, d);
  150.   r >>= GMP_NAIL_BITS;
  151. }
  152.       return r;
  153.     }
  154.   else
  155.     {
  156.       mp_limb_t  inv;
  157.       invert_limb (inv, d);
  158.       for (i = un - 1; i >= 0; i--)
  159. {
  160.   n0 = up[i] << GMP_NAIL_BITS;
  161.   udiv_qrnnd_preinv (dummy, r, r, n0, d, inv);
  162.   r >>= GMP_NAIL_BITS;
  163. }
  164.       return r;
  165.     }
  166. }
  167. mp_limb_t
  168. mpn_mod_1 (mp_srcptr ap, mp_size_t n, mp_limb_t b)
  169. {
  170.   ASSERT (n >= 0);
  171.   ASSERT (b != 0);
  172.   /* Should this be handled at all?  Rely on callers?  Note un==0 is currently
  173.      required by mpz/fdiv_r_ui.c and possibly other places.  */
  174.   if (n == 0)
  175.     return 0;
  176.   if (UNLIKELY ((b & GMP_NUMB_HIGHBIT) != 0))
  177.     {
  178.       if (BELOW_THRESHOLD (n, MOD_1N_TO_MOD_1_1_THRESHOLD))
  179. {
  180.   return mpn_mod_1_norm (ap, n, b);
  181. }
  182.       else
  183. {
  184.   mp_limb_t pre[4];
  185.   mpn_mod_1_1p_cps (pre, b);
  186.   return mpn_mod_1_1p (ap, n, b, pre);
  187. }
  188.     }
  189.   else
  190.     {
  191.       if (BELOW_THRESHOLD (n, MOD_1U_TO_MOD_1_1_THRESHOLD))
  192. {
  193.   return mpn_mod_1_unnorm (ap, n, b);
  194. }
  195.       else if (BELOW_THRESHOLD (n, MOD_1_1_TO_MOD_1_2_THRESHOLD))
  196. {
  197.   mp_limb_t pre[4];
  198.   mpn_mod_1_1p_cps (pre, b);
  199.   return mpn_mod_1_1p (ap, n, b << pre[1], pre);
  200. }
  201.       else if (BELOW_THRESHOLD (n, MOD_1_2_TO_MOD_1_4_THRESHOLD) || UNLIKELY (b > GMP_NUMB_MASK / 4))
  202. {
  203.   mp_limb_t pre[5];
  204.   mpn_mod_1s_2p_cps (pre, b);
  205.   return mpn_mod_1s_2p (ap, n, b << pre[1], pre);
  206. }
  207.       else
  208. {
  209.   mp_limb_t pre[7];
  210.   mpn_mod_1s_4p_cps (pre, b);
  211.   return mpn_mod_1s_4p (ap, n, b << pre[1], pre);
  212. }
  213.     }
  214. }