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

数学计算

开发平台:

Unix_Linux

  1. /* mpn_mullo_n -- multiply two n-limb numbers and return the low n limbs
  2.    of their products.
  3.    Contributed to the GNU project by Torbjorn Granlund and Marco Bodrato.
  4.    THIS IS (FOR NOW) AN INTERNAL FUNCTION.  IT IS ONLY SAFE TO REACH THIS
  5.    FUNCTION THROUGH DOCUMENTED INTERFACES.  IN FACT, IT IS ALMOST GUARANTEED
  6.    THAT IT'LL CHANGE OR DISAPPEAR IN A FUTURE GNU MP RELEASE.
  7. Copyright 2004, 2005, 2009, 2010 Free Software Foundation, Inc.
  8. This file is part of the GNU MP Library.
  9. The GNU MP Library is free software; you can redistribute it and/or modify
  10. it under the terms of the GNU Lesser General Public License as published by
  11. the Free Software Foundation; either version 3 of the License, or (at your
  12. option) any later version.
  13. The GNU MP Library is distributed in the hope that it will be useful, but
  14. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  15. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  16. License for more details.
  17. You should have received a copy of the GNU Lesser General Public License
  18. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  19. #include "gmp.h"
  20. #include "gmp-impl.h"
  21. #ifndef MULLO_BASECASE_THRESHOLD
  22. #define MULLO_BASECASE_THRESHOLD 0 /* never use mpn_mul_basecase */
  23. #endif
  24. #ifndef MULLO_DC_THRESHOLD
  25. #define MULLO_DC_THRESHOLD 3*MUL_TOOM22_THRESHOLD
  26. #endif
  27. #ifndef MULLO_MUL_N_THRESHOLD
  28. #define MULLO_MUL_N_THRESHOLD MUL_FFT_THRESHOLD
  29. #endif
  30. #if TUNE_PROGRAM_BUILD
  31. #define MAYBE_range_basecase 1
  32. #define MAYBE_range_toom22   1
  33. #else
  34. #define MAYBE_range_basecase                                           
  35.   ((MULLO_DC_THRESHOLD == 0 ? MULLO_BASECASE_THRESHOLD : MULLO_DC_THRESHOLD) < MUL_TOOM22_THRESHOLD*36/(36-11))
  36. #define MAYBE_range_toom22                                             
  37.   ((MULLO_DC_THRESHOLD == 0 ? MULLO_BASECASE_THRESHOLD : MULLO_DC_THRESHOLD) < MUL_TOOM33_THRESHOLD*36/(36-11) )
  38. #endif
  39. /*  THINK: The DC strategy uses different constants in different Toom's
  40.  ranges. Something smoother?
  41. */
  42. /*
  43.   Compute the least significant half of the product {xy,n}*{yp,n}, or
  44.   formally {rp,n} = {xy,n}*{yp,n} Mod (B^n).
  45.   Above the given threshold, the Divide and Conquer strategy is used.
  46.   The operands are split in two, and a full product plus two mullo
  47.   are used to obtain the final result. The more natural strategy is to
  48.   split in two halves, but this is far from optimal when a
  49.   sub-quadratic multiplication is used.
  50.   Mulders suggests an unbalanced split in favour of the full product,
  51.   split n = n1 + n2, where an = n1 <= n2 = (1-a)n; i.e. 0 < a <= 1/2.
  52.   To compute the value of a, we assume that the cost of mullo for a
  53.   given size ML(n) is a fraction of the cost of a full product with
  54.   same size M(n), and the cost M(n)=n^e for some exponent 1 < e <= 2;
  55.   then we can write:
  56.   ML(n) = 2*ML(an) + M((1-a)n) => k*M(n) = 2*k*M(n)*a^e + M(n)*(1-a)^e
  57.   Given a value for e, want to minimise the value of k, i.e. the
  58.   function k=(1-a)^e/(1-2*a^e).
  59.   With e=2, the exponent for schoolbook multiplication, the minimum is
  60.   given by the values a=1-a=1/2.
  61.   With e=log(3)/log(2), the exponent for Karatsuba (aka toom22),
  62.   Mulders compute (1-a) = 0.694... and we approximate a with 11/36.
  63.   Other possible approximations follow:
  64.   e=log(5)/log(3) [Toom-3] -> a ~= 9/40
  65.   e=log(7)/log(4) [Toom-4] -> a ~= 7/39
  66.   e=log(11)/log(6) [Toom-6] -> a ~= 1/8
  67.   e=log(15)/log(8) [Toom-8] -> a ~= 1/10
  68.   The values above where obtained with the following trivial commands
  69.   in the gp-pari shell:
  70. fun(e,a)=(1-a)^e/(1-2*a^e)
  71. mul(a,b,c)={local(m,x,p);if(b-c<1/10000,(b+c)/2,m=1;x=b;forstep(p=c,b,(b-c)/8,if(fun(a,p)<m,m=fun(a,p);x=p));mul(a,(b+x)/2,(c+x)/2))}
  72. contfracpnqn(contfrac(mul(log(2*2-1)/log(2),1/2,0),5))
  73. contfracpnqn(contfrac(mul(log(3*2-1)/log(3),1/2,0),5))
  74. contfracpnqn(contfrac(mul(log(4*2-1)/log(4),1/2,0),5))
  75. contfracpnqn(contfrac(mul(log(6*2-1)/log(6),1/2,0),3))
  76. contfracpnqn(contfrac(mul(log(8*2-1)/log(8),1/2,0),3))
  77.   ,
  78.   |
  79.   | 
  80.   +----,
  81.   |    |
  82.   |    |
  83.   |    |
  84.   |    | 
  85.   +----+--`
  86.   ^ n2 ^n1^
  87.   For an actual implementation, the assumption that M(n)=n^e is
  88.   incorrect, as a consequence also the assumption that ML(n)=k*M(n)
  89.   with a constant k is wrong.
  90.   But theory suggest us two things:
  91.   - the best the multiplication product is (lower e), the more k
  92.     approaches 1, and a approaches 0.
  93.   - A value for a smaller than optimal is probably less bad than a
  94.     bigger one: e.g. let e=log(3)/log(2), a=0.3058_ the optimal
  95.     value, and k(a)=0.808_ the mul/mullo speed ratio. We get
  96.     k(a+1/6)=0.929_ but k(a-1/6)=0.865_.
  97. */
  98. static mp_size_t
  99. mpn_mullo_n_itch (mp_size_t n)
  100. {
  101.   return 2*n;
  102. }
  103. /*
  104.     mpn_dc_mullo_n requires a scratch space of 2*n limbs at tp.
  105.     It accepts tp == rp.
  106. */
  107. static void
  108. mpn_dc_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n, mp_ptr tp)
  109. {
  110.   mp_size_t n2, n1;
  111.   ASSERT (n >= 2);
  112.   ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
  113.   ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
  114.   ASSERT (MPN_SAME_OR_SEPARATE2_P(rp, n, tp, 2*n));
  115.   /* Divide-and-conquer */
  116.   /* We need fractional approximation of the value 0 < a <= 1/2
  117.      giving the minimum in the function k=(1-a)^e/(1-2*a^e).
  118.   */
  119.   if (MAYBE_range_basecase && BELOW_THRESHOLD (n, MUL_TOOM22_THRESHOLD*36/(36-11)))
  120.     n1 = n >> 1;
  121.   else if (MAYBE_range_toom22 && BELOW_THRESHOLD (n, MUL_TOOM33_THRESHOLD*36/(36-11)))
  122.     n1 = n * 11 / (size_t) 36; /* n1 ~= n*(1-.694...) */
  123.   else if (BELOW_THRESHOLD (n, MUL_TOOM44_THRESHOLD*40/(40-9)))
  124.     n1 = n * 9 / (size_t) 40; /* n1 ~= n*(1-.775...) */
  125.   else if (BELOW_THRESHOLD (n, MUL_TOOM8H_THRESHOLD*10/9))
  126.     n1 = n * 7 / (size_t) 39; /* n1 ~= n*(1-.821...) */
  127.   /* n1 = n * 4 / (size_t) 31; // n1 ~= n*(1-.871...) [TOOM66] */
  128.   else
  129.     n1 = n / (size_t) 10; /* n1 ~= n*(1-.899...) [TOOM88] */
  130.   n2 = n - n1;
  131.   /* Split as x = x1 2^(n2 GMP_NUMB_BITS) + x0,
  132.       y = y1 2^(n2 GMP_NUMB_BITS) + y0 */
  133.   /* x0 * y0 */
  134.   mpn_mul_n (tp, xp, yp, n2);
  135.   MPN_COPY (rp, tp, n2);
  136.   /* x1 * y0 * 2^(n2 GMP_NUMB_BITS) */
  137.   if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
  138.     mpn_mul_basecase (tp + n, xp + n2, n1, yp, n1);
  139.   else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
  140.     mpn_mullo_basecase (tp + n, xp + n2, yp, n1);
  141.   else
  142.     mpn_dc_mullo_n (tp + n, xp + n2, yp, n1, tp + n);
  143.   mpn_add_n (rp + n2, tp + n2, tp + n, n1);
  144.   /* x0 * y1 * 2^(n2 GMP_NUMB_BITS) */
  145.   if (BELOW_THRESHOLD (n1, MULLO_BASECASE_THRESHOLD))
  146.     mpn_mul_basecase (tp + n, xp, n1, yp + n2, n1);
  147.   else if (BELOW_THRESHOLD (n1, MULLO_DC_THRESHOLD))
  148.     mpn_mullo_basecase (tp + n, xp, yp + n2, n1);
  149.   else
  150.     mpn_dc_mullo_n (tp + n, xp, yp + n2, n1, tp + n);
  151.   mpn_add_n (rp + n2, rp + n2, tp + n, n1);
  152. }
  153. /* Avoid zero allocations when MULLO_BASECASE_THRESHOLD is 0.  */
  154. #define MUL_BASECASE_ALLOC 
  155.  (MULLO_BASECASE_THRESHOLD_LIMIT == 0 ? 1 : 2*MULLO_BASECASE_THRESHOLD_LIMIT)
  156. /* FIXME: This function should accept a temporary area; dc_mullow_n
  157.    accepts a pointer tp, and handle the case tp == rp, do the same here.
  158.    Maybe recombine the two functions.
  159.    THINK: If mpn_mul_basecase is always faster than mpn_mullo_basecase
  160.   (typically thanks to mpn_addmul_2) should we unconditionally use
  161.   mpn_mul_n?
  162. */
  163. void
  164. mpn_mullo_n (mp_ptr rp, mp_srcptr xp, mp_srcptr yp, mp_size_t n)
  165. {
  166.   ASSERT (n >= 1);
  167.   ASSERT (! MPN_OVERLAP_P (rp, n, xp, n));
  168.   ASSERT (! MPN_OVERLAP_P (rp, n, yp, n));
  169.   if (BELOW_THRESHOLD (n, MULLO_BASECASE_THRESHOLD))
  170.     {
  171.       /* Allocate workspace of fixed size on stack: fast! */
  172.       mp_limb_t tp[MUL_BASECASE_ALLOC];
  173.       mpn_mul_basecase (tp, xp, n, yp, n);
  174.       MPN_COPY (rp, tp, n);
  175.     }
  176.   else if (BELOW_THRESHOLD (n, MULLO_DC_THRESHOLD))
  177.     {
  178.       mpn_mullo_basecase (rp, xp, yp, n);
  179.     }
  180.   else
  181.     {
  182.       mp_ptr tp;
  183.       TMP_DECL;
  184.       TMP_MARK;
  185.       tp = TMP_ALLOC_LIMBS (mpn_mullo_n_itch (n));
  186.       if (BELOW_THRESHOLD (n, MULLO_MUL_N_THRESHOLD))
  187. {
  188.   mpn_dc_mullo_n (rp, xp, yp, n, tp);
  189. }
  190.       else
  191. {
  192.   /* For really large operands, use plain mpn_mul_n but throw away upper n
  193.      limbs of result.  */
  194. #if !TUNE_PROGRAM_BUILD && (MULLO_MUL_N_THRESHOLD > MUL_FFT_THRESHOLD)
  195.   mpn_fft_mul (tp, xp, n, yp, n);
  196. #else
  197.   mpn_mul_n (tp, xp, yp, n);
  198. #endif
  199.   MPN_COPY (rp, tp, n);
  200. }
  201.       TMP_FREE;
  202.     }
  203. }