Home | History | Annotate | Line # | Download | only in ld80 
      1 /*-
      2  * SPDX-License-Identifier: BSD-3-Clause
      3  *
      4  * Copyright (c) 1985, 1993
      5  *	The Regents of the University of California.  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  * 3. Neither the name of the University nor the names of its contributors
     16  *    may be used to endorse or promote products derived from this software
     17  *    without specific prior written permission.
     18  *
     19  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
     20  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
     21  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
     22  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
     23  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     24  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
     25  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     26  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     27  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     28  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     29  * SUCH DAMAGE.
     30  */
     31 
     32 #include <sys/cdefs.h>
     33 /*
     34  * See bsdsrc/b_exp.c for implementation details.
     35  *
     36  * bsdrc/b_exp.c converted to long double by Steven G. Kargl.
     37  */
     38 
     39 #include "math_private.h"
     40 
     41 static const union ieee_ext_u
     42     p0u = LD80C(0xaaaaaaaaaaaaaaab,    -3,  1.66666666666666666671e-01L),
     43     p1u = LD80C(0xb60b60b60b60b59a,    -9, -2.77777777777777775377e-03L),
     44     p2u = LD80C(0x8ab355e008a3cfce,   -14,  6.61375661375629297465e-05L),
     45     p3u = LD80C(0xddebbc994b0c1376,   -20, -1.65343915327882529784e-06L),
     46     p4u = LD80C(0xb354784cb4ef4c41,   -25,  4.17535101591534118469e-08L),
     47     p5u = LD80C(0x913e8a718382ce75,   -30, -1.05679137034774806475e-09L),
     48     p6u = LD80C(0xe8f0042aa134502e,   -36,  2.64819349895429516863e-11L);
     49 #define	p1	(p0u.extu_ld)
     50 #define	p2	(p1u.extu_ld)
     51 #define	p3	(p2u.extu_ld)
     52 #define	p4	(p3u.extu_ld)
     53 #define	p5	(p4u.extu_ld)
     54 #define	p6	(p5u.extu_ld)
     55 #define	p7	(p6u.extu_ld)
     56 
     57 /*
     58  * lnhuge = (LDBL_MAX_EXP + 9) * log(2.)
     59  * lntiny = (LDBL_MIN_EXP - 64 - 10) * log(2.)
     60  * invln2 = 1 / log(2.)
     61  */
     62 static const union ieee_ext_u
     63 ln2hiu  = LD80C(0xb17217f700000000,  -1,  6.93147180369123816490e-01L),
     64 ln2lou  = LD80C(0xd1cf79abc9e3b398, -33,  1.90821492927058781614e-10L),
     65 lnhugeu = LD80C(0xb18b0c0330a8fad9,  13,  1.13627617309191834574e+04L),
     66 lntinyu = LD80C(0xb236f28a68bc3bd7,  13, -1.14057368561139000667e+04L),
     67 invln2u = LD80C(0xb8aa3b295c17f0bc,   0,  1.44269504088896340739e+00L);
     68 #define	ln2hi	(ln2hiu.extu_ld)
     69 #define ln2lo	(ln2lou.extu_ld)
     70 #define lnhuge	(lnhugeu.extu_ld)
     71 #define	lntiny	(lntinyu.extu_ld)
     72 #define	invln2	(invln2u.extu_ld)
     73 
     74 /* returns exp(r = x + c) for |c| < |x| with no overlap.  */
     75 
     76 static long double
     77 __exp__LD(long double x, long double c)
     78 {
     79 	long double hi, lo, z;
     80 	int k;
     81 
     82 	if (x != x)	/* x is NaN. */
     83 		return(x);
     84 
     85 	if (x <= lnhuge) {
     86 		if (x >= lntiny) {
     87 			/* argument reduction: x --> x - k*ln2 */
     88 			z = invln2 * x;
     89 			k = z + copysignl(0.5L, x);
     90 
     91 		    	/*
     92 			 * Express (x + c) - k * ln2 as hi - lo.
     93 			 * Let x = hi - lo rounded.
     94 			 */
     95 			hi = x - k * ln2hi;	/* Exact. */
     96 			lo = k * ln2lo - c;
     97 			x = hi - lo;
     98 
     99 			/* Return 2^k*[1+x+x*c/(2+c)]  */
    100 			z = x * x;
    101 			c = x - z * (p1 + z * (p2 + z * (p3 + z * (p4 +
    102 			    z * (p5 + z * (p6 + z * p7))))));
    103 			c = (x * c) / (2 - c);
    104 
    105 			return (ldexpl(1 + (hi - (lo - c)), k));
    106 		} else {
    107 			/* exp(-INF) is 0. exp(-big) underflows to 0.  */
    108 			return (isfinite(x) ? ldexpl(1., -5000) : 0);
    109 		}
    110 	} else
    111 		/* exp(INF) is INF, exp(+big#) overflows to INF */
    112 		return (isfinite(x) ? ldexpl(1., 5000) : x);
    113 }
    114