rts/gmp/mpz/addmul_ui.c

   1 /* mpz_addmul_ui(prodsum, multiplier, small_multiplicand) --
   2    Add MULTIPLICATOR times SMALL_MULTIPLICAND to PRODSUM.
   3
   4 Copyright (C) 1997, 2000 Free Software Foundation, Inc.
   5
   6 This file is part of the GNU MP Library.
   7
   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 2.1 of the License, or (at your
  11 option) any later version.
  12
  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
  18 You should have received a copy of the GNU Lesser General Public License
  19 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  20 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  21 MA 02111-1307, USA. */
  22
  23 #include "gmp.h"
  24 #include "gmp-impl.h"
  25
  26 static mp_limb_t mpn_neg1 _PROTO ((mp_ptr, mp_size_t));
  27
  28 #if 0
  29 #undef  MPN_NORMALIZE
  30 #define MPN_NORMALIZE(DST, NLIMBS) \
  31   do {                                                                  \
  32     while (--(NLIMBS) >= 0 && (DST)[NLIMBS] == 0)                       \
  33       ;                                                                 \
  34     (NLIMBS)++;                                                         \
  35   } while (0)
  36 #undef  MPN_NORMALIZE_NOT_ZERO
  37 #define MPN_NORMALIZE_NOT_ZERO(DST, NLIMBS) \
  38   do {                                                                  \
  39     while ((DST)[--(NLIMBS)] == 0)                                      \
  40       ;                                                                 \
  41     (NLIMBS)++;                                                         \
  42   } while (0)
  43 #endif
  44
  45 void
  46 #if __STDC__
  47 mpz_addmul_ui (mpz_ptr rz, mpz_srcptr az, unsigned long int bu)
  48 #else
  49 mpz_addmul_ui (rz, az, bu)
  50      mpz_ptr rz;
  51      mpz_srcptr az;
  52      unsigned long int bu;
  53 #endif
  54 {
  55   mp_size_t rn, an;
  56   mp_ptr rp, ap;
  57
  58   an = SIZ (az);
  59
  60   /* If either multiplier is zero, result is unaffected.  */
  61   if (bu == 0 || an == 0)
  62     return;
  63
  64   rn = SIZ (rz);
  65
  66   if (rn == 0)
  67     {
  68       mp_limb_t cy;
  69
  70       an = ABS (an);
  71       if (ALLOC (rz) <= an)
  72         _mpz_realloc (rz, an + 1);
  73       rp = PTR (rz);
  74       ap = PTR (az);
  75       cy = mpn_mul_1 (rp, ap, an, (mp_limb_t) bu);
  76       rp[an] = cy;
  77       an += cy != 0;
  78       SIZ (rz) = SIZ (az) >= 0 ? an : -an;
  79       return;
  80     }
  81
  82   if ((an ^ rn) >= 0)
  83     {
  84       /* Sign of operands are the same--really add.  */
  85       an = ABS (an);
  86       rn = ABS (rn);
  87       if (rn > an)
  88         {
  89           mp_limb_t cy;
  90           if (ALLOC (rz) <= rn)
  91             _mpz_realloc (rz, rn + 1);
  92           rp = PTR (rz);
  93           ap = PTR (az);
  94           cy = mpn_addmul_1 (rp, ap, an, (mp_limb_t) bu);
  95           cy = mpn_add_1 (rp + an, rp + an, rn - an, cy);
  96           rp[rn] = cy;
  97           rn += cy != 0;
  98           SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn;
  99           return;
 100         }
 101       else
 102         {
 103           mp_limb_t cy;
 104           if (ALLOC (rz) <= an)
 105             _mpz_realloc (rz, an + 1);
 106           rp = PTR (rz);
 107           ap = PTR (az);
 108           cy = mpn_addmul_1 (rp, ap, rn, (mp_limb_t) bu);
 109           if (an != rn)
 110             {
 111               mp_limb_t cy2;
 112               cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu);
 113               cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy);
 114             }
 115           rn = an;
 116           rp[rn] = cy;
 117           rn += cy != 0;
 118           SIZ (rz) = SIZ (rz) >= 0 ? rn : -rn;
 119           return;
 120         }
 121     }
 122   else
 123     {
 124       /* Sign of operands are different--actually subtract.  */
 125       an = ABS (an);
 126       rn = ABS (rn);
 127       if (rn > an)
 128         {
 129           mp_limb_t cy;
 130           rp = PTR (rz);
 131           ap = PTR (az);
 132           cy = mpn_submul_1 (rp, ap, an, (mp_limb_t) bu);
 133           cy = mpn_sub_1 (rp + an, rp + an, rn - an, cy);
 134           if (cy != 0)
 135             {
 136               mpn_neg1 (rp, rn);
 137               MPN_NORMALIZE_NOT_ZERO (rp, rn);
 138             }
 139           else
 140             {
 141               MPN_NORMALIZE (rp, rn);
 142               rn = -rn;
 143             }
 144
 145           SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn;
 146           return;
 147         }
 148       else
 149         {
 150           /* Tricky case.  We need to subtract an operand that might be larger
 151              than the minuend.  To avoid allocating temporary space, we compute
 152              a*b-r instead of r-a*b and then negate.  */
 153           mp_limb_t cy;
 154           if (ALLOC (rz) <= an)
 155             _mpz_realloc (rz, an + 1);
 156           rp = PTR (rz);
 157           ap = PTR (az);
 158           cy = mpn_submul_1 (rp, ap, rn, (mp_limb_t) bu);
 159           if (an != rn)
 160             {
 161               mp_limb_t cy2;
 162               cy -= mpn_neg1 (rp, rn);
 163               cy2 = mpn_mul_1 (rp + rn, ap + rn, an - rn, (mp_limb_t) bu);
 164               if (cy == ~(mp_limb_t) 0)
 165                 cy = cy2 - mpn_sub_1 (rp + rn, rp + rn, an - rn, (mp_limb_t) 1);
 166               else
 167                 cy = cy2 + mpn_add_1 (rp + rn, rp + rn, an - rn, cy);
 168               rp[an] = cy;
 169               rn = an + (cy != 0);
 170               rn -= rp[rn - 1] == 0;
 171             }
 172           else if (cy != 0)
 173             {
 174               cy -= mpn_neg1 (rp, rn);
 175               rp[an] = cy;
 176               rn = an + 1;
 177               MPN_NORMALIZE_NOT_ZERO (rp, rn);
 178             }
 179           else
 180             {
 181               rn = an;
 182               MPN_NORMALIZE (rp, rn);
 183               rn = -rn;
 184             }
 185
 186           SIZ (rz) = SIZ (rz) >= 0 ? -rn : rn;
 187           return;
 188         }
 189     }
 190 }
 191
 192 static mp_limb_t
 193 #if __STDC__
 194 mpn_neg1 (mp_ptr rp, mp_size_t rn)
 195 #else
 196 mpn_neg1 (rp, rn)
 197      mp_ptr rp;
 198      mp_size_t rn;
 199 #endif
 200 {
 201   mp_size_t i;
 202
 203   while (rn != 0 && rp[0] == 0)
 204     rp++, rn--;
 205
 206   if (rn != 0)
 207     {
 208       rp[0] = -rp[0];
 209       for (i = 1; i < rn; i++)
 210         rp[i] = ~rp[i];
 211       return 1;
 212     }
 213   return 0;
 214 }