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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_sqr_basecase -- Internal routine to square a natural number
  2.    of length n.
  3.    THIS IS AN INTERNAL FUNCTION WITH A MUTABLE INTERFACE.  IT IS ONLY
  4.    SAFE TO REACH THIS FUNCTION THROUGH DOCUMENTED INTERFACES.
  5. Copyright 1991, 1992, 1993, 1994, 1996, 1997, 2000, 2001, 2002, 2003, 2004,
  6. 2005, 2008 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. #include "longlong.h"
  21. #if HAVE_NATIVE_mpn_sqr_diagonal
  22. #define MPN_SQR_DIAGONAL(rp, up, n)
  23.   mpn_sqr_diagonal (rp, up, n)
  24. #else
  25. #define MPN_SQR_DIAGONAL(rp, up, n)
  26.   do {
  27.     mp_size_t _i;
  28.     for (_i = 0; _i < (n); _i++)
  29.       {
  30. mp_limb_t ul, lpl;
  31. ul = (up)[_i];
  32. umul_ppmm ((rp)[2 * _i + 1], lpl, ul, ul << GMP_NAIL_BITS);
  33. (rp)[2 * _i] = lpl >> GMP_NAIL_BITS;
  34.       }
  35.   } while (0)
  36. #endif
  37. #undef READY_WITH_mpn_sqr_basecase
  38. #if ! defined (READY_WITH_mpn_sqr_basecase) && HAVE_NATIVE_mpn_addmul_2s
  39. void
  40. mpn_sqr_basecase (mp_ptr rp, mp_srcptr up, mp_size_t n)
  41. {
  42.   mp_size_t i;
  43.   mp_limb_t tarr[2 * SQR_TOOM2_THRESHOLD];
  44.   mp_ptr tp = tarr;
  45.   mp_limb_t cy;
  46.   /* must fit 2*n limbs in tarr */
  47.   ASSERT (n <= SQR_TOOM2_THRESHOLD);
  48.   if ((n & 1) != 0)
  49.     {
  50.       if (n == 1)
  51. {
  52.   mp_limb_t ul, lpl;
  53.   ul = up[0];
  54.   umul_ppmm (rp[1], lpl, ul, ul << GMP_NAIL_BITS);
  55.   rp[0] = lpl >> GMP_NAIL_BITS;
  56.   return;
  57. }
  58.       MPN_ZERO (tp, n);
  59.       for (i = 0; i <= n - 2; i += 2)
  60. {
  61.   cy = mpn_addmul_2s (tp + 2 * i, up + i + 1, n - (i + 1), up + i);
  62.   tp[n + i] = cy;
  63. }
  64.     }
  65.   else
  66.     {
  67.       if (n == 2)
  68. {
  69.   rp[0] = 0;
  70.   rp[1] = 0;
  71.   rp[3] = mpn_addmul_2 (rp, up, 2, up);
  72.   return;
  73. }
  74.       MPN_ZERO (tp, n);
  75.       for (i = 0; i <= n - 4; i += 2)
  76. {
  77.   cy = mpn_addmul_2s (tp + 2 * i, up + i + 1, n - (i + 1), up + i);
  78.   tp[n + i] = cy;
  79. }
  80.       cy = mpn_addmul_1 (tp + 2 * n - 4, up + n - 1, 1, up[n - 2]);
  81.       tp[2 * n - 3] = cy;
  82.     }
  83.   MPN_SQR_DIAGONAL (rp, up, n);
  84. #if HAVE_NATIVE_mpn_addlsh1_n
  85.   cy = mpn_addlsh1_n (rp + 1, rp + 1, tp, 2 * n - 2);
  86. #else
  87.   cy = mpn_lshift (tp, tp, 2 * n - 2, 1);
  88.   cy += mpn_add_n (rp + 1, rp + 1, tp, 2 * n - 2);
  89. #endif
  90.   rp[2 * n - 1] += cy;
  91. }
  92. #define READY_WITH_mpn_sqr_basecase
  93. #endif
  94. #if ! defined (READY_WITH_mpn_sqr_basecase) && HAVE_NATIVE_mpn_addmul_2
  95. /* mpn_sqr_basecase using plain mpn_addmul_2.
  96.    This is tricky, since we have to let mpn_addmul_2 make some undesirable
  97.    multiplies, u[k]*u[k], that we would like to let mpn_sqr_diagonal handle.
  98.    This forces us to conditionally add or subtract the mpn_sqr_diagonal
  99.    results.  Examples of the product we form:
  100.    n = 4              n = 5 n = 6
  101.    u1u0 * u3u2u1      u1u0 * u4u3u2u1 u1u0 * u5u4u3u2u1
  102.    u2 * u3       u3u2 * u4u3 u3u2 * u5u4u3
  103. u4 * u5
  104.    add: u0 u2 u3      add: u0 u2 u4 add: u0 u2 u4 u5
  105.    sub: u1       sub: u1 u3 sub: u1 u3
  106. */
  107. void
  108. mpn_sqr_basecase (mp_ptr rp, mp_srcptr up, mp_size_t n)
  109. {
  110.   mp_size_t i;
  111.   mp_limb_t tarr[2 * SQR_TOOM2_THRESHOLD];
  112.   mp_ptr tp = tarr;
  113.   mp_limb_t cy;
  114.   /* must fit 2*n limbs in tarr */
  115.   ASSERT (n <= SQR_TOOM2_THRESHOLD);
  116.   if ((n & 1) != 0)
  117.     {
  118.       mp_limb_t x0, x1;
  119.       if (n == 1)
  120. {
  121.   mp_limb_t ul, lpl;
  122.   ul = up[0];
  123.   umul_ppmm (rp[1], lpl, ul, ul << GMP_NAIL_BITS);
  124.   rp[0] = lpl >> GMP_NAIL_BITS;
  125.   return;
  126. }
  127.       /* The code below doesn't like unnormalized operands.  Since such
  128.  operands are unusual, handle them with a dumb recursion.  */
  129.       if (up[n - 1] == 0)
  130. {
  131.   rp[2 * n - 2] = 0;
  132.   rp[2 * n - 1] = 0;
  133.   mpn_sqr_basecase (rp, up, n - 1);
  134.   return;
  135. }
  136.       MPN_ZERO (tp, n);
  137.       for (i = 0; i <= n - 2; i += 2)
  138. {
  139.   cy = mpn_addmul_2 (tp + 2 * i, up + i + 1, n - (i + 1), up + i);
  140.   tp[n + i] = cy;
  141. }
  142.       MPN_SQR_DIAGONAL (rp, up, n);
  143.       for (i = 2;; i += 4)
  144. {
  145.   x0 = rp[i + 0];
  146.   rp[i + 0] = (-x0) & GMP_NUMB_MASK;
  147.   x1 = rp[i + 1];
  148.   rp[i + 1] = (-x1 - (x0 != 0)) & GMP_NUMB_MASK;
  149.   __GMPN_SUB_1 (cy, rp + i + 2, rp + i + 2, 2, (x1 | x0) != 0);
  150.   if (i + 4 >= 2 * n)
  151.     break;
  152.   mpn_incr_u (rp + i + 4, cy);
  153. }
  154.     }
  155.   else
  156.     {
  157.       mp_limb_t x0, x1;
  158.       if (n == 2)
  159. {
  160.   rp[0] = 0;
  161.   rp[1] = 0;
  162.   rp[3] = mpn_addmul_2 (rp, up, 2, up);
  163.   return;
  164. }
  165.       /* The code below doesn't like unnormalized operands.  Since such
  166.  operands are unusual, handle them with a dumb recursion.  */
  167.       if (up[n - 1] == 0)
  168. {
  169.   rp[2 * n - 2] = 0;
  170.   rp[2 * n - 1] = 0;
  171.   mpn_sqr_basecase (rp, up, n - 1);
  172.   return;
  173. }
  174.       MPN_ZERO (tp, n);
  175.       for (i = 0; i <= n - 4; i += 2)
  176. {
  177.   cy = mpn_addmul_2 (tp + 2 * i, up + i + 1, n - (i + 1), up + i);
  178.   tp[n + i] = cy;
  179. }
  180.       cy = mpn_addmul_1 (tp + 2 * n - 4, up + n - 1, 1, up[n - 2]);
  181.       tp[2 * n - 3] = cy;
  182.       MPN_SQR_DIAGONAL (rp, up, n);
  183.       for (i = 2;; i += 4)
  184. {
  185.   x0 = rp[i + 0];
  186.   rp[i + 0] = (-x0) & GMP_NUMB_MASK;
  187.   x1 = rp[i + 1];
  188.   rp[i + 1] = (-x1 - (x0 != 0)) & GMP_NUMB_MASK;
  189.   if (i + 6 >= 2 * n)
  190.     break;
  191.   __GMPN_SUB_1 (cy, rp + i + 2, rp + i + 2, 2, (x1 | x0) != 0);
  192.   mpn_incr_u (rp + i + 4, cy);
  193. }
  194.       mpn_decr_u (rp + i + 2, (x1 | x0) != 0);
  195.     }
  196. #if HAVE_NATIVE_mpn_addlsh1_n
  197.   cy = mpn_addlsh1_n (rp + 1, rp + 1, tp, 2 * n - 2);
  198. #else
  199.   cy = mpn_lshift (tp, tp, 2 * n - 2, 1);
  200.   cy += mpn_add_n (rp + 1, rp + 1, tp, 2 * n - 2);
  201. #endif
  202.   rp[2 * n - 1] += cy;
  203. }
  204. #define READY_WITH_mpn_sqr_basecase
  205. #endif
  206. #if ! defined (READY_WITH_mpn_sqr_basecase)
  207. /* Default mpn_sqr_basecase using mpn_addmul_1.  */
  208. void
  209. mpn_sqr_basecase (mp_ptr rp, mp_srcptr up, mp_size_t n)
  210. {
  211.   mp_size_t i;
  212.   ASSERT (n >= 1);
  213.   ASSERT (! MPN_OVERLAP_P (rp, 2*n, up, n));
  214.   {
  215.     mp_limb_t ul, lpl;
  216.     ul = up[0];
  217.     umul_ppmm (rp[1], lpl, ul, ul << GMP_NAIL_BITS);
  218.     rp[0] = lpl >> GMP_NAIL_BITS;
  219.   }
  220.   if (n > 1)
  221.     {
  222.       mp_limb_t tarr[2 * SQR_TOOM2_THRESHOLD];
  223.       mp_ptr tp = tarr;
  224.       mp_limb_t cy;
  225.       /* must fit 2*n limbs in tarr */
  226.       ASSERT (n <= SQR_TOOM2_THRESHOLD);
  227.       cy = mpn_mul_1 (tp, up + 1, n - 1, up[0]);
  228.       tp[n - 1] = cy;
  229.       for (i = 2; i < n; i++)
  230. {
  231.   mp_limb_t cy;
  232.   cy = mpn_addmul_1 (tp + 2 * i - 2, up + i, n - i, up[i - 1]);
  233.   tp[n + i - 2] = cy;
  234. }
  235.       MPN_SQR_DIAGONAL (rp + 2, up + 1, n - 1);
  236.       {
  237. mp_limb_t cy;
  238. #if HAVE_NATIVE_mpn_addlsh1_n
  239. cy = mpn_addlsh1_n (rp + 1, rp + 1, tp, 2 * n - 2);
  240. #else
  241. cy = mpn_lshift (tp, tp, 2 * n - 2, 1);
  242. cy += mpn_add_n (rp + 1, rp + 1, tp, 2 * n - 2);
  243. #endif
  244. rp[2 * n - 1] += cy;
  245.       }
  246.     }
  247. }
  248. #endif