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

数学计算

开发平台:

Unix_Linux

  1. /* Interpolaton for the algorithm Toom-Cook 6.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, 2010 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 HAVE_NATIVE_mpn_sublsh_n
  21. #define DO_mpn_sublsh_n(dst,src,n,s,ws) mpn_sublsh_n(dst,dst,src,n,s)
  22. #else
  23. static mp_limb_t
  24. DO_mpn_sublsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws)
  25. {
  26. #if USE_MUL_1 && 0
  27.   return mpn_submul_1(dst,src,n,CNST_LIMB(1) <<(s));
  28. #else
  29.   mp_limb_t __cy;
  30.   __cy = mpn_lshift(ws,src,n,s);
  31.   return    __cy + mpn_sub_n(dst,dst,ws,n);
  32. #endif
  33. }
  34. #endif
  35. #if HAVE_NATIVE_mpn_addlsh_n
  36. #define DO_mpn_addlsh_n(dst,src,n,s,ws) mpn_addlsh_n(dst,dst,src,n,s)
  37. #else
  38. static mp_limb_t
  39. DO_mpn_addlsh_n(mp_ptr dst, mp_srcptr src, mp_size_t n, unsigned int s, mp_ptr ws)
  40. {
  41. #if USE_MUL_1 && 0
  42.   return mpn_addmul_1(dst,src,n,CNST_LIMB(1) <<(s));
  43. #else
  44.   mp_limb_t __cy;
  45.   __cy = mpn_lshift(ws,src,n,s);
  46.   return    __cy + mpn_add_n(dst,dst,ws,n);
  47. #endif
  48. }
  49. #endif
  50. #if HAVE_NATIVE_mpn_subrsh
  51. #define DO_mpn_subrsh(dst,nd,src,ns,s,ws) mpn_subrsh(dst,nd,src,ns,s)
  52. #else
  53. /* FIXME: This is not a correct definition, it assumes no carry */
  54. #define DO_mpn_subrsh(dst,nd,src,ns,s,ws)
  55. do {
  56.   mp_limb_t __cy;
  57.   MPN_DECR_U (dst, nd, src[0] >> s);
  58.   __cy = DO_mpn_sublsh_n (dst, src + 1, ns - 1, GMP_NUMB_BITS - s, ws);
  59.   MPN_DECR_U (dst + ns - 1, nd - ns + 1, __cy);
  60. } while (0)
  61. #endif
  62. #if GMP_NUMB_BITS < 21
  63. #error Not implemented: Both sublsh_n(,,,20) should be corrected.
  64. #endif
  65. #if GMP_NUMB_BITS < 16
  66. #error Not implemented: divexact_by42525 needs splitting.
  67. #endif
  68. #if GMP_NUMB_BITS < 12
  69. #error Not implemented: Hard to adapt...
  70. #endif
  71. /* FIXME: tuneup should decide the best variant */
  72. #ifndef AORSMUL_FASTER_AORS_AORSLSH
  73. #define AORSMUL_FASTER_AORS_AORSLSH 1
  74. #endif
  75. #ifndef AORSMUL_FASTER_AORS_2AORSLSH
  76. #define AORSMUL_FASTER_AORS_2AORSLSH 1
  77. #endif
  78. #ifndef AORSMUL_FASTER_2AORSLSH
  79. #define AORSMUL_FASTER_2AORSLSH 1
  80. #endif
  81. #ifndef AORSMUL_FASTER_3AORSLSH
  82. #define AORSMUL_FASTER_3AORSLSH 1
  83. #endif
  84. #define BINVERT_9 
  85.   ((((GMP_NUMB_MAX / 9) << (6 - GMP_NUMB_BITS % 6)) * 8 & GMP_NUMB_MAX) | 0x39)
  86. #define BINVERT_255 
  87.   (GMP_NUMB_MAX - ((GMP_NUMB_MAX / 255) << (8 - GMP_NUMB_BITS % 8)))
  88.   /* FIXME: find some more general expressions for 2835^-1, 42525^-1 */
  89. #if GMP_LIMB_BITS == 32
  90. #define BINVERT_2835  (GMP_NUMB_MASK & CNST_LIMB(0x53E3771B))
  91. #define BINVERT_42525 (GMP_NUMB_MASK & CNST_LIMB(0x9F314C35))
  92. #else
  93. #if GMP_LIMB_BITS == 64
  94. #define BINVERT_2835  (GMP_NUMB_MASK & CNST_LIMB(0x938CC70553E3771B))
  95. #define BINVERT_42525 (GMP_NUMB_MASK & CNST_LIMB(0xE7B40D449F314C35))
  96. #endif
  97. #endif
  98. #ifndef mpn_divexact_by255
  99. #if GMP_NUMB_BITS % 8 == 0
  100. #define mpn_divexact_by255(dst,src,size) 
  101.   (255 & 1 * mpn_bdiv_dbm1 (dst, src, size, __GMP_CAST (mp_limb_t, GMP_NUMB_MASK / 255)))
  102. #else
  103. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
  104. #define mpn_divexact_by255(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(255),BINVERT_255,0)
  105. #else
  106. #define mpn_divexact_by255(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(255))
  107. #endif
  108. #endif
  109. #endif
  110. #ifndef mpn_divexact_by9x4
  111. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1
  112. #define mpn_divexact_by9x4(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(9),BINVERT_9,2)
  113. #else
  114. #define mpn_divexact_by9x4(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(9)<<2)
  115. #endif
  116. #endif
  117. #ifndef mpn_divexact_by42525
  118. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_42525)
  119. #define mpn_divexact_by42525(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(42525),BINVERT_42525,0)
  120. #else
  121. #define mpn_divexact_by42525(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(42525))
  122. #endif
  123. #endif
  124. #ifndef mpn_divexact_by2835x4
  125. #if HAVE_NATIVE_mpn_pi1_bdiv_q_1 && defined(BINVERT_2835)
  126. #define mpn_divexact_by2835x4(dst,src,size) mpn_pi1_bdiv_q_1(dst,src,size,CNST_LIMB(2835),BINVERT_2835,2)
  127. #else
  128. #define mpn_divexact_by2835x4(dst,src,size) mpn_divexact_1(dst,src,size,CNST_LIMB(2835)<<2)
  129. #endif
  130. #endif
  131. /* Interpolation for Toom-6.5 (or Toom-6), using the evaluation
  132.    points: infinity(6.5 only), +-4, +-2, +-1, +-1/4, +-1/2, 0. More precisely,
  133.    we want to compute f(2^(GMP_NUMB_BITS * n)) for a polynomial f of
  134.    degree 11 (or 10), given the 12 (rsp. 11) values:
  135.      r0 = limit at infinity of f(x) / x^7,
  136.      r1 = f(4),f(-4),
  137.      r2 = f(2),f(-2),
  138.      r3 = f(1),f(-1),
  139.      r4 = f(1/4),f(-1/4),
  140.      r5 = f(1/2),f(-1/2),
  141.      r6 = f(0).
  142.    All couples of the form f(n),f(-n) must be already mixed with
  143.    toom_couple_handling(f(n),...,f(-n),...)
  144.    The result is stored in {pp, spt + 7*n (or 6*n)}.
  145.    At entry, r6 is stored at {pp, 2n},
  146.    r4 is stored at {pp + 3n, 3n + 1}.
  147.    r2 is stored at {pp + 7n, 3n + 1}.
  148.    r0 is stored at {pp +11n, spt}.
  149.    The other values are 3n+1 limbs each (with most significant limbs small).
  150.    Negative intermediate results are stored two-complemented.
  151.    Inputs are destroyed.
  152. */
  153. void
  154. mpn_toom_interpolate_12pts (mp_ptr pp, mp_ptr r1, mp_ptr r3, mp_ptr r5,
  155. mp_size_t n, mp_size_t spt, int half, mp_ptr wsi)
  156. {
  157.   mp_limb_t cy;
  158.   mp_size_t n3;
  159.   mp_size_t n3p1;
  160.   n3 = 3 * n;
  161.   n3p1 = n3 + 1;
  162. #define   r4    (pp + n3) /* 3n+1 */
  163. #define   r2    (pp + 7 * n) /* 3n+1 */
  164. #define   r0    (pp +11 * n) /* s+t <= 2*n */
  165.   /******************************* interpolation *****************************/
  166.   if (half != 0) {
  167.     cy = mpn_sub_n (r3, r3, r0, spt);
  168.     MPN_DECR_U (r3 + spt, n3p1 - spt, cy);
  169.     cy = DO_mpn_sublsh_n (r2, r0, spt, 10, wsi);
  170.     MPN_DECR_U (r2 + spt, n3p1 - spt, cy);
  171.     DO_mpn_subrsh(r5, n3p1, r0, spt, 2, wsi);
  172.     cy = DO_mpn_sublsh_n (r1, r0, spt, 20, wsi);
  173.     MPN_DECR_U (r1 + spt, n3p1 - spt, cy);
  174.     DO_mpn_subrsh(r4, n3p1, r0, spt, 4, wsi);
  175.   };
  176.   r4[n3] -= DO_mpn_sublsh_n (r4 + n, pp, 2 * n, 20, wsi);
  177.   DO_mpn_subrsh(r1 + n, 2 * n + 1, pp, 2 * n, 4, wsi);
  178. #if HAVE_NATIVE_mpn_add_n_sub_n
  179.   mpn_add_n_sub_n (r1, r4, r4, r1, n3p1);
  180. #else
  181.   ASSERT_NOCARRY(mpn_add_n (wsi, r1, r4, n3p1));
  182.   mpn_sub_n (r4, r4, r1, n3p1); /* can be negative */
  183.   MP_PTR_SWAP(r1, wsi);
  184. #endif
  185.   r5[n3] -= DO_mpn_sublsh_n (r5 + n, pp, 2 * n, 10, wsi);
  186.   DO_mpn_subrsh(r2 + n, 2 * n + 1, pp, 2 * n, 2, wsi);
  187. #if HAVE_NATIVE_mpn_add_n_sub_n
  188.   mpn_add_n_sub_n (r2, r5, r5, r2, n3p1);
  189. #else
  190.   mpn_sub_n (wsi, r5, r2, n3p1); /* can be negative */
  191.   ASSERT_NOCARRY(mpn_add_n (r2, r2, r5, n3p1));
  192.   MP_PTR_SWAP(r5, wsi);
  193. #endif
  194.   r3[n3] -= mpn_sub_n (r3+n, r3+n, pp, 2 * n);
  195. #if AORSMUL_FASTER_AORS_AORSLSH
  196.   mpn_submul_1 (r4, r5, n3p1, 257); /* can be negative */
  197. #else
  198.   mpn_sub_n (r4, r4, r5, n3p1); /* can be negative */
  199.   DO_mpn_sublsh_n (r4, r5, n3p1, 8, wsi); /* can be negative */
  200. #endif
  201.   /* A division by 2835x4 followsi. Warning: the operand can be negative! */
  202.   mpn_divexact_by2835x4(r4, r4, n3p1);
  203.   if ((r4[n3] & (GMP_NUMB_MAX << (GMP_NUMB_BITS-3))) != 0)
  204.     r4[n3] |= (GMP_NUMB_MAX << (GMP_NUMB_BITS-2));
  205. #if AORSMUL_FASTER_2AORSLSH
  206.   mpn_addmul_1 (r5, r4, n3p1, 60); /* can be negative */
  207. #else
  208.   DO_mpn_sublsh_n (r5, r4, n3p1, 2, wsi); /* can be negative */
  209.   DO_mpn_addlsh_n (r5, r4, n3p1, 6, wsi); /* can give a carry */
  210. #endif
  211.   mpn_divexact_by255(r5, r5, n3p1);
  212.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r2, r3, n3p1, 5, wsi));
  213. #if AORSMUL_FASTER_3AORSLSH
  214.   ASSERT_NOCARRY(mpn_submul_1 (r1, r2, n3p1, 100));
  215. #else
  216.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r1, r2, n3p1, 6, wsi));
  217.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r1, r2, n3p1, 5, wsi));
  218.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r1, r2, n3p1, 2, wsi));
  219. #endif
  220.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r1, r3, n3p1, 9, wsi));
  221.   mpn_divexact_by42525(r1, r1, n3p1);
  222. #if AORSMUL_FASTER_AORS_2AORSLSH
  223.   ASSERT_NOCARRY(mpn_submul_1 (r2, r1, n3p1, 225));
  224. #else
  225.   ASSERT_NOCARRY(mpn_sub_n (r2, r2, r1, n3p1));
  226.   ASSERT_NOCARRY(DO_mpn_addlsh_n (r2, r1, n3p1, 5, wsi));
  227.   ASSERT_NOCARRY(DO_mpn_sublsh_n (r2, r1, n3p1, 8, wsi));
  228. #endif
  229.   mpn_divexact_by9x4(r2, r2, n3p1);
  230.   ASSERT_NOCARRY(mpn_sub_n (r3, r3, r2, n3p1));
  231.   mpn_sub_n (r4, r2, r4, n3p1);
  232.   ASSERT_NOCARRY(mpn_rshift(r4, r4, n3p1, 1));
  233.   ASSERT_NOCARRY(mpn_sub_n (r2, r2, r4, n3p1));
  234.   mpn_add_n (r5, r5, r1, n3p1);
  235.   ASSERT_NOCARRY(mpn_rshift(r5, r5, n3p1, 1));
  236.   /* last interpolation steps... */
  237.   ASSERT_NOCARRY(mpn_sub_n (r3, r3, r1, n3p1));
  238.   ASSERT_NOCARRY(mpn_sub_n (r1, r1, r5, n3p1));
  239.   /* ... could be mixed with recomposition
  240. ||H-r5|M-r5|L-r5|   ||H-r1|M-r1|L-r1|
  241.   */
  242.   /***************************** recomposition *******************************/
  243.   /*
  244.     pp[] prior to operations:
  245.     |M r0|L r0|___||H r2|M r2|L r2|___||H r4|M r4|L r4|____|H_r6|L r6|pp
  246.     summation scheme for remaining operations:
  247.     |__12|n_11|n_10|n__9|n__8|n__7|n__6|n__5|n__4|n__3|n__2|n___|n___|pp
  248.     |M r0|L r0|___||H r2|M r2|L r2|___||H r4|M r4|L r4|____|H_r6|L r6|pp
  249. ||H r1|M r1|L r1|   ||H r3|M r3|L r3|   ||H_r5|M_r5|L_r5|
  250.   */
  251.   cy = mpn_add_n (pp + n, pp + n, r5, n);
  252.   cy = mpn_add_1 (pp + 2 * n, r5 + n, n, cy);
  253. #if HAVE_NATIVE_mpn_add_nc
  254.   cy = r5[n3] + mpn_add_nc(pp + n3, pp + n3, r5 + 2 * n, n, cy);
  255. #else
  256.   MPN_INCR_U (r5 + 2 * n, n + 1, cy);
  257.   cy = r5[n3] + mpn_add_n (pp + n3, pp + n3, r5 + 2 * n, n);
  258. #endif
  259.   MPN_INCR_U (pp + n3 + n, 2 * n + 1, cy);
  260.   pp[2 * n3]+= mpn_add_n (pp + 5 * n, pp + 5 * n, r3, n);
  261.   cy = mpn_add_1 (pp + 2 * n3, r3 + n, n, pp[2 * n3]);
  262. #if HAVE_NATIVE_mpn_add_nc
  263.   cy = r3[n3] + mpn_add_nc(pp + 7 * n, pp + 7 * n, r3 + 2 * n, n, cy);
  264. #else
  265.   MPN_INCR_U (r3 + 2 * n, n + 1, cy);
  266.   cy = r3[n3] + mpn_add_n (pp + 7 * n, pp + 7 * n, r3 + 2 * n, n);
  267. #endif
  268.   MPN_INCR_U (pp + 8 * n, 2 * n + 1, cy);
  269.   pp[10*n]+=mpn_add_n (pp + 9 * n, pp + 9 * n, r1, n);
  270.   if (half) {
  271.     cy = mpn_add_1 (pp + 10 * n, r1 + n, n, pp[10 * n]);
  272. #if HAVE_NATIVE_mpn_add_nc
  273.     if (LIKELY (spt > n)) {
  274.       cy = r1[n3] + mpn_add_nc(pp + 11 * n, pp + 11 * n, r1 + 2 * n, n, cy);
  275.       MPN_INCR_U (pp + 4 * n3, spt - n, cy);
  276.     } else {
  277.       ASSERT_NOCARRY(mpn_add_nc(pp + 11 * n, pp + 11 * n, r1 + 2 * n, spt, cy));
  278.     }
  279. #else
  280.     MPN_INCR_U (r1 + 2 * n, n + 1, cy);
  281.     if (LIKELY (spt > n)) {
  282.       cy = r1[n3] + mpn_add_n (pp + 11 * n, pp + 11 * n, r1 + 2 * n, n);
  283.       MPN_INCR_U (pp + 4 * n3, spt - n, cy);
  284.     } else {
  285.       ASSERT_NOCARRY(mpn_add_n (pp + 11 * n, pp + 11 * n, r1 + 2 * n, spt));
  286.     }
  287. #endif
  288.   } else {
  289.     ASSERT_NOCARRY(mpn_add_1 (pp + 10 * n, r1 + n, spt, pp[10 * n]));
  290.   }
  291. #undef   r0
  292. #undef   r2
  293. #undef   r4
  294. }