arch/parisc/math-emu/dfrem.c

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