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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom44_mul -- Multiply {ap,an} and {bp,bn} where an and bn are close in
  2.    size.  Or more accurately, bn <= an < (4/3)bn.
  3.    Contributed to the GNU project by Torbjorn Granlund 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 2006, 2007, 2008 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. /* Evaluate in: 0, +1, -1, +2, -2, 1/2, +inf
  22.   <-s--><--n--><--n--><--n-->
  23.    ____ ______ ______ ______
  24.   |_a3_|___a2_|___a1_|___a0_|
  25.    |b3_|___b2_|___b1_|___b0_|
  26.    <-t-><--n--><--n--><--n-->
  27.   v0  =   a0             *  b0              #    A(0)*B(0)
  28.   v1  = ( a0+ a1+ a2+ a3)*( b0+ b1+ b2+ b3) #    A(1)*B(1)      ah  <= 3   bh  <= 3
  29.   vm1 = ( a0- a1+ a2- a3)*( b0- b1+ b2- b3) #   A(-1)*B(-1)    |ah| <= 1  |bh| <= 1
  30.   v2  = ( a0+2a1+4a2+8a3)*( b0+2b1+4b2+8b3) #    A(2)*B(2)      ah  <= 14  bh  <= 14
  31.   vm2 = ( a0-2a1+4a2-8a3)*( b0-2b1+4b2-8b3) #    A(2)*B(2)      ah  <= 9  |bh| <= 9
  32.   vh  = (8a0+4a1+2a2+ a3)*(8b0+4b1+2b2+ b3) #  A(1/2)*B(1/2)    ah  <= 14  bh  <= 14
  33.   vinf=               a3 *          b2      #  A(inf)*B(inf)
  34. */
  35. #if TUNE_PROGRAM_BUILD
  36. #define MAYBE_mul_basecase 1
  37. #define MAYBE_mul_toom22   1
  38. #define MAYBE_mul_toom44   1
  39. #else
  40. #define MAYBE_mul_basecase
  41.   (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM22_THRESHOLD)
  42. #define MAYBE_mul_toom22
  43.   (MUL_TOOM44_THRESHOLD < 4 * MUL_TOOM33_THRESHOLD)
  44. #define MAYBE_mul_toom44
  45.   (MUL_FFT_THRESHOLD >= 4 * MUL_TOOM44_THRESHOLD)
  46. #endif
  47. #define TOOM44_MUL_N_REC(p, a, b, n, ws)
  48.   do {
  49.     if (MAYBE_mul_basecase
  50. && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD))
  51.       mpn_mul_basecase (p, a, n, b, n);
  52.     else if (MAYBE_mul_toom22
  53.      && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD))
  54.       mpn_toom22_mul (p, a, n, b, n, ws);
  55.     else if (! MAYBE_mul_toom44
  56.      || BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD))
  57.       mpn_toom33_mul (p, a, n, b, n, ws);
  58.     else
  59.       mpn_toom44_mul (p, a, n, b, n, ws);
  60.   } while (0)
  61. /* Use of scratch space. In the product area, we store
  62.       ___________________
  63.      |vinf|____|_v1_|_v0_|
  64.       s+t  2n-1 2n+1  2n
  65.    The other recursive products, vm1, v2, vm2, vh are stored in the
  66.    scratch area. When computing them, we use the product area for
  67.    intermediate values.
  68.    Next, we compute v1. We can store the intermediate factors at v0
  69.    and at vh + 2n + 2.
  70.    Finally, for v0 and vinf, factors are parts of the input operands,
  71.    and we need scratch space only for the recursive multiplication.
  72.    In all, if S(an) is the scratch need, the needed space is bounded by
  73.      S(an) <= 4 (2*ceil(an/4) + 1) + 1 + S(ceil(an/4) + 1)
  74.    which should give S(n) = 8 n/3 + c log(n) for some constant c.
  75. */
  76. void
  77. mpn_toom44_mul (mp_ptr pp,
  78. mp_srcptr ap, mp_size_t an,
  79. mp_srcptr bp, mp_size_t bn,
  80. mp_ptr scratch)
  81. {
  82.   mp_size_t n, s, t;
  83.   mp_limb_t cy;
  84.   enum toom7_flags flags;
  85. #define a0  ap
  86. #define a1  (ap + n)
  87. #define a2  (ap + 2*n)
  88. #define a3  (ap + 3*n)
  89. #define b0  bp
  90. #define b1  (bp + n)
  91. #define b2  (bp + 2*n)
  92. #define b3  (bp + 3*n)
  93.   ASSERT (an >= bn);
  94.   n = (an + 3) >> 2;
  95.   s = an - 3 * n;
  96.   t = bn - 3 * n;
  97.   ASSERT (0 < s && s <= n);
  98.   ASSERT (0 < t && t <= n);
  99.   ASSERT (s >= t);
  100.   /* NOTE: The multiplications to v2, vm2, vh and vm1 overwrites the
  101.    * following limb, so these must be computed in order, and we need a
  102.    * one limb gap to tp. */
  103. #define v0    pp /* 2n */
  104. #define v1    (pp + 2 * n) /* 2n+1 */
  105. #define vinf  (pp + 6 * n) /* s+t */
  106. #define v2    scratch /* 2n+1 */
  107. #define vm2   (scratch + 2 * n + 1) /* 2n+1 */
  108. #define vh    (scratch + 4 * n + 2) /* 2n+1 */
  109. #define vm1   (scratch + 6 * n + 3) /* 2n+1 */
  110. #define tp (scratch + 8*n + 5)
  111.   /* apx and bpx must not overlap with v1 */
  112. #define apx   pp /* n+1 */
  113. #define amx   (pp + n + 1) /* n+1 */
  114. #define bmx   (pp + 2*n + 2) /* n+1 */
  115. #define bpx   (pp + 4*n + 2) /* n+1 */
  116.   /* Total scratch need: 8*n + 5 + scratch for recursive calls. This
  117.      gives roughly 32 n/3 + log term. */
  118.   /* Compute apx = a0 + 2 a1 + 4 a2 + 8 a3 and amx = a0 - 2 a1 + 4 a2 - 8 a3.  */
  119.   flags = toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (apx, amx, ap, n, s, tp);
  120.   /* Compute bpx = b0 + 2 b1 + 4 b2 + 8 b3 and bmx = b0 - 2 b1 + 4 b2 - 8 b3.  */
  121.   flags ^= toom7_w1_neg & mpn_toom_eval_dgr3_pm2 (bpx, bmx, bp, n, t, tp);
  122.   TOOM44_MUL_N_REC (v2, apx, bpx, n + 1, tp); /* v2,  2n+1 limbs */
  123.   TOOM44_MUL_N_REC (vm2, amx, bmx, n + 1, tp); /* vm2,  2n+1 limbs */
  124.   /* Compute apx = 8 a0 + 4 a1 + 2 a2 + a3 = (((2*a0 + a1) * 2 + a2) * 2 + a3 */
  125. #if HAVE_NATIVE_mpn_addlsh1_n
  126.   cy = mpn_addlsh1_n (apx, a1, a0, n);
  127.   cy = 2*cy + mpn_addlsh1_n (apx, a2, apx, n);
  128.   if (s < n)
  129.     {
  130.       mp_limb_t cy2;
  131.       cy2 = mpn_addlsh1_n (apx, a3, apx, s);
  132.       apx[n] = 2*cy + mpn_lshift (apx + s, apx + s, n - s, 1);
  133.       MPN_INCR_U (apx + s, n+1-s, cy2);
  134.     }
  135.   else
  136.     apx[n] = 2*cy + mpn_addlsh1_n (apx, a3, apx, n);
  137. #else
  138.   cy = mpn_lshift (apx, a0, n, 1);
  139.   cy += mpn_add_n (apx, apx, a1, n);
  140.   cy = 2*cy + mpn_lshift (apx, apx, n, 1);
  141.   cy += mpn_add_n (apx, apx, a2, n);
  142.   cy = 2*cy + mpn_lshift (apx, apx, n, 1);
  143.   apx[n] = cy + mpn_add (apx, apx, n, a3, s);
  144. #endif
  145.   /* Compute bpx = 8 b0 + 4 b1 + 2 b2 + b3 = (((2*b0 + b1) * 2 + b2) * 2 + b3 */
  146. #if HAVE_NATIVE_mpn_addlsh1_n
  147.   cy = mpn_addlsh1_n (bpx, b1, b0, n);
  148.   cy = 2*cy + mpn_addlsh1_n (bpx, b2, bpx, n);
  149.   if (t < n)
  150.     {
  151.       mp_limb_t cy2;
  152.       cy2 = mpn_addlsh1_n (bpx, b3, bpx, t);
  153.       bpx[n] = 2*cy + mpn_lshift (bpx + t, bpx + t, n - t, 1);
  154.       MPN_INCR_U (bpx + t, n+1-t, cy2);
  155.     }
  156.   else
  157.     bpx[n] = 2*cy + mpn_addlsh1_n (bpx, b3, bpx, n);
  158. #else
  159.   cy = mpn_lshift (bpx, b0, n, 1);
  160.   cy += mpn_add_n (bpx, bpx, b1, n);
  161.   cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
  162.   cy += mpn_add_n (bpx, bpx, b2, n);
  163.   cy = 2*cy + mpn_lshift (bpx, bpx, n, 1);
  164.   bpx[n] = cy + mpn_add (bpx, bpx, n, b3, t);
  165. #endif
  166.   ASSERT (apx[n] < 15);
  167.   ASSERT (bpx[n] < 15);
  168.   TOOM44_MUL_N_REC (vh, apx, bpx, n + 1, tp); /* vh,  2n+1 limbs */
  169.   /* Compute apx = a0 + a1 + a2 + a3 and amx = a0 - a1 + a2 - a3.  */
  170.   flags |= toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (apx, amx, ap, n, s, tp);
  171.   /* Compute bpx = b0 + b1 + b2 + b3 bnd bmx = b0 - b1 + b2 - b3.  */
  172.   flags ^= toom7_w3_neg & mpn_toom_eval_dgr3_pm1 (bpx, bmx, bp, n, t, tp);
  173.   TOOM44_MUL_N_REC (vm1, amx, bmx, n + 1, tp); /* vm1,  2n+1 limbs */
  174.   /* Clobbers amx, bmx. */
  175.   TOOM44_MUL_N_REC (v1, apx, bpx, n + 1, tp); /* v1,  2n+1 limbs */
  176.   TOOM44_MUL_N_REC (v0, a0, b0, n, tp);
  177.   if (s > t)
  178.     mpn_mul (vinf, a3, s, b3, t);
  179.   else
  180.     TOOM44_MUL_N_REC (vinf, a3, b3, s, tp); /* vinf, s+t limbs */
  181.   mpn_toom_interpolate_7pts (pp, n, flags, vm2, vm1, v2, vh, s + t, tp);
  182. }