1/* Implementation of the squaring algorithm with Toom-Cook 6.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 9Copyright 2009 Free Software Foundation, Inc. 10 11This file is part of the GNU MP Library. 12 13The GNU MP Library is free software; you can redistribute it and/or modify 14it 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 20or 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 26or both in parallel, as here. 27 28The GNU MP Library is distributed in the hope that it will be useful, but 29WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 30or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 31for more details. 32 33You should have received copies of the GNU General Public License and the 34GNU Lesser General Public License along with the GNU MP Library. If not, 35see https://www.gnu.org/licenses/. */ 36 37 38#include "gmp-impl.h" 39 40 41#if GMP_NUMB_BITS < 21 42#error Not implemented. 43#endif 44 45 46#if TUNE_PROGRAM_BUILD 47#define MAYBE_sqr_basecase 1 48#define MAYBE_sqr_above_basecase 1 49#define MAYBE_sqr_toom2 1 50#define MAYBE_sqr_above_toom2 1 51#define MAYBE_sqr_toom3 1 52#define MAYBE_sqr_above_toom3 1 53#define MAYBE_sqr_above_toom4 1 54#else 55#ifdef SQR_TOOM8_THRESHOLD 56#define SQR_TOOM6_MAX ((SQR_TOOM8_THRESHOLD+6*2-1+5)/6) 57#else 58#define SQR_TOOM6_MAX \ 59 ((SQR_FFT_THRESHOLD <= MP_SIZE_T_MAX - (6*2-1+5)) ? \ 60 ((SQR_FFT_THRESHOLD+6*2-1+5)/6) \ 61 : MP_SIZE_T_MAX ) 62#endif 63#define MAYBE_sqr_basecase \ 64 (SQR_TOOM6_THRESHOLD < 6 * SQR_TOOM2_THRESHOLD) 65#define MAYBE_sqr_above_basecase \ 66 (SQR_TOOM6_MAX >= SQR_TOOM2_THRESHOLD) 67#define MAYBE_sqr_toom2 \ 68 (SQR_TOOM6_THRESHOLD < 6 * SQR_TOOM3_THRESHOLD) 69#define MAYBE_sqr_above_toom2 \ 70 (SQR_TOOM6_MAX >= SQR_TOOM3_THRESHOLD) 71#define MAYBE_sqr_toom3 \ 72 (SQR_TOOM6_THRESHOLD < 6 * SQR_TOOM4_THRESHOLD) 73#define MAYBE_sqr_above_toom3 \ 74 (SQR_TOOM6_MAX >= SQR_TOOM4_THRESHOLD) 75#define MAYBE_sqr_above_toom4 \ 76 (SQR_TOOM6_MAX >= SQR_TOOM6_THRESHOLD) 77#endif 78 79#define TOOM6_SQR_REC(p, a, n, ws) \ 80 do { \ 81 if (MAYBE_sqr_basecase && ( !MAYBE_sqr_above_basecase \ 82 || BELOW_THRESHOLD (n, SQR_TOOM2_THRESHOLD))) \ 83 mpn_sqr_basecase (p, a, n); \ 84 else if (MAYBE_sqr_toom2 && ( !MAYBE_sqr_above_toom2 \ 85 || BELOW_THRESHOLD (n, SQR_TOOM3_THRESHOLD))) \ 86 mpn_toom2_sqr (p, a, n, ws); \ 87 else if (MAYBE_sqr_toom3 && ( !MAYBE_sqr_above_toom3 \ 88 || BELOW_THRESHOLD (n, SQR_TOOM4_THRESHOLD))) \ 89 mpn_toom3_sqr (p, a, n, ws); \ 90 else if (! MAYBE_sqr_above_toom4 \ 91 || BELOW_THRESHOLD (n, SQR_TOOM6_THRESHOLD)) \ 92 mpn_toom4_sqr (p, a, n, ws); \ 93 else \ 94 mpn_toom6_sqr (p, a, n, ws); \ 95 } while (0) 96 97void 98mpn_toom6_sqr (mp_ptr pp, mp_srcptr ap, mp_size_t an, mp_ptr scratch) 99{ 100 mp_size_t n, s; 101 102 /***************************** decomposition *******************************/ 103 104 ASSERT( an >= 18 ); 105 106 n = 1 + (an - 1) / (size_t) 6; 107 108 s = an - 5 * n; 109 110 ASSERT (0 < s && s <= n); 111 112#define r4 (pp + 3 * n) /* 3n+1 */ 113#define r2 (pp + 7 * n) /* 3n+1 */ 114#define r0 (pp +11 * n) /* s+t <= 2*n */ 115#define r5 (scratch) /* 3n+1 */ 116#define r3 (scratch + 3 * n + 1) /* 3n+1 */ 117#define r1 (scratch + 6 * n + 2) /* 3n+1 */ 118#define v0 (pp + 7 * n) /* n+1 */ 119#define v2 (pp + 9 * n+2) /* n+1 */ 120#define wse (scratch + 9 * n + 3) /* 3n+1 */ 121 122 /* Alloc also 3n+1 limbs for ws... toom_interpolate_12pts may 123 need all of them, when DO_mpn_sublsh_n usea a scratch */ 124/* if (scratch== NULL) */ 125/* scratch = TMP_SALLOC_LIMBS (12 * n + 6); */ 126 127 /********************** evaluation and recursive calls *********************/ 128 /* $\pm1/2$ */ 129 mpn_toom_eval_pm2rexp (v2, v0, 5, ap, n, s, 1, pp); 130 TOOM6_SQR_REC(pp, v0, n + 1, wse); /* A(-1/2)*B(-1/2)*2^. */ 131 TOOM6_SQR_REC(r5, v2, n + 1, wse); /* A(+1/2)*B(+1/2)*2^. */ 132 mpn_toom_couple_handling (r5, 2 * n + 1, pp, 0, n, 1, 0); 133 134 /* $\pm1$ */ 135 mpn_toom_eval_pm1 (v2, v0, 5, ap, n, s, pp); 136 TOOM6_SQR_REC(pp, v0, n + 1, wse); /* A(-1)*B(-1) */ 137 TOOM6_SQR_REC(r3, v2, n + 1, wse); /* A(1)*B(1) */ 138 mpn_toom_couple_handling (r3, 2 * n + 1, pp, 0, n, 0, 0); 139 140 /* $\pm4$ */ 141 mpn_toom_eval_pm2exp (v2, v0, 5, ap, n, s, 2, pp); 142 TOOM6_SQR_REC(pp, v0, n + 1, wse); /* A(-4)*B(-4) */ 143 TOOM6_SQR_REC(r1, v2, n + 1, wse); /* A(+4)*B(+4) */ 144 mpn_toom_couple_handling (r1, 2 * n + 1, pp, 0, n, 2, 4); 145 146 /* $\pm1/4$ */ 147 mpn_toom_eval_pm2rexp (v2, v0, 5, ap, n, s, 2, pp); 148 TOOM6_SQR_REC(pp, v0, n + 1, wse); /* A(-1/4)*B(-1/4)*4^. */ 149 TOOM6_SQR_REC(r4, v2, n + 1, wse); /* A(+1/4)*B(+1/4)*4^. */ 150 mpn_toom_couple_handling (r4, 2 * n + 1, pp, 0, n, 2, 0); 151 152 /* $\pm2$ */ 153 mpn_toom_eval_pm2 (v2, v0, 5, ap, n, s, pp); 154 TOOM6_SQR_REC(pp, v0, n + 1, wse); /* A(-2)*B(-2) */ 155 TOOM6_SQR_REC(r2, v2, n + 1, wse); /* A(+2)*B(+2) */ 156 mpn_toom_couple_handling (r2, 2 * n + 1, pp, 0, n, 1, 2); 157 158#undef v0 159#undef v2 160 161 /* A(0)*B(0) */ 162 TOOM6_SQR_REC(pp, ap, n, wse); 163 164 mpn_toom_interpolate_12pts (pp, r1, r3, r5, n, 2 * s, 0, wse); 165 166#undef r0 167#undef r1 168#undef r2 169#undef r3 170#undef r4 171#undef r5 172 173} 174#undef TOOM6_SQR_REC 175#undef MAYBE_sqr_basecase 176#undef MAYBE_sqr_above_basecase 177#undef MAYBE_sqr_toom2 178#undef MAYBE_sqr_above_toom2 179#undef MAYBE_sqr_toom3 180#undef MAYBE_sqr_above_toom3 181#undef MAYBE_sqr_above_toom4 182