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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_toom32_mul -- Multiply {ap,an} and {bp,bn} where an is nominally 1.5
  2.    times as large as bn.  Or more accurately, bn < an < 3bn.
  3.    Contributed to the GNU project by Torbjorn Granlund.
  4.    Improvements by Marco Bodrato and Niels M鰈ler.
  5.    The idea of applying toom to unbalanced multiplication is due to Marco
  6.    Bodrato and Alberto Zanoni.
  7.    THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE.  IT IS ONLY
  8.    SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST
  9.    GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  10. Copyright 2006, 2007, 2008, 2009, 2010 Free Software Foundation, Inc.
  11. This file is part of the GNU MP Library.
  12. The GNU MP Library is free software; you can redistribute it and/or modify
  13. it under the terms of the GNU Lesser General Public License as published by
  14. the Free Software Foundation; either version 3 of the License, or (at your
  15. option) any later version.
  16. The GNU MP Library is distributed in the hope that it will be useful, but
  17. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  18. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  19. License for more details.
  20. You should have received a copy of the GNU Lesser General Public License
  21. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  22. #include "gmp.h"
  23. #include "gmp-impl.h"
  24. /* Evaluate in: -1, 0, +1, +inf
  25.   <-s-><--n--><--n-->
  26.    ___ ______ ______
  27.   |a2_|___a1_|___a0_|
  28. |_b1_|___b0_|
  29. <-t--><--n-->
  30.   v0  =  a0         * b0      #   A(0)*B(0)
  31.   v1  = (a0+ a1+ a2)*(b0+ b1) #   A(1)*B(1)      ah  <= 2  bh <= 1
  32.   vm1 = (a0- a1+ a2)*(b0- b1) #  A(-1)*B(-1)    |ah| <= 1  bh = 0
  33.   vinf=          a2 *     b1  # A(inf)*B(inf)
  34. */
  35. #define TOOM32_MUL_N_REC(p, a, b, n, ws)
  36.   do {
  37.     mpn_mul_n (p, a, b, n);
  38.   } while (0)
  39. void
  40. mpn_toom32_mul (mp_ptr pp,
  41. mp_srcptr ap, mp_size_t an,
  42. mp_srcptr bp, mp_size_t bn,
  43. mp_ptr scratch)
  44. {
  45.   mp_size_t n, s, t;
  46.   int vm1_neg;
  47.   mp_limb_t cy;
  48.   int hi;
  49.   mp_limb_t ap1_hi, bp1_hi;
  50. #define a0  ap
  51. #define a1  (ap + n)
  52. #define a2  (ap + 2 * n)
  53. #define b0  bp
  54. #define b1  (bp + n)
  55.   /* Required, to ensure that s + t >= n. */
  56.   ASSERT (bn + 2 <= an && an + 6 <= 3*bn);
  57.   n = 1 + (2 * an >= 3 * bn ? (an - 1) / (size_t) 3 : (bn - 1) >> 1);
  58.   s = an - 2 * n;
  59.   t = bn - n;
  60.   ASSERT (0 < s && s <= n);
  61.   ASSERT (0 < t && t <= n);
  62.   ASSERT (s + t >= n);
  63.   /* Product area of size an + bn = 3*n + s + t >= 4*n + 2. */
  64. #define ap1 (pp) /* n, most significant limb in ap1_hi */
  65. #define bp1 (pp + n) /* n, most significant bit in bp1_hi */
  66. #define am1 (pp + 2*n) /* n, most significant bit in hi */
  67. #define bm1 (pp + 3*n) /* n */
  68. #define v1 (scratch) /* 2n + 1 */
  69. #define vm1 (pp) /* 2n + 1 */
  70. #define scratch_out (scratch + 2*n + 1) /* Currently unused. */
  71.   /* Scratch need: 2*n + 1 + scratch for the recursive multiplications. */
  72.   /* FIXME: Keep v1[2*n] and vm1[2*n] in scalar variables? */
  73.   /* Compute ap1 = a0 + a1 + a3, am1 = a0 - a1 + a3 */
  74.   ap1_hi = mpn_add (ap1, a0, n, a2, s);
  75. #if HAVE_NATIVE_mpn_add_n_sub_n
  76.   if (ap1_hi == 0 && mpn_cmp (ap1, a1, n) < 0)
  77.     {
  78.       ap1_hi = mpn_add_n_sub_n (ap1, am1, a1, ap1, n) >> 1;
  79.       hi = 0;
  80.       vm1_neg = 1;
  81.     }
  82.   else
  83.     {
  84.       cy = mpn_add_n_sub_n (ap1, am1, ap1, a1, n);
  85.       hi = ap1_hi - (cy & 1);
  86.       ap1_hi += (cy >> 1);
  87.       vm1_neg = 0;
  88.     }
  89. #else
  90.   if (ap1_hi == 0 && mpn_cmp (ap1, a1, n) < 0)
  91.     {
  92.       ASSERT_NOCARRY (mpn_sub_n (am1, a1, ap1, n));
  93.       hi = 0;
  94.       vm1_neg = 1;
  95.     }
  96.   else
  97.     {
  98.       hi = ap1_hi - mpn_sub_n (am1, ap1, a1, n);
  99.       vm1_neg = 0;
  100.     }
  101.   ap1_hi += mpn_add_n (ap1, ap1, a1, n);
  102. #endif
  103.   /* Compute bp1 = b0 + b1 and bm1 = b0 - b1. */
  104.   if (t == n)
  105.     {
  106. #if HAVE_NATIVE_mpn_add_n_sub_n
  107.       if (mpn_cmp (b0, b1, n) < 0)
  108. {
  109.   cy = mpn_add_n_sub_n (bp1, bm1, b1, b0, n);
  110.   vm1_neg ^= 1;
  111. }
  112.       else
  113. {
  114.   cy = mpn_add_n_sub_n (bp1, bm1, b0, b1, n);
  115. }
  116.       bp1_hi = cy >> 1;
  117. #else
  118.       bp1_hi = mpn_add_n (bp1, b0, b1, n);
  119.       if (mpn_cmp (b0, b1, n) < 0)
  120. {
  121.   ASSERT_NOCARRY (mpn_sub_n (bm1, b1, b0, n));
  122.   vm1_neg ^= 1;
  123. }
  124.       else
  125. {
  126.   ASSERT_NOCARRY (mpn_sub_n (bm1, b0, b1, n));
  127. }
  128. #endif
  129.     }
  130.   else
  131.     {
  132.       /* FIXME: Should still use mpn_add_n_sub_n for the main part. */
  133.       bp1_hi = mpn_add (bp1, b0, n, b1, t);
  134.       if (mpn_zero_p (b0 + t, n - t) && mpn_cmp (b0, b1, t) < 0)
  135. {
  136.   ASSERT_NOCARRY (mpn_sub_n (bm1, b1, b0, t));
  137.   MPN_ZERO (bm1 + t, n - t);
  138.   vm1_neg ^= 1;
  139. }
  140.       else
  141. {
  142.   ASSERT_NOCARRY (mpn_sub (bm1, b0, n, b1, t));
  143. }
  144.     }
  145.   TOOM32_MUL_N_REC (v1, ap1, bp1, n, scratch_out);
  146.   if (ap1_hi == 1)
  147.     {
  148.       cy = bp1_hi + mpn_add_n (v1 + n, v1 + n, bp1, n);
  149.     }
  150.   else if (ap1_hi == 2)
  151.     {
  152. #if HAVE_NATIVE_mpn_addlsh1_n
  153.       cy = 2 * bp1_hi + mpn_addlsh1_n (v1 + n, v1 + n, bp1, n);
  154. #else
  155.       cy = 2 * bp1_hi + mpn_addmul_1 (v1 + n, bp1, n, CNST_LIMB(2));
  156. #endif
  157.     }
  158.   else
  159.     cy = 0;
  160.   if (bp1_hi != 0)
  161.     cy += mpn_add_n (v1 + n, v1 + n, ap1, n);
  162.   v1[2 * n] = cy;
  163.   TOOM32_MUL_N_REC (vm1, am1, bm1, n, scratch_out);
  164.   if (hi)
  165.     hi = mpn_add_n (vm1+n, vm1+n, bm1, n);
  166.   vm1[2*n] = hi;
  167.   /* v1 <-- (v1 + vm1) / 2 = x0 + x2 */
  168.   if (vm1_neg)
  169.     {
  170. #if HAVE_NATIVE_mpn_rsh1sub_n
  171.       mpn_rsh1sub_n (v1, v1, vm1, 2*n+1);
  172. #else
  173.       mpn_sub_n (v1, v1, vm1, 2*n+1);
  174.       ASSERT_NOCARRY (mpn_rshift (v1, v1, 2*n+1, 1));
  175. #endif
  176.     }
  177.   else
  178.     {
  179. #if HAVE_NATIVE_mpn_rsh1add_n
  180.       mpn_rsh1add_n (v1, v1, vm1, 2*n+1);
  181. #else
  182.       mpn_add_n (v1, v1, vm1, 2*n+1);
  183.       ASSERT_NOCARRY (mpn_rshift (v1, v1, 2*n+1, 1));
  184. #endif
  185.     }
  186.   /* We get x1 + x3 = (x0 + x2) - (x0 - x1 + x2 - x3), and hence
  187.      y = x1 + x3 + (x0 + x2) * B
  188.        = (x0 + x2) * B + (x0 + x2) - vm1.
  189.      y is 3*n + 1 limbs, y = y0 + y1 B + y2 B^2. We store them as
  190.      follows: y0 at scratch, y1 at pp + 2*n, and y2 at scratch + n
  191.      (already in place, except for carry propagation).
  192.      We thus add
  193.    B^3  B^2   B    1
  194.     |    |    |    |
  195.    +-----+----+
  196.  + |  x0 + x2 |
  197.    +----+-----+----+
  198.  +      |  x0 + x2 |
  199. +----------+
  200.  -      |  vm1     |
  201.  --+----++----+----+-
  202.    | y2  | y1 | y0 |
  203.    +-----+----+----+
  204.   Since we store y0 at the same location as the low half of x0 + x2, we
  205.   need to do the middle sum first. */
  206.   hi = vm1[2*n];
  207.   cy = mpn_add_n (pp + 2*n, v1, v1 + n, n);
  208.   MPN_INCR_U (v1 + n, n + 1, cy + v1[2*n]);
  209.   /* FIXME: Can we get rid of this second vm1_neg conditional by
  210.      swapping the location of +1 and -1 values? */
  211.   if (vm1_neg)
  212.     {
  213.       cy = mpn_add_n (v1, v1, vm1, n);
  214.       hi += mpn_add_nc (pp + 2*n, pp + 2*n, vm1 + n, n, cy);
  215.       MPN_INCR_U (v1 + n, n+1, hi);
  216.     }
  217.   else
  218.     {
  219.       cy = mpn_sub_n (v1, v1, vm1, n);
  220.       hi += mpn_sub_nc (pp + 2*n, pp + 2*n, vm1 + n, n, cy);
  221.       MPN_DECR_U (v1 + n, n+1, hi);
  222.     }
  223.   TOOM32_MUL_N_REC (pp, a0, b0, n, scratch_out);
  224.   /* vinf, s+t limbs.  Use mpn_mul for now, to handle unbalanced operands */
  225.   if (s > t)  mpn_mul (pp+3*n, a2, s, b1, t);
  226.   else        mpn_mul (pp+3*n, b1, t, a2, s);
  227.   /* Remaining interpolation.
  228.      y * B + x0 + x3 B^3 - x0 B^2 - x3 B
  229.      = (x1 + x3) B + (x0 + x2) B^2 + x0 + x3 B^3 - x0 B^2 - x3 B
  230.      = y0 B + y1 B^2 + y3 B^3 + Lx0 + H x0 B
  231.        + L x3 B^3 + H x3 B^4 - Lx0 B^2 - H x0 B^3 - L x3 B - H x3 B^2
  232.      = L x0 + (y0 + H x0 - L x3) B + (y1 - L x0 - H x3) B^2
  233.        + (y2 - (H x0 - L x3)) B^3 + H x3 B^4
  234.   B^4       B^3       B^2        B         1
  235.  |         |         |         |         |         |
  236.    +-------+                   +---------+---------+
  237.    |  Hx3  |                   | Hx0-Lx3 |    Lx0  |
  238.    +------+----------+---------+---------+---------+
  239.   |    y2    |  y1     |   y0    |
  240.   ++---------+---------+---------+
  241.   -| Hx0-Lx3 | - Lx0   |
  242.    +---------+---------+
  243.       | - Hx3  |
  244.       +--------+
  245.     We must take into account the carry from Hx0 - Lx3.
  246.   */
  247.   cy = mpn_sub_n (pp + n, pp + n, pp+3*n, n);
  248.   hi = scratch[2*n] + cy;
  249.   cy = mpn_sub_nc (pp + 2*n, pp + 2*n, pp, n, cy);
  250.   hi -= mpn_sub_nc (pp + 3*n, scratch + n, pp + n, n, cy);
  251.   hi += mpn_add (pp + n, pp + n, 3*n, scratch, n);
  252.   /* FIXME: Is support for s + t == n needed? */
  253.   if (LIKELY (s + t > n))
  254.     {
  255.       hi -= mpn_sub (pp + 2*n, pp + 2*n, 2*n, pp + 4*n, s+t-n);
  256.       if (hi < 0)
  257. MPN_DECR_U (pp + 4*n, s+t-n, -hi);
  258.       else
  259. MPN_INCR_U (pp + 4*n, s+t-n, hi);
  260.     }
  261.   else
  262.     ASSERT (hi == 0);
  263. }