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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom_interpolate_7pts -- Interpolate for toom44, 53, 62.
  2.    Contributed to the GNU project by Niels M鰈ler.
  3.    Improvements by Marco Bodrato.
  4.    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
  5.    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  6.    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  7. Copyright 2006, 2007, 2009 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. #define BINVERT_3 MODLIMB_INVERSE_3
  22. #define BINVERT_9 
  23.   ((((GMP_NUMB_MAX / 9) << (6 - GMP_NUMB_BITS % 6)) * 8 & GMP_NUMB_MAX) | 0x39)
  24. #define BINVERT_15 
  25.   ((((GMP_NUMB_MAX >> (GMP_NUMB_BITS % 4)) / 15) * 14 * 16 & GMP_NUMB_MAX) + 15))
  26. /* For the various mpn_divexact_byN here, fall back to using either
  27.    mpn_pi1_bdiv_q_1 or mpn_divexact_1.  The former has less overhead and is
  28.    many faster if it is native.  For now, since mpn_divexact_1 is native on
  29.    several platforms where mpn_pi1_bdiv_q_1 does not yet exist, do not use
  30.    mpn_pi1_bdiv_q_1 unconditionally.  FIXME.  */
  31. /* For odd divisors, mpn_divexact_1 works fine with two's complement. */
  32. #ifndef mpn_divexact_by3
  33. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
  34. #define mpn_divexact_by3(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,3,BINVERT_3,0)
  35. #else
  36. #define mpn_divexact_by3(dst,src,size) mpn_divexact_1(dst,src,size,3)
  37. #endif
  38. #endif
  39. #ifndef mpn_divexact_by9
  40. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
  41. #define mpn_divexact_by9(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,9,BINVERT_9,0)
  42. #else
  43. #define mpn_divexact_by9(dst,src,size) mpn_divexact_1(dst,src,size,9)
  44. #endif
  45. #endif
  46. #ifndef mpn_divexact_by15
  47. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
  48. #define mpn_divexact_by15(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,15,BINVERT_15,0)
  49. #else
  50. #define mpn_divexact_by15(dst,src,size) mpn_divexact_1(dst,src,size,15)
  51. #endif
  52. #endif
  53. /* Interpolation for toom4, using the evaluation points 0, infinity,
  54.    1, -1, 2, -2, 1/2. More precisely, we want to compute
  55.    f(2^(GMP_NUMB_BITS * n)) for a polynomial f of degree 6, given the
  56.    seven values
  57.      w0 = f(0),
  58.      w1 = f(-2),
  59.      w2 = f(1),
  60.      w3 = f(-1),
  61.      w4 = f(2)
  62.      w5 = 64 * f(1/2)
  63.      w6 = limit at infinity of f(x) / x^6,
  64.    The result is 6*n + w6n limbs. At entry, w0 is stored at {rp, 2n },
  65.    w2 is stored at { rp + 2n, 2n+1 }, and w6 is stored at { rp + 6n,
  66.    w6n }. The other values are 2n + 1 limbs each (with most
  67.    significant limbs small). f(-1) and f(-1/2) may be negative, signs
  68.    determined by the flag bits. Inputs are destroyed.
  69.    Needs (2*n + 1) limbs of temporary storage.
  70. */
  71. void
  72. mpn_toom_interpolate_7pts (mp_ptr rp, mp_size_t n, enum toom7_flags flags,
  73.    mp_ptr w1, mp_ptr w3, mp_ptr w4, mp_ptr w5,
  74.    mp_size_t w6n, mp_ptr tp)
  75. {
  76.   mp_size_t m;
  77.   mp_limb_t cy;
  78.   m = 2*n + 1;
  79. #define w0 rp
  80. #define w2 (rp + 2*n)
  81. #define w6 (rp + 6*n)
  82.   ASSERT (w6n > 0);
  83.   ASSERT (w6n <= 2*n);
  84.   /* Using formulas similar to Marco Bodrato's
  85.      W5 = W5 + W4
  86.      W1 =(W4 - W1)/2
  87.      W4 = W4 - W0
  88.      W4 =(W4 - W1)/4 - W6*16
  89.      W3 =(W2 - W3)/2
  90.      W2 = W2 - W3
  91.      W5 = W5 - W2*65      May be negative.
  92.      W2 = W2 - W6 - W0
  93.      W5 =(W5 + W2*45)/2   Now >= 0 again.
  94.      W4 =(W4 - W2)/3
  95.      W2 = W2 - W4
  96.      W1 = W5 - W1         May be negative.
  97.      W5 =(W5 - W3*8)/9
  98.      W3 = W3 - W5
  99.      W1 =(W1/15 + W5)/2   Now >= 0 again.
  100.      W5 = W5 - W1
  101.      where W0 = f(0), W1 = f(-2), W2 = f(1), W3 = f(-1),
  102.    W4 = f(2), W5 = f(1/2), W6 = f(oo),
  103.      Note that most intermediate results are positive; the ones that
  104.      may be negative are represented in two's complement. We must
  105.      never shift right a value that may be negative, since that would
  106.      invalidate the sign bit. On the other hand, divexact by odd
  107.      numbers work fine with two's complement.
  108.   */
  109.   mpn_add_n (w5, w5, w4, m);
  110.   if (flags & toom7_w1_neg)
  111.     {
  112. #ifdef HAVE_NATIVE_mpn_rsh1add_n
  113.       mpn_rsh1add_n (w1, w1, w4, m);
  114. #else
  115.       mpn_add_n (w1, w1, w4, m);  ASSERT (!(w1[0] & 1));
  116.       mpn_rshift (w1, w1, m, 1);
  117. #endif
  118.     }
  119.   else
  120.     {
  121. #ifdef HAVE_NATIVE_mpn_rsh1sub_n
  122.       mpn_rsh1sub_n (w1, w4, w1, m);
  123. #else
  124.       mpn_sub_n (w1, w4, w1, m);  ASSERT (!(w1[0] & 1));
  125.       mpn_rshift (w1, w1, m, 1);
  126. #endif
  127.     }
  128.   mpn_sub (w4, w4, m, w0, 2*n);
  129.   mpn_sub_n (w4, w4, w1, m);  ASSERT (!(w4[0] & 3));
  130.   mpn_rshift (w4, w4, m, 2); /* w4>=0 */
  131.   tp[w6n] = mpn_lshift (tp, w6, w6n, 4);
  132.   mpn_sub (w4, w4, m, tp, w6n+1);
  133.   if (flags & toom7_w3_neg)
  134.     {
  135. #ifdef HAVE_NATIVE_mpn_rsh1add_n
  136.       mpn_rsh1add_n (w3, w3, w2, m);
  137. #else
  138.       mpn_add_n (w3, w3, w2, m);  ASSERT (!(w3[0] & 1));
  139.       mpn_rshift (w3, w3, m, 1);
  140. #endif
  141.     }
  142.   else
  143.     {
  144. #ifdef HAVE_NATIVE_mpn_rsh1sub_n
  145.       mpn_rsh1sub_n (w3, w2, w3, m);
  146. #else
  147.       mpn_sub_n (w3, w2, w3, m);  ASSERT (!(w3[0] & 1));
  148.       mpn_rshift (w3, w3, m, 1);
  149. #endif
  150.     }
  151.   mpn_sub_n (w2, w2, w3, m);
  152.   mpn_submul_1 (w5, w2, m, 65);
  153.   mpn_sub (w2, w2, m, w6, w6n);
  154.   mpn_sub (w2, w2, m, w0, 2*n);
  155.   mpn_addmul_1 (w5, w2, m, 45);  ASSERT (!(w5[0] & 1));
  156.   mpn_rshift (w5, w5, m, 1);
  157.   mpn_sub_n (w4, w4, w2, m);
  158.   mpn_divexact_by3 (w4, w4, m);
  159.   mpn_sub_n (w2, w2, w4, m);
  160.   mpn_sub_n (w1, w5, w1, m);
  161.   mpn_lshift (tp, w3, m, 3);
  162.   mpn_sub_n (w5, w5, tp, m);
  163.   mpn_divexact_by9 (w5, w5, m);
  164.   mpn_sub_n (w3, w3, w5, m);
  165.   mpn_divexact_by15 (w1, w1, m);
  166.   mpn_add_n (w1, w1, w5, m);  ASSERT (!(w1[0] & 1));
  167.   mpn_rshift (w1, w1, m, 1); /* w1>=0 now */
  168.   mpn_sub_n (w5, w5, w1, m);
  169.   /* These bounds are valid for the 4x4 polynomial product of toom44,
  170.    * and they are conservative for toom53 and toom62. */
  171.   ASSERT (w1[2*n] < 2);
  172.   ASSERT (w2[2*n] < 3);
  173.   ASSERT (w3[2*n] < 4);
  174.   ASSERT (w4[2*n] < 3);
  175.   ASSERT (w5[2*n] < 2);
  176.   /* Addition chain. Note carries and the 2n'th limbs that need to be
  177.    * added in.
  178.    *
  179.    * Special care is needed for w2[2n] and the corresponding carry,
  180.    * since the "simple" way of adding it all together would overwrite
  181.    * the limb at wp[2*n] and rp[4*n] (same location) with the sum of
  182.    * the high half of w3 and the low half of w4.
  183.    *
  184.    *         7    6    5    4    3    2    1    0
  185.    *    |    |    |    |    |    |    |    |    |
  186.    *                  ||w3 (2n+1)|
  187.    *             ||w4 (2n+1)|
  188.    *        ||w5 (2n+1)|        ||w1 (2n+1)|
  189.    *  + | w6 (w6n)|        ||w2 (2n+1)| w0 (2n) |  (share storage with r)
  190.    *  -----------------------------------------------
  191.    *  r |    |    |    |    |    |    |    |    |
  192.    *        c7   c6   c5   c4   c3                 Carries to propagate
  193.    */
  194.   cy = mpn_add_n (rp + n, rp + n, w1, m);
  195.   MPN_INCR_U (w2 + n + 1, n , cy);
  196.   cy = mpn_add_n (rp + 3*n, rp + 3*n, w3, n);
  197.   MPN_INCR_U (w3 + n, n + 1, w2[2*n] + cy);
  198.   cy = mpn_add_n (rp + 4*n, w3 + n, w4, n);
  199.   MPN_INCR_U (w4 + n, n + 1, w3[2*n] + cy);
  200.   cy = mpn_add_n (rp + 5*n, w4 + n, w5, n);
  201.   MPN_INCR_U (w5 + n, n + 1, w4[2*n] + cy);
  202.   if (w6n > n + 1)
  203.     ASSERT_NOCARRY (mpn_add (rp + 6*n, rp + 6*n, w6n, w5 + n, n + 1));
  204.   else
  205.     {
  206.       ASSERT_NOCARRY (mpn_add_n (rp + 6*n, rp + 6*n, w5 + n, w6n));
  207. #if WANT_ASSERT
  208.       {
  209. mp_size_t i;
  210. for (i = w6n; i <= n; i++)
  211.   ASSERT (w5[n + i] == 0);
  212.       }
  213. #endif
  214.     }
  215. }