source/libs/gmp/gmp-src/mpz/aorsmul.c

   1 /* mpz_addmul, mpz_submul -- add or subtract multiple.
   2
   3 Copyright 2001, 2004, 2005, 2012 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 either:
   9
  10   * the GNU Lesser General Public License as published by the Free
  11     Software Foundation; either version 3 of the License, or (at your
  12     option) any later version.
  13
  14 or
  15
  16   * the GNU General Public License as published by the Free Software
  17     Foundation; either version 2 of the License, or (at your option) any
  18     later version.
  19
  20 or both in parallel, as here.
  21
  22 The GNU MP Library is distributed in the hope that it will be useful, but
  23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  24 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  25 for more details.
  26
  27 You should have received copies of the GNU General Public License and the
  28 GNU Lesser General Public License along with the GNU MP Library.  If not,
  29 see https://www.gnu.org/licenses/.  */
  30
  31 #include "gmp.h"
  32 #include "gmp-impl.h"
  33
  34
  35 /* expecting x and y both with non-zero high limbs */
  36 #define mpn_cmp_twosizes_lt(xp,xsize, yp,ysize)                 \
  37   ((xsize) < (ysize)                                            \
  38    || ((xsize) == (ysize) && mpn_cmp (xp, yp, xsize) < 0))
  39
  40
  41 /* sub>=0 means an addmul w += x*y, sub<0 means a submul w -= x*y.
  42
  43    The signs of w, x and y are fully accounted for by each flipping "sub".
  44
  45    The sign of w is retained for the result, unless the absolute value
  46    submul underflows, in which case it flips.  */
  47
  48 static void __gmpz_aorsmul (REGPARM_3_1 (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub)) REGPARM_ATTR (1);
  49 #define mpz_aorsmul(w,x,y,sub)  __gmpz_aorsmul (REGPARM_3_1 (w, x, y, sub))
  50
  51 REGPARM_ATTR (1) static void
  52 mpz_aorsmul (mpz_ptr w, mpz_srcptr x, mpz_srcptr y, mp_size_t sub)
  53 {
  54   mp_size_t  xsize, ysize, tsize, wsize, wsize_signed;
  55   mp_ptr     wp, tp;
  56   mp_limb_t  c, high;
  57   TMP_DECL;
  58
  59   /* w unaffected if x==0 or y==0 */
  60   xsize = SIZ(x);
  61   ysize = SIZ(y);
  62   if (xsize == 0 || ysize == 0)
  63     return;
  64
  65   /* make x the bigger of the two */
  66   if (ABS(ysize) > ABS(xsize))
  67     {
  68       MPZ_SRCPTR_SWAP (x, y);
  69       MP_SIZE_T_SWAP (xsize, ysize);
  70     }
  71
  72   sub ^= ysize;
  73   ysize = ABS(ysize);
  74
  75   /* use mpn_addmul_1/mpn_submul_1 if possible */
  76   if (ysize == 1)
  77     {
  78       mpz_aorsmul_1 (w, x, PTR(y)[0], sub);
  79       return;
  80     }
  81
  82   sub ^= xsize;
  83   xsize = ABS(xsize);
  84
  85   wsize_signed = SIZ(w);
  86   sub ^= wsize_signed;
  87   wsize = ABS(wsize_signed);
  88
  89   tsize = xsize + ysize;
  90   wp = MPZ_REALLOC (w, MAX (wsize, tsize) + 1);
  91
  92   if (wsize_signed == 0)
  93     {
  94       /* Nothing to add to, just set w=x*y.  No w==x or w==y overlap here,
  95          since we know x,y!=0 but w==0.  */
  96       high = mpn_mul (wp, PTR(x),xsize, PTR(y),ysize);
  97       tsize -= (high == 0);
  98       SIZ(w) = (sub >= 0 ? tsize : -tsize);
  99       return;
 100     }
 101
 102   TMP_MARK;
 103   tp = TMP_ALLOC_LIMBS (tsize);
 104
 105   high = mpn_mul (tp, PTR(x),xsize, PTR(y),ysize);
 106   tsize -= (high == 0);
 107   ASSERT (tp[tsize-1] != 0);
 108   if (sub >= 0)
 109     {
 110       mp_srcptr up    = wp;
 111       mp_size_t usize = wsize;
 112
 113       if (usize < tsize)
 114         {
 115           up    = tp;
 116           usize = tsize;
 117           tp    = wp;
 118           tsize = wsize;
 119
 120           wsize = usize;
 121         }
 122
 123       c = mpn_add (wp, up,usize, tp,tsize);
 124       wp[wsize] = c;
 125       wsize += (c != 0);
 126     }
 127   else
 128     {
 129       mp_srcptr up    = wp;
 130       mp_size_t usize = wsize;
 131
 132       if (mpn_cmp_twosizes_lt (up,usize, tp,tsize))
 133         {
 134           up    = tp;
 135           usize = tsize;
 136           tp    = wp;
 137           tsize = wsize;
 138
 139           wsize = usize;
 140           wsize_signed = -wsize_signed;
 141         }
 142
 143       ASSERT_NOCARRY (mpn_sub (wp, up,usize, tp,tsize));
 144       wsize = usize;
 145       MPN_NORMALIZE (wp, wsize);
 146     }
 147
 148   SIZ(w) = (wsize_signed >= 0 ? wsize : -wsize);
 149
 150   TMP_FREE;
 151 }
 152
 153
 154 void
 155 mpz_addmul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
 156 {
 157   mpz_aorsmul (w, u, v, (mp_size_t) 0);
 158 }
 159
 160 void
 161 mpz_submul (mpz_ptr w, mpz_srcptr u, mpz_srcptr v)
 162 {
 163   mpz_aorsmul (w, u, v, (mp_size_t) -1);
 164 }