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