[BACK]Return to ldexp_ieee754.c CVS log [TXT][DIR] Up to [cvs.NetBSD.org] / src / lib / libc / gen

Annotation of src/lib/libc/gen/ldexp_ieee754.c, Revision 1.2

1.2     ! kleink      1: /*     $NetBSD: ldexp_ieee754.c,v 1.1 2003/05/12 15:15:16 kleink Exp $ */
1.1       kleink      2:
                      3: /*-
                      4:  * Copyright (c) 1999 The NetBSD Foundation, Inc.
                      5:  * All rights reserved.
                      6:  *
                      7:  * This code is derived from software contributed to The NetBSD Foundation
                      8:  * by Charles M. Hannum.
                      9:  *
                     10:  * Redistribution and use in source and binary forms, with or without
                     11:  * modification, are permitted provided that the following conditions
                     12:  * are met:
                     13:  * 1. Redistributions of source code must retain the above copyright
                     14:  *    notice, this list of conditions and the following disclaimer.
                     15:  * 2. Redistributions in binary form must reproduce the above copyright
                     16:  *    notice, this list of conditions and the following disclaimer in the
                     17:  *    documentation and/or other materials provided with the distribution.
                     18:  * 3. All advertising materials mentioning features or use of this software
                     19:  *    must display the following acknowledgement:
                     20:  *        This product includes software developed by the NetBSD
                     21:  *        Foundation, Inc. and its contributors.
                     22:  * 4. Neither the name of The NetBSD Foundation nor the names of its
                     23:  *    contributors may be used to endorse or promote products derived
                     24:  *    from this software without specific prior written permission.
                     25:  *
                     26:  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
                     27:  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
                     28:  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
                     29:  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
                     30:  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
                     31:  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
                     32:  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
                     33:  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
                     34:  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
                     35:  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
                     36:  * POSSIBILITY OF SUCH DAMAGE.
                     37:  */
                     38:
                     39: #include <sys/cdefs.h>
                     40: #if defined(LIBC_SCCS) && !defined(lint)
1.2     ! kleink     41: __RCSID("$NetBSD: ldexp_ieee754.c,v 1.1 2003/05/12 15:15:16 kleink Exp $");
1.1       kleink     42: #endif /* LIBC_SCCS and not lint */
                     43:
                     44: #include <sys/types.h>
                     45: #include <machine/ieee.h>
                     46: #include <errno.h>
                     47: #include <math.h>
                     48:
                     49: /*
                     50:  * Multiply the given value by 2^expon.
                     51:  */
                     52: double
                     53: ldexp(double val, int expon)
                     54: {
                     55:        int oldexp, newexp;
1.2     ! kleink     56:        union ieee_double_u u, mul;
1.1       kleink     57:
1.2     ! kleink     58:        u.dlbu_d = val;
        !            59:        oldexp = u.dblu_dbl.dbl_exp;
1.1       kleink     60:
                     61:        /*
                     62:         * If input is zero, Inf or NaN, just return it.
                     63:         */
1.2     ! kleink     64:        if (u.dblu_d == 0.0 || oldexp == DBL_EXP_INFNAN)
1.1       kleink     65:                return (val);
                     66:
                     67:        if (oldexp == 0) {
                     68:                /*
                     69:                 * u.v is denormal.  We must adjust it so that the exponent
                     70:                 * arithmetic below will work.
                     71:                 */
                     72:                if (expon <= DBL_EXP_BIAS) {
                     73:                        /*
                     74:                         * Optimization: if the scaling can be done in a single
                     75:                         * multiply, or underflows, just do it now.
                     76:                         */
                     77:                        if (expon <= -DBL_FRACBITS) {
                     78:                                errno = ERANGE;
                     79:                                return (0.0);
                     80:                        }
1.2     ! kleink     81:                        mul.dblu_d = 0.0;
        !            82:                        mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
        !            83:                        u.dblu_d *= mul.dblu_d;
        !            84:                        if (u.dblu_d == 0.0) {
1.1       kleink     85:                                errno = ERANGE;
                     86:                                return (0.0);
                     87:                        }
1.2     ! kleink     88:                        return (u.dblu_d);
1.1       kleink     89:                } else {
                     90:                        /*
                     91:                         * We know that expon is very large, and therefore the
                     92:                         * result cannot be denormal (though it may be Inf).
                     93:                         * Shift u.v by just enough to make it normal.
                     94:                         */
1.2     ! kleink     95:                        mul.dblu_d = 0.0;
        !            96:                        mul.dblu_dbl.dbl_exp = DBL_FRACBITS + DBL_EXP_BIAS;
        !            97:                        u.dblu_d *= mul.dblu_d;
1.1       kleink     98:                        expon -= DBL_FRACBITS;
1.2     ! kleink     99:                        oldexp = u.dblu_dbl.dbl_exp;
1.1       kleink    100:                }
                    101:        }
                    102:
                    103:        /*
                    104:         * u.v is now normalized and oldexp has been adjusted if necessary.
                    105:         * Calculate the new exponent and check for underflow and overflow.
                    106:         */
                    107:        newexp = oldexp + expon;
                    108:
                    109:        if (newexp <= 0) {
                    110:                /*
                    111:                 * The output number is either denormal or underflows (see
                    112:                 * comments in machine/ieee.h).
                    113:                 */
                    114:                if (newexp <= -DBL_FRACBITS) {
                    115:                        errno = ERANGE;
                    116:                        return (0.0);
                    117:                }
                    118:                /*
                    119:                 * Denormalize the result.  We do this with a multiply.  If
                    120:                 * expon is very large, it won't fit in a double, so we have
                    121:                 * to adjust the exponent first.  This is safe because we know
                    122:                 * that u.v is normal at this point.
                    123:                 */
                    124:                if (expon <= -DBL_EXP_BIAS) {
1.2     ! kleink    125:                        u.dblu_dbl.dbl_exp = 1;
1.1       kleink    126:                        expon += oldexp - 1;
                    127:                }
1.2     ! kleink    128:                mul.dblu_d = 0.0;
        !           129:                mul.dblu_dbl.dbl_exp = expon + DBL_EXP_BIAS;
        !           130:                u.dblu_d *= mul.dblu_d;
        !           131:                return (u.dblu_d);
1.1       kleink    132:        } else if (newexp >= DBL_EXP_INFNAN) {
                    133:                /*
                    134:                 * The result overflowed; return +/-Inf.
                    135:                 */
1.2     ! kleink    136:                u.dblu_dbl.dbl_exp = DBL_EXP_INFNAN;
        !           137:                u.dblu_dbl.dbl_frach = 0;
        !           138:                u.dblu_dbl.dbl_fracl = 0;
1.1       kleink    139:                errno = ERANGE;
1.2     ! kleink    140:                return (u.dblu_d);
1.1       kleink    141:        } else {
                    142:                /*
                    143:                 * The result is normal; just replace the old exponent with the
                    144:                 * new one.
                    145:                 */
1.2     ! kleink    146:                u.dblu_dbl.dbl_exp = newexp;
        !           147:                return (u.dblu_d);
1.1       kleink    148:        }
                    149: }

CVSweb <webmaster@jp.NetBSD.org>