/src/gmp-6.2.1/mpn/mul.c

Source (jump to first uncovered line)
/* 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-03-09 06:52

Line	Count	Source (jump to first uncovered line)
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	535k	#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	18.3k	#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	2.64M	{
118	2.64M	ASSERT (un >= vn);
119	2.64M	ASSERT (vn >= 1);
120	2.64M	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, up, un));
121	2.64M	ASSERT (! MPN_OVERLAP_P (prodp, un+vn, vp, vn));
122
123	2.64M	if (BELOW_THRESHOLD (un, MUL_TOOM22_THRESHOLD))
124	1.51M	{
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	1.51M	mpn_mul_basecase (prodp, up, un, vp, vn);
130	1.51M	}
131	1.13M	else if (un == vn)
132	467k	{
133	467k	mpn_mul_n (prodp, up, vp, un);
134	467k	}
135	662k	else if (vn < MUL_TOOM22_THRESHOLD)
136	267k	{ /* 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	267k	if (un <= MUL_BASECASE_MAX_UN
141	267k	#if HAVE_NATIVE_mpn_mul_2
142	267k	\|\| vn <= 2
143		#else
144		\|\| vn == 1
145		#endif
146	267k	)
147	267k	mpn_mul_basecase (prodp, up, un, vp, vn);
148	0	else
149	0	{
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	0	mp_limb_t tp[MUL_TOOM22_THRESHOLD_LIMIT];
175	0	mp_limb_t cy;
176	0	ASSERT (MUL_TOOM22_THRESHOLD <= MUL_TOOM22_THRESHOLD_LIMIT);
177
178	0	mpn_mul_basecase (prodp, up, MUL_BASECASE_MAX_UN, vp, vn);
179	0	prodp += MUL_BASECASE_MAX_UN;
180	0	MPN_COPY (tp, prodp, vn); /* preserve high triangle */
181	0	up += MUL_BASECASE_MAX_UN;
182	0	un -= MUL_BASECASE_MAX_UN;
183	0	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	0	if (un > vn)
194	0	{
195	0	mpn_mul_basecase (prodp, up, un, vp, vn);
196	0	}
197	0	else
198	0	{
199	0	ASSERT (un > 0);
200	0	mpn_mul_basecase (prodp, vp, vn, up, un);
201	0	}
202	0	cy = mpn_add_n (prodp, prodp, tp, vn); /* add back preserved triangle */
203	0	mpn_incr_u (prodp + vn, cy);
204	0	}
205	267k	}
206	394k	else if (BELOW_THRESHOLD (vn, MUL_TOOM33_THRESHOLD))
207	333k	{
208		/* Use ToomX2 variants */
209	333k	mp_ptr scratch;
210	333k	TMP_SDECL; TMP_SMARK;
211
212	333k	#define ITCH_TOOMX2 (9 * vn / 2 + GMP_NUMB_BITS * 2)
213	333k	scratch = TMP_SALLOC_LIMBS (ITCH_TOOMX2);
214	333k	ASSERT (mpn_toom22_mul_itch ((5vn-1)/4, vn) <= ITCH_TOOMX2); / 5vn/2+ */
215	333k	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX2); / 7vn/6+ */
216	333k	ASSERT (mpn_toom42_mul_itch (3 * vn - 1, vn) <= ITCH_TOOMX2); /* 9vn/2+ */
217	333k	#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	333k	if (un >= 3 * vn)
224	1.38k	{
225	1.38k	mp_limb_t cy;
226	1.38k	mp_ptr ws;
227
228		/* The maximum ws usage is for the mpn_mul result. */
229	1.38k	ws = TMP_SALLOC_LIMBS (4 * vn);
230
231	1.38k	mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
232	1.38k	un -= 2 * vn;
233	1.38k	up += 2 * vn;
234	1.38k	prodp += 2 * vn;
235
236	3.79k	while (un >= 3 * vn)
237	2.41k	{
238	2.41k	mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
239	2.41k	un -= 2 * vn;
240	2.41k	up += 2 * vn;
241	2.41k	cy = mpn_add_n (prodp, prodp, ws, vn);
242	2.41k	MPN_COPY (prodp + vn, ws + vn, 2 * vn);
243	2.41k	mpn_incr_u (prodp + vn, cy);
244	2.41k	prodp += 2 * vn;
245	2.41k	}
246
247		/* vn <= un < 3vn */
248
249	1.38k	if (4 * un < 5 * vn)
250	83	mpn_toom22_mul (ws, up, un, vp, vn, scratch);
251	1.30k	else if (4 * un < 7 * vn)
252	403	mpn_toom32_mul (ws, up, un, vp, vn, scratch);
253	898	else
254	898	mpn_toom42_mul (ws, up, un, vp, vn, scratch);
255
256	1.38k	cy = mpn_add_n (prodp, prodp, ws, vn);
257	1.38k	MPN_COPY (prodp + vn, ws + vn, un);
258	1.38k	mpn_incr_u (prodp + vn, cy);
259	1.38k	}
260	332k	else
261	332k	{
262	332k	if (4 * un < 5 * vn)
263	323k	mpn_toom22_mul (prodp, up, un, vp, vn, scratch);
264	8.74k	else if (4 * un < 7 * vn)
265	8.31k	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
266	433	else
267	433	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
268	332k	}
269	333k	TMP_SFREE;
270	333k	}
271	61.4k	else if (BELOW_THRESHOLD ((un + vn) >> 1, MUL_FFT_THRESHOLD) \|\|
272	61.4k	BELOW_THRESHOLD (3 * vn, MUL_FFT_THRESHOLD))
273	61.4k	{
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	61.4k	if (BELOW_THRESHOLD (vn, MUL_TOOM44_THRESHOLD) \|\| !TOOM44_OK (un, vn))
279	46.0k	{
280		/* Use ToomX3 variants */
281	46.0k	mp_ptr scratch;
282	46.0k	TMP_DECL; TMP_MARK;
283
284	46.0k	#define ITCH_TOOMX3 (4 * vn + GMP_NUMB_BITS)
285	46.0k	scratch = TMP_ALLOC_LIMBS (ITCH_TOOMX3);
286	46.0k	ASSERT (mpn_toom33_mul_itch ((7vn-1)/6, vn) <= ITCH_TOOMX3); / 7vn/2+ */
287	46.0k	ASSERT (mpn_toom43_mul_itch ((3vn-1)/2, vn) <= ITCH_TOOMX3); / 9vn/4+ */
288	46.0k	ASSERT (mpn_toom32_mul_itch ((7vn-1)/4, vn) <= ITCH_TOOMX3); / 7vn/6+ */
289	46.0k	ASSERT (mpn_toom53_mul_itch ((11vn-1)/6, vn) <= ITCH_TOOMX3); / 11vn/3+ */
290	46.0k	ASSERT (mpn_toom42_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
291	46.0k	ASSERT (mpn_toom63_mul_itch ((5vn-1)/2, vn) <= ITCH_TOOMX3); / 15vn/4+ */
292	46.0k	#undef ITCH_TOOMX3
293
294	46.0k	if (2 * un >= 5 * vn)
295	2.89k	{
296	2.89k	mp_limb_t cy;
297	2.89k	mp_ptr ws;
298
299		/* The maximum ws usage is for the mpn_mul result. */
300	2.89k	ws = TMP_ALLOC_LIMBS (7 * vn >> 1);
301
302	2.89k	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
303	749	mpn_toom42_mul (prodp, up, 2 * vn, vp, vn, scratch);
304	2.14k	else
305	2.14k	mpn_toom63_mul (prodp, up, 2 * vn, vp, vn, scratch);
306	2.89k	un -= 2 * vn;
307	2.89k	up += 2 * vn;
308	2.89k	prodp += 2 * vn;
309
310	10.9k	while (2 * un >= 5 * vn) /* un >= 2.5vn */
311	8.04k	{
312	8.04k	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
313	2.23k	mpn_toom42_mul (ws, up, 2 * vn, vp, vn, scratch);
314	5.80k	else
315	5.80k	mpn_toom63_mul (ws, up, 2 * vn, vp, vn, scratch);
316	8.04k	un -= 2 * vn;
317	8.04k	up += 2 * vn;
318	8.04k	cy = mpn_add_n (prodp, prodp, ws, vn);
319	8.04k	MPN_COPY (prodp + vn, ws + vn, 2 * vn);
320	8.04k	mpn_incr_u (prodp + vn, cy);
321	8.04k	prodp += 2 * vn;
322	8.04k	}
323
324		/* vn / 2 <= un < 2.5vn */
325
326	2.89k	if (un < vn)
327	668	mpn_mul (ws, vp, vn, up, un);
328	2.22k	else
329	2.22k	mpn_mul (ws, up, un, vp, vn);
330
331	2.89k	cy = mpn_add_n (prodp, prodp, ws, vn);
332	2.89k	MPN_COPY (prodp + vn, ws + vn, un);
333	2.89k	mpn_incr_u (prodp + vn, cy);
334	2.89k	}
335	43.1k	else
336	43.1k	{
337	43.1k	if (6 * un < 7 * vn)
338	38.4k	mpn_toom33_mul (prodp, up, un, vp, vn, scratch);
339	4.78k	else if (2 * un < 3 * vn)
340	2.99k	{
341	2.99k	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM43_THRESHOLD))
342	413	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
343	2.58k	else
344	2.58k	mpn_toom43_mul (prodp, up, un, vp, vn, scratch);
345	2.99k	}
346	1.79k	else if (6 * un < 11 * vn)
347	772	{
348	772	if (4 * un < 7 * vn)
349	560	{
350	560	if (BELOW_THRESHOLD (vn, MUL_TOOM32_TO_TOOM53_THRESHOLD))
351	225	mpn_toom32_mul (prodp, up, un, vp, vn, scratch);
352	335	else
353	335	mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
354	560	}
355	212	else
356	212	{
357	212	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM53_THRESHOLD))
358	95	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
359	117	else
360	117	mpn_toom53_mul (prodp, up, un, vp, vn, scratch);
361	212	}
362	772	}
363	1.01k	else
364	1.01k	{
365	1.01k	if (BELOW_THRESHOLD (vn, MUL_TOOM42_TO_TOOM63_THRESHOLD))
366	128	mpn_toom42_mul (prodp, up, un, vp, vn, scratch);
367	890	else
368	890	mpn_toom63_mul (prodp, up, un, vp, vn, scratch);
369	1.01k	}
370	43.1k	}
371	46.0k	TMP_FREE;
372	46.0k	}
373	15.3k	else
374	15.3k	{
375	15.3k	mp_ptr scratch;
376	15.3k	TMP_DECL; TMP_MARK;
377
378	15.3k	if (BELOW_THRESHOLD (vn, MUL_TOOM6H_THRESHOLD))
379	4.49k	{
380	4.49k	scratch = TMP_SALLOC_LIMBS (mpn_toom44_mul_itch (un, vn));
381	4.49k	mpn_toom44_mul (prodp, up, un, vp, vn, scratch);
382	4.49k	}
383	10.8k	else if (BELOW_THRESHOLD (vn, MUL_TOOM8H_THRESHOLD))
384	4.28k	{
385	4.28k	scratch = TMP_SALLOC_LIMBS (mpn_toom6h_mul_itch (un, vn));
386	4.28k	mpn_toom6h_mul (prodp, up, un, vp, vn, scratch);
387	4.28k	}
388	6.60k	else
389	6.60k	{
390	6.60k	scratch = TMP_ALLOC_LIMBS (mpn_toom8h_mul_itch (un, vn));
391	6.60k	mpn_toom8h_mul (prodp, up, un, vp, vn, scratch);
392	6.60k	}
393	15.3k	TMP_FREE;
394	15.3k	}
395	61.4k	}
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	2.64M	return prodp[un + vn - 1]; /* historic */
441	2.64M	}