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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom_interpolate_5pts -- Interpolate for toom3, 33, 42.
  2.    Contributed to the GNU project by Robert Harley.
  3.    Improvements by Paul Zimmermann and 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 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2009 Free Software
  8. Foundation, Inc.
  9. This file is part of the GNU MP Library.
  10. The GNU MP Library is free software; you can redistribute it and/or modify it
  11. under the terms of the GNU Lesser General Public License as published by the
  12. Free Software Foundation; either version 3 of the License, or (at your
  13. option) any later version.
  14. The GNU MP Library is distributed in the hope that it will be useful, but
  15. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  16. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License
  17. for more details.
  18. You should have received a copy of the GNU Lesser General Public License
  19. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  20. #include "gmp.h"
  21. #include "gmp-impl.h"
  22. void
  23. mpn_toom_interpolate_5pts (mp_ptr c, mp_ptr v2, mp_ptr vm1,
  24.    mp_size_t k, mp_size_t twor, int sa,
  25.    mp_limb_t vinf0)
  26. {
  27.   mp_limb_t cy, saved;
  28.   mp_size_t twok;
  29.   mp_size_t kk1;
  30.   mp_ptr c1, v1, c3, vinf;
  31.   twok = k + k;
  32.   kk1 = twok + 1;
  33.   c1 = c  + k;
  34.   v1 = c1 + k;
  35.   c3 = v1 + k;
  36.   vinf = c3 + k;
  37. #define v0 (c)
  38.   /* (1) v2 <- v2-vm1 < v2+|vm1|,       (16 8 4 2 1) - (1 -1 1 -1  1) =
  39.      thus 0 <= v2 < 50*B^(2k) < 2^6*B^(2k)             (15 9 3  3  0)
  40.   */
  41.   if (sa)
  42.     ASSERT_NOCARRY (mpn_add_n (v2, v2, vm1, kk1));
  43.   else
  44.     ASSERT_NOCARRY (mpn_sub_n (v2, v2, vm1, kk1));
  45.   /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
  46.        v0       v1       hi(vinf)       |vm1|     v2-vm1      EMPTY */
  47.   ASSERT_NOCARRY (mpn_divexact_by3 (v2, v2, kk1));    /* v2 <- v2 / 3 */
  48.       /* (5 3 1 1 0)*/
  49.   /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
  50.        v0       v1      hi(vinf)       |vm1|     (v2-vm1)/3    EMPTY */
  51.   /* (2) vm1 <- tm1 := (v1 - vm1) / 2  [(1 1 1 1 1) - (1 -1 1 -1 1)] / 2 =
  52.      tm1 >= 0                                         (0  1 0  1 0)
  53.      No carry comes out from {v1, kk1} +/- {vm1, kk1},
  54.      and the division by two is exact.
  55.      If (sa!=0) the sign of vm1 is negative */
  56.   if (sa)
  57.     {
  58. #ifdef HAVE_NATIVE_mpn_rsh1add_n
  59.       mpn_rsh1add_n (vm1, v1, vm1, kk1);
  60. #else
  61.       ASSERT_NOCARRY (mpn_add_n (vm1, v1, vm1, kk1));
  62.       ASSERT_NOCARRY (mpn_rshift (vm1, vm1, kk1, 1));
  63. #endif
  64.     }
  65.   else
  66.     {
  67. #ifdef HAVE_NATIVE_mpn_rsh1sub_n
  68.       mpn_rsh1sub_n (vm1, v1, vm1, kk1);
  69. #else
  70.       ASSERT_NOCARRY (mpn_sub_n (vm1, v1, vm1, kk1));
  71.       ASSERT_NOCARRY (mpn_rshift (vm1, vm1, kk1, 1));
  72. #endif
  73.     }
  74.   /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
  75.        v0       v1        hi(vinf)       tm1     (v2-vm1)/3    EMPTY */
  76.   /* (3) v1 <- t1 := v1 - v0    (1 1 1 1 1) - (0 0 0 0 1) = (1 1 1 1 0)
  77.      t1 >= 0
  78.   */
  79.   vinf[0] -= mpn_sub_n (v1, v1, c, twok);
  80.   /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
  81.        v0     v1-v0        hi(vinf)       tm1     (v2-vm1)/3    EMPTY */
  82.   /* (4) v2 <- t2 := ((v2-vm1)/3-t1)/2 = (v2-vm1-3*t1)/6
  83.      t2 >= 0                  [(5 3 1 1 0) - (1 1 1 1 0)]/2 = (2 1 0 0 0)
  84.   */
  85. #ifdef HAVE_NATIVE_mpn_rsh1sub_n
  86.   mpn_rsh1sub_n (v2, v2, v1, kk1);
  87. #else
  88.   ASSERT_NOCARRY (mpn_sub_n (v2, v2, v1, kk1));
  89.   ASSERT_NOCARRY (mpn_rshift (v2, v2, kk1, 1));
  90. #endif
  91.   /* {c,2k} {c+2k,2k+1} {c+4k+1,2r-1} {t,2k+1} {t+2k+1,2k+1} {t+4k+2,2r}
  92.        v0     v1-v0        hi(vinf)     tm1    (v2-vm1-3t1)/6    EMPTY */
  93.   /* (5) v1 <- t1-tm1           (1 1 1 1 0) - (0 1 0 1 0) = (1 0 1 0 0)
  94.      result is v1 >= 0
  95.   */
  96.   ASSERT_NOCARRY (mpn_sub_n (v1, v1, vm1, kk1));
  97.   /* We do not need to read the value in vm1, so we add it in {c+k, ...} */
  98.   cy = mpn_add_n (c1, c1, vm1, kk1);
  99.   MPN_INCR_U (c3 + 1, twor + k - 1, cy); /* 2n-(3k+1) = 2r+k-1 */
  100.   /* Memory allocated for vm1 is now free, it can be recycled ...*/
  101.   /* (6) v2 <- v2 - 2*vinf,     (2 1 0 0 0) - 2*(1 0 0 0 0) = (0 1 0 0 0)
  102.      result is v2 >= 0 */
  103.   saved = vinf[0];       /* Remember v1's highest byte (will be overwritten). */
  104.   vinf[0] = vinf0;       /* Set the right value for vinf0                     */
  105. #ifdef HAVE_NATIVE_mpn_sublsh1_n
  106.   cy = mpn_sublsh1_n (v2, v2, vinf, twor);
  107. #else
  108.   /* Overwrite unused vm1 */
  109.   cy = mpn_lshift (vm1, vinf, twor, 1);
  110.   cy += mpn_sub_n (v2, v2, vm1, twor);
  111. #endif
  112.   MPN_DECR_U (v2 + twor, kk1 - twor, cy);
  113.   /* Current matrix is
  114.      [1 0 0 0 0; vinf
  115.       0 1 0 0 0; v2
  116.       1 0 1 0 0; v1
  117.       0 1 0 1 0; vm1
  118.       0 0 0 0 1] v0
  119.      Some vaues already are in-place (we added vm1 in the correct position)
  120.      | vinf|  v1 |  v0 |
  121.       | vm1 |
  122.      One still is in a separated area
  123. | +v2 |
  124.      We have to compute v1-=vinf; vm1 -= v2,
  125.    |-vinf|
  126.       | -v2 |
  127.      Carefully reordering operations we can avoid to compute twice the sum
  128.      of the high half of v2 plus the low half of vinf.
  129.   */
  130.   /* Add the high half of t2 in {vinf} */
  131.   if ( LIKELY(twor > k + 1) ) { /* This is the expected flow  */
  132.     cy = mpn_add_n (vinf, vinf, v2 + k, k + 1);
  133.     MPN_INCR_U (c3 + kk1, twor - k - 1, cy); /* 2n-(5k+1) = 2r-k-1 */
  134.   } else { /* triggered only by very unbalanced cases like
  135.       (k+k+(k-2))x(k+k+1) , should be handled by toom32 */
  136.     ASSERT_NOCARRY (mpn_add_n (vinf, vinf, v2 + k, twor));
  137.   }
  138.   /* (7) v1 <- v1 - vinf,       (1 0 1 0 0) - (1 0 0 0 0) = (0 0 1 0 0)
  139.      result is >= 0 */
  140.   /* Side effect: we also subtracted (high half) vm1 -= v2 */
  141.   cy = mpn_sub_n (v1, v1, vinf, twor);          /* vinf is at most twor long.  */
  142.   vinf0 = vinf[0];                     /* Save again the right value for vinf0 */
  143.   vinf[0] = saved;
  144.   MPN_DECR_U (v1 + twor, kk1 - twor, cy);       /* Treat the last bytes.       */
  145.   /* (8) vm1 <- vm1-v2          (0 1 0 1 0) - (0 1 0 0 0) = (0 0 0 1 0)
  146.      Operate only on the low half.
  147.   */
  148.   cy = mpn_sub_n (c1, c1, v2, k);
  149.   MPN_DECR_U (v1, kk1, cy);
  150.   /********************* Beginning the final phase **********************/
  151.   /* Most of the recomposition was done */
  152.   /* add t2 in {c+3k, ...}, but only the low half */
  153.   cy = mpn_add_n (c3, c3, v2, k);
  154.   vinf[0] += cy;
  155.   ASSERT(vinf[0] >= cy); /* No carry */
  156.   MPN_INCR_U (vinf, twor, vinf0); /* Add vinf0, propagate carry. */
  157. #undef v0
  158. }