/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-02-14 06:49

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