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

数学计算

开发平台:

Unix_Linux

  1. /* sqrmod_bnm1.c -- squaring 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. /* Input is {ap,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. static void
  27. mpn_bc_sqrmod_bnm1 (mp_ptr rp, mp_srcptr ap, mp_size_t rn, mp_ptr tp)
  28. {
  29.   mp_limb_t cy;
  30.   ASSERT (0 < rn);
  31.   mpn_sqr (tp, ap, rn);
  32.   cy = mpn_add_n (rp, tp, tp + rn, rn);
  33.   /* If cy == 1, then the value of rp is at most B^rn - 2, so there can
  34.    * be no overflow when adding in the carry. */
  35.   MPN_INCR_U (rp, rn, cy);
  36. }
  37. /* Input is {ap,rn+1}; output is {rp,rn+1}, in
  38.    semi-normalised representation, computation is mod B^rn + 1. Needs
  39.    a scratch area of 2rn + 2 limbs at tp; tp == rp is allowed.
  40.    Output is normalised. */
  41. static void
  42. mpn_bc_sqrmod_bnp1 (mp_ptr rp, mp_srcptr ap, mp_size_t rn, mp_ptr tp)
  43. {
  44.   mp_limb_t cy;
  45.   ASSERT (0 < rn);
  46.   mpn_sqr (tp, ap, rn + 1);
  47.   ASSERT (tp[2*rn+1] == 0);
  48.   ASSERT (tp[2*rn] < GMP_NUMB_MAX);
  49.   cy = tp[2*rn] + mpn_sub_n (rp, tp, tp+rn, rn);
  50.   rp[rn] = 0;
  51.   MPN_INCR_U (rp, rn+1, cy );
  52. }
  53. /* Computes {rp,MIN(rn,2an)} <- {ap,an}^2 Mod(B^rn-1)
  54.  *
  55.  * The result is expected to be ZERO if and only if the operand
  56.  * already is. Otherwise the class [0] Mod(B^rn-1) is represented by
  57.  * B^rn-1.
  58.  * It should not be a problem if sqrmod_bnm1 is used to
  59.  * compute the full square with an <= 2*rn, because this condition
  60.  * implies (B^an-1)^2 < (B^rn-1) .
  61.  *
  62.  * Requires rn/4 < an <= rn
  63.  * Scratch need: rn/2 + (need for recursive call OR rn + 3). This gives
  64.  *
  65.  * S(n) <= rn/2 + MAX (rn + 4, S(n/2)) <= 3/2 rn + 4
  66.  */
  67. void
  68. mpn_sqrmod_bnm1 (mp_ptr rp, mp_size_t rn, mp_srcptr ap, mp_size_t an, mp_ptr tp)
  69. {
  70.   ASSERT (0 < an);
  71.   ASSERT (an <= rn);
  72.   if ((rn & 1) != 0 || BELOW_THRESHOLD (rn, SQRMOD_BNM1_THRESHOLD))
  73.     {
  74.       if (UNLIKELY (an < rn))
  75. {
  76.   if (UNLIKELY (2*an <= rn))
  77.     {
  78.       mpn_sqr (rp, ap, an);
  79.     }
  80.   else
  81.     {
  82.       mp_limb_t cy;
  83.       mpn_sqr (tp, ap, an);
  84.       cy = mpn_add (rp, tp, rn, tp + rn, 2*an - rn);
  85.       MPN_INCR_U (rp, rn, cy);
  86.     }
  87. }
  88.       else
  89. mpn_bc_sqrmod_bnm1 (rp, ap, rn, tp);
  90.     }
  91.   else
  92.     {
  93.       mp_size_t n;
  94.       mp_limb_t cy;
  95.       mp_limb_t hi;
  96.       n = rn >> 1;
  97.       ASSERT (2*an > n);
  98.       /* Compute xm = a^2 mod (B^n - 1), xp = a^2 mod (B^n + 1)
  99.  and crt together as
  100.  x = -xp * B^n + (B^n + 1) * [ (xp + xm)/2 mod (B^n-1)]
  101.       */
  102. #define a0 ap
  103. #define a1 (ap + n)
  104. #define xp  tp /* 2n + 2 */
  105.       /* am1  maybe in {xp, n} */
  106. #define sp1 (tp + 2*n + 2)
  107.       /* ap1  maybe in {sp1, n + 1} */
  108.       {
  109. mp_srcptr am1;
  110. mp_size_t anm;
  111. mp_ptr so;
  112. if (LIKELY (an > n))
  113.   {
  114.     so = xp + n;
  115.     am1 = xp;
  116.     cy = mpn_add (xp, a0, n, a1, an - n);
  117.     MPN_INCR_U (xp, n, cy);
  118.     anm = n;
  119.   }
  120. else
  121.   {
  122.     so = xp;
  123.     am1 = a0;
  124.     anm = an;
  125.   }
  126. mpn_sqrmod_bnm1 (rp, n, am1, anm, so);
  127.       }
  128.       {
  129. int       k;
  130. mp_srcptr ap1;
  131. mp_size_t anp;
  132. if (LIKELY (an > n)) {
  133.   ap1 = sp1;
  134.   cy = mpn_sub (sp1, a0, n, a1, an - n);
  135.   sp1[n] = 0;
  136.   MPN_INCR_U (sp1, n + 1, cy);
  137.   anp = n + ap1[n];
  138. } else {
  139.   ap1 = a0;
  140.   anp = an;
  141. }
  142. if (BELOW_THRESHOLD (n, MUL_FFT_MODF_THRESHOLD))
  143.   k=0;
  144. else
  145.   {
  146.     int mask;
  147.     k = mpn_fft_best_k (n, 1);
  148.     mask = (1<<k) -1;
  149.     while (n & mask) {k--; mask >>=1;};
  150.   }
  151. if (k >= FFT_FIRST_K)
  152.   xp[n] = mpn_mul_fft (xp, n, ap1, anp, ap1, anp, k);
  153. else if (UNLIKELY (ap1 == a0))
  154.   {
  155.     ASSERT (anp <= n);
  156.     ASSERT (2*anp > n);
  157.     mpn_sqr (xp, a0, an);
  158.     anp = 2*an - n;
  159.     cy = mpn_sub (xp, xp, n, xp + n, anp);
  160.     xp[n] = 0;
  161.     MPN_INCR_U (xp, n+1, cy);
  162.   }
  163. else
  164.   mpn_bc_sqrmod_bnp1 (xp, ap1, n, xp);
  165.       }
  166.       /* Here the CRT recomposition begins.
  167.  xm <- (xp + xm)/2 = (xp + xm)B^n/2 mod (B^n-1)
  168.  Division by 2 is a bitwise rotation.
  169.  Assumes xp normalised mod (B^n+1).
  170.  The residue class [0] is represented by [B^n-1]; except when
  171.  both input are ZERO.
  172.       */
  173. #if HAVE_NATIVE_mpn_rsh1add_n || HAVE_NATIVE_mpn_rsh1add_nc
  174. #if HAVE_NATIVE_mpn_rsh1add_nc
  175.       cy = mpn_rsh1add_nc(rp, rp, xp, n, xp[n]); /* B^n = 1 */
  176.       hi = cy << (GMP_NUMB_BITS - 1);
  177.       cy = 0;
  178.       /* next update of rp[n-1] will set cy = 1 only if rp[n-1]+=hi
  179.  overflows, i.e. a further increment will not overflow again. */
  180. #else /* ! _nc */
  181.       cy = xp[n] + mpn_rsh1add_n(rp, rp, xp, n); /* B^n = 1 */
  182.       hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
  183.       cy >>= 1;
  184.       /* cy = 1 only if xp[n] = 1 i.e. {xp,n} = ZERO, this implies that
  185.  the rsh1add was a simple rshift: the top bit is 0. cy=1 => hi=0. */
  186. #endif
  187. #if GMP_NAIL_BITS == 0
  188.       add_ssaaaa(cy, rp[n-1], cy, rp[n-1], 0, hi);
  189. #else
  190.       cy += (hi & rp[n-1]) >> (GMP_NUMB_BITS-1);
  191.       rp[n-1] ^= hi;
  192. #endif
  193. #else /* ! HAVE_NATIVE_mpn_rsh1add_n */
  194. #if HAVE_NATIVE_mpn_add_nc
  195.       cy = mpn_add_nc(rp, rp, xp, n, xp[n]);
  196. #else /* ! _nc */
  197.       cy = xp[n] + mpn_add_n(rp, rp, xp, n); /* xp[n] == 1 implies {xp,n} == ZERO */
  198. #endif
  199.       cy += (rp[0]&1);
  200.       mpn_rshift(rp, rp, n, 1);
  201.       ASSERT (cy <= 2);
  202.       hi = (cy<<(GMP_NUMB_BITS-1))&GMP_NUMB_MASK; /* (cy&1) << ... */
  203.       cy >>= 1;
  204.       /* We can have cy != 0 only if hi = 0... */
  205.       ASSERT ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0);
  206.       rp[n-1] |= hi;
  207.       /* ... rp[n-1] + cy can not overflow, the following INCR is correct. */
  208. #endif
  209.       ASSERT (cy <= 1);
  210.       /* Next increment can not overflow, read the previous comments about cy. */
  211.       ASSERT ((cy == 0) || ((rp[n-1] & GMP_NUMB_HIGHBIT) == 0));
  212.       MPN_INCR_U(rp, n, cy);
  213.       /* Compute the highest half:
  214.  ([(xp + xm)/2 mod (B^n-1)] - xp ) * B^n
  215.        */
  216.       if (UNLIKELY (2*an < rn))
  217. {
  218.   /* Note that in this case, the only way the result can equal
  219.      zero mod B^{rn} - 1 is if the input is zero, and
  220.      then the output of both the recursive calls and this CRT
  221.      reconstruction is zero, not B^{rn} - 1. */
  222.   cy = mpn_sub_n (rp + n, rp, xp, 2*an - n);
  223.   /* FIXME: This subtraction of the high parts is not really
  224.      necessary, we do it to get the carry out, and for sanity
  225.      checking. */
  226.   cy = xp[n] + mpn_sub_nc (xp + 2*an - n, rp + 2*an - n,
  227.    xp + 2*an - n, rn - 2*an, cy);
  228.   ASSERT (mpn_zero_p (xp + 2*an - n+1, rn - 1 - 2*an));
  229.   cy = mpn_sub_1 (rp, rp, 2*an, cy);
  230.   ASSERT (cy == (xp + 2*an - n)[0]);
  231. }
  232.       else
  233. {
  234.   cy = xp[n] + mpn_sub_n (rp + n, rp, xp, n);
  235.   /* cy = 1 only if {xp,n+1} is not ZERO, i.e. {rp,n} is not ZERO.
  236.      DECR will affect _at most_ the lowest n limbs. */
  237.   MPN_DECR_U (rp, 2*n, cy);
  238. }
  239. #undef a0
  240. #undef a1
  241. #undef xp
  242. #undef sp1
  243.     }
  244. }
  245. mp_size_t
  246. mpn_sqrmod_bnm1_next_size (mp_size_t n)
  247. {
  248.   mp_size_t nh;
  249.   if (BELOW_THRESHOLD (n,     SQRMOD_BNM1_THRESHOLD))
  250.     return n;
  251.   if (BELOW_THRESHOLD (n, 4 * (SQRMOD_BNM1_THRESHOLD - 1) + 1))
  252.     return (n + (2-1)) & (-2);
  253.   if (BELOW_THRESHOLD (n, 8 * (SQRMOD_BNM1_THRESHOLD - 1) + 1))
  254.     return (n + (4-1)) & (-4);
  255.   nh = (n + 1) >> 1;
  256.   if (BELOW_THRESHOLD (nh, SQR_FFT_MODF_THRESHOLD))
  257.     return (n + (8-1)) & (-8);
  258.   return 2 * mpn_fft_next_size (nh, mpn_fft_best_k (nh, 1));
  259. }