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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom_interpolate_6pts -- Interpolate for toom43, 52
  2.    Contributed to the GNU project by Marco Bodrato.
  3.    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
  4.    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  5.    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  6. Copyright 2009, 2010 Free 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. /* For odd divisors, mpn_divexact_1 works fine with two's complement. */
  21. #ifndef mpn_divexact_by3
  22. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && MODLIMB_INVERSE_3
  23. #define mpn_divexact_by3(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,3,MODLIMB_INVERSE_3,0)
  24. #else
  25. #define mpn_divexact_by3(dst,src,size) mpn_divexact_1(dst,src,size,3)
  26. #endif
  27. #endif
  28. /* Interpolation for Toom-3.5, using the evaluation points: infinity,
  29.    1, -1, 2, -2. More precisely, we want to compute
  30.    f(2^(GMP_NUMB_BITS * n)) for a polynomial f of degree 5, given the
  31.    six values
  32.      w5 = f(0),
  33.      w4 = f(-1),
  34.      w3 = f(1)
  35.      w2 = f(-2),
  36.      w1 = f(2),
  37.      w0 = limit at infinity of f(x) / x^5,
  38.    The result is stored in {pp, 5*n + w0n}. At entry, w5 is stored at
  39.    {pp, 2n}, w3 is stored at {pp + 2n, 2n+1}, and w0 is stored at
  40.    {pp + 5n, w0n}. The other values are 2n + 1 limbs each (with most
  41.    significant limbs small). f(-1) and f(-2) may be negative, signs
  42.    determined by the flag bits. All intermediate results are positive.
  43.    Inputs are destroyed.
  44.    Interpolation sequence was taken from the paper: "Integer and
  45.    Polynomial Multiplication: Towards Optimal Toom-Cook Matrices".
  46.    Some slight variations were introduced: adaptation to "gmp
  47.    instruction set", and a final saving of an operation by interlacing
  48.    interpolation and recomposition phases.
  49. */
  50. void
  51. mpn_toom_interpolate_6pts (mp_ptr pp, mp_size_t n, enum toom6_flags flags,
  52.    mp_ptr w4, mp_ptr w2, mp_ptr w1,
  53.    mp_size_t w0n)
  54. {
  55.   mp_limb_t cy;
  56.   /* cy6 can be stored in w1[2*n], cy4 in w4[0], embankment in w2[0] */
  57.   mp_limb_t cy4, cy6, embankment;
  58.   ASSERT( n > 0 );
  59.   ASSERT( 2*n >= w0n && w0n > 0 );
  60. #define w5  pp /* 2n   */
  61. #define w3  (pp + 2 * n) /* 2n+1 */
  62. #define w0  (pp + 5 * n) /* w0n  */
  63.   /* Interpolate with sequence:
  64.      W2 =(W1 - W2)>>2
  65.      W1 =(W1 - W5)>>1
  66.      W1 =(W1 - W2)>>1
  67.      W4 =(W3 - W4)>>1
  68.      W2 =(W2 - W4)/3
  69.      W3 = W3 - W4 - W5
  70.      W1 =(W1 - W3)/3
  71.      // Last steps are mixed with recomposition...
  72.      W2 = W2 - W0<<2
  73.      W4 = W4 - W2
  74.      W3 = W3 - W1
  75.      W2 = W2 - W0
  76.   */
  77.   /* W2 =(W1 - W2)>>2 */
  78.   if (flags & toom6_vm2_neg)
  79.     mpn_add_n (w2, w1, w2, 2 * n + 1);
  80.   else
  81.     mpn_sub_n (w2, w1, w2, 2 * n + 1);
  82.   mpn_rshift (w2, w2, 2 * n + 1, 2);
  83.   /* W1 =(W1 - W5)>>1 */
  84.   w1[2*n] -= mpn_sub_n (w1, w1, w5, 2*n);
  85.   mpn_rshift (w1, w1, 2 * n + 1, 1);
  86.   /* W1 =(W1 - W2)>>1 */
  87. #if HAVE_NATIVE_mpn_rsh1sub_n
  88.   mpn_rsh1sub_n (w1, w1, w2, 2 * n + 1);
  89. #else
  90.   mpn_sub_n (w1, w1, w2, 2 * n + 1);
  91.   mpn_rshift (w1, w1, 2 * n + 1, 1);
  92. #endif
  93.   /* W4 =(W3 - W4)>>1 */
  94.   if (flags & toom6_vm1_neg)
  95.     {
  96. #if HAVE_NATIVE_mpn_rsh1add_n
  97.       mpn_rsh1add_n (w4, w3, w4, 2 * n + 1);
  98. #else
  99.       mpn_add_n (w4, w3, w4, 2 * n + 1);
  100.       mpn_rshift (w4, w4, 2 * n + 1, 1);
  101. #endif
  102.     }
  103.   else
  104.     {
  105. #if HAVE_NATIVE_mpn_rsh1sub_n
  106.       mpn_rsh1sub_n (w4, w3, w4, 2 * n + 1);
  107. #else
  108.       mpn_sub_n (w4, w3, w4, 2 * n + 1);
  109.       mpn_rshift (w4, w4, 2 * n + 1, 1);
  110. #endif
  111.     }
  112.   /* W2 =(W2 - W4)/3 */
  113.   mpn_sub_n (w2, w2, w4, 2 * n + 1);
  114.   mpn_divexact_by3 (w2, w2, 2 * n + 1);
  115.   /* W3 = W3 - W4 - W5 */
  116.   mpn_sub_n (w3, w3, w4, 2 * n + 1);
  117.   w3[2 * n] -= mpn_sub_n (w3, w3, w5, 2 * n);
  118.   /* W1 =(W1 - W3)/3 */
  119.   mpn_sub_n (w1, w1, w3, 2 * n + 1);
  120.   mpn_divexact_by3 (w1, w1, 2 * n + 1);
  121.   /*
  122.     [1 0 0 0 0 0;
  123.      0 1 0 0 0 0;
  124.      1 0 1 0 0 0;
  125.      0 1 0 1 0 0;
  126.      1 0 1 0 1 0;
  127.      0 0 0 0 0 1]
  128.     pp[] prior to operations:
  129.      |_H w0__|_L w0__|______||_H w3__|_L w3__|_H w5__|_L w5__|
  130.     summation scheme for remaining operations:
  131.      |______________5|n_____4|n_____3|n_____2|n______|n______|pp
  132.      |_H w0__|_L w0__|______||_H w3__|_L w3__|_H w5__|_L w5__|
  133.     || H w4  | L w4  |
  134.     || H w2  | L w2  |
  135.     || H w1  | L w1  |
  136.     ||-H w1  |-L w1  |
  137.      |-H w0  |-L w0 ||-H w2  |-L w2  |
  138.   */
  139.   cy = mpn_add_n (pp + n, pp + n, w4, 2 * n + 1);
  140.   MPN_INCR_U (pp + 3 * n + 1, n, cy);
  141.   /* W2 -= W0<<2 */
  142. #if HAVE_NATIVE_mpn_sublsh_n || HAVE_NATIVE_mpn_sublsh2_n
  143. #if HAVE_NATIVE_mpn_sublsh2_n
  144.   cy = mpn_sublsh2_n(w2, w0, w0n);
  145. #else
  146.   cy = mpn_sublsh_n(w2, w0, w0n, 2);
  147. #endif
  148. #else
  149.   /* {W4,2*n+1} is now free and can be overwritten. */
  150.   cy = mpn_lshift(w4, w0, w0n, 2);
  151.   cy+= mpn_sub_n(w2, w2, w4, w0n);
  152. #endif
  153.   MPN_DECR_U (w2 + w0n, 2 * n + 1 - w0n, cy);
  154.   /* W4L = W4L - W2L */
  155.   cy = mpn_sub_n (pp + n, pp + n, w2, n);
  156.   MPN_DECR_U (w3, 2 * n + 1, cy);
  157.   /* W3H = W3H + W2L */
  158.   cy4 = w3[2 * n] + mpn_add_n (pp + 3 * n, pp + 3 * n, w2, n);
  159.   /* W1L + W2H */
  160.   cy = w2[2 * n] + mpn_add_n (pp + 4 * n, w1, w2 + n, n);
  161.   MPN_INCR_U (w1 + n, n + 1, cy);
  162.   /* W0 = W0 + W1H */
  163.   if (LIKELY (w0n > n))
  164.     cy6 = w1[2 * n] + mpn_add_n (w0, w0, w1 + n, n);
  165.   else
  166.     cy6 = mpn_add_n (w0, w0, w1 + n, w0n);
  167.   /*
  168.     summation scheme for the next operation:
  169.      |...____5|n_____4|n_____3|n_____2|n______|n______|pp
  170.      |...w0___|_w1_w2_|_H w3__|_L w3__|_H w5__|_L w5__|
  171.      ...-w0___|-w1_w2 |
  172.   */
  173.   /* if(LIKELY(w0n>n)) the two operands below DO overlap! */
  174.   cy = mpn_sub_n (pp + 2 * n, pp + 2 * n, pp + 4 * n, n + w0n);
  175.   /* embankment is a "dirty trick" to avoid carry/borrow propagation
  176.      beyond allocated memory */
  177.   embankment = w0[w0n - 1] - 1;
  178.   w0[w0n - 1] = 1;
  179.   if (LIKELY (w0n > n)) {
  180.     if ( LIKELY(cy4 > cy6) )
  181.       MPN_INCR_U (pp + 4 * n, w0n + n, cy4 - cy6);
  182.     else
  183.       MPN_DECR_U (pp + 4 * n, w0n + n, cy6 - cy4);
  184.     MPN_DECR_U (pp + 3 * n + w0n, 2 * n, cy);
  185.     MPN_INCR_U (w0 + n, w0n - n, cy6);
  186.   } else {
  187.     MPN_INCR_U (pp + 4 * n, w0n + n, cy4);
  188.     MPN_DECR_U (pp + 3 * n + w0n, 2 * n, cy + cy6);
  189.   }
  190.   w0[w0n - 1] += embankment;
  191. #undef w5
  192. #undef w3
  193. #undef w0
  194. }