toom8_sqr.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:7k
源码类别:

数学计算

开发平台:

Unix_Linux

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