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

数学计算

开发平台:

Unix_Linux

  1. /* mpz_addmul_ui, mpz_submul_ui - add or subtract small multiple.
  2.    THE mpz_aorsmul_1 FUNCTION IN THIS FILE IS FOR INTERNAL USE ONLY AND IS
  3.    ALMOST CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR
  4.    COMPLETELY IN FUTURE GNU MP RELEASES.
  5. Copyright 2001, 2002, 2004, 2005 Free Software Foundation, Inc.
  6. This file is part of the GNU MP Library.
  7. The GNU MP Library is free software; you can redistribute it and/or modify
  8. it under the terms of the GNU Lesser General Public License as published by
  9. the Free Software Foundation; either version 3 of the License, or (at your
  10. option) any later version.
  11. The GNU MP Library is distributed in the hope that it will be useful, but
  12. WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  13. or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  14. License for more details.
  15. You should have received a copy of the GNU Lesser General Public License
  16. along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
  17. #include "gmp.h"
  18. #include "gmp-impl.h"
  19. #if HAVE_NATIVE_mpn_mul_1c
  20. #define MPN_MUL_1C(cout, dst, src, size, n, cin)        
  21.   do {                                                  
  22.     (cout) = mpn_mul_1c (dst, src, size, n, cin);       
  23.   } while (0)
  24. #else
  25. #define MPN_MUL_1C(cout, dst, src, size, n, cin)        
  26.   do {                                                  
  27.     mp_limb_t __cy;                                     
  28.     __cy = mpn_mul_1 (dst, src, size, n);               
  29.     (cout) = __cy + mpn_add_1 (dst, dst, size, cin);    
  30.   } while (0)
  31. #endif
  32. /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
  33.    All that's needed to account for negative w or x is to flip "sub".
  34.    The final w will retain its sign, unless an underflow occurs in a submul
  35.    of absolute values, in which case it's flipped.
  36.    If x has more limbs than w, then mpn_submul_1 followed by mpn_com is
  37.    used.  The alternative would be mpn_mul_1 into temporary space followed
  38.    by mpn_sub_n.  Avoiding temporary space seem good, and submul+com stands
  39.    a chance of being faster since it involves only one set of carry
  40.    propagations, not two.  Note that doing an addmul_1 with a
  41.    twos-complement negative y doesn't work, because it effectively adds an
  42.    extra x * 2^GMP_LIMB_BITS.  */
  43. REGPARM_ATTR(1) void
  44. mpz_aorsmul_1 (mpz_ptr w, mpz_srcptr x, mp_limb_t y, mp_size_t sub)
  45. {
  46.   mp_size_t  xsize, wsize, wsize_signed, new_wsize, min_size, dsize;
  47.   mp_srcptr  xp;
  48.   mp_ptr     wp;
  49.   mp_limb_t  cy;
  50.   /* w unaffected if x==0 or y==0 */
  51.   xsize = SIZ (x);
  52.   if (xsize == 0 || y == 0)
  53.     return;
  54.   sub ^= xsize;
  55.   xsize = ABS (xsize);
  56.   wsize_signed = SIZ (w);
  57.   if (wsize_signed == 0)
  58.     {
  59.       /* nothing to add to, just set x*y, "sub" gives the sign */
  60.       MPZ_REALLOC (w, xsize+1);
  61.       wp = PTR (w);
  62.       cy = mpn_mul_1 (wp, PTR(x), xsize, y);
  63.       wp[xsize] = cy;
  64.       xsize += (cy != 0);
  65.       SIZ (w) = (sub >= 0 ? xsize : -xsize);
  66.       return;
  67.     }
  68.   sub ^= wsize_signed;
  69.   wsize = ABS (wsize_signed);
  70.   new_wsize = MAX (wsize, xsize);
  71.   MPZ_REALLOC (w, new_wsize+1);
  72.   wp = PTR (w);
  73.   xp = PTR (x);
  74.   min_size = MIN (wsize, xsize);
  75.   if (sub >= 0)
  76.     {
  77.       /* addmul of absolute values */
  78.       cy = mpn_addmul_1 (wp, xp, min_size, y);
  79.       wp += min_size;
  80.       xp += min_size;
  81.       dsize = xsize - wsize;
  82. #if HAVE_NATIVE_mpn_mul_1c
  83.       if (dsize > 0)
  84.         cy = mpn_mul_1c (wp, xp, dsize, y, cy);
  85.       else if (dsize < 0)
  86.         {
  87.           dsize = -dsize;
  88.           cy = mpn_add_1 (wp, wp, dsize, cy);
  89.         }
  90. #else
  91.       if (dsize != 0)
  92.         {
  93.           mp_limb_t  cy2;
  94.           if (dsize > 0)
  95.             cy2 = mpn_mul_1 (wp, xp, dsize, y);
  96.           else
  97.             {
  98.               dsize = -dsize;
  99.               cy2 = 0;
  100.             }
  101.           cy = cy2 + mpn_add_1 (wp, wp, dsize, cy);
  102.         }
  103. #endif
  104.       wp[dsize] = cy;
  105.       new_wsize += (cy != 0);
  106.     }
  107.   else
  108.     {
  109.       /* submul of absolute values */
  110.       cy = mpn_submul_1 (wp, xp, min_size, y);
  111.       if (wsize >= xsize)
  112.         {
  113.           /* if w bigger than x, then propagate borrow through it */
  114.           if (wsize != xsize)
  115.             cy = mpn_sub_1 (wp+xsize, wp+xsize, wsize-xsize, cy);
  116.           if (cy != 0)
  117.             {
  118.               /* Borrow out of w, take twos complement negative to get
  119.                  absolute value, flip sign of w.  */
  120.               wp[new_wsize] = ~-cy;  /* extra limb is 0-cy */
  121.               mpn_com (wp, wp, new_wsize);
  122.               new_wsize++;
  123.               MPN_INCR_U (wp, new_wsize, CNST_LIMB(1));
  124.               wsize_signed = -wsize_signed;
  125.             }
  126.         }
  127.       else /* wsize < xsize */
  128.         {
  129.           /* x bigger than w, so want x*y-w.  Submul has given w-x*y, so
  130.              take twos complement and use an mpn_mul_1 for the rest.  */
  131.           mp_limb_t  cy2;
  132.           /* -(-cy*b^n + w-x*y) = (cy-1)*b^n + ~(w-x*y) + 1 */
  133.           mpn_com (wp, wp, wsize);
  134.           cy += mpn_add_1 (wp, wp, wsize, CNST_LIMB(1));
  135.           cy -= 1;
  136.           /* If cy-1 == -1 then hold that -1 for latter.  mpn_submul_1 never
  137.              returns cy==MP_LIMB_T_MAX so that value always indicates a -1. */
  138.           cy2 = (cy == MP_LIMB_T_MAX);
  139.           cy += cy2;
  140.           MPN_MUL_1C (cy, wp+wsize, xp+wsize, xsize-wsize, y, cy);
  141.           wp[new_wsize] = cy;
  142.           new_wsize += (cy != 0);
  143.           /* Apply any -1 from above.  The value at wp+wsize is non-zero
  144.              because y!=0 and the high limb of x will be non-zero.  */
  145.           if (cy2)
  146.             MPN_DECR_U (wp+wsize, new_wsize-wsize, CNST_LIMB(1));
  147.           wsize_signed = -wsize_signed;
  148.         }
  149.       /* submul can produce high zero limbs due to cancellation, both when w
  150.          has more limbs or x has more  */
  151.       MPN_NORMALIZE (wp, new_wsize);
  152.     }
  153.   SIZ (w) = (wsize_signed >= 0 ? new_wsize : -new_wsize);
  154.   ASSERT (new_wsize == 0 || PTR(w)[new_wsize-1] != 0);
  155. }
  156. void
  157. mpz_addmul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
  158. {
  159. #if BITS_PER_ULONG > GMP_NUMB_BITS
  160.   if (UNLIKELY (y > GMP_NUMB_MAX && SIZ(x) != 0))
  161.     {
  162.       mpz_t t;
  163.       mp_ptr tp;
  164.       mp_size_t xn;
  165.       TMP_DECL;
  166.       TMP_MARK;
  167.       xn = SIZ (x);
  168.       MPZ_TMP_INIT (t, ABS (xn) + 1);
  169.       tp = PTR (t);
  170.       tp[0] = 0;
  171.       MPN_COPY (tp + 1, PTR(x), ABS (xn));
  172.       SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
  173.       mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) 0);
  174.       PTR(t) = tp + 1;
  175.       SIZ(t) = xn;
  176.       mpz_aorsmul_1 (w, t, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) 0);
  177.       TMP_FREE;
  178.       return;
  179.     }
  180. #endif
  181.   mpz_aorsmul_1 (w, x, (mp_limb_t) y, (mp_size_t) 0);
  182. }
  183. void
  184. mpz_submul_ui (mpz_ptr w, mpz_srcptr x, unsigned long y)
  185. {
  186. #if BITS_PER_ULONG > GMP_NUMB_BITS
  187.   if (y > GMP_NUMB_MAX && SIZ(x) != 0)
  188.     {
  189.       mpz_t t;
  190.       mp_ptr tp;
  191.       mp_size_t xn;
  192.       TMP_DECL;
  193.       TMP_MARK;
  194.       xn = SIZ (x);
  195.       MPZ_TMP_INIT (t, ABS (xn) + 1);
  196.       tp = PTR (t);
  197.       tp[0] = 0;
  198.       MPN_COPY (tp + 1, PTR(x), ABS (xn));
  199.       SIZ(t) = xn >= 0 ? xn + 1 : xn - 1;
  200.       mpz_aorsmul_1 (w, t, (mp_limb_t) y >> GMP_NUMB_BITS, (mp_size_t) -1);
  201.       PTR(t) = tp + 1;
  202.       SIZ(t) = xn;
  203.       mpz_aorsmul_1 (w, t, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) -1);
  204.       TMP_FREE;
  205.       return;
  206.     }
  207. #endif
  208.   mpz_aorsmul_1 (w, x, (mp_limb_t) y & GMP_NUMB_MASK, (mp_size_t) -1);
  209. }