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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom4_sqr -- Square {ap,an}.
  2.    Contributed to the GNU project by Torbjorn Granlund and 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 2006, 2007, 2008, 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. /* Evaluate in: -1, -1/2, 0, +1/2, +1, +2, +inf
  21.   <-s--><--n--><--n--><--n-->
  22.    ____ ______ ______ ______
  23.   |_a3_|___a2_|___a1_|___a0_|
  24.   v0  =   a0             ^2 #    A(0)^2
  25.   v1  = ( a0+ a1+ a2+ a3)^2 #    A(1)^2   ah  <= 3
  26.   vm1 = ( a0- a1+ a2- a3)^2 #   A(-1)^2  |ah| <= 1
  27.   v2  = ( a0+2a1+4a2+8a3)^2 #    A(2)^2   ah  <= 14
  28.   vh  = (8a0+4a1+2a2+ a3)^2 #  A(1/2)^2   ah  <= 14
  29.   vmh = (8a0-4a1+2a2- a3)^2 # A(-1/2)^2  -4<=ah<=9
  30.   vinf=               a3 ^2 #  A(inf)^2
  31. */
  32. #if TUNE_PROGRAM_BUILD
  33. #define MAYBE_sqr_basecase 1
  34. #define MAYBE_sqr_toom2   1
  35. #define MAYBE_sqr_toom4   1
  36. #else
  37. #define MAYBE_sqr_basecase
  38.   (SQR_TOOM4_THRESHOLD < 4 * SQR_TOOM2_THRESHOLD)
  39. #define MAYBE_sqr_toom2
  40.   (SQR_TOOM4_THRESHOLD < 4 * SQR_TOOM3_THRESHOLD)
  41. #define MAYBE_sqr_toom4
  42.   (SQR_FFT_THRESHOLD >= 4 * SQR_TOOM4_THRESHOLD)
  43. #endif
  44. #define TOOM4_SQR_REC(p, a, n, ws)
  45.   do {
  46.     if (MAYBE_sqr_basecase
  47. && BELOW_THRESHOLD (n, SQR_TOOM2_THRESHOLD))
  48.       mpn_sqr_basecase (p, a, n);
  49.     else if (MAYBE_sqr_toom2
  50.      && BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD))
  51.       mpn_toom2_sqr (p, a, n, ws);
  52.     else if (! MAYBE_sqr_toom4
  53.      || BELOW_THRESHOLD (n, SQR_TOOM4_THRESHOLD))
  54.       mpn_toom3_sqr (p, a, n, ws);
  55.     else
  56.       mpn_toom4_sqr (p, a, n, ws);
  57.   } while (0)
  58. void
  59. mpn_toom4_sqr (mp_ptr pp,
  60.        mp_srcptr ap, mp_size_t an,
  61.        mp_ptr scratch)
  62. {
  63.   mp_size_t n, s;
  64.   mp_limb_t cy;
  65. #define a0  ap
  66. #define a1  (ap + n)
  67. #define a2  (ap + 2*n)
  68. #define a3  (ap + 3*n)
  69.   n = (an + 3) >> 2;
  70.   s = an - 3 * n;
  71.   ASSERT (0 < s && s <= n);
  72.   /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
  73.    * following limb, so these must be computed in order, and we need a
  74.    * one limb gap to tp. */
  75. #define v0    pp /* 2n */
  76. #define v1    (pp + 2 * n) /* 2n+1 */
  77. #define vinf  (pp + 6 * n) /* s+t */
  78. #define v2    scratch /* 2n+1 */
  79. #define vm2   (scratch + 2 * n + 1) /* 2n+1 */
  80. #define vh    (scratch + 4 * n + 2) /* 2n+1 */
  81. #define vm1   (scratch + 6 * n + 3) /* 2n+1 */
  82. #define tp (scratch + 8*n + 5)
  83.   /* No overlap with v1 */
  84. #define apx   pp /* n+1 */
  85. #define amx   (pp + 4*n + 2) /* n+1 */
  86.   /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
  87.      gives roughly 32 n/3 + log term. */
  88.   /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3.  */
  89.   mpn_toom_eval_dgr3_pm2 (apx, amx, ap, n, s, tp);
  90.   TOOM4_SQR_REC (v2, apx, n + 1, tp); /* v2,  2n+1 limbs */
  91.   TOOM4_SQR_REC (vm2, amx, n + 1, tp); /* vm2,  2n+1 limbs */
  92.   /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
  93. #if HAVE_NATIVE_mpn_addlsh1_n
  94.   cy = mpn_addlsh1_n (apx, a1, a0, n);
  95.   cy = 2*cy + mpn_addlsh1_n (apx, a2, apx, n);
  96.   if (s < n)
  97.     {
  98.       mp_limb_t cy2;
  99.       cy2 = mpn_addlsh1_n (apx, a3, apx, s);
  100.       apx[n] = 2*cy + mpn_lshift (apx + s, apx + s, n - s, 1);
  101.       MPN_INCR_U (apx + s, n+1-s, cy2);
  102.     }
  103.   else
  104.     apx[n] = 2*cy + mpn_addlsh1_n (apx, a3, apx, n);
  105. #else
  106.   cy = mpn_lshift (apx, a0, n, 1);
  107.   cy += mpn_add_n (apx, apx, a1, n);
  108.   cy = 2*cy + mpn_lshift (apx, apx, n, 1);
  109.   cy += mpn_add_n (apx, apx, a2, n);
  110.   cy = 2*cy + mpn_lshift (apx, apx, n, 1);
  111.   apx[n] = cy + mpn_add (apx, apx, n, a3, s);
  112. #endif
  113.   ASSERT (apx[n] < 15);
  114.   TOOM4_SQR_REC (vh, apx, n + 1, tp); /* vh,  2n+1 limbs */
  115.   /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3.  */
  116.   mpn_toom_eval_dgr3_pm1 (apx, amx, ap, n, s, tp);
  117.   TOOM4_SQR_REC (v1, apx, n + 1, tp); /* v1,  2n+1 limbs */
  118.   TOOM4_SQR_REC (vm1, amx, n + 1, tp); /* vm1,  2n+1 limbs */
  119.   TOOM4_SQR_REC (v0, a0, n, tp);
  120.   TOOM4_SQR_REC (vinf, a3, s, tp); /* vinf, 2s limbs */
  121.   mpn_toom_interpolate_7pts (pp, n, 0, vm2, vm1, v2, vh, 2*s, tp);
  122. }