sys/arch/hppa/spmath/dfrem.c

   1 /*      $NetBSD: dfrem.c,v 1.3 2005/12/11 12:17:40 christos Exp $       */
   2
   3 /*      $OpenBSD: dfrem.c,v 1.4 2001/03/29 03:58:17 mickey Exp $        */
   4
   5 /*
   6  * Copyright 1996 1995 by Open Software Foundation, Inc.
   7  *              All Rights Reserved
   8  *
   9  * Permission to use, copy, modify, and distribute this software and
  10  * its documentation for any purpose and without fee is hereby granted,
  11  * provided that the above copyright notice appears in all copies and
  12  * that both the copyright notice and this permission notice appear in
  13  * supporting documentation.
  14  *
  15  * OSF DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE
  16  * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  17  * FOR A PARTICULAR PURPOSE.
  18  *
  19  * IN NO EVENT SHALL OSF BE LIABLE FOR ANY SPECIAL, INDIRECT, OR
  20  * CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
  21  * LOSS OF USE, DATA OR PROFITS, WHETHER IN ACTION OF CONTRACT,
  22  * NEGLIGENCE, OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION
  23  * WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
  24  *
  25  */
  26 /*
  27  * pmk1.1
  28  */
  29 /*
  30  * (c) Copyright 1986 HEWLETT-PACKARD COMPANY
  31  *
  32  * To anyone who acknowledges that this file is provided "AS IS"
  33  * without any express or implied warranty:
  34  *     permission to use, copy, modify, and distribute this file
  35  * for any purpose is hereby granted without fee, provided that
  36  * the above copyright notice and this notice appears in all
  37  * copies, and that the name of Hewlett-Packard Company not be
  38  * used in advertising or publicity pertaining to distribution
  39  * of the software without specific, written prior permission.
  40  * Hewlett-Packard Company makes no representations about the
  41  * suitability of this software for any purpose.
  42  */
  43
  44 #include <sys/cdefs.h>
  45 __KERNEL_RCSID(0, "$NetBSD: dfrem.c,v 1.3 2005/12/11 12:17:40 christos Exp $");
  46
  47 #include "../spmath/float.h"
  48 #include "../spmath/dbl_float.h"
  49
  50 /*
  51  *  Double Precision Floating-point Remainder
  52  */
  53 int
  54 dbl_frem(srcptr1,srcptr2,dstptr,status)
  55
  56 dbl_floating_point *srcptr1, *srcptr2, *dstptr;
  57 unsigned int *status;
  58 {
  59         register unsigned int opnd1p1, opnd1p2, opnd2p1, opnd2p2;
  60         register unsigned int resultp1, resultp2;
  61         register int opnd1_exponent, opnd2_exponent, dest_exponent, stepcount;
  62         register int roundup = false;
  63
  64         Dbl_copyfromptr(srcptr1,opnd1p1,opnd1p2);
  65         Dbl_copyfromptr(srcptr2,opnd2p1,opnd2p2);
  66         /*
  67          * check first operand for NaN's or infinity
  68          */
  69         if ((opnd1_exponent = Dbl_exponent(opnd1p1)) == DBL_INFINITY_EXPONENT) {
  70                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
  71                         if (Dbl_isnotnan(opnd2p1,opnd2p2)) {
  72                                 /* invalid since first operand is infinity */
  73                                 if (Is_invalidtrap_enabled())
  74                                         return(INVALIDEXCEPTION);
  75                                 Set_invalidflag();
  76                                 Dbl_makequietnan(resultp1,resultp2);
  77                                 Dbl_copytoptr(resultp1,resultp2,dstptr);
  78                                 return(NOEXCEPTION);
  79                         }
  80                 }
  81                 else {
  82                         /*
  83                          * is NaN; signaling or quiet?
  84                          */
  85                         if (Dbl_isone_signaling(opnd1p1)) {
  86                                 /* trap if INVALIDTRAP enabled */
  87                                 if (Is_invalidtrap_enabled())
  88                                         return(INVALIDEXCEPTION);
  89                                 /* make NaN quiet */
  90                                 Set_invalidflag();
  91                                 Dbl_set_quiet(opnd1p1);
  92                         }
  93                         /*
  94                          * is second operand a signaling NaN?
  95                          */
  96                         else if (Dbl_is_signalingnan(opnd2p1)) {
  97                                 /* trap if INVALIDTRAP enabled */
  98                                 if (Is_invalidtrap_enabled())
  99                                         return(INVALIDEXCEPTION);
 100                                 /* make NaN quiet */
 101                                 Set_invalidflag();
 102                                 Dbl_set_quiet(opnd2p1);
 103                                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
 104                                 return(NOEXCEPTION);
 105                         }
 106                         /*
 107                          * return quiet NaN
 108                          */
 109                         Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
 110                         return(NOEXCEPTION);
 111                 }
 112         }
 113         /*
 114          * check second operand for NaN's or infinity
 115          */
 116         if ((opnd2_exponent = Dbl_exponent(opnd2p1)) == DBL_INFINITY_EXPONENT) {
 117                 if (Dbl_iszero_mantissa(opnd2p1,opnd2p2)) {
 118                         /*
 119                          * return first operand
 120                          */
 121                         Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
 122                         return(NOEXCEPTION);
 123                 }
 124                 /*
 125                  * is NaN; signaling or quiet?
 126                  */
 127                 if (Dbl_isone_signaling(opnd2p1)) {
 128                         /* trap if INVALIDTRAP enabled */
 129                         if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 130                         /* make NaN quiet */
 131                         Set_invalidflag();
 132                         Dbl_set_quiet(opnd2p1);
 133                 }
 134                 /*
 135                  * return quiet NaN
 136                  */
 137                 Dbl_copytoptr(opnd2p1,opnd2p2,dstptr);
 138                 return(NOEXCEPTION);
 139         }
 140         /*
 141          * check second operand for zero
 142          */
 143         if (Dbl_iszero_exponentmantissa(opnd2p1,opnd2p2)) {
 144                 /* invalid since second operand is zero */
 145                 if (Is_invalidtrap_enabled()) return(INVALIDEXCEPTION);
 146                 Set_invalidflag();
 147                 Dbl_makequietnan(resultp1,resultp2);
 148                 Dbl_copytoptr(resultp1,resultp2,dstptr);
 149                 return(NOEXCEPTION);
 150         }
 151
 152         /*
 153          * get sign of result
 154          */
 155         resultp1 = opnd1p1;
 156
 157         /*
 158          * check for denormalized operands
 159          */
 160         if (opnd1_exponent == 0) {
 161                 /* check for zero */
 162                 if (Dbl_iszero_mantissa(opnd1p1,opnd1p2)) {
 163                         Dbl_copytoptr(opnd1p1,opnd1p2,dstptr);
 164                         return(NOEXCEPTION);
 165                 }
 166                 /* normalize, then continue */
 167                 opnd1_exponent = 1;
 168                 Dbl_normalize(opnd1p1,opnd1p2,opnd1_exponent);
 169         }
 170         else {
 171                 Dbl_clear_signexponent_set_hidden(opnd1p1);
 172         }
 173         if (opnd2_exponent == 0) {
 174                 /* normalize, then continue */
 175                 opnd2_exponent = 1;
 176                 Dbl_normalize(opnd2p1,opnd2p2,opnd2_exponent);
 177         }
 178         else {
 179                 Dbl_clear_signexponent_set_hidden(opnd2p1);
 180         }
 181
 182         /* find result exponent and divide step loop count */
 183         dest_exponent = opnd2_exponent - 1;
 184         stepcount = opnd1_exponent - opnd2_exponent;
 185
 186         /*
 187          * check for opnd1/opnd2 < 1
 188          */
 189         if (stepcount < 0) {
 190                 /*
 191                  * check for opnd1/opnd2 > 1/2
 192                  *
 193                  * In this case n will round to 1, so
 194                  *    r = opnd1 - opnd2
 195                  */
 196                 if (stepcount == -1 &&
 197                     Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
 198                         /* set sign */
 199                         Dbl_allp1(resultp1) = ~Dbl_allp1(resultp1);
 200                         /* align opnd2 with opnd1 */
 201                         Dbl_leftshiftby1(opnd2p1,opnd2p2);
 202                         Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,
 203                          opnd2p1,opnd2p2);
 204                         /* now normalize */
 205                         while (Dbl_iszero_hidden(opnd2p1)) {
 206                                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
 207                                 dest_exponent--;
 208                         }
 209                         Dbl_set_exponentmantissa(resultp1,resultp2,opnd2p1,opnd2p2);
 210                         goto testforunderflow;
 211                 }
 212                 /*
 213                  * opnd1/opnd2 <= 1/2
 214                  *
 215                  * In this case n will round to zero, so
 216                  *    r = opnd1
 217                  */
 218                 Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
 219                 dest_exponent = opnd1_exponent;
 220                 goto testforunderflow;
 221         }
 222
 223         /*
 224          * Generate result
 225          *
 226          * Do iterative subtract until remainder is less than operand 2.
 227          */
 228         while (stepcount-- > 0 && (Dbl_allp1(opnd1p1) || Dbl_allp2(opnd1p2))) {
 229                 if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
 230                         Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
 231                 }
 232                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
 233         }
 234         /*
 235          * Do last subtract, then determine which way to round if remainder
 236          * is exactly 1/2 of opnd2
 237          */
 238         if (Dbl_isnotlessthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
 239                 Dbl_subtract(opnd1p1,opnd1p2,opnd2p1,opnd2p2,opnd1p1,opnd1p2);
 240                 roundup = true;
 241         }
 242         if (stepcount > 0 || Dbl_iszero(opnd1p1,opnd1p2)) {
 243                 /* division is exact, remainder is zero */
 244                 Dbl_setzero_exponentmantissa(resultp1,resultp2);
 245                 Dbl_copytoptr(resultp1,resultp2,dstptr);
 246                 return(NOEXCEPTION);
 247         }
 248
 249         /*
 250          * Check for cases where opnd1/opnd2 < n
 251          *
 252          * In this case the result's sign will be opposite that of
 253          * opnd1.  The mantissa also needs some correction.
 254          */
 255         Dbl_leftshiftby1(opnd1p1,opnd1p2);
 256         if (Dbl_isgreaterthan(opnd1p1,opnd1p2,opnd2p1,opnd2p2)) {
 257                 Dbl_invert_sign(resultp1);
 258                 Dbl_leftshiftby1(opnd2p1,opnd2p2);
 259                 Dbl_subtract(opnd2p1,opnd2p2,opnd1p1,opnd1p2,opnd1p1,opnd1p2);
 260         }
 261         /* check for remainder being exactly 1/2 of opnd2 */
 262         else if (Dbl_isequal(opnd1p1,opnd1p2,opnd2p1,opnd2p2) && roundup) {
 263                 Dbl_invert_sign(resultp1);
 264         }
 265
 266         /* normalize result's mantissa */
 267         while (Dbl_iszero_hidden(opnd1p1)) {
 268                 dest_exponent--;
 269                 Dbl_leftshiftby1(opnd1p1,opnd1p2);
 270         }
 271         Dbl_set_exponentmantissa(resultp1,resultp2,opnd1p1,opnd1p2);
 272
 273         /*
 274          * Test for underflow
 275          */
 276     testforunderflow:
 277         if (dest_exponent <= 0) {
 278                 /* trap if UNDERFLOWTRAP enabled */
 279                 if (Is_underflowtrap_enabled()) {
 280                         /*
 281                          * Adjust bias of result
 282                          */
 283                         Dbl_setwrapped_exponent(resultp1,dest_exponent,unfl);
 284                         /* frem is always exact */
 285                         Dbl_copytoptr(resultp1,resultp2,dstptr);
 286                         return(UNDERFLOWEXCEPTION);
 287                 }
 288                 /*
 289                  * denormalize result or set to signed zero
 290                  */
 291                 if (dest_exponent >= (1 - DBL_P)) {
 292                         Dbl_rightshift_exponentmantissa(resultp1,resultp2,
 293                          1-dest_exponent);
 294                 }
 295                 else {
 296                         Dbl_setzero_exponentmantissa(resultp1,resultp2);
 297                 }
 298         }
 299         else Dbl_set_exponent(resultp1,dest_exponent);
 300         Dbl_copytoptr(resultp1,resultp2,dstptr);
 301         return(NOEXCEPTION);
 302 }