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 } 172