LXR linux/arch/parisc/math-emu/sfrem.c

   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/sfrem.c               $Revision: 1.1 $
  26 *
  27 *  Purpose:
  28 *      Single Precision Floating-point Remainder
  29 *
  30 *  External Interfaces:
  31 *      sgl_frem(srcptr1,srcptr2,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
  43#include "float.h"
  44#include "sgl_float.h"
  45
  46/*
  47 *  Single Precision Floating-point Remainder
  48 */
  49
  50int
  51sgl_frem (sgl_floating_point * srcptr1, sgl_floating_point * srcptr2,
  52          sgl_floating_point * dstptr, unsigned int *status)
  53{
  54        register unsigned int opnd1, opnd2, result;
  55        register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
  56        register boolean roundup = FALSE;
  57
  58        opnd1 = *srcptr1;
  59        opnd2 = *srcptr2;
  60        /*
  61         * check first operand for NaN's or infinity
  62         */
  63        if ((opnd1_exponent = Sgl_exponent(opnd1)) == SGL_INFINITY_EXPONENT) {
  64                if (Sgl_iszero_mantissa(opnd1)) {
  65                        if (Sgl_isnotnan(opnd2)) {
  66                                /* invalid since first operand is infinity */
  67                                if (Is_invalidtrap_enabled()) 
  68                                        return(INVALIDEXCEPTION);
  69                                Set_invalidflag();
  70                                Sgl_makequietnan(result);
  71                                *dstptr = result;
  72                                return(NOEXCEPTION);
  73                        }
  74                }
  75                else {
  76                        /*
  77                         * is NaN; signaling or quiet?
  78                         */
  79                        if (Sgl_isone_signaling(opnd1)) {
  80                                /* trap if INVALIDTRAP enabled */
  81                                if (Is_invalidtrap_enabled()) 
  82                                        return(INVALIDEXCEPTION);
  83                                /* make NaN quiet */
  84                                Set_invalidflag();
  85                                Sgl_set_quiet(opnd1);
  86                        }
  87                        /* 
  88                         * is second operand a signaling NaN? 
  89                         */
  90                        else if (Sgl_is_signalingnan(opnd2)) {
  91                                /* trap if INVALIDTRAP enabled */
  92                                if (Is_invalidtrap_enabled()) 
  93                                        return(INVALIDEXCEPTION);
  94                                /* make NaN quiet */
  95                                Set_invalidflag();
  96                                Sgl_set_quiet(opnd2);
  97                                *dstptr = opnd2;
  98                                return(NOEXCEPTION);
  99                        }
 100                        /*
 101                         * return quiet NaN
 102                         */
 103                        *dstptr = opnd1;
 104                        return(NOEXCEPTION);
 105                }
 106        } 
 107        /*
 108         * check second operand for NaN's or infinity
 109         */
 110        if ((opnd2_exponent = Sgl_exponent(opnd2)) == SGL_INFINITY_EXPONENT) {
 111                if (Sgl_iszero_mantissa(opnd2)) {
 112                        /*
 113                         * return first operand
 114                         */
 115                        *dstptr = opnd1;
 116                        return(NOEXCEPTION);
 117                }
 118                /*
 119                 * is NaN; signaling or quiet?
 120                 */
 121                if (Sgl_isone_signaling(opnd2)) {
 122                        /* trap if INVALIDTRAP enabled */
 123                        if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 124                        /* make NaN quiet */
 125                        Set_invalidflag();
 126                        Sgl_set_quiet(opnd2);
 127                }
 128                /*
 129                 * return quiet NaN
 130                 */
 131                *dstptr = opnd2;
 132                return(NOEXCEPTION);
 133        }
 134        /*
 135         * check second operand for zero
 136         */
 137        if (Sgl_iszero_exponentmantissa(opnd2)) {
 138                /* invalid since second operand is zero */
 139                if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 140                Set_invalidflag();
 141                Sgl_makequietnan(result);
 142                *dstptr = result;
 143                return(NOEXCEPTION);
 144        }
 145
 146        /* 
 147         * get sign of result
 148         */
 149        result = opnd1;  
 150
 151        /* 
 152         * check for denormalized operands
 153         */
 154        if (opnd1_exponent == 0) {
 155                /* check for zero */
 156                if (Sgl_iszero_mantissa(opnd1)) {
 157                        *dstptr = opnd1;
 158                        return(NOEXCEPTION);
 159                }
 160                /* normalize, then continue */
 161                opnd1_exponent = 1;
 162                Sgl_normalize(opnd1,opnd1_exponent);
 163        }
 164        else {
 165                Sgl_clear_signexponent_set_hidden(opnd1);
 166        }
 167        if (opnd2_exponent == 0) {
 168                /* normalize, then continue */
 169                opnd2_exponent = 1;
 170                Sgl_normalize(opnd2,opnd2_exponent);
 171        }
 172        else {
 173                Sgl_clear_signexponent_set_hidden(opnd2);
 174        }
 175
 176        /* find result exponent and divide step loop count */
 177        dest_exponent = opnd2_exponent - 1;
 178        stepcount = opnd1_exponent - opnd2_exponent;
 179
 180        /*
 181         * check for opnd1/opnd2 < 1
 182         */
 183        if (stepcount < 0) {
 184                /*
 185                 * check for opnd1/opnd2 > 1/2
 186                 *
 187                 * In this case n will round to 1, so 
 188                 *    r = opnd1 - opnd2 
 189                 */
 190                if (stepcount == -1 && Sgl_isgreaterthan(opnd1,opnd2)) {
 191                        Sgl_all(result) = ~Sgl_all(result);   /* set sign */
 192                        /* align opnd2 with opnd1 */
 193                        Sgl_leftshiftby1(opnd2); 
 194                        Sgl_subtract(opnd2,opnd1,opnd2);
 195                        /* now normalize */
 196                        while (Sgl_iszero_hidden(opnd2)) {
 197                                Sgl_leftshiftby1(opnd2);
 198                                dest_exponent--;
 199                        }
 200                        Sgl_set_exponentmantissa(result,opnd2);
 201                        goto testforunderflow;
 202                }
 203                /*
 204                 * opnd1/opnd2 <= 1/2
 205                 *
 206                 * In this case n will round to zero, so 
 207                 *    r = opnd1
 208                 */
 209                Sgl_set_exponentmantissa(result,opnd1);
 210                dest_exponent = opnd1_exponent;
 211                goto testforunderflow;
 212        }
 213
 214        /*
 215         * Generate result
 216         *
 217         * Do iterative subtract until remainder is less than operand 2.
 218         */
 219        while (stepcount-- > 0 && Sgl_all(opnd1)) {
 220                if (Sgl_isnotlessthan(opnd1,opnd2))
 221                        Sgl_subtract(opnd1,opnd2,opnd1);
 222                Sgl_leftshiftby1(opnd1);
 223        }
 224        /*
 225         * Do last subtract, then determine which way to round if remainder 
 226         * is exactly 1/2 of opnd2 
 227         */
 228        if (Sgl_isnotlessthan(opnd1,opnd2)) {
 229                Sgl_subtract(opnd1,opnd2,opnd1);
 230                roundup = TRUE;
 231        }
 232        if (stepcount > 0 || Sgl_iszero(opnd1)) {
 233                /* division is exact, remainder is zero */
 234                Sgl_setzero_exponentmantissa(result);
 235                *dstptr = result;
 236                return(NOEXCEPTION);
 237        }
 238
 239        /* 
 240         * Check for cases where opnd1/opnd2 < n 
 241         *
 242         * In this case the result's sign will be opposite that of
 243         * opnd1.  The mantissa also needs some correction.
 244         */
 245        Sgl_leftshiftby1(opnd1);
 246        if (Sgl_isgreaterthan(opnd1,opnd2)) {
 247                Sgl_invert_sign(result);
 248                Sgl_subtract((opnd2<<1),opnd1,opnd1);
 249        }
 250        /* check for remainder being exactly 1/2 of opnd2 */
 251        else if (Sgl_isequal(opnd1,opnd2) && roundup) { 
 252                Sgl_invert_sign(result);
 253        }
 254
 255        /* normalize result's mantissa */
 256        while (Sgl_iszero_hidden(opnd1)) {
 257                dest_exponent--;
 258                Sgl_leftshiftby1(opnd1);
 259        }
 260        Sgl_set_exponentmantissa(result,opnd1);
 261
 262        /* 
 263         * Test for underflow
 264         */
 265    testforunderflow:
 266        if (dest_exponent <= 0) {
 267                /* trap if UNDERFLOWTRAP enabled */
 268                if (Is_underflowtrap_enabled()) {
 269                        /*
 270                         * Adjust bias of result
 271                         */
 272                        Sgl_setwrapped_exponent(result,dest_exponent,unfl);
 273                        *dstptr = result;
 274                        /* frem is always exact */
 275                        return(UNDERFLOWEXCEPTION);
 276                }
 277                /*
 278                 * denormalize result or set to signed zero
 279                 */
 280                if (dest_exponent >= (1 - SGL_P)) {
 281                        Sgl_rightshift_exponentmantissa(result,1-dest_exponent);
 282                }
 283                else {
 284                        Sgl_setzero_exponentmantissa(result);
 285                }
 286        }
 287        else Sgl_set_exponent(result,dest_exponent);
 288        *dstptr = result;
 289        return(NOEXCEPTION);
 290}
 291