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

     1  /* mpf_add -- Add two floats.
     2  
     3  Copyright 1993, 1994, 1996, 2000, 2001, 2005 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  void
    35  mpf_add (mpf_ptr r, mpf_srcptr u, mpf_srcptr v)
    36  {
    37    mp_srcptr up, vp;
    38    mp_ptr rp, tp;
    39    mp_size_t usize, vsize, rsize;
    40    mp_size_t prec;
    41    mp_exp_t uexp;
    42    mp_size_t ediff;
    43    mp_limb_t cy;
    44    int negate;
    45    TMP_DECL;
    46  
    47    usize = u->_mp_size;
    48    vsize = v->_mp_size;
    49  
    50    /* Handle special cases that don't work in generic code below.  */
    51    if (usize == 0)
    52      {
    53      set_r_v_maybe:
    54        if (r != v)
    55          mpf_set (r, v);
    56        return;
    57      }
    58    if (vsize == 0)
    59      {
    60        v = u;
    61        goto set_r_v_maybe;
    62      }
    63  
    64    /* If signs of U and V are different, perform subtraction.  */
    65    if ((usize ^ vsize) < 0)
    66      {
    67        __mpf_struct v_negated;
    68        v_negated._mp_size = -vsize;
    69        v_negated._mp_exp = v->_mp_exp;
    70        v_negated._mp_d = v->_mp_d;
    71        mpf_sub (r, u, &v_negated);
    72        return;
    73      }
    74  
    75    TMP_MARK;
    76  
    77    /* Signs are now known to be the same.  */
    78    negate = usize < 0;
    79  
    80    /* Make U be the operand with the largest exponent.  */
    81    if (u->_mp_exp < v->_mp_exp)
    82      {
    83        mpf_srcptr t;
    84        t = u; u = v; v = t;
    85        usize = u->_mp_size;
    86        vsize = v->_mp_size;
    87      }
    88  
    89    usize = ABS (usize);
    90    vsize = ABS (vsize);
    91    up = u->_mp_d;
    92    vp = v->_mp_d;
    93    rp = r->_mp_d;
    94    prec = r->_mp_prec;
    95    uexp = u->_mp_exp;
    96    ediff = u->_mp_exp - v->_mp_exp;
    97  
    98    /* If U extends beyond PREC, ignore the part that does.  */
    99    if (usize > prec)
   100      {
   101        up += usize - prec;
   102        usize = prec;
   103      }
   104  
   105    /* If V extends beyond PREC, ignore the part that does.
   106       Note that this may make vsize negative.  */
   107    if (vsize + ediff > prec)
   108      {
   109        vp += vsize + ediff - prec;
   110        vsize = prec - ediff;
   111      }
   112  
   113  #if 0
   114    /* Locate the least significant non-zero limb in (the needed parts
   115       of) U and V, to simplify the code below.  */
   116    while (up[0] == 0)
   117      up++, usize--;
   118    while (vp[0] == 0)
   119      vp++, vsize--;
   120  #endif
   121  
   122    /* Allocate temp space for the result.  Allocate
   123       just vsize + ediff later???  */
   124    tp = TMP_ALLOC_LIMBS (prec);
   125  
   126    if (ediff >= prec)
   127      {
   128        /* V completely cancelled.  */
   129        if (rp != up)
   130  	MPN_COPY_INCR (rp, up, usize);
   131        rsize = usize;
   132      }
   133    else
   134      {
   135        /* uuuu     |  uuuu     |  uuuu     |  uuuu     |  uuuu    */
   136        /* vvvvvvv  |  vv       |    vvvvv  |    v      |       vv */
   137  
   138        if (usize > ediff)
   139  	{
   140  	  /* U and V partially overlaps.  */
   141  	  if (vsize + ediff <= usize)
   142  	    {
   143  	      /* uuuu     */
   144  	      /*   v      */
   145  	      mp_size_t size;
   146  	      size = usize - ediff - vsize;
   147  	      MPN_COPY (tp, up, size);
   148  	      cy = mpn_add (tp + size, up + size, usize - size, vp, vsize);
   149  	      rsize = usize;
   150  	    }
   151  	  else
   152  	    {
   153  	      /* uuuu     */
   154  	      /*   vvvvv  */
   155  	      mp_size_t size;
   156  	      size = vsize + ediff - usize;
   157  	      MPN_COPY (tp, vp, size);
   158  	      cy = mpn_add (tp + size, up, usize, vp + size, usize - ediff);
   159  	      rsize = vsize + ediff;
   160  	    }
   161  	}
   162        else
   163  	{
   164  	  /* uuuu     */
   165  	  /*      vv  */
   166  	  mp_size_t size;
   167  	  size = vsize + ediff - usize;
   168  	  MPN_COPY (tp, vp, vsize);
   169  	  MPN_ZERO (tp + vsize, ediff - usize);
   170  	  MPN_COPY (tp + size, up, usize);
   171  	  cy = 0;
   172  	  rsize = size + usize;
   173  	}
   174  
   175        MPN_COPY (rp, tp, rsize);
   176        rp[rsize] = cy;
   177        rsize += cy;
   178        uexp += cy;
   179      }
   180  
   181    r->_mp_size = negate ? -rsize : rsize;
   182    r->_mp_exp = uexp;
   183    TMP_FREE;
   184  }