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

数学计算

开发平台:

Unix_Linux

  1. /* Implementation of the algorithm for Toom-Cook 4.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. /* Stores |{ap,n}-{bp,n}| in {rp,n}, returns the sign. */
  21. static int
  22. abs_sub_n (mp_ptr rp, mp_srcptr ap, mp_srcptr bp, mp_size_t n)
  23. {
  24.   mp_limb_t  x, y;
  25.   while (--n >= 0)
  26.     {
  27.       x = ap[n];
  28.       y = bp[n];
  29.       if (x != y)
  30. {
  31.   n++;
  32.   if (x > y)
  33.     {
  34.       mpn_sub_n (rp, ap, bp, n);
  35.       return 0;
  36.     }
  37.   else
  38.     {
  39.       mpn_sub_n (rp, bp, ap, n);
  40.       return ~0;
  41.     }
  42. }
  43.       rp[n] = 0;
  44.     }
  45.   return 0;
  46. }
  47. static int
  48. abs_sub_add_n (mp_ptr rm, mp_ptr rp, mp_srcptr rs, mp_size_t n) {
  49.   int result;
  50.   result = abs_sub_n (rm, rp, rs, n);
  51.   ASSERT_NOCARRY(mpn_add_n (rp, rp, rs, n));
  52.   return result;
  53. }
  54. /* Toom-4.5, the splitting 6x3 unbalanced version.
  55.    Evaluate in: infinity, +4, -4, +2, -2, +1, -1, 0.
  56.   <--s-><--n--><--n--><--n--><--n--><--n-->
  57.    ____ ______ ______ ______ ______ ______
  58.   |_a5_|__a4__|__a3__|__a2__|__a1__|__a0__|
  59. |b2_|__b1__|__b0__|
  60. <-t-><--n--><--n-->
  61. */
  62. #define TOOM_63_MUL_N_REC(p, a, b, n, ws)
  63.   do { mpn_mul_n (p, a, b, n);
  64.   } while (0)
  65. #define TOOM_63_MUL_REC(p, a, na, b, nb, ws)
  66.   do { mpn_mul (p, a, na, b, nb);
  67.   } while (0)
  68. void
  69. mpn_toom63_mul (mp_ptr pp,
  70. mp_srcptr ap, mp_size_t an,
  71. mp_srcptr bp, mp_size_t bn, mp_ptr scratch)
  72. {
  73.   mp_size_t n, s, t;
  74.   mp_limb_t cy;
  75.   int sign;
  76.   /***************************** decomposition *******************************/
  77. #define a5  (ap + 5 * n)
  78. #define b0  (bp + 0 * n)
  79. #define b1  (bp + 1 * n)
  80. #define b2  (bp + 2 * n)
  81.   ASSERT (an >= bn);
  82.   n = 1 + (an >= 2 * bn ? (an - 1) / (size_t) 6 : (bn - 1) / (size_t) 3);
  83.   s = an - 5 * n;
  84.   t = bn - 2 * n;
  85.   ASSERT (0 < s && s <= n);
  86.   ASSERT (0 < t && t <= n);
  87.   /* WARNING! it assumes s+t>n */
  88.   ASSERT ( s + t > n );
  89.   ASSERT ( s + t > 4);
  90.   /* WARNING! it assumes n>1 */
  91.   ASSERT ( n > 2);
  92. #define   r8    pp /* 2n   */
  93. #define   r7    scratch /* 3n+1 */
  94. #define   r5    (pp + 3*n) /* 3n+1 */
  95. #define   v0    (pp + 3*n) /* n+1 */
  96. #define   v1    (pp + 4*n+1) /* n+1 */
  97. #define   v2    (pp + 5*n+2) /* n+1 */
  98. #define   v3    (pp + 6*n+3) /* n+1 */
  99. #define   r3    (scratch + 3 * n + 1) /* 3n+1 */
  100. #define   r1    (pp + 7*n) /* s+t <= 2*n */
  101. #define   ws    (scratch + 6 * n + 2) /* ??? */
  102.   /* Alloc also 3n+1 limbs for ws... mpn_toom_interpolate_8pts may
  103.      need all of them, when DO_mpn_sublsh_n usea a scratch  */
  104. /*   if (scratch == NULL) scratch = TMP_SALLOC_LIMBS (9 * n + 3); */
  105.   /********************** evaluation and recursive calls *********************/
  106.   /* $pm4$ */
  107.   sign = mpn_toom_eval_pm2exp (v2, v0, 5, ap, n, s, 2, pp);
  108.   pp[n] = mpn_lshift (pp, b1, n, 2); /* 4b1 */
  109.   /* FIXME: use addlsh */
  110.   v3[t] = mpn_lshift (v3, b2, t, 4);/* 16b2 */
  111.   if ( n == t )
  112.     v3[n]+= mpn_add_n (v3, v3, b0, n); /* 16b2+b0 */
  113.   else
  114.     v3[n] = mpn_add (v3, b0, n, v3, t+1); /* 16b2+b0 */
  115.   sign ^= abs_sub_add_n (v1, v3, pp, n + 1);
  116.   TOOM_63_MUL_N_REC(pp, v0, v1, n + 1, ws); /* A(-4)*B(-4) */
  117.   TOOM_63_MUL_N_REC(r3, v2, v3, n + 1, ws); /* A(+4)*B(+4) */
  118.   mpn_toom_couple_handling (r3, 2*n+1, pp, sign, n, 2, 4);
  119.   /* $pm1$ */
  120.   sign = mpn_toom_eval_pm1 (v2, v0, 5, ap, n, s,    pp);
  121.   /* Compute bs1 and bsm1. Code taken from toom33 */
  122.   cy = mpn_add (ws, b0, n, b2, t);
  123. #if HAVE_NATIVE_mpn_add_n_sub_n
  124.   if (cy == 0 && mpn_cmp (ws, b1, n) < 0)
  125.     {
  126.       cy = mpn_add_n_sub_n (v3, v1, b1, ws, n);
  127.       v3[n] = 0;
  128.       v1[n] = 0;
  129.       sign = ~sign;
  130.     }
  131.   else
  132.     {
  133.       mp_limb_t cy2;
  134.       cy2 = mpn_add_n_sub_n (v3, v1, ws, b1, n);
  135.       w3[n] = cy + (cy2 >> 1);
  136.       w1[n] = cy - (cy & 1);
  137.     }
  138. #else
  139.   v3[n] = cy + mpn_add_n (v3, ws, b1, n);
  140.   if (cy == 0 && mpn_cmp (ws, b1, n) < 0)
  141.     {
  142.       mpn_sub_n (v1, b1, ws, n);
  143.       v1[n] = 0;
  144.       sign = ~sign;
  145.     }
  146.   else
  147.     {
  148.       cy -= mpn_sub_n (v1, ws, b1, n);
  149.       v1[n] = cy;
  150.     }
  151. #endif
  152.   TOOM_63_MUL_N_REC(pp, v0, v1, n + 1, ws); /* A(-1)*B(-1) */
  153.   TOOM_63_MUL_N_REC(r7, v2, v3, n + 1, ws); /* A(1)*B(1) */
  154.   mpn_toom_couple_handling (r7, 2*n+1, pp, sign, n, 0, 0);
  155.   /* $pm2$ */
  156.   sign = mpn_toom_eval_pm2 (v2, v0, 5, ap, n, s, pp);
  157.   pp[n] = mpn_lshift (pp, b1, n, 1); /* 2b1 */
  158.   /* FIXME: use addlsh or addlsh2 */
  159.   v3[t] = mpn_lshift (v3, b2, t, 2);/* 4b2 */
  160.   if ( n == t )
  161.     v3[n]+= mpn_add_n (v3, v3, b0, n); /* 4b2+b0 */
  162.   else
  163.     v3[n] = mpn_add (v3, b0, n, v3, t+1); /* 4b2+b0 */
  164.   sign ^= abs_sub_add_n (v1, v3, pp, n + 1);
  165.   TOOM_63_MUL_N_REC(pp, v0, v1, n + 1, ws); /* A(-2)*B(-2) */
  166.   TOOM_63_MUL_N_REC(r5, v2, v3, n + 1, ws); /* A(+2)*B(+2) */
  167.   mpn_toom_couple_handling (r5, 2*n+1, pp, sign, n, 1, 2);
  168.   /* A(0)*B(0) */
  169.   TOOM_63_MUL_N_REC(pp, ap, bp, n, ws);
  170.   /* Infinity */
  171.   if (s > t) {
  172.     TOOM_63_MUL_REC(r1, a5, s, b2, t, ws);
  173.   } else {
  174.     TOOM_63_MUL_REC(r1, b2, t, a5, s, ws);
  175.   };
  176.   mpn_toom_interpolate_8pts (pp, n, r3, r7, s + t, ws);
  177. #undef a5
  178. #undef b0
  179. #undef b1
  180. #undef b2
  181. #undef r1
  182. #undef r3
  183. #undef r5
  184. #undef v0
  185. #undef v1
  186. #undef v2
  187. #undef v3
  188. #undef r7
  189. #undef r8
  190. #undef ws
  191. }