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