/src/gmp/mpn/mul.c

Source
/* mpn_mul -- Multiply two natural numbers.

   Contributed to the GNU project by Torbjorn Granlund.

Copyright 1991, 1993, 1994, 1996, 1997, 1999-2003, 2005-2007, 2009, 2010, 2012,
2014, 2019 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 2 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the GNU MP Library.  If not,
see https://www.gnu.org/licenses/.  */

#include "gmp-impl.h"


#ifndef MUL_BASECASE_MAX_UN
#define MUL_BASECASE_MAX_UN 500
#endif

/* Areas where the different toom algorithms can be called (extracted
   from the t-toom*.c files, and ignoring small constant offsets):

   1/6  1/5 1/4 4/13 1/3 3/8 2/5 5/11 1/2 3/5 2/3 3/4 4/5   1 vn/un
                                        4/7              6/7
               6/11
                                       |--------------------| toom22 (small)
                                                           || toom22 (large)
                                                       |xxxx| toom22 called
                      |-------------------------------------| toom32
                                         |xxxxxxxxxxxxxxxx| | toom32 called
                                               |------------| toom33
                                                          |x| toom33 called
             |---------------------------------|            | toom42
                |xxxxxxxxxxxxxxxxxxxxxxxx|            | toom42 called
                                       |--------------------| toom43
                                               |xxxxxxxxxx|   toom43 called
         |-----------------------------|                      toom52 (unused)
                                                   |--------| toom44
               |xxxxxxxx| toom44 called
                              |--------------------|        | toom53
                                        |xxxxxx|              toom53 called
    |-------------------------|                               toom62 (unused)
                                           |----------------| toom54 (unused)
                      |--------------------|                  toom63
                        |xxxxxxxxx|                   | toom63 called
                          |---------------------------------| toom6h
               |xxxxxxxx| toom6h called
                                  |-------------------------| toom8h (32 bit)
                 |------------------------------------------| toom8h (64 bit)
               |xxxxxxxx| toom8h called
*/

#define TOOM33_OK(an,bn) (6 + 2 * an < 3 * bn)
#define TOOM44_OK(an,bn) (12 + 3 * an < 4 * bn)

/* Multiply the natural numbers u (pointed to by UP, with UN limbs) and v
   (pointed to by VP, with VN limbs), and store the result at PRODP.  The
   result is UN + VN limbs.  Return the most significant limb of the result.

   NOTE: The space pointed to by PRODP is overwritten before finished with U
   and V, so overlap is an error.

   Argument constraints:
   1. UN >= VN.
   2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from
      the multiplier and the multiplicand.  */

/*
  * The cutoff lines in the toomX2 and toomX3 code are now exactly between the
    ideal lines of the surrounding algorithms.  Is that optimal?

  * The toomX3 code now uses a structure similar to the one of toomX2, except
    that it loops longer in the unbalanced case.  The result is that the
    remaining area might have un < vn.  Should we fix the toomX2 code in a
    similar way?

  * The toomX3 code is used for the largest non-FFT unbalanced operands.  It
    therefore calls mpn_mul recursively for certain cases.

  * Allocate static temp space using THRESHOLD variables (except for toom44
    when !WANT_FFT).  That way, we can typically have no TMP_ALLOC at all.

  * We sort ToomX2 algorithms together, assuming the toom22, toom32, toom42
    have the same vn threshold.  This is not true, we should actually use
    mul_basecase for slightly larger operands for toom32 than for toom22, and
    even larger for toom42.

  * That problem is even more prevalent for toomX3.  We therefore use special
    THRESHOLD variables there.
*/

mp_limb_t
mpn_mul (mp_ptr prodp,
   mp_srcptr up, mp_size_t un,
   mp_srcptr vp, mp_size_t vn)
{
  ASSERT (un >= vn);
  ASSERT (vn >= 1);
  ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
  ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));

  if (BELOW_THRESHOLD (un, MUL_TOOM22_THRESHOLD))
    {
      /* When un (and thus vn) is below the toom22 range, do mul_basecase.
   Test un and not vn here not to thwart the un >> vn code below.
   This special case is not necessary, but cuts the overhead for the
   smallest operands. */
      mpn_mul_basecase (prodp, up, un, vp, vn);
    }
  else if (un == vn)
    {
      mpn_mul_n (prodp, up, vp, un);
    }
  else if (vn < MUL_TOOM22_THRESHOLD)
    { /* plain schoolbook multiplication */

      /* Unless un is very large, or else if have an applicable mpn_mul_N,
   perform basecase multiply directly.  */
      if (un <= MUL_BASECASE_MAX_UN
#if HAVE_NATIVE_mpn_mul_2
    || vn <= 2
#else
    || vn == 1
#endif
    )
  mpn_mul_basecase (prodp, up, un, vp, vn);
      else
  {
    /* We have un >> MUL_BASECASE_MAX_UN > vn.  For better memory
       locality, split up[] into MUL_BASECASE_MAX_UN pieces and multiply
       these pieces with the vp[] operand.  After each such partial
       multiplication (but the last) we copy the most significant vn
       limbs into a temporary buffer since that part would otherwise be
       overwritten by the next multiplication.  After the next
       multiplication, we add it back.  This illustrates the situation:

                                                    -->vn<--
                                                      |  |<------- un ------->|
                                                         _____________________|
                                                        X                    /|
                                                      /XX__________________/  |
                                    _____________________                     |
                                   X                    /                     |
                                 /XX__________________/                       |
               _____________________                                          |
              /                    /                                          |
            /____________________/                                            |
      ==================================================================

      The parts marked with X are the parts whose sums are copied into
      the temporary buffer.  */

    mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
    mp_limb_t cy;
    ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);

    mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
    prodp += MUL_BASECASE_MAX_UN;
    MPN_COPY (tp, prodp, vn);   /* preserve high triangle */
    up += MUL_BASECASE_MAX_UN;
    un -= MUL_BASECASE_MAX_UN;
    while (un > MUL_BASECASE_MAX_UN)
      {
        mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
        cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
        mpn_incr_u (prodp + vn, cy);
        prodp += MUL_BASECASE_MAX_UN;
        MPN_COPY (tp, prodp, vn);   /* preserve high triangle */
        up += MUL_BASECASE_MAX_UN;
        un -= MUL_BASECASE_MAX_UN;
      }
    if (un > vn)
      {
        mpn_mul_basecase (prodp, up, un, vp, vn);
      }
    else
      {
        ASSERT (un > 0);
        mpn_mul_basecase (prodp, vp, vn, up, un);
      }
    cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
    mpn_incr_u (prodp + vn, cy);
  }
    }
  else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
    {
      /* Use ToomX2 variants */
      mp_ptr scratch;
      TMP_SDECL; TMP_SMARK;

#define ITCH_TOOMX2 (9 * vn / 2 + GMP_NUMB_BITS * 2)
      scratch = TMP_SALLOC_LIMBS (ITCH_TOOMX2);
      ASSERT (mpn_toom22_mul_itch ((5*vn-1)/4, vn) <= ITCH_TOOMX2); /* 5vn/2+ */
      ASSERT (mpn_toom32_mul_itch ((7*vn-1)/4, vn) <= ITCH_TOOMX2); /* 7vn/6+ */
      ASSERT (mpn_toom42_mul_itch (3 * vn - 1, vn) <= ITCH_TOOMX2); /* 9vn/2+ */
#undef ITCH_TOOMX2

      /* FIXME: This condition (repeated in the loop below) leaves from a vn*vn
   square to a (3vn-1)*vn rectangle.  Leaving such a rectangle is hardly
   wise; we would get better balance by slightly moving the bound.  We
   will sometimes end up with un < vn, like in the X3 arm below.  */
      if (un >= 3 * vn)
  {
    mp_limb_t cy;
    mp_ptr ws;

    /* The maximum ws usage is for the mpn_mul result.  */
    ws = TMP_SALLOC_LIMBS (4 * vn);

    mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
    un -= 2 * vn;
    up += 2 * vn;
    prodp += 2 * vn;

    while (un >= 3 * vn)
      {
        mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
        un -= 2 * vn;
        up += 2 * vn;
        cy = mpn_add_n (prodp, prodp, ws, vn);
        MPN_COPY (prodp + vn, ws + vn, 2 * vn);
        mpn_incr_u (prodp + vn, cy);
        prodp += 2 * vn;
      }

    /* vn <= un < 3vn */

    if (4 * un < 5 * vn)
      mpn_toom22_mul (ws, up, un, vp, vn, scratch);
    else if (4 * un < 7 * vn)
      mpn_toom32_mul (ws, up, un, vp, vn, scratch);
    else
      mpn_toom42_mul (ws, up, un, vp, vn, scratch);

    cy = mpn_add_n (prodp, prodp, ws, vn);
    MPN_COPY (prodp + vn, ws + vn, un);
    mpn_incr_u (prodp + vn, cy);
  }
      else
  {
    if (4 * un < 5 * vn)
      mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
    else if (4 * un < 7 * vn)
      mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
    else
      mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
  }
      TMP_SFREE;
    }
  else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) ||
     BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
    {
      /* Handle the largest operands that are not in the FFT range.  The 2nd
   condition makes very unbalanced operands avoid the FFT code (except
   perhaps as coefficient products of the Toom code.  */

      if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) || !TOOM44_OK (un, vn))
  {
    /* Use ToomX3 variants */
    mp_ptr scratch;
    TMP_DECL; TMP_MARK;

#define ITCH_TOOMX3 (4 * vn + GMP_NUMB_BITS)
    scratch = TMP_ALLOC_LIMBS (ITCH_TOOMX3);
    ASSERT (mpn_toom33_mul_itch ((7*vn-1)/6, vn) <= ITCH_TOOMX3); /* 7vn/2+ */
    ASSERT (mpn_toom43_mul_itch ((3*vn-1)/2, vn) <= ITCH_TOOMX3); /* 9vn/4+ */
    ASSERT (mpn_toom32_mul_itch ((7*vn-1)/4, vn) <= ITCH_TOOMX3); /* 7vn/6+ */
    ASSERT (mpn_toom53_mul_itch ((11*vn-1)/6, vn) <= ITCH_TOOMX3); /* 11vn/3+ */
    ASSERT (mpn_toom42_mul_itch ((5*vn-1)/2, vn) <= ITCH_TOOMX3); /* 15vn/4+ */
    ASSERT (mpn_toom63_mul_itch ((5*vn-1)/2, vn) <= ITCH_TOOMX3); /* 15vn/4+ */
#undef ITCH_TOOMX3

    if (2 * un >= 5 * vn)
      {
        mp_limb_t cy;
        mp_ptr ws;

        /* The maximum ws usage is for the mpn_mul result.  */
        ws = TMP_ALLOC_LIMBS (7 * vn >> 1);

        if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
    mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
        else
    mpn_toom63_mul (prodp, up, 2 * vn, vp, vn, scratch);
        un -= 2 * vn;
        up += 2 * vn;
        prodp += 2 * vn;

        while (2 * un >= 5 * vn) /* un >= 2.5vn */
    {
      if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
        mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
      else
        mpn_toom63_mul (ws, up, 2 * vn, vp, vn, scratch);
      un -= 2 * vn;
      up += 2 * vn;
      cy = mpn_add_n (prodp, prodp, ws, vn);
      MPN_COPY (prodp + vn, ws + vn, 2 * vn);
      mpn_incr_u (prodp + vn, cy);
      prodp += 2 * vn;
    }

        /* vn / 2 <= un < 2.5vn */

        if (un < vn)
    mpn_mul (ws, vp, vn, up, un);
        else
    mpn_mul (ws, up, un, vp, vn);

        cy = mpn_add_n (prodp, prodp, ws, vn);
        MPN_COPY (prodp + vn, ws + vn, un);
        mpn_incr_u (prodp + vn, cy);
      }
    else
      {
        if (6 * un < 7 * vn)
    mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
        else if (2 * un < 3 * vn)
    {
      if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
        mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
      else
        mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
    }
        else if (6 * un < 11 * vn)
    {
      if (4 * un < 7 * vn)
        {
          if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
      mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
          else
      mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
        }
      else
        {
          if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM53_THRESHOLD))
      mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
          else
      mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
        }
    }
        else
    {
      if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
        mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
      else
        mpn_toom63_mul (prodp, up, un, vp, vn, scratch);
    }
      }
    TMP_FREE;
  }
      else
  {
    mp_ptr scratch;
    TMP_DECL; TMP_MARK;

    if (BELOW_THRESHOLD (vn, MUL_TOOM6H_THRESHOLD))
      {
        scratch = TMP_SALLOC_LIMBS (mpn_toom44_mul_itch (un, vn));
        mpn_toom44_mul (prodp, up, un, vp, vn, scratch);
      }
    else if (BELOW_THRESHOLD (vn, MUL_TOOM8H_THRESHOLD))
      {
        scratch = TMP_SALLOC_LIMBS (mpn_toom6h_mul_itch (un, vn));
        mpn_toom6h_mul (prodp, up, un, vp, vn, scratch);
      }
    else
      {
        scratch = TMP_ALLOC_LIMBS (mpn_toom8h_mul_itch (un, vn));
        mpn_toom8h_mul (prodp, up, un, vp, vn, scratch);
      }
    TMP_FREE;
  }
    }
  else
    {
      if (un >= 8 * vn)
  {
    mp_limb_t cy;
    mp_ptr ws;
    TMP_DECL; TMP_MARK;

    /* The maximum ws usage is for the mpn_mul result.  */
    ws = TMP_BALLOC_LIMBS (9 * vn >> 1);

    mpn_fft_mul (prodp, up, 3 * vn, vp, vn);
    un -= 3 * vn;
    up += 3 * vn;
    prodp += 3 * vn;

    while (2 * un >= 7 * vn) /* un >= 3.5vn  */
      {
        mpn_fft_mul (ws, up, 3 * vn, vp, vn);
        un -= 3 * vn;
        up += 3 * vn;
        cy = mpn_add_n (prodp, prodp, ws, vn);
        MPN_COPY (prodp + vn, ws + vn, 3 * vn);
        mpn_incr_u (prodp + vn, cy);
        prodp += 3 * vn;
      }

    /* vn / 2 <= un < 3.5vn */

    if (un < vn)
      mpn_mul (ws, vp, vn, up, un);
    else
      mpn_mul (ws, up, un, vp, vn);

    cy = mpn_add_n (prodp, prodp, ws, vn);
    MPN_COPY (prodp + vn, ws + vn, un);
    mpn_incr_u (prodp + vn, cy);

    TMP_FREE;
  }
      else
  mpn_fft_mul (prodp, up, un, vp, vn);
    }

  return prodp[un + vn - 1];  /* historic */
}

Coverage Report

Created: 2025-12-31 06:37

Line	Count	Source
1		/* mpn_mul -- Multiply two natural numbers.
2
3		Contributed to the GNU project by Torbjorn Granlund.
4
5		Copyright 1991, 1993, 1994, 1996, 1997, 1999-2003, 2005-2007, 2009, 2010, 2012,
6		2014, 2019 Free Software Foundation, Inc.
7
8		This file is part of the GNU MP Library.
9
10		The GNU MP Library is free software; you can redistribute it and/or modify
11		it under the terms of either:
12
13		* the GNU Lesser General Public License as published by the Free
14		Software Foundation; either version 3 of the License, or (at your
15		option) any later version.
16
17		or
18
19		* the GNU General Public License as published by the Free Software
20		Foundation; either version 2 of the License, or (at your option) any
21		later version.
22
23		or both in parallel, as here.
24
25		The GNU MP Library is distributed in the hope that it will be useful, but
26		WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27		or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
28		for more details.
29
30		You should have received copies of the GNU General Public License and the
31		GNU Lesser General Public License along with the GNU MP Library. If not,
32		see https://www.gnu.org/licenses/. */
33
34		#include "gmp-impl.h"
35
36
37		#ifndef MUL_BASECASE_MAX_UN
38	29.0k	#define MUL_BASECASE_MAX_UN 500
39		#endif
40
41		/* Areas where the different toom algorithms can be called (extracted
42		from the t-toom*.c files, and ignoring small constant offsets):
43
44		1/6 1/5 1/4 4/13 1/3 3/8 2/5 5/11 1/2 3/5 2/3 3/4 4/5 1 vn/un
45		4/7 6/7
46		6/11
47		\|--------------------\| toom22 (small)
48		\|\| toom22 (large)
49		\|xxxx\| toom22 called
50		\|-------------------------------------\| toom32
51		\|xxxxxxxxxxxxxxxx\| \| toom32 called
52		\|------------\| toom33
53		\|x\| toom33 called
54		\|---------------------------------\| \| toom42
55		\|xxxxxxxxxxxxxxxxxxxxxxxx\| \| toom42 called
56		\|--------------------\| toom43
57		\|xxxxxxxxxx\| toom43 called
58		\|-----------------------------\| toom52 (unused)
59		\|--------\| toom44
60		\|xxxxxxxx\| toom44 called
61		\|--------------------\| \| toom53
62		\|xxxxxx\| toom53 called
63		\|-------------------------\| toom62 (unused)
64		\|----------------\| toom54 (unused)
65		\|--------------------\| toom63
66		\|xxxxxxxxx\| \| toom63 called
67		\|---------------------------------\| toom6h
68		\|xxxxxxxx\| toom6h called
69		\|-------------------------\| toom8h (32 bit)
70		\|------------------------------------------\| toom8h (64 bit)
71		\|xxxxxxxx\| toom8h called
72		*/
73
74		#define TOOM33_OK(an,bn) (6 + 2 * an < 3 * bn)
75	0	#define TOOM44_OK(an,bn) (12 + 3 * an < 4 * bn)
76
77		/* Multiply the natural numbers u (pointed to by UP, with UN limbs) and v
78		(pointed to by VP, with VN limbs), and store the result at PRODP. The
79		result is UN + VN limbs. Return the most significant limb of the result.
80
81		NOTE: The space pointed to by PRODP is overwritten before finished with U
82		and V, so overlap is an error.
83
84		Argument constraints:
85		1. UN >= VN.
86		2. PRODP != UP and PRODP != VP, i.e. the destination must be distinct from
87		the multiplier and the multiplicand. */
88
89		/*
90		* The cutoff lines in the toomX2 and toomX3 code are now exactly between the
91		ideal lines of the surrounding algorithms. Is that optimal?
92
93		* The toomX3 code now uses a structure similar to the one of toomX2, except
94		that it loops longer in the unbalanced case. The result is that the
95		remaining area might have un < vn. Should we fix the toomX2 code in a
96		similar way?
97
98		* The toomX3 code is used for the largest non-FFT unbalanced operands. It
99		therefore calls mpn_mul recursively for certain cases.
100
101		* Allocate static temp space using THRESHOLD variables (except for toom44
102		when !WANT_FFT). That way, we can typically have no TMP_ALLOC at all.
103
104		* We sort ToomX2 algorithms together, assuming the toom22, toom32, toom42
105		have the same vn threshold. This is not true, we should actually use
106		mul_basecase for slightly larger operands for toom32 than for toom22, and
107		even larger for toom42.
108
109		* That problem is even more prevalent for toomX3. We therefore use special
110		THRESHOLD variables there.
111		*/
112
113		mp_limb_t
114		mpn_mul (mp_ptr prodp,
115		mp_srcptr up, mp_size_t un,
116		mp_srcptr vp, mp_size_t vn)
117	129k	{
118	129k	ASSERT (un >= vn);
119	129k	ASSERT (vn >= 1);
120	129k	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
121	129k	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));
122
123	129k	if (BELOW_THRESHOLD (un, MUL_TOOM22_THRESHOLD))
124	95.0k	{
125		/* When un (and thus vn) is below the toom22 range, do mul_basecase.
126		Test un and not vn here not to thwart the un >> vn code below.
127		This special case is not necessary, but cuts the overhead for the
128		smallest operands. */
129	95.0k	mpn_mul_basecase (prodp, up, un, vp, vn);
130	95.0k	}
131	34.2k	else if (un == vn)
132	4.68k	{
133	4.68k	mpn_mul_n (prodp, up, vp, un);
134	4.68k	}
135	29.5k	else if (vn < MUL_TOOM22_THRESHOLD)
136	14.4k	{ /* plain schoolbook multiplication */
137
138		/* Unless un is very large, or else if have an applicable mpn_mul_N,
139		perform basecase multiply directly. */
140	14.4k	if (un <= MUL_BASECASE_MAX_UN
141		#if HAVE_NATIVE_mpn_mul_2
142		\|\| vn <= 2
143		#else
144	31	\|\| vn == 1
145	14.4k	#endif
146	14.4k	)
147	14.4k	mpn_mul_basecase (prodp, up, un, vp, vn);
148	31	else
149	31	{
150		/* We have un >> MUL_BASECASE_MAX_UN > vn. For better memory
151		locality, split up[] into MUL_BASECASE_MAX_UN pieces and multiply
152		these pieces with the vp[] operand. After each such partial
153		multiplication (but the last) we copy the most significant vn
154		limbs into a temporary buffer since that part would otherwise be
155		overwritten by the next multiplication. After the next
156		multiplication, we add it back. This illustrates the situation:
157
158		-->vn<--
159		\| \|<------- un ------->\|
160		_____________________\|
161		X /\|
162		/XX__________________/ \|
163		_____________________ \|
164		X / \|
165		/XX__________________/ \|
166		_____________________ \|
167		/ / \|
168		/____________________/ \|
169		==================================================================
170
171		The parts marked with X are the parts whose sums are copied into
172		the temporary buffer. */
173
174	31	mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
175	31	mp_limb_t cy;
176	31	ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);
177
178	31	mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
179	31	prodp += MUL_BASECASE_MAX_UN;
180	31	MPN_COPY (tp, prodp, vn); /* preserve high triangle */
181	31	up += MUL_BASECASE_MAX_UN;
182	31	un -= MUL_BASECASE_MAX_UN;
183	31	while (un > MUL_BASECASE_MAX_UN)
184	0	{
185	0	mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
186	0	cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
187	0	mpn_incr_u (prodp + vn, cy);
188	0	prodp += MUL_BASECASE_MAX_UN;
189	0	MPN_COPY (tp, prodp, vn); /* preserve high triangle */
190	0	up += MUL_BASECASE_MAX_UN;
191	0	un -= MUL_BASECASE_MAX_UN;
192	0	}
193	31	if (un > vn)
194	23	{
195	23	mpn_mul_basecase (prodp, up, un, vp, vn);
196	23	}
197	8	else
198	8	{
199	8	ASSERT (un > 0);
200	8	mpn_mul_basecase (prodp, vp, vn, up, un);
201	8	}
202	31	cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
203	31	mpn_incr_u (prodp + vn, cy);
204	31	}
205	14.4k	}
206	15.1k	else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
207	14.5k	{
208		/* Use ToomX2 variants */
209	14.5k	mp_ptr scratch;
210	14.5k	TMP_SDECL; TMP_SMARK;
211
212	14.5k	#define ITCH_TOOMX2 (9 * vn / 2 + GMP_NUMB_BITS * 2)
213	14.5k	scratch = TMP_SALLOC_LIMBS (ITCH_TOOMX2);
214	14.5k	ASSERT (mpn_toom22_mul_itch ((5vn-1)/4, vn) <= ITCH_TOOMX2); / 5vn/2+ */
215	14.5k	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX2); / 7vn/6+ */
216	14.5k	ASSERT (mpn_toom42_mul_itch (3 * vn - 1, vn) <= ITCH_TOOMX2); /* 9vn/2+ */
217	14.5k	#undef ITCH_TOOMX2
218
219		/* FIXME: This condition (repeated in the loop below) leaves from a vn*vn
220		square to a (3vn-1)*vn rectangle. Leaving such a rectangle is hardly
221		wise; we would get better balance by slightly moving the bound. We
222		will sometimes end up with un < vn, like in the X3 arm below. */
223	14.5k	if (un >= 3 * vn)
224	1.00k	{
225	1.00k	mp_limb_t cy;
226	1.00k	mp_ptr ws;
227
228		/* The maximum ws usage is for the mpn_mul result. */
229	1.00k	ws = TMP_SALLOC_LIMBS (4 * vn);
230
231	1.00k	mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
232	1.00k	un -= 2 * vn;
233	1.00k	up += 2 * vn;
234	1.00k	prodp += 2 * vn;
235
236	1.37k	while (un >= 3 * vn)
237	377	{
238	377	mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
239	377	un -= 2 * vn;
240	377	up += 2 * vn;
241	377	cy = mpn_add_n (prodp, prodp, ws, vn);
242	377	MPN_COPY (prodp + vn, ws + vn, 2 * vn);
243	377	mpn_incr_u (prodp + vn, cy);
244	377	prodp += 2 * vn;
245	377	}
246
247		/* vn <= un < 3vn */
248
249	1.00k	if (4 * un < 5 * vn)
250	386	mpn_toom22_mul (ws, up, un, vp, vn, scratch);
251	614	else if (4 * un < 7 * vn)
252	313	mpn_toom32_mul (ws, up, un, vp, vn, scratch);
253	301	else
254	301	mpn_toom42_mul (ws, up, un, vp, vn, scratch);
255
256	1.00k	cy = mpn_add_n (prodp, prodp, ws, vn);
257	1.00k	MPN_COPY (prodp + vn, ws + vn, un);
258	1.00k	mpn_incr_u (prodp + vn, cy);
259	1.00k	}
260	13.5k	else
261	13.5k	{
262	13.5k	if (4 * un < 5 * vn)
263	6.27k	mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
264	7.25k	else if (4 * un < 7 * vn)
265	5.22k	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
266	2.02k	else
267	2.02k	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
268	13.5k	}
269	14.5k	TMP_SFREE;
270	14.5k	}
271	614	else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) \|\|
272	0	BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
273	614	{
274		/* Handle the largest operands that are not in the FFT range. The 2nd
275		condition makes very unbalanced operands avoid the FFT code (except
276		perhaps as coefficient products of the Toom code. */
277
278	614	if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) \|\| !TOOM44_OK (un, vn))
279	614	{
280		/* Use ToomX3 variants */
281	614	mp_ptr scratch;
282	614	TMP_DECL; TMP_MARK;
283
284	614	#define ITCH_TOOMX3 (4 * vn + GMP_NUMB_BITS)
285	614	scratch = TMP_ALLOC_LIMBS (ITCH_TOOMX3);
286	614	ASSERT (mpn_toom33_mul_itch ((7vn-1)/6, vn) <= ITCH_TOOMX3); / 7vn/2+ */
287	614	ASSERT (mpn_toom43_mul_itch ((3vn-1)/2, vn) <= ITCH_TOOMX3); / 9vn/4+ */
288	614	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX3); / 7vn/6+ */
289	614	ASSERT (mpn_toom53_mul_itch ((11vn-1)/6, vn) <= ITCH_TOOMX3); / 11vn/3+ */
290	614	ASSERT (mpn_toom42_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
291	614	ASSERT (mpn_toom63_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
292	614	#undef ITCH_TOOMX3
293
294	614	if (2 * un >= 5 * vn)
295	0	{
296	0	mp_limb_t cy;
297	0	mp_ptr ws;
298
299		/* The maximum ws usage is for the mpn_mul result. */
300	0	ws = TMP_ALLOC_LIMBS (7 * vn >> 1);
301
302	0	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
303	0	mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
304	0	else
305	0	mpn_toom63_mul (prodp, up, 2 * vn, vp, vn, scratch);
306	0	un -= 2 * vn;
307	0	up += 2 * vn;
308	0	prodp += 2 * vn;
309
310	0	while (2 * un >= 5 * vn) /* un >= 2.5vn */
311	0	{
312	0	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
313	0	mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
314	0	else
315	0	mpn_toom63_mul (ws, up, 2 * vn, vp, vn, scratch);
316	0	un -= 2 * vn;
317	0	up += 2 * vn;
318	0	cy = mpn_add_n (prodp, prodp, ws, vn);
319	0	MPN_COPY (prodp + vn, ws + vn, 2 * vn);
320	0	mpn_incr_u (prodp + vn, cy);
321	0	prodp += 2 * vn;
322	0	}
323
324		/* vn / 2 <= un < 2.5vn */
325
326	0	if (un < vn)
327	0	mpn_mul (ws, vp, vn, up, un);
328	0	else
329	0	mpn_mul (ws, up, un, vp, vn);
330
331	0	cy = mpn_add_n (prodp, prodp, ws, vn);
332	0	MPN_COPY (prodp + vn, ws + vn, un);
333	0	mpn_incr_u (prodp + vn, cy);
334	0	}
335	614	else
336	614	{
337	614	if (6 * un < 7 * vn)
338	62	mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
339	552	else if (2 * un < 3 * vn)
340	508	{
341	508	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
342	0	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
343	508	else
344	508	mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
345	508	}
346	44	else if (6 * un < 11 * vn)
347	44	{
348	44	if (4 * un < 7 * vn)
349	44	{
350	44	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
351	44	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
352	0	else
353	0	mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
354	44	}
355	0	else
356	0	{
357	0	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM53_THRESHOLD))
358	0	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
359	0	else
360	0	mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
361	0	}
362	44	}
363	0	else
364	0	{
365	0	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
366	0	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
367	0	else
368	0	mpn_toom63_mul (prodp, up, un, vp, vn, scratch);
369	0	}
370	614	}
371	614	TMP_FREE;
372	614	}
373	0	else
374	0	{
375	0	mp_ptr scratch;
376	0	TMP_DECL; TMP_MARK;
377
378	0	if (BELOW_THRESHOLD (vn, MUL_TOOM6H_THRESHOLD))
379	0	{
380	0	scratch = TMP_SALLOC_LIMBS (mpn_toom44_mul_itch (un, vn));
381	0	mpn_toom44_mul (prodp, up, un, vp, vn, scratch);
382	0	}
383	0	else if (BELOW_THRESHOLD (vn, MUL_TOOM8H_THRESHOLD))
384	0	{
385	0	scratch = TMP_SALLOC_LIMBS (mpn_toom6h_mul_itch (un, vn));
386	0	mpn_toom6h_mul (prodp, up, un, vp, vn, scratch);
387	0	}
388	0	else
389	0	{
390	0	scratch = TMP_ALLOC_LIMBS (mpn_toom8h_mul_itch (un, vn));
391	0	mpn_toom8h_mul (prodp, up, un, vp, vn, scratch);
392	0	}
393	0	TMP_FREE;
394	0	}
395	614	}
396	0	else
397	0	{
398	0	if (un >= 8 * vn)
399	0	{
400	0	mp_limb_t cy;
401	0	mp_ptr ws;
402	0	TMP_DECL; TMP_MARK;
403
404		/* The maximum ws usage is for the mpn_mul result. */
405	0	ws = TMP_BALLOC_LIMBS (9 * vn >> 1);
406
407	0	mpn_fft_mul (prodp, up, 3 * vn, vp, vn);
408	0	un -= 3 * vn;
409	0	up += 3 * vn;
410	0	prodp += 3 * vn;
411
412	0	while (2 * un >= 7 * vn) /* un >= 3.5vn */
413	0	{
414	0	mpn_fft_mul (ws, up, 3 * vn, vp, vn);
415	0	un -= 3 * vn;
416	0	up += 3 * vn;
417	0	cy = mpn_add_n (prodp, prodp, ws, vn);
418	0	MPN_COPY (prodp + vn, ws + vn, 3 * vn);
419	0	mpn_incr_u (prodp + vn, cy);
420	0	prodp += 3 * vn;
421	0	}
422
423		/* vn / 2 <= un < 3.5vn */
424
425	0	if (un < vn)
426	0	mpn_mul (ws, vp, vn, up, un);
427	0	else
428	0	mpn_mul (ws, up, un, vp, vn);
429
430	0	cy = mpn_add_n (prodp, prodp, ws, vn);
431	0	MPN_COPY (prodp + vn, ws + vn, un);
432	0	mpn_incr_u (prodp + vn, cy);
433
434	0	TMP_FREE;
435	0	}
436	0	else
437	0	mpn_fft_mul (prodp, up, un, vp, vn);
438	0	}
439
440	129k	return prodp[un + vn - 1]; /* historic */
441	129k	}