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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_mul -- Multiply two natural numbers.
  2.    Contributed to the GNU project by Torbjorn Granlund.
  3. Copyright 1991, 1993, 1994, 1996, 1997, 1999, 2000, 2001, 2002, 2003, 2005,
  4. 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
  5. This file is part of the GNU MP Library.
  6. The GNU MP Library is free software; you can redistribute it and/or modify
  7. it under the terms of the GNU Lesser General Public License as published by
  8. the Free Software Foundation; either version 3 of the License, or (at your
  9. option) any later version.
  10. The GNU MP Library is distributed in the hope that it will be useful, but
  11. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  13. License for more details.
  14. You should have received a copy of the GNU Lesser General Public License
  15. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  16. #include "gmp.h"
  17. #include "gmp-impl.h"
  18. #ifndef MUL_BASECASE_MAX_UN
  19. #define MUL_BASECASE_MAX_UN 500
  20. #endif
  21. #define TOOM33_OK(an,bn) (6 + 2 * an < 3 * bn)
  22. #define TOOM44_OK(an,bn) (12 + 3 * an < 4 * bn)
  23. /* Multiply the natural numbers u (pointed to by UP, with UN limbs) and v
  24.    (pointed to by VP, with VN limbs), and store the result at PRODP.  The
  25.    result is UN + VN limbs.  Return the most significant limb of the result.
  26.    NOTE: The space pointed to by PRODP is overwritten before finished with U
  27.    and V, so overlap is an error.
  28.    Argument constraints:
  29.    1. UN >= VN.
  30.    2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from
  31.       the multiplier and the multiplicand.  */
  32. /*
  33.   * The cutoff lines in the toomX2 and toomX3 code are now exactly between the
  34.     ideal lines of the surrounding algorithms.  Is that optimal?
  35.   * The toomX3 code now uses a structure similar to the one of toomX2, except
  36.     that it loops longer in the unbalanced case.  The result is that the
  37.     remaining area might have un < vn.  Should we fix the toomX2 code in a
  38.     similar way?
  39.   * The toomX3 code is used for the largest non-FFT unbalanced operands.  It
  40.     therefore calls mpn_mul recursively for certain cases.
  41.   * Allocate static temp space using THRESHOLD variables (except for toom44
  42.     when !WANT_FFT).  That way, we can typically have no TMP_ALLOC at all.
  43.   * We sort ToomX2 algorithms together, assuming the toom22, toom32, toom42
  44.     have the same vn threshold.  This is not true, we should actually use
  45.     mul_basecase for slightly larger operands for toom32 than for toom22, and
  46.     even larger for toom42.
  47.   * That problem is even more prevalent for toomX3.  We therefore use special
  48.     THRESHOLD variables there.
  49.   * Is our ITCH allocation correct?
  50. */
  51. #define ITCH (16*vn + 100)
  52. mp_limb_t
  53. mpn_mul (mp_ptr prodp,
  54.  mp_srcptr up, mp_size_t un,
  55.  mp_srcptr vp, mp_size_t vn)
  56. {
  57.   ASSERT (un >= vn);
  58.   ASSERT (vn >= 1);
  59.   ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
  60.   ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));
  61.   if (un == vn)
  62.     {
  63.       if (up == vp)
  64. mpn_sqr (prodp, up, un);
  65.       else
  66. mpn_mul_n (prodp, up, vp, un);
  67.     }
  68.   else if (vn < MUL_TOOM22_THRESHOLD)
  69.     { /* plain schoolbook multiplication */
  70.       /* Unless un is very large, or else if have an applicable mpn_mul_N,
  71.  perform basecase multiply directly.  */
  72.       if (un <= MUL_BASECASE_MAX_UN
  73. #if HAVE_NATIVE_mpn_mul_2
  74.   || vn <= 2
  75. #else
  76.   || vn == 1
  77. #endif
  78.   )
  79. mpn_mul_basecase (prodp, up, un, vp, vn);
  80.       else
  81. {
  82.   /* We have un >> MUL_BASECASE_MAX_UN > vn.  For better memory
  83.      locality, split up[] into MUL_BASECASE_MAX_UN pieces and multiply
  84.      these pieces with the vp[] operand.  After each such partial
  85.      multiplication (but the last) we copy the most significant vn
  86.      limbs into a temporary buffer since that part would otherwise be
  87.      overwritten by the next multiplication.  After the next
  88.      multiplication, we add it back.  This illustrates the situation:
  89.                                                     -->vn<--
  90.                                                       |  |<------- un ------->|
  91.                                                          _____________________|
  92.                                                         X                    /|
  93.                                                       /XX__________________/  |
  94.                                     _____________________                     |
  95.                                    X                    /                     |
  96.                                  /XX__________________/                       |
  97.                _____________________                                          |
  98.               /                    /                                          |
  99.             /____________________/                                            |
  100.     ==================================================================
  101.     The parts marked with X are the parts whose sums are copied into
  102.     the temporary buffer.  */
  103.   mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
  104.   mp_limb_t cy;
  105.   ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);
  106.   mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
  107.   prodp += MUL_BASECASE_MAX_UN;
  108.   MPN_COPY (tp, prodp, vn); /* preserve high triangle */
  109.   up += MUL_BASECASE_MAX_UN;
  110.   un -= MUL_BASECASE_MAX_UN;
  111.   while (un > MUL_BASECASE_MAX_UN)
  112.     {
  113.       mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
  114.       cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
  115.       mpn_incr_u (prodp + vn, cy);
  116.       prodp += MUL_BASECASE_MAX_UN;
  117.       MPN_COPY (tp, prodp, vn); /* preserve high triangle */
  118.       up += MUL_BASECASE_MAX_UN;
  119.       un -= MUL_BASECASE_MAX_UN;
  120.     }
  121.   if (un > vn)
  122.     {
  123.       mpn_mul_basecase (prodp, up, un, vp, vn);
  124.     }
  125.   else
  126.     {
  127.       ASSERT (un > 0);
  128.       mpn_mul_basecase (prodp, vp, vn, up, un);
  129.     }
  130.   cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
  131.   mpn_incr_u (prodp + vn, cy);
  132. }
  133.     }
  134.   else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
  135.     {
  136.       /* Use ToomX2 variants */
  137.       mp_ptr scratch;
  138.       TMP_SDECL; TMP_SMARK;
  139.       scratch = TMP_SALLOC_LIMBS (ITCH);
  140.       /* FIXME: This condition (repeated in the loop below) leaves from a vn*vn
  141.  square to a (3vn-1)*vn rectangle.  Leaving such a rectangle is hardly
  142.  wise; we would get better balance by slightly moving the bound.  We
  143.  will sometimes end up with un < vn, like the the X3 arm below.  */
  144.       if (un >= 3 * vn)
  145. {
  146.   mp_limb_t cy;
  147.   mp_ptr ws;
  148.   /* The maximum ws usage is for the mpn_mul result.  */
  149.   ws = TMP_SALLOC_LIMBS (4 * vn);
  150.   mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
  151.   un -= 2 * vn;
  152.   up += 2 * vn;
  153.   prodp += 2 * vn;
  154.   while (un >= 3 * vn)
  155.     {
  156.       mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
  157.       un -= 2 * vn;
  158.       up += 2 * vn;
  159.       cy = mpn_add_n (prodp, prodp, ws, vn);
  160.       MPN_COPY (prodp + vn, ws + vn, 2 * vn);
  161.       mpn_incr_u (prodp + vn, cy);
  162.       prodp += 2 * vn;
  163.     }
  164.   /* vn <= un < 3vn */
  165.   if (4 * un < 5 * vn)
  166.     mpn_toom22_mul (ws, up, un, vp, vn, scratch);
  167.   else if (4 * un < 7 * vn)
  168.     mpn_toom32_mul (ws, up, un, vp, vn, scratch);
  169.   else
  170.     mpn_toom42_mul (ws, up, un, vp, vn, scratch);
  171.   cy = mpn_add_n (prodp, prodp, ws, vn);
  172.   MPN_COPY (prodp + vn, ws + vn, un);
  173.   mpn_incr_u (prodp + vn, cy);
  174. }
  175.       else
  176. {
  177.   if (4 * un < 5 * vn)
  178.     mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
  179.   else if (4 * un < 7 * vn)
  180.     mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
  181.   else
  182.     mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
  183. }
  184.       TMP_SFREE;
  185.     }
  186.   else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) ||
  187.    BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
  188.     {
  189.       /* Handle the largest operands that are not in the FFT range.  The 2nd
  190.  condition makes very unbalanced operands avoid the FFT code (except
  191.  perhaps as coefficient products of the Toom code.  */
  192.       if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) || !TOOM44_OK (un, vn))
  193. {
  194.   /* Use ToomX3 variants */
  195.   mp_ptr scratch;
  196.   TMP_SDECL; TMP_SMARK;
  197.   scratch = TMP_SALLOC_LIMBS (ITCH);
  198.   if (2 * un >= 5 * vn)
  199.     {
  200.       mp_limb_t cy;
  201.       mp_ptr ws;
  202.       /* The maximum ws usage is for the mpn_mul result.  */
  203.       ws = TMP_SALLOC_LIMBS (7 * vn >> 1);
  204.       if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
  205. mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
  206.       else
  207. mpn_toom63_mul (prodp, up, 2 * vn, vp, vn, scratch);
  208.       un -= 2 * vn;
  209.       up += 2 * vn;
  210.       prodp += 2 * vn;
  211.       while (2 * un >= 5 * vn) /* un >= 2.5vn */
  212. {
  213.   if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
  214.     mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
  215.   else
  216.     mpn_toom63_mul (ws, up, 2 * vn, vp, vn, scratch);
  217.   un -= 2 * vn;
  218.   up += 2 * vn;
  219.   cy = mpn_add_n (prodp, prodp, ws, vn);
  220.   MPN_COPY (prodp + vn, ws + vn, 2 * vn);
  221.   mpn_incr_u (prodp + vn, cy);
  222.   prodp += 2 * vn;
  223. }
  224.       /* vn / 2 <= un < 2.5vn */
  225.       if (un < vn)
  226. mpn_mul (ws, vp, vn, up, un);
  227.       else
  228. mpn_mul (ws, up, un, vp, vn);
  229.       cy = mpn_add_n (prodp, prodp, ws, vn);
  230.       MPN_COPY (prodp + vn, ws + vn, un);
  231.       mpn_incr_u (prodp + vn, cy);
  232.     }
  233.   else
  234.     {
  235.       if (6 * un < 7 * vn)
  236. mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
  237.       else if (2 * un < 3 * vn)
  238. {
  239.   if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
  240.     mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
  241.   else
  242.     mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
  243. }
  244.       else if (6 * un < 11 * vn)
  245. {
  246.   if (4 * un < 7 * vn)
  247.     {
  248.       if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
  249. mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
  250.       else
  251. mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
  252.     }
  253.   else
  254.     {
  255.       if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM53_THRESHOLD))
  256. mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
  257.       else
  258. mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
  259.     }
  260. }
  261.       else
  262. {
  263.   if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
  264.     mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
  265.   else
  266.     mpn_toom63_mul (prodp, up, un, vp, vn, scratch);
  267. }
  268.     }
  269.   TMP_SFREE;
  270. }
  271.       else
  272. {
  273.   mp_ptr scratch;
  274.   TMP_DECL; TMP_MARK;
  275.   if (BELOW_THRESHOLD (vn, MUL_TOOM6H_THRESHOLD))
  276.     {
  277.       scratch = TMP_ALLOC_LIMBS (mpn_toom44_mul_itch (un, vn));
  278.       mpn_toom44_mul (prodp, up, un, vp, vn, scratch);
  279.     }
  280.   else if (BELOW_THRESHOLD (vn, MUL_TOOM8H_THRESHOLD))
  281.     {
  282.       scratch = TMP_ALLOC_LIMBS (mpn_toom6h_mul_itch (un, vn));
  283.       mpn_toom6h_mul (prodp, up, un, vp, vn, scratch);
  284.     }
  285.   else
  286.     {
  287.       scratch = TMP_ALLOC_LIMBS (mpn_toom8h_mul_itch (un, vn));
  288.       mpn_toom8h_mul (prodp, up, un, vp, vn, scratch);
  289.     }
  290.   TMP_FREE;
  291. }
  292.     }
  293.   else
  294.     {
  295.       if (un >= 8 * vn)
  296. {
  297.   mp_limb_t cy;
  298.   mp_ptr ws;
  299.   TMP_DECL; TMP_MARK;
  300.   /* The maximum ws usage is for the mpn_mul result.  */
  301.   ws = TMP_BALLOC_LIMBS (9 * vn >> 1);
  302.   mpn_fft_mul (prodp, up, 3 * vn, vp, vn);
  303.   un -= 3 * vn;
  304.   up += 3 * vn;
  305.   prodp += 3 * vn;
  306.   while (2 * un >= 7 * vn) /* un >= 3.5vn  */
  307.     {
  308.       mpn_fft_mul (ws, up, 3 * vn, vp, vn);
  309.       un -= 3 * vn;
  310.       up += 3 * vn;
  311.       cy = mpn_add_n (prodp, prodp, ws, vn);
  312.       MPN_COPY (prodp + vn, ws + vn, 3 * vn);
  313.       mpn_incr_u (prodp + vn, cy);
  314.       prodp += 3 * vn;
  315.     }
  316.   /* vn / 2 <= un < 3.5vn */
  317.   if (un < vn)
  318.     mpn_mul (ws, vp, vn, up, un);
  319.   else
  320.     mpn_mul (ws, up, un, vp, vn);
  321.   cy = mpn_add_n (prodp, prodp, ws, vn);
  322.   MPN_COPY (prodp + vn, ws + vn, un);
  323.   mpn_incr_u (prodp + vn, cy);
  324.   TMP_FREE;
  325. }
  326.       else
  327. mpn_fft_mul (prodp, up, un, vp, vn);
  328.     }
  329.   return prodp[un + vn - 1]; /* historic */
  330. }