linux/arch/parisc/math-emu/sfadd.c
<<
>>
Prefs
   1// SPDX-License-Identifier: GPL-2.0-or-later
   2/*
   3 * Linux/PA-RISC Project (http://www.parisc-linux.org/)
   4 *
   5 * Floating-point emulation code
   6 *  Copyright (C) 2001 Hewlett-Packard (Paul Bame) <bame@debian.org>
   7 */
   8/*
   9 * BEGIN_DESC
  10 *
  11 *  File:
  12 *      @(#)    pa/spmath/sfadd.c               $Revision: 1.1 $
  13 *
  14 *  Purpose:
  15 *      Single_add: add two single precision values.
  16 *
  17 *  External Interfaces:
  18 *      sgl_fadd(leftptr, rightptr, dstptr, status)
  19 *
  20 *  Internal Interfaces:
  21 *
  22 *  Theory:
  23 *      <<please update with a overview of the operation of this file>>
  24 *
  25 * END_DESC
  26*/
  27
  28
  29#include "float.h"
  30#include "sgl_float.h"
  31
  32/*
  33 * Single_add: add two single precision values.
  34 */
  35int
  36sgl_fadd(
  37    sgl_floating_point *leftptr,
  38    sgl_floating_point *rightptr,
  39    sgl_floating_point *dstptr,
  40    unsigned int *status)
  41    {
  42    register unsigned int left, right, result, extent;
  43    register unsigned int signless_upper_left, signless_upper_right, save;
  44    
  45    
  46    register int result_exponent, right_exponent, diff_exponent;
  47    register int sign_save, jumpsize;
  48    register boolean inexact = FALSE;
  49    register boolean underflowtrap;
  50        
  51    /* Create local copies of the numbers */
  52    left = *leftptr;
  53    right = *rightptr;
  54
  55    /* A zero "save" helps discover equal operands (for later),  *
  56     * and is used in swapping operands (if needed).             */
  57    Sgl_xortointp1(left,right,/*to*/save);
  58
  59    /*
  60     * check first operand for NaN's or infinity
  61     */
  62    if ((result_exponent = Sgl_exponent(left)) == SGL_INFINITY_EXPONENT)
  63        {
  64        if (Sgl_iszero_mantissa(left)) 
  65            {
  66            if (Sgl_isnotnan(right)) 
  67                {
  68                if (Sgl_isinfinity(right) && save!=0) 
  69                    {
  70                    /* 
  71                     * invalid since operands are opposite signed infinity's
  72                     */
  73                    if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  74                    Set_invalidflag();
  75                    Sgl_makequietnan(result);
  76                    *dstptr = result;
  77                    return(NOEXCEPTION);
  78                    }
  79                /*
  80                 * return infinity
  81                 */
  82                *dstptr = left;
  83                return(NOEXCEPTION);
  84                }
  85            }
  86        else 
  87            {
  88            /*
  89             * is NaN; signaling or quiet?
  90             */
  91            if (Sgl_isone_signaling(left)) 
  92                {
  93                /* trap if INVALIDTRAP enabled */
  94                if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
  95                /* make NaN quiet */
  96                Set_invalidflag();
  97                Sgl_set_quiet(left);
  98                }
  99            /* 
 100             * is second operand a signaling NaN? 
 101             */
 102            else if (Sgl_is_signalingnan(right)) 
 103                {
 104                /* trap if INVALIDTRAP enabled */
 105                if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 106                /* make NaN quiet */
 107                Set_invalidflag();
 108                Sgl_set_quiet(right);
 109                *dstptr = right;
 110                return(NOEXCEPTION);
 111                }
 112            /*
 113             * return quiet NaN
 114             */
 115            *dstptr = left;
 116            return(NOEXCEPTION);
 117            }
 118        } /* End left NaN or Infinity processing */
 119    /*
 120     * check second operand for NaN's or infinity
 121     */
 122    if (Sgl_isinfinity_exponent(right)) 
 123        {
 124        if (Sgl_iszero_mantissa(right)) 
 125            {
 126            /* return infinity */
 127            *dstptr = right;
 128            return(NOEXCEPTION);
 129            }
 130        /*
 131         * is NaN; signaling or quiet?
 132         */
 133        if (Sgl_isone_signaling(right)) 
 134            {
 135            /* trap if INVALIDTRAP enabled */
 136            if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 137            /* make NaN quiet */
 138            Set_invalidflag();
 139            Sgl_set_quiet(right);
 140            }
 141        /*
 142         * return quiet NaN
 143         */
 144        *dstptr = right;
 145        return(NOEXCEPTION);
 146        } /* End right NaN or Infinity processing */
 147
 148    /* Invariant: Must be dealing with finite numbers */
 149
 150    /* Compare operands by removing the sign */
 151    Sgl_copytoint_exponentmantissa(left,signless_upper_left);
 152    Sgl_copytoint_exponentmantissa(right,signless_upper_right);
 153
 154    /* sign difference selects add or sub operation. */
 155    if(Sgl_ismagnitudeless(signless_upper_left,signless_upper_right))
 156        {
 157        /* Set the left operand to the larger one by XOR swap *
 158         *  First finish the first word using "save"          */
 159        Sgl_xorfromintp1(save,right,/*to*/right);
 160        Sgl_xorfromintp1(save,left,/*to*/left);
 161        result_exponent = Sgl_exponent(left);
 162        }
 163    /* Invariant:  left is not smaller than right. */ 
 164
 165    if((right_exponent = Sgl_exponent(right)) == 0)
 166        {
 167        /* Denormalized operands.  First look for zeroes */
 168        if(Sgl_iszero_mantissa(right)) 
 169            {
 170            /* right is zero */
 171            if(Sgl_iszero_exponentmantissa(left))
 172                {
 173                /* Both operands are zeros */
 174                if(Is_rounding_mode(ROUNDMINUS))
 175                    {
 176                    Sgl_or_signs(left,/*with*/right);
 177                    }
 178                else
 179                    {
 180                    Sgl_and_signs(left,/*with*/right);
 181                    }
 182                }
 183            else 
 184                {
 185                /* Left is not a zero and must be the result.  Trapped
 186                 * underflows are signaled if left is denormalized.  Result
 187                 * is always exact. */
 188                if( (result_exponent == 0) && Is_underflowtrap_enabled() )
 189                    {
 190                    /* need to normalize results mantissa */
 191                    sign_save = Sgl_signextendedsign(left);
 192                    Sgl_leftshiftby1(left);
 193                    Sgl_normalize(left,result_exponent);
 194                    Sgl_set_sign(left,/*using*/sign_save);
 195                    Sgl_setwrapped_exponent(left,result_exponent,unfl);
 196                    *dstptr = left;
 197                    return(UNDERFLOWEXCEPTION);
 198                    }
 199                }
 200            *dstptr = left;
 201            return(NOEXCEPTION);
 202            }
 203
 204        /* Neither are zeroes */
 205        Sgl_clear_sign(right);  /* Exponent is already cleared */
 206        if(result_exponent == 0 )
 207            {
 208            /* Both operands are denormalized.  The result must be exact
 209             * and is simply calculated.  A sum could become normalized and a
 210             * difference could cancel to a true zero. */
 211            if( (/*signed*/int) save < 0 )
 212                {
 213                Sgl_subtract(left,/*minus*/right,/*into*/result);
 214                if(Sgl_iszero_mantissa(result))
 215                    {
 216                    if(Is_rounding_mode(ROUNDMINUS))
 217                        {
 218                        Sgl_setone_sign(result);
 219                        }
 220                    else
 221                        {
 222                        Sgl_setzero_sign(result);
 223                        }
 224                    *dstptr = result;
 225                    return(NOEXCEPTION);
 226                    }
 227                }
 228            else
 229                {
 230                Sgl_addition(left,right,/*into*/result);
 231                if(Sgl_isone_hidden(result))
 232                    {
 233                    *dstptr = result;
 234                    return(NOEXCEPTION);
 235                    }
 236                }
 237            if(Is_underflowtrap_enabled())
 238                {
 239                /* need to normalize result */
 240                sign_save = Sgl_signextendedsign(result);
 241                Sgl_leftshiftby1(result);
 242                Sgl_normalize(result,result_exponent);
 243                Sgl_set_sign(result,/*using*/sign_save);
 244                Sgl_setwrapped_exponent(result,result_exponent,unfl);
 245                *dstptr = result;
 246                return(UNDERFLOWEXCEPTION);
 247                }
 248            *dstptr = result;
 249            return(NOEXCEPTION);
 250            }
 251        right_exponent = 1;     /* Set exponent to reflect different bias
 252                                 * with denomalized numbers. */
 253        }
 254    else
 255        {
 256        Sgl_clear_signexponent_set_hidden(right);
 257        }
 258    Sgl_clear_exponent_set_hidden(left);
 259    diff_exponent = result_exponent - right_exponent;
 260
 261    /* 
 262     * Special case alignment of operands that would force alignment 
 263     * beyond the extent of the extension.  A further optimization
 264     * could special case this but only reduces the path length for this
 265     * infrequent case.
 266     */
 267    if(diff_exponent > SGL_THRESHOLD)
 268        {
 269        diff_exponent = SGL_THRESHOLD;
 270        }
 271    
 272    /* Align right operand by shifting to right */
 273    Sgl_right_align(/*operand*/right,/*shifted by*/diff_exponent,
 274     /*and lower to*/extent);
 275
 276    /* Treat sum and difference of the operands separately. */
 277    if( (/*signed*/int) save < 0 )
 278        {
 279        /*
 280         * Difference of the two operands.  Their can be no overflow.  A
 281         * borrow can occur out of the hidden bit and force a post
 282         * normalization phase.
 283         */
 284        Sgl_subtract_withextension(left,/*minus*/right,/*with*/extent,/*into*/result);
 285        if(Sgl_iszero_hidden(result))
 286            {
 287            /* Handle normalization */
 288            /* A straightforward algorithm would now shift the result
 289             * and extension left until the hidden bit becomes one.  Not
 290             * all of the extension bits need participate in the shift.
 291             * Only the two most significant bits (round and guard) are
 292             * needed.  If only a single shift is needed then the guard
 293             * bit becomes a significant low order bit and the extension
 294             * must participate in the rounding.  If more than a single 
 295             * shift is needed, then all bits to the right of the guard 
 296             * bit are zeros, and the guard bit may or may not be zero. */
 297            sign_save = Sgl_signextendedsign(result);
 298            Sgl_leftshiftby1_withextent(result,extent,result);
 299
 300            /* Need to check for a zero result.  The sign and exponent
 301             * fields have already been zeroed.  The more efficient test
 302             * of the full object can be used.
 303             */
 304            if(Sgl_iszero(result))
 305                /* Must have been "x-x" or "x+(-x)". */
 306                {
 307                if(Is_rounding_mode(ROUNDMINUS)) Sgl_setone_sign(result);
 308                *dstptr = result;
 309                return(NOEXCEPTION);
 310                }
 311            result_exponent--;
 312            /* Look to see if normalization is finished. */
 313            if(Sgl_isone_hidden(result))
 314                {
 315                if(result_exponent==0)
 316                    {
 317                    /* Denormalized, exponent should be zero.  Left operand *
 318                     * was normalized, so extent (guard, round) was zero    */
 319                    goto underflow;
 320                    }
 321                else
 322                    {
 323                    /* No further normalization is needed. */
 324                    Sgl_set_sign(result,/*using*/sign_save);
 325                    Ext_leftshiftby1(extent);
 326                    goto round;
 327                    }
 328                }
 329
 330            /* Check for denormalized, exponent should be zero.  Left    * 
 331             * operand was normalized, so extent (guard, round) was zero */
 332            if(!(underflowtrap = Is_underflowtrap_enabled()) &&
 333               result_exponent==0) goto underflow;
 334
 335            /* Shift extension to complete one bit of normalization and
 336             * update exponent. */
 337            Ext_leftshiftby1(extent);
 338
 339            /* Discover first one bit to determine shift amount.  Use a
 340             * modified binary search.  We have already shifted the result
 341             * one position right and still not found a one so the remainder
 342             * of the extension must be zero and simplifies rounding. */
 343            /* Scan bytes */
 344            while(Sgl_iszero_hiddenhigh7mantissa(result))
 345                {
 346                Sgl_leftshiftby8(result);
 347                if((result_exponent -= 8) <= 0  && !underflowtrap)
 348                    goto underflow;
 349                }
 350            /* Now narrow it down to the nibble */
 351            if(Sgl_iszero_hiddenhigh3mantissa(result))
 352                {
 353                /* The lower nibble contains the normalizing one */
 354                Sgl_leftshiftby4(result);
 355                if((result_exponent -= 4) <= 0 && !underflowtrap)
 356                    goto underflow;
 357                }
 358            /* Select case were first bit is set (already normalized)
 359             * otherwise select the proper shift. */
 360            if((jumpsize = Sgl_hiddenhigh3mantissa(result)) > 7)
 361                {
 362                /* Already normalized */
 363                if(result_exponent <= 0) goto underflow;
 364                Sgl_set_sign(result,/*using*/sign_save);
 365                Sgl_set_exponent(result,/*using*/result_exponent);
 366                *dstptr = result;
 367                return(NOEXCEPTION);
 368                }
 369            Sgl_sethigh4bits(result,/*using*/sign_save);
 370            switch(jumpsize) 
 371                {
 372                case 1:
 373                    {
 374                    Sgl_leftshiftby3(result);
 375                    result_exponent -= 3;
 376                    break;
 377                    }
 378                case 2:
 379                case 3:
 380                    {
 381                    Sgl_leftshiftby2(result);
 382                    result_exponent -= 2;
 383                    break;
 384                    }
 385                case 4:
 386                case 5:
 387                case 6:
 388                case 7:
 389                    {
 390                    Sgl_leftshiftby1(result);
 391                    result_exponent -= 1;
 392                    break;
 393                    }
 394                }
 395            if(result_exponent > 0) 
 396                {
 397                Sgl_set_exponent(result,/*using*/result_exponent);
 398                *dstptr = result;
 399                return(NOEXCEPTION); /* Sign bit is already set */
 400                }
 401            /* Fixup potential underflows */
 402          underflow:
 403            if(Is_underflowtrap_enabled())
 404                {
 405                Sgl_set_sign(result,sign_save);
 406                Sgl_setwrapped_exponent(result,result_exponent,unfl);
 407                *dstptr = result;
 408                /* inexact = FALSE; */
 409                return(UNDERFLOWEXCEPTION);
 410                }
 411            /* 
 412             * Since we cannot get an inexact denormalized result,
 413             * we can now return.
 414             */
 415            Sgl_right_align(result,/*by*/(1-result_exponent),extent);
 416            Sgl_clear_signexponent(result);
 417            Sgl_set_sign(result,sign_save);
 418            *dstptr = result;
 419            return(NOEXCEPTION);
 420            } /* end if(hidden...)... */
 421        /* Fall through and round */
 422        } /* end if(save < 0)... */
 423    else 
 424        {
 425        /* Add magnitudes */
 426        Sgl_addition(left,right,/*to*/result);
 427        if(Sgl_isone_hiddenoverflow(result))
 428            {
 429            /* Prenormalization required. */
 430            Sgl_rightshiftby1_withextent(result,extent,extent);
 431            Sgl_arithrightshiftby1(result);
 432            result_exponent++;
 433            } /* end if hiddenoverflow... */
 434        } /* end else ...add magnitudes... */
 435    
 436    /* Round the result.  If the extension is all zeros,then the result is
 437     * exact.  Otherwise round in the correct direction.  No underflow is
 438     * possible. If a postnormalization is necessary, then the mantissa is
 439     * all zeros so no shift is needed. */
 440  round:
 441    if(Ext_isnotzero(extent))
 442        {
 443        inexact = TRUE;
 444        switch(Rounding_mode())
 445            {
 446            case ROUNDNEAREST: /* The default. */
 447            if(Ext_isone_sign(extent))
 448                {
 449                /* at least 1/2 ulp */
 450                if(Ext_isnotzero_lower(extent)  ||
 451                  Sgl_isone_lowmantissa(result))
 452                    {
 453                    /* either exactly half way and odd or more than 1/2ulp */
 454                    Sgl_increment(result);
 455                    }
 456                }
 457            break;
 458
 459            case ROUNDPLUS:
 460            if(Sgl_iszero_sign(result))
 461                {
 462                /* Round up positive results */
 463                Sgl_increment(result);
 464                }
 465            break;
 466            
 467            case ROUNDMINUS:
 468            if(Sgl_isone_sign(result))
 469                {
 470                /* Round down negative results */
 471                Sgl_increment(result);
 472                }
 473            
 474            case ROUNDZERO:;
 475            /* truncate is simple */
 476            } /* end switch... */
 477        if(Sgl_isone_hiddenoverflow(result)) result_exponent++;
 478        }
 479    if(result_exponent == SGL_INFINITY_EXPONENT)
 480        {
 481        /* Overflow */
 482        if(Is_overflowtrap_enabled())
 483            {
 484            Sgl_setwrapped_exponent(result,result_exponent,ovfl);
 485            *dstptr = result;
 486            if (inexact)
 487                if (Is_inexacttrap_enabled())
 488                    return(OVERFLOWEXCEPTION | INEXACTEXCEPTION);
 489                else Set_inexactflag();
 490            return(OVERFLOWEXCEPTION);
 491            }
 492        else
 493            {
 494            Set_overflowflag();
 495            inexact = TRUE;
 496            Sgl_setoverflow(result);
 497            }
 498        }
 499    else Sgl_set_exponent(result,result_exponent);
 500    *dstptr = result;
 501    if(inexact) 
 502        if(Is_inexacttrap_enabled()) return(INEXACTEXCEPTION);
 503        else Set_inexactflag();
 504    return(NOEXCEPTION);
 505    }
 506