home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 4: GNU Archives / Linux Cubed Series 4 - GNU Archives.iso / gnu / glibc-1.09 / glibc-1 / glibc-1.09.1 / sysdeps / generic / mul_n.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-03  |  11.7 KB  |  421 lines
  1. /* __mpn_mul_n -- Multiply two natural numbers of length n.
  2.  
  3. Copyright (C) 1991, 1992, 1993, 1994 Free Software Foundation, Inc.
  4.  
  5. This file is part of the GNU MP Library.
  6.  
  7. The GNU MP Library is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU Library General Public License as published by
  9. the Free Software Foundation; either version 2 of the License, or (at your
  10. option) any later version.
  11.  
  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 Library General Public
  15. License for more details.
  16.  
  17. You should have received a copy of the GNU Library General Public License
  18. along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  19. the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
  20.  
  21. #include "gmp.h"
  22. #include "gmp-impl.h"
  23.  
  24. /* Multiply the natural numbers u (pointed to by UP) and v (pointed to by VP),
  25.    both with SIZE limbs, and store the result at PRODP.  2 * SIZE limbs are
  26.    always stored.  Return the most significant limb.
  27.  
  28.    Argument constraints:
  29.    1. PRODP != UP and PRODP != VP, i.e. the destination
  30.       must be distinct from the multiplier and the multiplicand.  */
  31.  
  32. /* If KARATSUBA_THRESHOLD is not already defined, define it to a
  33.    value which is good on most machines.  */
  34. #ifndef KARATSUBA_THRESHOLD
  35. #define KARATSUBA_THRESHOLD 32
  36. #endif
  37.  
  38. /* The code can't handle KARATSUBA_THRESHOLD smaller than 2.  */
  39. #if KARATSUBA_THRESHOLD < 2
  40. #undef KARATSUBA_THRESHOLD
  41. #define KARATSUBA_THRESHOLD 2
  42. #endif
  43.  
  44. void
  45. #if __STDC__
  46. ____mpn_mul_n (mp_ptr, mp_srcptr, mp_srcptr, mp_size_t, mp_ptr);
  47. #else
  48. ____mpn_mul_n ();
  49. #endif
  50.  
  51. /* Handle simple cases with traditional multiplication.
  52.  
  53.    This is the most critical code of multiplication.  All multiplies rely
  54.    on this, both small and huge.  Small ones arrive here immediately.  Huge
  55.    ones arrive here as this is the base case for Karatsuba's recursive
  56.    algorithm below.  */
  57.  
  58. void
  59. #if __STDC__
  60. ____mpn_mul_n_basecase (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
  61. #else
  62. ____mpn_mul_n_basecase (prodp, up, vp, size)
  63.      mp_ptr prodp;
  64.      mp_srcptr up;
  65.      mp_srcptr vp;
  66.      mp_size_t size;
  67. #endif
  68. {
  69.   mp_size_t i;
  70.   mp_limb cy_limb;
  71.   mp_limb v_limb;
  72.  
  73.   /* Multiply by the first limb in V separately, as the result can be
  74.      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
  75.   v_limb = vp[0];
  76.   if (v_limb <= 1)
  77.     {
  78.       if (v_limb == 1)
  79.     MPN_COPY (prodp, up, size);
  80.       else
  81.     MPN_ZERO (prodp, size);
  82.       cy_limb = 0;
  83.     }
  84.   else
  85.     cy_limb = __mpn_mul_1 (prodp, up, size, v_limb);
  86.  
  87.   prodp[size] = cy_limb;
  88.   prodp++;
  89.  
  90.   /* For each iteration in the outer loop, multiply one limb from
  91.      U with one limb from V, and add it to PROD.  */
  92.   for (i = 1; i < size; i++)
  93.     {
  94.       v_limb = vp[i];
  95.       if (v_limb <= 1)
  96.     {
  97.       cy_limb = 0;
  98.       if (v_limb == 1)
  99.         cy_limb = __mpn_add_n (prodp, prodp, up, size);
  100.     }
  101.       else
  102.     cy_limb = __mpn_addmul_1 (prodp, up, size, v_limb);
  103.  
  104.       prodp[size] = cy_limb;
  105.       prodp++;
  106.     }
  107. }
  108.  
  109. void
  110. #if __STDC__
  111. ____mpn_mul_n (mp_ptr prodp,
  112.          mp_srcptr up, mp_srcptr vp, mp_size_t size, mp_ptr tspace)
  113. #else
  114. ____mpn_mul_n (prodp, up, vp, size, tspace)
  115.      mp_ptr prodp;
  116.      mp_srcptr up;
  117.      mp_srcptr vp;
  118.      mp_size_t size;
  119.      mp_ptr tspace;
  120. #endif
  121. {
  122.   if ((size & 1) != 0)
  123.     {
  124.       /* The size is odd, the code code below doesn't handle that.
  125.      Multiply the least significant (size - 1) limbs with a recursive
  126.      call, and handle the most significant limb of S1 and S2
  127.      separately.  */
  128.       /* A slightly faster way to do this would be to make the Karatsuba
  129.      code below behave as if the size were even, and let it check for
  130.      odd size in the end.  I.e., in essence move this code to the end.
  131.      Doing so would save us a recursive call, and potentially make the
  132.      stack grow a lot less.  */
  133.  
  134.       mp_size_t esize = size - 1;    /* even size */
  135.       mp_limb cy_limb;
  136.  
  137.       MPN_MUL_N_RECURSE (prodp, up, vp, esize, tspace);
  138.       cy_limb = __mpn_addmul_1 (prodp + esize, up, esize, vp[esize]);
  139.       prodp[esize + esize] = cy_limb;
  140.       cy_limb = __mpn_addmul_1 (prodp + esize, vp, size, up[esize]);
  141.  
  142.       prodp[esize + size] = cy_limb;
  143.     }
  144.   else
  145.     {
  146.       /* Anatolij Alekseevich Karatsuba's divide-and-conquer algorithm.
  147.  
  148.      Split U in two pieces, U1 and U0, such that
  149.      U = U0 + U1*(B**n),
  150.      and V in V1 and V0, such that
  151.      V = V0 + V1*(B**n).
  152.  
  153.      UV is then computed recursively using the identity
  154.  
  155.         2n   n          n                     n
  156.      UV = (B  + B )U V  +  B (U -U )(V -V )  +  (B + 1)U V
  157.             1 1        1  0   0  1              0 0
  158.  
  159.      Where B = 2**BITS_PER_MP_LIMB.  */
  160.  
  161.       mp_size_t hsize = size >> 1;
  162.       mp_limb cy;
  163.       int negflg;
  164.  
  165.       /*** Product H.     ________________  ________________
  166.             |_____U1 x V1____||____U0 x V0_____|  */
  167.       /* Put result in upper part of PROD and pass low part of TSPACE
  168.      as new TSPACE.  */
  169.       MPN_MUL_N_RECURSE (prodp + size, up + hsize, vp + hsize, hsize, tspace);
  170.  
  171.       /*** Product M.     ________________
  172.             |_(U1-U0)(V0-V1)_|  */
  173.       if (__mpn_cmp (up + hsize, up, hsize) >= 0)
  174.     {
  175.       __mpn_sub_n (prodp, up + hsize, up, hsize);
  176.       negflg = 0;
  177.     }
  178.       else
  179.     {
  180.       __mpn_sub_n (prodp, up, up + hsize, hsize);
  181.       negflg = 1;
  182.     }
  183.       if (__mpn_cmp (vp + hsize, vp, hsize) >= 0)
  184.     {
  185.       __mpn_sub_n (prodp + hsize, vp + hsize, vp, hsize);
  186.       negflg ^= 1;
  187.     }
  188.       else
  189.     {
  190.       __mpn_sub_n (prodp + hsize, vp, vp + hsize, hsize);
  191.       /* No change of NEGFLG.  */
  192.     }
  193.       /* Read temporary operands from low part of PROD.
  194.      Put result in low part of TSPACE using upper part of TSPACE
  195.      as new TSPACE.  */
  196.       MPN_MUL_N_RECURSE (tspace, prodp, prodp + hsize, hsize, tspace + size);
  197.  
  198.       /*** Add/copy product H.  */
  199.       MPN_COPY (prodp + hsize, prodp + size, hsize);
  200.       cy = __mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
  201.  
  202.       /*** Add product M (if NEGFLG M is a negative number).  */
  203.       if (negflg)
  204.     cy -= __mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
  205.       else
  206.     cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
  207.  
  208.       /*** Product L.     ________________  ________________
  209.             |________________||____U0 x V0_____|  */
  210.       /* Read temporary operands from low part of PROD.
  211.      Put result in low part of TSPACE using upper part of TSPACE
  212.      as new TSPACE.  */
  213.       MPN_MUL_N_RECURSE (tspace, up, vp, hsize, tspace + size);
  214.  
  215.       /*** Add/copy Product L (twice).  */
  216.  
  217.       cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
  218.       if (cy)
  219.     {
  220.       if (cy > 0)
  221.         __mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
  222.       else
  223.         {
  224.           __mpn_sub_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
  225.           abort ();
  226.         }
  227.     }
  228.  
  229.       MPN_COPY (prodp, tspace, hsize);
  230.       cy = __mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
  231.       if (cy)
  232.     __mpn_add_1 (prodp + size, prodp + size, size, 1);
  233.     }
  234. }
  235.  
  236. void
  237. #if __STDC__
  238. ____mpn_sqr_n_basecase (mp_ptr prodp, mp_srcptr up, mp_size_t size)
  239. #else
  240. ____mpn_sqr_n_basecase (prodp, up, size)
  241.      mp_ptr prodp;
  242.      mp_srcptr up;
  243.      mp_size_t size;
  244. #endif
  245. {
  246.   mp_size_t i;
  247.   mp_limb cy_limb;
  248.   mp_limb v_limb;
  249.  
  250.   /* Multiply by the first limb in V separately, as the result can be
  251.      stored (not added) to PROD.  We also avoid a loop for zeroing.  */
  252.   v_limb = up[0];
  253.   if (v_limb <= 1)
  254.     {
  255.       if (v_limb == 1)
  256.     MPN_COPY (prodp, up, size);
  257.       else
  258.     MPN_ZERO (prodp, size);
  259.       cy_limb = 0;
  260.     }
  261.   else
  262.     cy_limb = __mpn_mul_1 (prodp, up, size, v_limb);
  263.  
  264.   prodp[size] = cy_limb;
  265.   prodp++;
  266.  
  267.   /* For each iteration in the outer loop, multiply one limb from
  268.      U with one limb from V, and add it to PROD.  */
  269.   for (i = 1; i < size; i++)
  270.     {
  271.       v_limb = up[i];
  272.       if (v_limb <= 1)
  273.     {
  274.       cy_limb = 0;
  275.       if (v_limb == 1)
  276.         cy_limb = __mpn_add_n (prodp, prodp, up, size);
  277.     }
  278.       else
  279.     cy_limb = __mpn_addmul_1 (prodp, up, size, v_limb);
  280.  
  281.       prodp[size] = cy_limb;
  282.       prodp++;
  283.     }
  284. }
  285.  
  286. void
  287. #if __STDC__
  288. ____mpn_sqr_n (mp_ptr prodp,
  289.          mp_srcptr up, mp_size_t size, mp_ptr tspace)
  290. #else
  291. ____mpn_sqr_n (prodp, up, size, tspace)
  292.      mp_ptr prodp;
  293.      mp_srcptr up;
  294.      mp_size_t size;
  295.      mp_ptr tspace;
  296. #endif
  297. {
  298.   if ((size & 1) != 0)
  299.     {
  300.       /* The size is odd, the code code below doesn't handle that.
  301.      Multiply the least significant (size - 1) limbs with a recursive
  302.      call, and handle the most significant limb of S1 and S2
  303.      separately.  */
  304.       /* A slightly faster way to do this would be to make the Karatsuba
  305.      code below behave as if the size were even, and let it check for
  306.      odd size in the end.  I.e., in essence move this code to the end.
  307.      Doing so would save us a recursive call, and potentially make the
  308.      stack grow a lot less.  */
  309.  
  310.       mp_size_t esize = size - 1;    /* even size */
  311.       mp_limb cy_limb;
  312.  
  313.       MPN_SQR_N_RECURSE (prodp, up, esize, tspace);
  314.       cy_limb = __mpn_addmul_1 (prodp + esize, up, esize, up[esize]);
  315.       prodp[esize + esize] = cy_limb;
  316.       cy_limb = __mpn_addmul_1 (prodp + esize, up, size, up[esize]);
  317.  
  318.       prodp[esize + size] = cy_limb;
  319.     }
  320.   else
  321.     {
  322.       mp_size_t hsize = size >> 1;
  323.       mp_limb cy;
  324.  
  325.       /*** Product H.     ________________  ________________
  326.             |_____U1 x U1____||____U0 x U0_____|  */
  327.       /* Put result in upper part of PROD and pass low part of TSPACE
  328.      as new TSPACE.  */
  329.       MPN_SQR_N_RECURSE (prodp + size, up + hsize, hsize, tspace);
  330.  
  331.       /*** Product M.     ________________
  332.             |_(U1-U0)(U0-U1)_|  */
  333.       if (__mpn_cmp (up + hsize, up, hsize) >= 0)
  334.     {
  335.       __mpn_sub_n (prodp, up + hsize, up, hsize);
  336.     }
  337.       else
  338.     {
  339.       __mpn_sub_n (prodp, up, up + hsize, hsize);
  340.     }
  341.  
  342.       /* Read temporary operands from low part of PROD.
  343.      Put result in low part of TSPACE using upper part of TSPACE
  344.      as new TSPACE.  */
  345.       MPN_SQR_N_RECURSE (tspace, prodp, hsize, tspace + size);
  346.  
  347.       /*** Add/copy product H.  */
  348.       MPN_COPY (prodp + hsize, prodp + size, hsize);
  349.       cy = __mpn_add_n (prodp + size, prodp + size, prodp + size + hsize, hsize);
  350.  
  351.       /*** Add product M (if NEGFLG M is a negative number).  */
  352.       cy -= __mpn_sub_n (prodp + hsize, prodp + hsize, tspace, size);
  353.  
  354.       /*** Product L.     ________________  ________________
  355.             |________________||____U0 x U0_____|  */
  356.       /* Read temporary operands from low part of PROD.
  357.      Put result in low part of TSPACE using upper part of TSPACE
  358.      as new TSPACE.  */
  359.       MPN_SQR_N_RECURSE (tspace, up, hsize, tspace + size);
  360.  
  361.       /*** Add/copy Product L (twice).  */
  362.  
  363.       cy += __mpn_add_n (prodp + hsize, prodp + hsize, tspace, size);
  364.       if (cy)
  365.     {
  366.       if (cy > 0)
  367.         __mpn_add_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
  368.       else
  369.         {
  370.           __mpn_sub_1 (prodp + hsize + size, prodp + hsize + size, hsize, cy);
  371.           abort ();
  372.         }
  373.     }
  374.  
  375.       MPN_COPY (prodp, tspace, hsize);
  376.       cy = __mpn_add_n (prodp + hsize, prodp + hsize, tspace + hsize, hsize);
  377.       if (cy)
  378.     __mpn_add_1 (prodp + size, prodp + size, size, 1);
  379.     }
  380. }
  381.  
  382. /* This should be made into an inline function in gmp.h.  */
  383. inline void
  384. #if __STDC__
  385. __mpn_mul_n (mp_ptr prodp, mp_srcptr up, mp_srcptr vp, mp_size_t size)
  386. #else
  387. __mpn_mul_n (prodp, up, vp, size)
  388.      mp_ptr prodp;
  389.      mp_srcptr up;
  390.      mp_srcptr vp;
  391.      mp_size_t size;
  392. #endif
  393. {
  394.   if (up == vp)
  395.     {
  396.       if (size < KARATSUBA_THRESHOLD)
  397.     {
  398.       ____mpn_sqr_n_basecase (prodp, up, size);
  399.     }
  400.       else
  401.     {
  402.       mp_ptr tspace;
  403.       tspace = (mp_ptr) alloca (2 * size * BYTES_PER_MP_LIMB);
  404.       ____mpn_sqr_n (prodp, up, size, tspace);
  405.     }
  406.     }
  407.   else
  408.     {
  409.       if (size < KARATSUBA_THRESHOLD)
  410.     {
  411.       ____mpn_mul_n_basecase (prodp, up, vp, size);
  412.     }
  413.       else
  414.     {
  415.       mp_ptr tspace;
  416.       tspace = (mp_ptr) alloca (2 * size * BYTES_PER_MP_LIMB);
  417.       ____mpn_mul_n (prodp, up, vp, size, tspace);
  418.     }
  419.     }
  420. }
  421.