github.com/aergoio/aergo@v1.3.1/libtool/src/gmp-6.1.2/mpf/mul_ui.c (about)

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