/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: 2026-05-16 06:52

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	26.5k	#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	114k	{
118	114k	ASSERT (un >= vn);
119	114k	ASSERT (vn >= 1);
120	114k	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
121	114k	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));
122
123	114k	if (BELOW_THRESHOLD (un, MUL_TOOM22_THRESHOLD))
124	83.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	83.0k	mpn_mul_basecase (prodp, up, un, vp, vn);
130	83.0k	}
131	31.0k	else if (un == vn)
132	4.22k	{
133	4.22k	mpn_mul_n (prodp, up, vp, un);
134	4.22k	}
135	26.7k	else if (vn < MUL_TOOM22_THRESHOLD)
136	13.1k	{ /* 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	13.1k	if (un <= MUL_BASECASE_MAX_UN
141		#if HAVE_NATIVE_mpn_mul_2
142		\|\| vn <= 2
143		#else
144	41	\|\| vn == 1
145	13.1k	#endif
146	13.1k	)
147	13.1k	mpn_mul_basecase (prodp, up, un, vp, vn);
148	41	else
149	41	{
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	41	mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
175	41	mp_limb_t cy;
176	41	ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);
177
178	41	mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
179	41	prodp += MUL_BASECASE_MAX_UN;
180	41	MPN_COPY (tp, prodp, vn); /* preserve high triangle */
181	41	up += MUL_BASECASE_MAX_UN;
182	41	un -= MUL_BASECASE_MAX_UN;
183	41	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	41	if (un > vn)
194	35	{
195	35	mpn_mul_basecase (prodp, up, un, vp, vn);
196	35	}
197	6	else
198	6	{
199	6	ASSERT (un > 0);
200	6	mpn_mul_basecase (prodp, vp, vn, up, un);
201	6	}
202	41	cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
203	41	mpn_incr_u (prodp + vn, cy);
204	41	}
205	13.1k	}
206	13.5k	else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
207	13.0k	{
208		/* Use ToomX2 variants */
209	13.0k	mp_ptr scratch;
210	13.0k	TMP_SDECL; TMP_SMARK;
211
212	13.0k	#define ITCH_TOOMX2 (9 * vn / 2 + GMP_NUMB_BITS * 2)
213	13.0k	scratch = TMP_SALLOC_LIMBS (ITCH_TOOMX2);
214	13.0k	ASSERT (mpn_toom22_mul_itch ((5vn-1)/4, vn) <= ITCH_TOOMX2); / 5vn/2+ */
215	13.0k	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX2); / 7vn/6+ */
216	13.0k	ASSERT (mpn_toom42_mul_itch (3 * vn - 1, vn) <= ITCH_TOOMX2); /* 9vn/2+ */
217	13.0k	#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	13.0k	if (un >= 3 * vn)
224	905	{
225	905	mp_limb_t cy;
226	905	mp_ptr ws;
227
228		/* The maximum ws usage is for the mpn_mul result. */
229	905	ws = TMP_SALLOC_LIMBS (4 * vn);
230
231	905	mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
232	905	un -= 2 * vn;
233	905	up += 2 * vn;
234	905	prodp += 2 * vn;
235
236	1.21k	while (un >= 3 * vn)
237	307	{
238	307	mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
239	307	un -= 2 * vn;
240	307	up += 2 * vn;
241	307	cy = mpn_add_n (prodp, prodp, ws, vn);
242	307	MPN_COPY (prodp + vn, ws + vn, 2 * vn);
243	307	mpn_incr_u (prodp + vn, cy);
244	307	prodp += 2 * vn;
245	307	}
246
247		/* vn <= un < 3vn */
248
249	905	if (4 * un < 5 * vn)
250	349	mpn_toom22_mul (ws, up, un, vp, vn, scratch);
251	556	else if (4 * un < 7 * vn)
252	285	mpn_toom32_mul (ws, up, un, vp, vn, scratch);
253	271	else
254	271	mpn_toom42_mul (ws, up, un, vp, vn, scratch);
255
256	905	cy = mpn_add_n (prodp, prodp, ws, vn);
257	905	MPN_COPY (prodp + vn, ws + vn, un);
258	905	mpn_incr_u (prodp + vn, cy);
259	905	}
260	12.1k	else
261	12.1k	{
262	12.1k	if (4 * un < 5 * vn)
263	5.65k	mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
264	6.49k	else if (4 * un < 7 * vn)
265	4.63k	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
266	1.85k	else
267	1.85k	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
268	12.1k	}
269	13.0k	TMP_SFREE;
270	13.0k	}
271	548	else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) \|\|
272	0	BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
273	548	{
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	548	if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) \|\| !TOOM44_OK (un, vn))
279	548	{
280		/* Use ToomX3 variants */
281	548	mp_ptr scratch;
282	548	TMP_DECL; TMP_MARK;
283
284	548	#define ITCH_TOOMX3 (4 * vn + GMP_NUMB_BITS)
285	548	scratch = TMP_ALLOC_LIMBS (ITCH_TOOMX3);
286	548	ASSERT (mpn_toom33_mul_itch ((7vn-1)/6, vn) <= ITCH_TOOMX3); / 7vn/2+ */
287	548	ASSERT (mpn_toom43_mul_itch ((3vn-1)/2, vn) <= ITCH_TOOMX3); / 9vn/4+ */
288	548	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX3); / 7vn/6+ */
289	548	ASSERT (mpn_toom53_mul_itch ((11vn-1)/6, vn) <= ITCH_TOOMX3); / 11vn/3+ */
290	548	ASSERT (mpn_toom42_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
291	548	ASSERT (mpn_toom63_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
292	548	#undef ITCH_TOOMX3
293
294	548	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	548	else
336	548	{
337	548	if (6 * un < 7 * vn)
338	62	mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
339	486	else if (2 * un < 3 * vn)
340	446	{
341	446	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
342	0	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
343	446	else
344	446	mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
345	446	}
346	40	else if (6 * un < 11 * vn)
347	40	{
348	40	if (4 * un < 7 * vn)
349	40	{
350	40	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
351	40	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
352	0	else
353	0	mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
354	40	}
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	40	}
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	548	}
371	548	TMP_FREE;
372	548	}
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	548	}
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	114k	return prodp[un + vn - 1]; /* historic */
441	114k	}