LXR linux/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
  37int
  38dbl_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}
 285