Return to s_rintl.c CVS log

Up to [cvs.NetBSD.org] / src / lib / libm / src

Annotation of src/lib/libm/src/s_rintl.c, Revision 1.4.4.3

1.4.4.2   tls         1: /*     $NetBSD$        */
                      2:
                      3: /*-
                      4:  * Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
                      5:  * All rights reserved.
                      6:  *
                      7:  * Redistribution and use in source and binary forms, with or without
                      8:  * modification, are permitted provided that the following conditions
                      9:  * are met:
                     10:  * 1. Redistributions of source code must retain the above copyright
                     11:  *    notice, this list of conditions and the following disclaimer.
                     12:  * 2. Redistributions in binary form must reproduce the above copyright
                     13:  *    notice, this list of conditions and the following disclaimer in the
                     14:  *    documentation and/or other materials provided with the distribution.
                     15:  *
                     16:  * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
                     17:  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
                     18:  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
                     19:  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
                     20:  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
                     21:  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
                     22:  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
                     23:  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
                     24:  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
                     25:  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
                     26:  * SUCH DAMAGE.
                     27:  */
                     28:
                     29: #include <sys/cdefs.h>
                     30: #if 0
                     31: __FBSDID("$FreeBSD: src/lib/msun/src/s_rintl.c,v 1.5 2008/02/22 11:59:05 bde Exp $");
                     32: #else
                     33: __RCSID("$NetBSD$");
                     34: #endif
                     35:
                     36: #include <float.h>
                     37: #include <machine/ieee.h>
                     38:
                     39: #include "math.h"
                     40: #include "math_private.h"
                     41:
                     42: #ifdef __HAVE_LONG_DOUBLE
                     43: static const float
                     44: shift[2] = {
                     45: #if EXT_FRACBITS == 64
                     46:        0x1.0p63, -0x1.0p63
                     47: #elif EXT_FRACBITS == 113
                     48:        0x1.0p112, -0x1.0p112
1.4.4.3 ! tls        49: #elif EXT_FRACBITS == 112
        !            50:        0x1.0p111, -0x1.0p111
1.4.4.2   tls        51: #else
                     52: #error "Unsupported long double format"
                     53: #endif
                     54: };
                     55: static const float zero[2] = { 0.0, -0.0 };
                     56:
                     57: long double
                     58: rintl(long double x)
                     59: {
                     60:        union ieee_ext_u u;
                     61:        uint32_t expsign;
                     62:        int ex, sign;
                     63:
                     64:        u.extu_ld = x;
                     65:        u.extu_ext.ext_frach &= ~0x80000000;
                     66:        expsign = u.extu_ext.ext_sign;
                     67:        ex = expsign & 0x7fff;
                     68:
                     69:        if (ex >= EXT_EXP_BIAS + EXT_FRACBITS - 1) {
                     70:                if (ex == EXT_EXP_BIAS + EXT_FRACBITS)
                     71:                        return (x + x); /* Inf, NaN, or unsupported format */
                     72:                return (x);             /* finite and already an integer */
                     73:        }
                     74:        sign = expsign >> 15;
                     75:
                     76:        /*
                     77:         * The following code assumes that intermediate results are
                     78:         * evaluated in long double precision. If they are evaluated in
                     79:         * greater precision, double rounding may occur, and if they are
                     80:         * evaluated in less precision (as on i386), results will be
                     81:         * wildly incorrect.
                     82:         */
                     83:        x += shift[sign];
                     84:        x -= shift[sign];
                     85:
                     86:        /*
                     87:         * If the result is +-0, then it must have the same sign as x, but
                     88:         * the above calculation doesn't always give this.  Fix up the sign.
                     89:         */
                     90:        if (ex < EXT_EXP_BIAS && x == 0.0L)
                     91:                return (zero[sign]);
                     92:
                     93:        return (x);
                     94: }
                     95: #endif