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

数学计算

开发平台:

Unix_Linux

  1. /* mulmod_bnm1.c -- multiplication mod B^n-1.
  2.    Contributed to the GNU project by Niels M鰈ler, Torbjorn Granlund and
  3.    Marco Bodrato.
  4.    THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES.  IT IS ONLY
  5.    SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  6.    GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  7. Copyright 2009, 2010 Free Software Foundation, Inc.
  8. This file is part of the GNU MP Library.
  9. The GNU MP Library is free software; you can redistribute it and/or modify
  10. it under the terms of the GNU Lesser General Public License as published by
  11. the Free Software Foundation; either version 3 of the License, or (at your
  12. option) any later version.
  13. The GNU MP Library is distributed in the hope that it will be useful, but
  14. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16. License for more details.
  17. You should have received a copy of the GNU Lesser General Public License
  18. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19. #include "gmp.h"
  20. #include "gmp-impl.h"
  21. #include "longlong.h"
  22. /* Inputs are {ap,rn} and {bp,rn}; output is {rp,rn}, computation is
  23.    mod B^rn - 1, and values are semi-normalised; zero is represented
  24.    as either 0 or B^n - 1.  Needs a scratch of 2rn limbs at tp.
  25.    tp==rp is allowed. */
  26. void
  27. mpn_bc_mulmod_bnm1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn,
  28.     mp_ptr tp)
  29. {
  30.   mp_limb_t cy;
  31.   ASSERT (0 < rn);
  32.   mpn_mul_n (tp, ap, bp, rn);
  33.   cy = mpn_add_n (rp, tp, tp + rn, rn);
  34.   /* If cy == 1, then the value of rp is at most B^rn - 2, so there can
  35.    * be no overflow when adding in the carry. */
  36.   MPN_INCR_U (rp, rn, cy);
  37. }
  38. /* Inputs are {ap,rn+1} and {bp,rn+1}; output is {rp,rn+1}, in
  39.    semi-normalised representation, computation is mod B^rn + 1. Needs
  40.    a scratch area of 2rn + 2 limbs at tp; tp == rp is allowed.
  41.    Output is normalised. */
  42. static void
  43. mpn_bc_mulmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t rn,
  44.     mp_ptr tp)
  45. {
  46.   mp_limb_t cy;
  47.   ASSERT (0 < rn);
  48.   mpn_mul_n (tp, ap, bp, rn + 1);
  49.   ASSERT (tp[2*rn+1] == 0);
  50.   ASSERT (tp[2*rn] < GMP_NUMB_MAX);
  51.   cy = tp[2*rn] + mpn_sub_n (rp, tp, tp+rn, rn);
  52.   rp[rn] = 0;
  53.   MPN_INCR_U (rp, rn+1, cy );
  54. }
  55. /* Computes {rp,MIN(rn,an+bn)} <- {ap,an}*{bp,bn} Mod(B^rn-1)
  56.  *
  57.  * The result is expected to be ZERO if and only if one of the operand
  58.  * already is. Otherwise the class [0] Mod(B^rn-1) is represented by
  59.  * B^rn-1. This should not be a problem if mulmod_bnm1 is used to
  60.  * combine results and obtain a natural number when one knows in
  61.  * advance that the final value is less than (B^rn-1).
  62.  * Moreover it should not be a problem if mulmod_bnm1 is used to
  63.  * compute the full product with an+bn <= rn, because this condition
  64.  * implies (B^an-1)(B^bn-1) < (B^rn-1) .
  65.  *
  66.  * Requires 0 < bn <= an <= rn and an + bn > rn/2
  67.  * Scratch need: rn + (need for recursive call OR rn + 4). This gives
  68.  *
  69.  * S(n) <= rn + MAX (rn + 4, S(n/2)) <= 2rn + 4
  70.  */
  71. void
  72. mpn_mulmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an, mp_srcptr bp, mp_size_t bn, mp_ptr tp)
  73. {
  74.   ASSERT (0 < bn);
  75.   ASSERT (bn <= an);
  76.   ASSERT (an <= rn);
  77.   if ((rn & 1) != 0 || BELOW_THRESHOLD (rn, MULMOD_BNM1_THRESHOLD))
  78.     {
  79.       if (UNLIKELY (bn < rn))
  80. {
  81.   if (UNLIKELY (an + bn <= rn))
  82.     {
  83.       mpn_mul (rp, ap, an, bp, bn);
  84.     }
  85.   else
  86.     {
  87.       mp_limb_t cy;
  88.       mpn_mul (tp, ap, an, bp, bn);
  89.       cy = mpn_add (rp, tp, rn, tp + rn, an + bn - rn);
  90.       MPN_INCR_U (rp, rn, cy);
  91.     }
  92. }
  93.       else
  94. mpn_bc_mulmod_bnm1 (rp, ap, bp, rn, tp);
  95.     }
  96.   else
  97.     {
  98.       mp_size_t n;
  99.       mp_limb_t cy;
  100.       mp_limb_t hi;
  101.       n = rn >> 1;
  102.       /* We need at least an + bn >= n, to be able to fit one of the
  103.  recursive products at rp. Requiring strict inequality makes
  104.  the coded slightly simpler. If desired, we could avoid this
  105.  restriction by initially halving rn as long as rn is even and
  106.  an + bn <= rn/2. */
  107.       ASSERT (an + bn > n);
  108.       /* Compute xm = a*b mod (B^n - 1), xp = a*b mod (B^n + 1)
  109.  and crt together as
  110.  x = -xp * B^n + (B^n + 1) * [ (xp + xm)/2 mod (B^n-1)]
  111.       */
  112. #define a0 ap
  113. #define a1 (ap + n)
  114. #define b0 bp
  115. #define b1 (bp + n)
  116. #define xp  tp /* 2n + 2 */
  117.       /* am1  maybe in {xp, n} */
  118.       /* bm1  maybe in {xp + n, n} */
  119. #define sp1 (tp + 2*n + 2)
  120.       /* ap1  maybe in {sp1, n + 1} */
  121.       /* bp1  maybe in {sp1 + n + 1, n + 1} */
  122.       {
  123. mp_srcptr am1, bm1;
  124. mp_size_t anm, bnm;
  125. mp_ptr so;
  126. if (LIKELY (an > n))
  127.   {
  128.     am1 = xp;
  129.     cy = mpn_add (xp, a0, n, a1, an - n);
  130.     MPN_INCR_U (xp, n, cy);
  131.     anm = n;
  132.     if (LIKELY (bn > n))
  133.       {
  134. bm1 = xp + n;
  135. cy = mpn_add (xp + n, b0, n, b1, bn - n);
  136. MPN_INCR_U (xp + n, n, cy);
  137. bnm = n;
  138. so = xp + 2*n;
  139.       }
  140.     else
  141.       {
  142. so = xp + n;
  143. bm1 = b0;
  144. bnm = bn;
  145.       }
  146.   }
  147. else
  148.   {
  149.     so = xp;
  150.     am1 = a0;
  151.     anm = an;
  152.     bm1 = b0;
  153.     bnm = bn;
  154.   }
  155. mpn_mulmod_bnm1 (rp, n, am1, anm, bm1, bnm, so);
  156.       }
  157.       {
  158. int       k;
  159. mp_srcptr ap1, bp1;
  160. mp_size_t anp, bnp;
  161. if (LIKELY (an > n)) {
  162.   ap1 = sp1;
  163.   cy = mpn_sub (sp1, a0, n, a1, an - n);
  164.   sp1[n] = 0;
  165.   MPN_INCR_U (sp1, n + 1, cy);
  166.   anp = n + ap1[n];
  167. } else {
  168.   ap1 = a0;
  169.   anp = an;
  170. }
  171. if (LIKELY (bn > n)) {
  172.   bp1 = sp1 + n + 1;
  173.   cy = mpn_sub (sp1 + n + 1, b0, n, b1, bn - n);
  174.   sp1[2*n+1] = 0;
  175.   MPN_INCR_U (sp1 + n + 1, n + 1, cy);
  176.   bnp = n + bp1[n];
  177. } else {
  178.   bp1 = b0;
  179.   bnp = bn;
  180. }
  181. if (BELOW_THRESHOLD (n, MUL_FFT_MODF_THRESHOLD))
  182.   k=0;
  183. else
  184.   {
  185.     int mask;
  186.     k = mpn_fft_best_k (n, 0);
  187.     mask = (1<<k) -1;
  188.     while (n & mask) {k--; mask >>=1;};
  189.   }
  190. if (k >= FFT_FIRST_K)
  191.   xp[n] = mpn_mul_fft (xp, n, ap1, anp, bp1, bnp, k);
  192. else if (UNLIKELY (bp1 == b0))
  193.   {
  194.     ASSERT (anp + bnp <= 2*n+1);
  195.     ASSERT (anp + bnp > n);
  196.     ASSERT (anp >= bnp);
  197.     mpn_mul (xp, ap1, anp, bp1, bnp);
  198.     anp = anp + bnp - n;
  199.     ASSERT (anp <= n || xp[2*n]==0);
  200.     anp-= anp > n;
  201.     cy = mpn_sub (xp, xp, n, xp + n, anp);
  202.     xp[n] = 0;
  203.     MPN_INCR_U (xp, n+1, cy);
  204.   }
  205. else
  206.   mpn_bc_mulmod_bnp1 (xp, ap1, bp1, n, xp);
  207.       }
  208.       /* Here the CRT recomposition begins.
  209.  xm <- (xp + xm)/2 = (xp + xm)B^n/2 mod (B^n-1)
  210.  Division by 2 is a bitwise rotation.
  211.  Assumes xp normalised mod (B^n+1).
  212.  The residue class [0] is represented by [B^n-1]; except when
  213.  both input are ZERO.
  214.       */
  215. #if HAVE_NATIVE_mpn_rsh1add_n || HAVE_NATIVE_mpn_rsh1add_nc
  216. #if HAVE_NATIVE_mpn_rsh1add_nc
  217.       cy = mpn_rsh1add_nc(rp, rp, xp, n, xp[n]); /* B^n = 1 */
  218.       hi = cy << (GMP_NUMB_BITS - 1);
  219.       cy = 0;
  220.       /* next update of rp[n-1] will set cy = 1 only if rp[n-1]+=hi
  221.  overflows, i.e. a further increment will not overflow again. */
  222. #else /* ! _nc */
  223.       cy = xp[n] + mpn_rsh1add_n(rp, rp, xp, n); /* B^n = 1 */
  224.       hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
  225.       cy >>= 1;
  226.       /* cy = 1 only if xp[n] = 1 i.e. {xp,n} = ZERO, this implies that
  227.  the rsh1add was a simple rshift: the top bit is 0. cy=1 => hi=0. */
  228. #endif
  229. #if GMP_NAIL_BITS == 0
  230.       add_ssaaaa(cy, rp[n-1], cy, rp[n-1], 0, hi);
  231. #else
  232.       cy += (hi & rp[n-1]) >> (GMP_NUMB_BITS-1);
  233.       rp[n-1] ^= hi;
  234. #endif
  235. #else /* ! HAVE_NATIVE_mpn_rsh1add_n */
  236. #if HAVE_NATIVE_mpn_add_nc
  237.       cy = mpn_add_nc(rp, rp, xp, n, xp[n]);
  238. #else /* ! _nc */
  239.       cy = xp[n] + mpn_add_n(rp, rp, xp, n); /* xp[n] == 1 implies {xp,n} == ZERO */
  240. #endif
  241.       cy += (rp[0]&1);
  242.       mpn_rshift(rp, rp, n, 1);
  243.       ASSERT (cy <= 2);
  244.       hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
  245.       cy >>= 1;
  246.       /* We can have cy != 0 only if hi = 0... */
  247.       ASSERT ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0);
  248.       rp[n-1] |= hi;
  249.       /* ... rp[n-1] + cy can not overflow, the following INCR is correct. */
  250. #endif
  251.       ASSERT (cy <= 1);
  252.       /* Next increment can not overflow, read the previous comments about cy. */
  253.       ASSERT ((cy == 0) || ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0));
  254.       MPN_INCR_U(rp, n, cy);
  255.       /* Compute the highest half:
  256.  ([(xp + xm)/2 mod (B^n-1)] - xp ) * B^n
  257.        */
  258.       if (UNLIKELY (an + bn < rn))
  259. {
  260.   /* Note that in this case, the only way the result can equal
  261.      zero mod B^{rn} - 1 is if one of the inputs is zero, and
  262.      then the output of both the recursive calls and this CRT
  263.      reconstruction is zero, not B^{rn} - 1. Which is good,
  264.      since the latter representation doesn't fit in the output
  265.      area.*/
  266.   cy = mpn_sub_n (rp + n, rp, xp, an + bn - n);
  267.   /* FIXME: This subtraction of the high parts is not really
  268.      necessary, we do it to get the carry out, and for sanity
  269.      checking. */
  270.   cy = xp[n] + mpn_sub_nc (xp + an + bn - n, rp + an + bn - n,
  271.    xp + an + bn - n, rn - (an + bn), cy);
  272.   ASSERT (an + bn == rn - 1 ||
  273.   mpn_zero_p (xp + an + bn - n + 1, rn - 1 - (an + bn)));
  274.   cy = mpn_sub_1 (rp, rp, an + bn, cy);
  275.   ASSERT (cy == (xp + an + bn - n)[0]);
  276. }
  277.       else
  278. {
  279.   cy = xp[n] + mpn_sub_n (rp + n, rp, xp, n);
  280.   /* cy = 1 only if {xp,n+1} is not ZERO, i.e. {rp,n} is not ZERO.
  281.      DECR will affect _at most_ the lowest n limbs. */
  282.   MPN_DECR_U (rp, 2*n, cy);
  283. }
  284. #undef a0
  285. #undef a1
  286. #undef b0
  287. #undef b1
  288. #undef xp
  289. #undef sp1
  290.     }
  291. }
  292. mp_size_t
  293. mpn_mulmod_bnm1_next_size (mp_size_t n)
  294. {
  295.   mp_size_t nh;
  296.   if (BELOW_THRESHOLD (n,     MULMOD_BNM1_THRESHOLD))
  297.     return n;
  298.   if (BELOW_THRESHOLD (n, 4 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
  299.     return (n + (2-1)) & (-2);
  300.   if (BELOW_THRESHOLD (n, 8 * (MULMOD_BNM1_THRESHOLD - 1) + 1))
  301.     return (n + (4-1)) & (-4);
  302.   nh = (n + 1) >> 1;
  303.   if (BELOW_THRESHOLD (nh, MUL_FFT_MODF_THRESHOLD))
  304.     return (n + (8-1)) & (-8);
  305.   return 2 * mpn_fft_next_size (nh, mpn_fft_best_k (nh, 0));
  306. }