contrib/gmp/mpf/mul_ui.c

   1 /* mpf_mul_ui -- Multiply a float and an unsigned integer.
   2
   3 Copyright 1993, 1994, 1996, 2001, 2003, 2004 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 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
  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 Lesser General Public
  15 License for more details.
  16
  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
  20 #include "gmp.h"
  21 #include "gmp-impl.h"
  22 #include "longlong.h"
  23
  24
  25 /* The core operation is a multiply of PREC(r) limbs from u by v, producing
  26    either PREC(r) or PREC(r)+1 result limbs.  If u is shorter than PREC(r),
  27    then we take only as much as it has.  If u is longer we incorporate a
  28    carry from the lower limbs.
  29
  30    If u has just 1 extra limb, then the carry to add is high(up[0]*v).  That
  31    is of course what mpn_mul_1 would do if it was called with PREC(r)+1
  32    limbs of input.
  33
  34    If u has more than 1 extra limb, then there can be a further carry bit
  35    out of lower uncalculated limbs (the way the low of one product adds to
  36    the high of the product below it).  This is of course what an mpn_mul_1
  37    would do if it was called with the full u operand.  But we instead work
  38    downwards explicitly, until a carry occurs or until a value other than
  39    GMP_NUMB_MAX occurs (that being the only value a carry bit can propagate
  40    across).
  41
  42    The carry determination normally requires two umul_ppmm's, only rarely
  43    will GMP_NUMB_MAX occur and require further products.
  44
  45    The carry limb is conveniently added into the mul_1 using mpn_mul_1c when
  46    that function exists, otherwise a subsequent mpn_add_1 is needed.
  47
  48    Clearly when mpn_mul_1c is used the carry must be calculated first.  But
  49    this is also the case when add_1 is used, since if r==u and ABSIZ(r) >
  50    PREC(r) then the mpn_mul_1 overwrites the low part of the input.
  51
  52    A reuse r==u with size > prec can occur from a size PREC(r)+1 in the
  53    usual way, or it can occur from an mpf_set_prec_raw leaving a bigger
  54    sized value.  In both cases we can end up calling mpn_mul_1 with
  55    overlapping src and dst regions, but this will be with dst < src and such
  56    an overlap is permitted.
  57
  58    Not done:
  59
  60    No attempt is made to determine in advance whether the result will be
  61    PREC(r) or PREC(r)+1 limbs.  If it's going to be PREC(r)+1 then we could
  62    take one less limb from u and generate just PREC(r), that of course
  63    satisfying application requested precision.  But any test counting bits
  64    or forming the high product would almost certainly take longer than the
  65    incremental cost of an extra limb in mpn_mul_1.
  66
  67    Enhancements:
  68
  69    Repeated mpf_mul_ui's with an even v will accumulate low zero bits on the
  70    result, leaving low zero limbs after a while, which it might be nice to
  71    strip to save work in subsequent operations.  Calculating the low limb
  72    explicitly would let us direct mpn_mul_1 to put the balance at rp when
  73    the low is zero (instead of normally rp+1).  But it's not clear whether
  74    this would be worthwhile.  Explicit code for the low limb will probably
  75    be slower than having it done in mpn_mul_1, so we need to consider how
  76    often a zero will be stripped and how much that's likely to save
  77    later.  */
  78
  79 void
  80 mpf_mul_ui (mpf_ptr r, mpf_srcptr u, unsigned long int v)
  81 {
  82   mp_srcptr up;
  83   mp_size_t usize;
  84   mp_size_t size;
  85   mp_size_t prec, excess;
  86   mp_limb_t cy_limb, vl, cbit, cin;
  87   mp_ptr rp;
  88
  89   usize = u->_mp_size;
  90   if (UNLIKELY (v == 0) || UNLIKELY (usize == 0))
  91     {
  92       r->_mp_size = 0;
  93       r->_mp_exp = 0;
  94       return;
  95     }
  96
  97 #if BITS_PER_ULONG > GMP_NUMB_BITS  /* avoid warnings about shift amount */
  98   if (v > GMP_NUMB_MAX)
  99     {
 100       mpf_t     vf;
 101       mp_limb_t vp[2];
 102       vp[0] = v & GMP_NUMB_MASK;
 103       vp[1] = v >> GMP_NUMB_BITS;
 104       PTR(vf) = vp;
 105       SIZ(vf) = 2;
 106       ASSERT_CODE (PREC(vf) = 2);
 107       EXP(vf) = 2;
 108       mpf_mul (r, u, vf);
 109       return;
 110     }
 111 #endif
 112
 113   size = ABS (usize);
 114   prec = r->_mp_prec;
 115   up = u->_mp_d;
 116   vl = v;
 117   excess = size - prec;
 118   cin = 0;
 119
 120   if (excess > 0)
 121     {
 122       /* up is bigger than desired rp, shorten it to prec limbs and
 123          determine a carry-in */
 124
 125       mp_limb_t  vl_shifted = vl << GMP_NAIL_BITS;
 126       mp_limb_t  hi, lo, next_lo, sum;
 127       mp_size_t  i;
 128
 129       /* high limb of top product */
 130       i = excess - 1;
 131       umul_ppmm (cin, lo, up[i], vl_shifted);
 132
 133       /* and carry bit out of products below that, if any */
 134       for (;;)
 135         {
 136           i--;
 137           if (i < 0)
 138             break;
 139
 140           umul_ppmm (hi, next_lo, up[i], vl_shifted);
 141           lo >>= GMP_NAIL_BITS;
 142           ADDC_LIMB (cbit, sum, hi, lo);
 143           cin += cbit;
 144           lo = next_lo;
 145
 146           /* Continue only if the sum is GMP_NUMB_MAX.  GMP_NUMB_MAX is the
 147              only value a carry from below can propagate across.  If we've
 148              just seen the carry out (ie. cbit!=0) then sum!=GMP_NUMB_MAX,
 149              so this test stops us for that case too.  */
 150           if (LIKELY (sum != GMP_NUMB_MAX))
 151             break;
 152         }
 153
 154       up += excess;
 155       size = prec;
 156     }
 157
 158   rp = r->_mp_d;
 159 #if HAVE_NATIVE_mpn_mul_1c
 160   cy_limb = mpn_mul_1c (rp, up, size, vl, cin);
 161 #else
 162   cy_limb = mpn_mul_1 (rp, up, size, vl);
 163   __GMPN_ADD_1 (cbit, rp, rp, size, cin);
 164   cy_limb += cbit;
 165 #endif
 166   rp[size] = cy_limb;
 167   cy_limb = cy_limb != 0;
 168   r->_mp_exp = u->_mp_exp + cy_limb;
 169   size += cy_limb;
 170   r->_mp_size = usize >= 0 ? size : -size;
 171 }