1 /* Implementation of the squaring algorithm with Toom-Cook 8.5-way.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 THE FUNCTION IN THIS FILE IS INTERNAL WITH A MUTABLE INTERFACE. IT IS ONLY
6 SAFE TO REACH IT THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
7 GUARANTEED THAT IT WILL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
8
9 Copyright 2009, 2012 Free Software Foundation, Inc.
10
11 This file is part of the GNU MP Library.
12
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of either:
15
16 * the GNU Lesser General Public License as published by the Free
17 Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
19
20 or
21
22 * the GNU General Public License as published by the Free Software
23 Foundation; either version 2 of the License, or (at your option) any
24 later version.
25
26 or both in parallel, as here.
27
28 The GNU MP Library is distributed in the hope that it will be useful, but
29 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
30 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
31 for more details.
32
33 You should have received copies of the GNU General Public License and the
34 GNU Lesser General Public License along with the GNU MP Library. If not,
35 see https://www.gnu.org/licenses/. */
36
37
38 #include "gmp-impl.h"
39
40 #if GMP_NUMB_BITS < 29
41 #error Not implemented.
42 #endif
43
44 #if GMP_NUMB_BITS < 43
45 #define BIT_CORRECTION 1
46 #define CORRECTION_BITS GMP_NUMB_BITS
47 #else
48 #define BIT_CORRECTION 0
49 #define CORRECTION_BITS 0
50 #endif
51
52 #ifndef SQR_TOOM8_THRESHOLD
53 #define SQR_TOOM8_THRESHOLD MUL_TOOM8H_THRESHOLD
54 #endif
55
56 #ifndef SQR_TOOM6_THRESHOLD
57 #define SQR_TOOM6_THRESHOLD MUL_TOOM6H_THRESHOLD
58 #endif
59
60 #if TUNE_PROGRAM_BUILD
61 #define MAYBE_sqr_basecase 1
62 #define MAYBE_sqr_above_basecase 1
63 #define MAYBE_sqr_toom2 1
64 #define MAYBE_sqr_above_toom2 1
65 #define MAYBE_sqr_toom3 1
66 #define MAYBE_sqr_above_toom3 1
67 #define MAYBE_sqr_toom4 1
68 #define MAYBE_sqr_above_toom4 1
69 #define MAYBE_sqr_above_toom6 1
70 #else
71 #define SQR_TOOM8_MAX \
72 ((SQR_FFT_THRESHOLD <= MP_SIZE_T_MAX - (8*2-1+7)) ? \
73 ((SQR_FFT_THRESHOLD+8*2-1+7)/8) \
74 : MP_SIZE_T_MAX )
75 #define MAYBE_sqr_basecase \
76 (SQR_TOOM8_THRESHOLD < 8 * SQR_TOOM2_THRESHOLD)
77 #define MAYBE_sqr_above_basecase \
78 (SQR_TOOM8_MAX >= SQR_TOOM2_THRESHOLD)
79 #define MAYBE_sqr_toom2 \
80 (SQR_TOOM8_THRESHOLD < 8 * SQR_TOOM3_THRESHOLD)
81 #define MAYBE_sqr_above_toom2 \
82 (SQR_TOOM8_MAX >= SQR_TOOM3_THRESHOLD)
83 #define MAYBE_sqr_toom3 \
84 (SQR_TOOM8_THRESHOLD < 8 * SQR_TOOM4_THRESHOLD)
85 #define MAYBE_sqr_above_toom3 \
86 (SQR_TOOM8_MAX >= SQR_TOOM4_THRESHOLD)
87 #define MAYBE_sqr_toom4 \
88 (SQR_TOOM8_THRESHOLD < 8 * SQR_TOOM6_THRESHOLD)
89 #define MAYBE_sqr_above_toom4 \
90 (SQR_TOOM8_MAX >= SQR_TOOM6_THRESHOLD)
91 #define MAYBE_sqr_above_toom6 \
92 (SQR_TOOM8_MAX >= SQR_TOOM8_THRESHOLD)
93 #endif
94
95 #define TOOM8_SQR_REC(p, a, f, p2, a2, n, ws) \
96 do { \
97 if (MAYBE_sqr_basecase && ( !MAYBE_sqr_above_basecase \
98 || BELOW_THRESHOLD (n, SQR_TOOM2_THRESHOLD))) { \
99 mpn_sqr_basecase (p, a, n); \
100 if (f) mpn_sqr_basecase (p2, a2, n); \
101 } else if (MAYBE_sqr_toom2 && ( !MAYBE_sqr_above_toom2 \
102 || BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD))) { \
103 mpn_toom2_sqr (p, a, n, ws); \
104 if (f) mpn_toom2_sqr (p2, a2, n, ws); \
105 } else if (MAYBE_sqr_toom3 && ( !MAYBE_sqr_above_toom3 \
106 || BELOW_THRESHOLD (n, SQR_TOOM4_THRESHOLD))) { \
107 mpn_toom3_sqr (p, a, n, ws); \
108 if (f) mpn_toom3_sqr (p2, a2, n, ws); \
109 } else if (MAYBE_sqr_toom4 && ( !MAYBE_sqr_above_toom4 \
110 || BELOW_THRESHOLD (n, SQR_TOOM6_THRESHOLD))) { \
111 mpn_toom4_sqr (p, a, n, ws); \
112 if (f) mpn_toom4_sqr (p2, a2, n, ws); \
113 } else if (! MAYBE_sqr_above_toom6 \
114 || BELOW_THRESHOLD (n, SQR_TOOM8_THRESHOLD)) { \
115 mpn_toom6_sqr (p, a, n, ws); \
116 if (f) mpn_toom6_sqr (p2, a2, n, ws); \
117 } else { \
118 mpn_toom8_sqr (p, a, n, ws); \
119 if (f) mpn_toom8_sqr (p2, a2, n, ws); \
120 } \
121 } while (0)
122
123 void
124 mpn_toom8_sqr (mp_ptr pp, mp_srcptr ap, mp_size_t an, mp_ptr scratch)
125 {
126 mp_size_t n, s;
127
128 /***************************** decomposition *******************************/
129
130 ASSERT ( an >= 40 );
131
132 n = 1 + ((an - 1)>>3);
133
134 s = an - 7 * n;
135
136 ASSERT (0 < s && s <= n);
137 ASSERT ( s + s > 3 );
138
139 #define r6 (pp + 3 * n) /* 3n+1 */
140 #define r4 (pp + 7 * n) /* 3n+1 */
141 #define r2 (pp +11 * n) /* 3n+1 */
142 #define r0 (pp +15 * n) /* s+t <= 2*n */
143 #define r7 (scratch) /* 3n+1 */
144 #define r5 (scratch + 3 * n + 1) /* 3n+1 */
145 #define r3 (scratch + 6 * n + 2) /* 3n+1 */
146 #define r1 (scratch + 9 * n + 3) /* 3n+1 */
147 #define v0 (pp +11 * n) /* n+1 */
148 #define v2 (pp +13 * n+2) /* n+1 */
149 #define wse (scratch +12 * n + 4) /* 3n+1 */
150
151 /* Alloc also 3n+1 limbs for ws... toom_interpolate_16pts may
152 need all of them, when DO_mpn_sublsh_n usea a scratch */
153 /* if (scratch == NULL) */
154 /* scratch = TMP_SALLOC_LIMBS (30 * n + 6); */
155
156 /********************** evaluation and recursive calls *********************/
157 /* $\pm1/8$ */
158 mpn_toom_eval_pm2rexp (v2, v0, 7, ap, n, s, 3, pp);
159 /* A(-1/8)*B(-1/8)*8^. */ /* A(+1/8)*B(+1/8)*8^. */
160 TOOM8_SQR_REC(pp, v0, 2, r7, v2, n + 1, wse);
161 mpn_toom_couple_handling (r7, 2 * n + 1 + BIT_CORRECTION, pp, 0, n, 3, 0);
162
163 /* $\pm1/4$ */
164 mpn_toom_eval_pm2rexp (v2, v0, 7, ap, n, s, 2, pp);
165 /* A(-1/4)*B(-1/4)*4^. */ /* A(+1/4)*B(+1/4)*4^. */
166 TOOM8_SQR_REC(pp, v0, 2, r5, v2, n + 1, wse);
167 mpn_toom_couple_handling (r5, 2 * n + 1, pp, 0, n, 2, 0);
168
169 /* $\pm2$ */
170 mpn_toom_eval_pm2 (v2, v0, 7, ap, n, s, pp);
171 /* A(-2)*B(-2) */ /* A(+2)*B(+2) */
172 TOOM8_SQR_REC(pp, v0, 2, r3, v2, n + 1, wse);
173 mpn_toom_couple_handling (r3, 2 * n + 1, pp, 0, n, 1, 2);
174
175 /* $\pm8$ */
176 mpn_toom_eval_pm2exp (v2, v0, 7, ap, n, s, 3, pp);
177 /* A(-8)*B(-8) */ /* A(+8)*B(+8) */
178 TOOM8_SQR_REC(pp, v0, 2, r1, v2, n + 1, wse);
179 mpn_toom_couple_handling (r1, 2 * n + 1 + BIT_CORRECTION, pp, 0, n, 3, 6);
180
181 /* $\pm1/2$ */
182 mpn_toom_eval_pm2rexp (v2, v0, 7, ap, n, s, 1, pp);
183 /* A(-1/2)*B(-1/2)*2^. */ /* A(+1/2)*B(+1/2)*2^. */
184 TOOM8_SQR_REC(pp, v0, 2, r6, v2, n + 1, wse);
185 mpn_toom_couple_handling (r6, 2 * n + 1, pp, 0, n, 1, 0);
186
187 /* $\pm1$ */
188 mpn_toom_eval_pm1 (v2, v0, 7, ap, n, s, pp);
189 /* A(-1)*B(-1) */ /* A(1)*B(1) */
190 TOOM8_SQR_REC(pp, v0, 2, r4, v2, n + 1, wse);
191 mpn_toom_couple_handling (r4, 2 * n + 1, pp, 0, n, 0, 0);
192
193 /* $\pm4$ */
194 mpn_toom_eval_pm2exp (v2, v0, 7, ap, n, s, 2, pp);
195 /* A(-4)*B(-4) */ /* A(+4)*B(+4) */
196 TOOM8_SQR_REC(pp, v0, 2, r2, v2, n + 1, wse);
197 mpn_toom_couple_handling (r2, 2 * n + 1, pp, 0, n, 2, 4);
198
199 #undef v0
200 #undef v2
201
202 /* A(0)*B(0) */
203 TOOM8_SQR_REC(pp, ap, 0, pp, ap, n, wse);
204
205 mpn_toom_interpolate_16pts (pp, r1, r3, r5, r7, n, 2 * s, 0, wse);
206
207 #undef r0
208 #undef r1
209 #undef r2
210 #undef r3
211 #undef r4
212 #undef r5
213 #undef r6
214 #undef wse
215
216 }
217
218 #undef TOOM8_SQR_REC
219 #undef MAYBE_sqr_basecase
220 #undef MAYBE_sqr_above_basecase
221 #undef MAYBE_sqr_toom2
222 #undef MAYBE_sqr_above_toom2
223 #undef MAYBE_sqr_toom3
224 #undef MAYBE_sqr_above_toom3
225 #undef MAYBE_sqr_above_toom4