Minix/minix
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Synchronize on NetBSD-CVS (2013/12/1 12:00:00 UTC)

Browse Source 
- Fix for possible unset uid/gid in toproto
 - Fix for default mtree style
 - Update libelf
 - Importing libexecinfo
 - Resynchronize GCC, mpc, gmp, mpfr
 - build.sh: Replace params with show-params.
     This has been done as the make target has been renamed in the same
     way, while a new target named params has been added. This new
     target generates a file containing all the parameters, instead of
     printing it on the console.
 - Update test48 with new etc/services (Fix by Ben Gras <ben@minix3.org)
     get getservbyport() out of the inner loop

Change-Id: Ie6ad5226fa2621ff9f0dee8782ea48f9443d2091
This commit is contained in:
Lionel Sambuc
2013-12-06 12:04:52 +01:00
parent ff10274392
commit 84d9c625bf
4655 changed files with 379317 additions and 151059 deletions
Download Patch File Download Diff File Expand all files Collapse all files

137
lib/libm/src/b_exp.c Normal file
View File

@@ -0,0 +1,137 @@
/* $NetBSD: b_exp.c,v 1.1 2012/05/05 17:54:14 christos Exp $ */
/*
* Copyright (c) 1985, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/* @(#)exp.c 8.1 (Berkeley) 6/4/93 */
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: release/9.0.0/lib/msun/bsdsrc/b_exp.c 176449 2008-02-22 02:26:51Z das $");
#else
__RCSID("$NetBSD: b_exp.c,v 1.1 2012/05/05 17:54:14 christos Exp $");
#endif
/* EXP(X)
* RETURN THE EXPONENTIAL OF X
* DOUBLE PRECISION (IEEE 53 bits, VAX D FORMAT 56 BITS)
* CODED IN C BY K.C. NG, 1/19/85;
* REVISED BY K.C. NG on 2/6/85, 2/15/85, 3/7/85, 3/24/85, 4/16/85, 6/14/86.
*
* Required system supported functions:
* scalb(x,n)
* copysign(x,y)
* finite(x)
*
* Method:
* 1. Argument Reduction: given the input x, find r and integer k such
* that
* x = k*ln2 + r, |r| <= 0.5*ln2 .
* r will be represented as r := z+c for better accuracy.
*
* 2. Compute exp(r) by
*
* exp(r) = 1 + r + r*R1/(2-R1),
* where
* R1 = x - x^2*(p1+x^2*(p2+x^2*(p3+x^2*(p4+p5*x^2)))).
*
* 3. exp(x) = 2^k * exp(r) .
*
* Special cases:
* exp(INF) is INF, exp(NaN) is NaN;
* exp(-INF)= 0;
* for finite argument, only exp(0)=1 is exact.
*
* Accuracy:
* exp(x) returns the exponential of x nearly rounded. In a test run
* with 1,156,000 random arguments on a VAX, the maximum observed
* error was 0.869 ulps (units in the last place).
*/
#include "math.h"
#include "math_private.h"
static const double p1 = 0x1.555555555553ep-3;
static const double p2 = -0x1.6c16c16bebd93p-9;
static const double p3 = 0x1.1566aaf25de2cp-14;
static const double p4 = -0x1.bbd41c5d26bf1p-20;
static const double p5 = 0x1.6376972bea4d0p-25;
static const double ln2hi = 0x1.62e42fee00000p-1;
static const double ln2lo = 0x1.a39ef35793c76p-33;
static const double lnhuge = 0x1.6602b15b7ecf2p9;
static const double lntiny = -0x1.77af8ebeae354p9;
static const double invln2 = 0x1.71547652b82fep0;
/* returns exp(r = x + c) for |c| < |x| with no overlap. */
double
__exp__D(double x, double c)
{
double z,hi,lo;
int k;
if (x != x) /* x is NaN */
return(x);
if ( x <= lnhuge ) {
if ( x >= lntiny ) {
/* argument reduction : x --> x - k*ln2 */
z = invln2*x;
k = z + copysign(.5, x);
/* express (x+c)-k*ln2 as hi-lo and let x=hi-lo rounded */
hi=(x-k*ln2hi); /* Exact. */
x= hi - (lo = k*ln2lo-c);
/* return 2^k*[1+x+x*c/(2+c)] */
z=x*x;
c= x - z*(p1+z*(p2+z*(p3+z*(p4+z*p5))));
c = (x*c)/(2.0-c);
return scalb(1.+(hi-(lo - c)), k);
}
/* end of x > lntiny */
else
/* exp(-big#) underflows to zero */
if(finite(x)) return(scalb(1.0,-5000));
/* exp(-INF) is zero */
else return(0.0);
}
/* end of x < lnhuge */
else
/* exp(INF) is INF, exp(+big#) overflows to INF */
return( finite(x) ? scalb(1.0,5000) : x);
}

406
lib/libm/src/b_log.c Normal file
View File

@@ -0,0 +1,406 @@
/* $NetBSD: b_log.c,v 1.1 2012/05/05 17:54:14 christos Exp $ */
/*
* Copyright (c) 1992, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/* @(#)log.c 8.2 (Berkeley) 11/30/93 */
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: release/9.0.0/lib/msun/bsdsrc/b_log.c 176449 2008-02-22 02:26:51Z das $");
#else
__RCSID("$NetBSD: b_log.c,v 1.1 2012/05/05 17:54:14 christos Exp $");
#endif
#include <errno.h>
#include "math.h"
#include "math_private.h"
/* Table-driven natural logarithm.
*
* This code was derived, with minor modifications, from:
* Peter Tang, "Table-Driven Implementation of the
* Logarithm in IEEE Floating-Point arithmetic." ACM Trans.
* Math Software, vol 16. no 4, pp 378-400, Dec 1990).
*
* Calculates log(2^m*F*(1+f/F)), |f/j| <= 1/256,
* where F = j/128 for j an integer in [0, 128].
*
* log(2^m) = log2_hi*m + log2_tail*m
* since m is an integer, the dominant term is exact.
* m has at most 10 digits (for subnormal numbers),
* and log2_hi has 11 trailing zero bits.
*
* log(F) = logF_hi[j] + logF_lo[j] is in tabular form in log_table.h
* logF_hi[] + 512 is exact.
*
* log(1+f/F) = 2*f/(2*F + f) + 1/12 * (2*f/(2*F + f))**3 + ...
* the leading term is calculated to extra precision in two
* parts, the larger of which adds exactly to the dominant
* m and F terms.
* There are two cases:
* 1. when m, j are non-zero (m | j), use absolute
* precision for the leading term.
* 2. when m = j = 0, |1-x| < 1/256, and log(x) ~= (x-1).
* In this case, use a relative precision of 24 bits.
* (This is done differently in the original paper)
*
* Special cases:
* 0 return signalling -Inf
* neg return signalling NaN
* +Inf return +Inf
*/
#define N 128
/* Table of log(Fj) = logF_head[j] + logF_tail[j], for Fj = 1+j/128.
* Used for generation of extend precision logarithms.
* The constant 35184372088832 is 2^45, so the divide is exact.
* It ensures correct reading of logF_head, even for inaccurate
* decimal-to-binary conversion routines. (Everybody gets the
* right answer for integers less than 2^53.)
* Values for log(F) were generated using error < 10^-57 absolute
* with the bc -l package.
*/
static const double A1 = .08333333333333178827;
static const double A2 = .01250000000377174923;
static const double A3 = .002232139987919447809;
static const double A4 = .0004348877777076145742;
static const double logF_head[N+1] = {
0.,
.007782140442060381246,
.015504186535963526694,
.023167059281547608406,
.030771658666765233647,
.038318864302141264488,
.045809536031242714670,
.053244514518837604555,
.060624621816486978786,
.067950661908525944454,
.075223421237524235039,
.082443669210988446138,
.089612158689760690322,
.096729626458454731618,
.103796793681567578460,
.110814366340264314203,
.117783035656430001836,
.124703478501032805070,
.131576357788617315236,
.138402322859292326029,
.145182009844575077295,
.151916042025732167530,
.158605030176659056451,
.165249572895390883786,
.171850256926518341060,
.178407657472689606947,
.184922338493834104156,
.191394852999565046047,
.197825743329758552135,
.204215541428766300668,
.210564769107350002741,
.216873938300523150246,
.223143551314024080056,
.229374101064877322642,
.235566071312860003672,
.241719936886966024758,
.247836163904594286577,
.253915209980732470285,
.259957524436686071567,
.265963548496984003577,
.271933715484010463114,
.277868451003087102435,
.283768173130738432519,
.289633292582948342896,
.295464212893421063199,
.301261330578199704177,
.307025035294827830512,
.312755710004239517729,
.318453731118097493890,
.324119468654316733591,
.329753286372579168528,
.335355541920762334484,
.340926586970454081892,
.346466767346100823488,
.351976423156884266063,
.357455888922231679316,
.362905493689140712376,
.368325561158599157352,
.373716409793814818840,
.379078352934811846353,
.384411698910298582632,
.389716751140440464951,
.394993808240542421117,
.400243164127459749579,
.405465108107819105498,
.410659924985338875558,
.415827895143593195825,
.420969294644237379543,
.426084395310681429691,
.431173464818130014464,
.436236766774527495726,
.441274560805140936281,
.446287102628048160113,
.451274644139630254358,
.456237433481874177232,
.461175715122408291790,
.466089729924533457960,
.470979715219073113985,
.475845904869856894947,
.480688529345570714212,
.485507815781602403149,
.490303988045525329653,
.495077266798034543171,
.499827869556611403822,
.504556010751912253908,
.509261901790523552335,
.513945751101346104405,
.518607764208354637958,
.523248143765158602036,
.527867089620485785417,
.532464798869114019908,
.537041465897345915436,
.541597282432121573947,
.546132437597407260909,
.550647117952394182793,
.555141507540611200965,
.559615787935399566777,
.564070138285387656651,
.568504735352689749561,
.572919753562018740922,
.577315365035246941260,
.581691739635061821900,
.586049045003164792433,
.590387446602107957005,
.594707107746216934174,
.599008189645246602594,
.603290851438941899687,
.607555250224322662688,
.611801541106615331955,
.616029877215623855590,
.620240409751204424537,
.624433288012369303032,
.628608659422752680256,
.632766669570628437213,
.636907462236194987781,
.641031179420679109171,
.645137961373620782978,
.649227946625615004450,
.653301272011958644725,
.657358072709030238911,
.661398482245203922502,
.665422632544505177065,
.669430653942981734871,
.673422675212350441142,
.677398823590920073911,
.681359224807238206267,
.685304003098281100392,
.689233281238557538017,
.693147180560117703862
};
static const double logF_tail[N+1] = {
0.,
-.00000000000000543229938420049,
.00000000000000172745674997061,
-.00000000000001323017818229233,
-.00000000000001154527628289872,
-.00000000000000466529469958300,
.00000000000005148849572685810,
-.00000000000002532168943117445,
-.00000000000005213620639136504,
-.00000000000001819506003016881,
.00000000000006329065958724544,
.00000000000008614512936087814,
-.00000000000007355770219435028,
.00000000000009638067658552277,
.00000000000007598636597194141,
.00000000000002579999128306990,
-.00000000000004654729747598444,
-.00000000000007556920687451336,
.00000000000010195735223708472,
-.00000000000017319034406422306,
-.00000000000007718001336828098,
.00000000000010980754099855238,
-.00000000000002047235780046195,
-.00000000000008372091099235912,
.00000000000014088127937111135,
.00000000000012869017157588257,
.00000000000017788850778198106,
.00000000000006440856150696891,
.00000000000016132822667240822,
-.00000000000007540916511956188,
-.00000000000000036507188831790,
.00000000000009120937249914984,
.00000000000018567570959796010,
-.00000000000003149265065191483,
-.00000000000009309459495196889,
.00000000000017914338601329117,
-.00000000000001302979717330866,
.00000000000023097385217586939,
.00000000000023999540484211737,
.00000000000015393776174455408,
-.00000000000036870428315837678,
.00000000000036920375082080089,
-.00000000000009383417223663699,
.00000000000009433398189512690,
.00000000000041481318704258568,
-.00000000000003792316480209314,
.00000000000008403156304792424,
-.00000000000034262934348285429,
.00000000000043712191957429145,
-.00000000000010475750058776541,
-.00000000000011118671389559323,
.00000000000037549577257259853,
.00000000000013912841212197565,
.00000000000010775743037572640,
.00000000000029391859187648000,
-.00000000000042790509060060774,
.00000000000022774076114039555,
.00000000000010849569622967912,
-.00000000000023073801945705758,
.00000000000015761203773969435,
.00000000000003345710269544082,
-.00000000000041525158063436123,
.00000000000032655698896907146,
-.00000000000044704265010452446,
.00000000000034527647952039772,
-.00000000000007048962392109746,
.00000000000011776978751369214,
-.00000000000010774341461609578,
.00000000000021863343293215910,
.00000000000024132639491333131,
.00000000000039057462209830700,
-.00000000000026570679203560751,
.00000000000037135141919592021,
-.00000000000017166921336082431,
-.00000000000028658285157914353,
-.00000000000023812542263446809,
.00000000000006576659768580062,
-.00000000000028210143846181267,
.00000000000010701931762114254,
.00000000000018119346366441110,
.00000000000009840465278232627,
-.00000000000033149150282752542,
-.00000000000018302857356041668,
-.00000000000016207400156744949,
.00000000000048303314949553201,
-.00000000000071560553172382115,
.00000000000088821239518571855,
-.00000000000030900580513238244,
-.00000000000061076551972851496,
.00000000000035659969663347830,
.00000000000035782396591276383,
-.00000000000046226087001544578,
.00000000000062279762917225156,
.00000000000072838947272065741,
.00000000000026809646615211673,
-.00000000000010960825046059278,
.00000000000002311949383800537,
-.00000000000058469058005299247,
-.00000000000002103748251144494,
-.00000000000023323182945587408,
-.00000000000042333694288141916,
-.00000000000043933937969737844,
.00000000000041341647073835565,
.00000000000006841763641591466,
.00000000000047585534004430641,
.00000000000083679678674757695,
-.00000000000085763734646658640,
.00000000000021913281229340092,
-.00000000000062242842536431148,
-.00000000000010983594325438430,
.00000000000065310431377633651,
-.00000000000047580199021710769,
-.00000000000037854251265457040,
.00000000000040939233218678664,
.00000000000087424383914858291,
.00000000000025218188456842882,
-.00000000000003608131360422557,
-.00000000000050518555924280902,
.00000000000078699403323355317,
-.00000000000067020876961949060,
.00000000000016108575753932458,
.00000000000058527188436251509,
-.00000000000035246757297904791,
-.00000000000018372084495629058,
.00000000000088606689813494916,
.00000000000066486268071468700,
.00000000000063831615170646519,
.00000000000025144230728376072,
-.00000000000017239444525614834
};
/*
* Extra precision variant, returning struct {double a, b;};
* log(x) = a+b to 63 bits, with a rounded to 26 bits.
*/
struct Double
__log__D(double x)
{
int m, j;
double F, f, g, q, u, v, u2;
volatile double u1;
struct Double r;
/* Argument reduction: 1 <= g < 2; x/2^m = g; */
/* y = F*(1 + f/F) for |f| <= 2^-8 */
m = logb(x);
g = ldexp(x, -m);
if (m == -1022) {
j = logb(g), m += j;
g = ldexp(g, -j);
}
j = N*(g-1) + .5;
F = (1.0/N) * j + 1;
f = g - F;
g = 1/(2*F+f);
u = 2*f*g;
v = u*u;
q = u*v*(A1 + v*(A2 + v*(A3 + v*A4)));
if (m | j)
u1 = u + 513, u1 -= 513;
else
u1 = u, TRUNC(u1);
u2 = (2.0*(f - F*u1) - u1*f) * g;
u1 += m*logF_head[N] + logF_head[j];
u2 += logF_tail[j]; u2 += q;
u2 += logF_tail[N]*m;
r.a = u1 + u2; /* Only difference is here */
TRUNC(r.a);
r.b = (u1 - r.a) + u2;
return (r);
}

317
lib/libm/src/b_tgamma.c Normal file
View File

@@ -0,0 +1,317 @@
/* $NetBSD: b_tgamma.c,v 1.1 2012/05/05 17:54:14 christos Exp $ */
/*-
* Copyright (c) 1992, 1993
* The Regents of the University of California. All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the University of
* California, Berkeley and its contributors.
* 4. Neither the name of the University nor the names of its contributors
* may be used to endorse or promote products derived from this software
* without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/* @(#)gamma.c 8.1 (Berkeley) 6/4/93 */
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: release/9.0.0/lib/msun/bsdsrc/b_tgamma.c 176449 2008-02-22 02:26:51Z das $");
#else
__RCSID("$NetBSD: b_tgamma.c,v 1.1 2012/05/05 17:54:14 christos Exp $");
#endif
/*
* This code by P. McIlroy, Oct 1992;
*
* The financial support of UUNET Communications Services is greatfully
* acknowledged.
*/
#include "math.h"
#include "math_private.h"
/* METHOD:
* x < 0: Use reflection formula, G(x) = pi/(sin(pi*x)*x*G(x))
* At negative integers, return NaN and raise invalid.
*
* x < 6.5:
* Use argument reduction G(x+1) = xG(x) to reach the
* range [1.066124,2.066124]. Use a rational
* approximation centered at the minimum (x0+1) to
* ensure monotonicity.
*
* x >= 6.5: Use the asymptotic approximation (Stirling's formula)
* adjusted for equal-ripples:
*
* log(G(x)) ~= (x-.5)*(log(x)-1) + .5(log(2*pi)-1) + 1/x*P(1/(x*x))
*
* Keep extra precision in multiplying (x-.5)(log(x)-1), to
* avoid premature round-off.
*
* Special values:
* -Inf: return NaN and raise invalid;
* negative integer: return NaN and raise invalid;
* other x ~< 177.79: return +-0 and raise underflow;
* +-0: return +-Inf and raise divide-by-zero;
* finite x ~> 171.63: return +Inf and raise overflow;
* +Inf: return +Inf;
* NaN: return NaN.
*
* Accuracy: tgamma(x) is accurate to within
* x > 0: error provably < 0.9ulp.
* Maximum observed in 1,000,000 trials was .87ulp.
* x < 0:
* Maximum observed error < 4ulp in 1,000,000 trials.
*/
static double neg_gam(double);
static double small_gam(double);
static double smaller_gam(double);
static struct Double large_gam(double);
static struct Double ratfun_gam(double, double);
/*
* Rational approximation, A0 + x*x*P(x)/Q(x), on the interval
* [1.066.., 2.066..] accurate to 4.25e-19.
*/
#define LEFT -.3955078125 /* left boundary for rat. approx */
#define x0 .461632144968362356785 /* xmin - 1 */
#define a0_hi 0.88560319441088874992
#define a0_lo -.00000000000000004996427036469019695
#define P0 6.21389571821820863029017800727e-01
#define P1 2.65757198651533466104979197553e-01
#define P2 5.53859446429917461063308081748e-03
#define P3 1.38456698304096573887145282811e-03
#define P4 2.40659950032711365819348969808e-03
#define Q0 1.45019531250000000000000000000e+00
#define Q1 1.06258521948016171343454061571e+00
#define Q2 -2.07474561943859936441469926649e-01
#define Q3 -1.46734131782005422506287573015e-01
#define Q4 3.07878176156175520361557573779e-02
#define Q5 5.12449347980666221336054633184e-03
#define Q6 -1.76012741431666995019222898833e-03
#define Q7 9.35021023573788935372153030556e-05
#define Q8 6.13275507472443958924745652239e-06
/*
* Constants for large x approximation (x in [6, Inf])
* (Accurate to 2.8*10^-19 absolute)
*/
#define lns2pi_hi 0.418945312500000
#define lns2pi_lo -.000006779295327258219670263595
#define Pa0 8.33333333333333148296162562474e-02
#define Pa1 -2.77777777774548123579378966497e-03
#define Pa2 7.93650778754435631476282786423e-04
#define Pa3 -5.95235082566672847950717262222e-04
#define Pa4 8.41428560346653702135821806252e-04
#define Pa5 -1.89773526463879200348872089421e-03
#define Pa6 5.69394463439411649408050664078e-03
#define Pa7 -1.44705562421428915453880392761e-02
static const double zero = 0., one = 1.0, tiny = 1e-300;
double
tgamma(double x)
{
struct Double u;
if (x >= 6) {
if(x > 171.63)
return (x / zero);
u = large_gam(x);
return(__exp__D(u.a, u.b));
} else if (x >= 1.0 + LEFT + x0)
return (small_gam(x));
else if (x > 1.e-17)
return (smaller_gam(x));
else if (x > -1.e-17) {
if (x != 0.0)
u.a = one - tiny; /* raise inexact */
return (one/x);
} else if (!finite(x))
return (x - x); /* x is NaN or -Inf */
else
return (neg_gam(x));
}
/*
* Accurate to max(ulp(1/128) absolute, 2^-66 relative) error.
*/
static struct Double
large_gam(double x)
{
double z, p;
struct Double t, u, v;
z = one/(x*x);
p = Pa0+z*(Pa1+z*(Pa2+z*(Pa3+z*(Pa4+z*(Pa5+z*(Pa6+z*Pa7))))));
p = p/x;
u = __log__D(x);
u.a -= one;
v.a = (x -= .5);
TRUNC(v.a);
v.b = x - v.a;
t.a = v.a*u.a; /* t = (x-.5)*(log(x)-1) */
t.b = v.b*u.a + x*u.b;
/* return t.a + t.b + lns2pi_hi + lns2pi_lo + p */
t.b += lns2pi_lo; t.b += p;
u.a = lns2pi_hi + t.b; u.a += t.a;
u.b = t.a - u.a;
u.b += lns2pi_hi; u.b += t.b;
return (u);
}
/*
* Good to < 1 ulp. (provably .90 ulp; .87 ulp on 1,000,000 runs.)
* It also has correct monotonicity.
*/
static double
small_gam(double x)
{
double y, ym1, t;
struct Double yy, r;
y = x - one;
ym1 = y - one;
if (y <= 1.0 + (LEFT + x0)) {
yy = ratfun_gam(y - x0, 0);
return (yy.a + yy.b);
}
r.a = y;
TRUNC(r.a);
yy.a = r.a - one;
y = ym1;
yy.b = r.b = y - yy.a;
/* Argument reduction: G(x+1) = x*G(x) */
for (ym1 = y-one; ym1 > LEFT + x0; y = ym1--, yy.a--) {
t = r.a*yy.a;
r.b = r.a*yy.b + y*r.b;
r.a = t;
TRUNC(r.a);
r.b += (t - r.a);
}
/* Return r*tgamma(y). */
yy = ratfun_gam(y - x0, 0);
y = r.b*(yy.a + yy.b) + r.a*yy.b;
y += yy.a*r.a;
return (y);
}
/*
* Good on (0, 1+x0+LEFT]. Accurate to 1ulp.
*/
static double
smaller_gam(double x)
{
double t, d;
struct Double r, xx;
if (x < x0 + LEFT) {
t = x, TRUNC(t);
d = (t+x)*(x-t);
t *= t;
xx.a = (t + x), TRUNC(xx.a);
xx.b = x - xx.a; xx.b += t; xx.b += d;
t = (one-x0); t += x;
d = (one-x0); d -= t; d += x;
x = xx.a + xx.b;
} else {
xx.a = x, TRUNC(xx.a);
xx.b = x - xx.a;
t = x - x0;
d = (-x0 -t); d += x;
}
r = ratfun_gam(t, d);
d = r.a/x, TRUNC(d);
r.a -= d*xx.a; r.a -= d*xx.b; r.a += r.b;
return (d + r.a/x);
}
/*
* returns (z+c)^2 * P(z)/Q(z) + a0
*/
static struct Double
ratfun_gam(double z, double c)
{
double p, q;
struct Double r, t;
q = Q0 +z*(Q1+z*(Q2+z*(Q3+z*(Q4+z*(Q5+z*(Q6+z*(Q7+z*Q8)))))));
p = P0 + z*(P1 + z*(P2 + z*(P3 + z*P4)));
/* return r.a + r.b = a0 + (z+c)^2*p/q, with r.a truncated to 26 bits. */
p = p/q;
t.a = z, TRUNC(t.a); /* t ~= z + c */
t.b = (z - t.a) + c;
t.b *= (t.a + z);
q = (t.a *= t.a); /* t = (z+c)^2 */
TRUNC(t.a);
t.b += (q - t.a);
r.a = p, TRUNC(r.a); /* r = P/Q */
r.b = p - r.a;
t.b = t.b*p + t.a*r.b + a0_lo;
t.a *= r.a; /* t = (z+c)^2*(P/Q) */
r.a = t.a + a0_hi, TRUNC(r.a);
r.b = ((a0_hi-r.a) + t.a) + t.b;
return (r); /* r = a0 + t */
}
static double
neg_gam(double x)
{
int sgn = 1;
struct Double lg, lsine;
double y, z;
y = ceil(x);
if (y == x) /* Negative integer. */
return ((x - x) / zero);
z = y - x;
if (z > 0.5)
z = one - z;
y = 0.5 * y;
if (y == ceil(y))
sgn = -1;
if (z < .25)
z = sin(M_PI*z);
else
z = cos(M_PI*(0.5-z));
/* Special case: G(1-x) = Inf; G(x) may be nonzero. */
if (x < -170) {
if (x < -190)
return ((double)sgn*tiny*tiny);
y = one - x; /* exact: 128 < |x| < 255 */
lg = large_gam(y);
lsine = __log__D(M_PI/z); /* = TRUNC(log(u)) + small */
lg.a -= lsine.a; /* exact (opposite signs) */
lg.b -= lsine.b;
y = -(lg.a + lg.b);
z = (y + lg.a) + lg.b;
y = __exp__D(y, z);
if (sgn < 0) y = -y;
return (y);
}
y = one-x;
if (one-y == x)
y = tgamma(y);
else /* 1-x is inexact */
y = -x*tgamma(-x);
if (sgn < 0) y = -y;
return (M_PI / (y*z));
}

6
lib/libm/src/e_fmod.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $");
__RCSID("$NetBSD: e_fmod.c,v 1.12 2013/11/19 14:02:59 joerg Exp $");
#endif
/*
@@ -26,6 +26,10 @@ __RCSID("$NetBSD: e_fmod.c,v 1.11 2002/05/26 22:01:49 wiz Exp $");
static const double one = 1.0, Zero[] = {0.0, -0.0,};
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(__ieee754_fmodl, __ieee754_fmod)
#endif
double
__ieee754_fmod(double x, double y)
{

158
lib/libm/src/e_fmodl.c Normal file
View File

@@ -0,0 +1,158 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: e_fmodl.c,v 1.2 2013/11/14 15:25:22 martin Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/e_fmodl.c 181063 2008-07-31 20:09:47Z das $");
#endif
#include "namespace.h"
#include <float.h>
#include <stdint.h>
#include "math.h"
#include "math_private.h"
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#define BIAS (LDBL_MAX_EXP - 1)
#if EXT_FRACLBITS > 32
typedef uint64_t manl_t;
#else
typedef uint32_t manl_t;
#endif
#if EXT_FRACHBITS > 32
typedef uint64_t manh_t;
#else
typedef uint32_t manh_t;
#endif
/*
* These macros add and remove an explicit integer bit in front of the
* fractional mantissa, if the architecture doesn't have such a bit by
* default already.
*/
#ifdef LDBL_IMPLICIT_NBIT
#define SET_NBIT(hx) ((hx) | (1ULL << EXT_FRACHBITS))
#define HFRAC_BITS EXT_FRACHBITS
#define LDBL_NBIT 0
#else
#define SET_NBIT(hx) (hx)
#define HFRAC_BITS (EXT_FRACHBITS - 1)
#endif
#define MANL_SHIFT (EXT_FRACLBITS - 1)
static const long double one = 1.0, Zero[] = {0.0, -0.0,};
/*
* fmodl(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*
* Assumptions:
* - The low part of the mantissa fits in a manl_t exactly.
* - The high part of the mantissa fits in an int64_t with enough room
* for an explicit integer bit in front of the fractional bits.
*/
long double
__ieee754_fmodl(long double x, long double y)
{
union ieee_ext_u ux = { .extu_ld = x, };
union ieee_ext_u uy = { .extu_ld = y, };
int64_t hx,hz; /* We need a carry bit even if EXT_FRACHBITS is 32. */
manh_t hy;
manl_t lx,ly,lz;
int ix,iy,n,sx;
sx = ux.extu_sign;
/* purge off exception values */
if((uy.extu_exp|uy.extu_frach|uy.extu_fracl)==0 || /* y=0 */
(ux.extu_exp == BIAS + LDBL_MAX_EXP) || /* or x not finite */
(uy.extu_exp == BIAS + LDBL_MAX_EXP &&
((uy.extu_frach&~LDBL_NBIT)|uy.extu_fracl)!=0)) /* or y is NaN */
return (x*y)/(x*y);
if(ux.extu_exp<=uy.extu_exp) {
if((ux.extu_exp<uy.extu_exp) ||
(ux.extu_frach<=uy.extu_frach &&
(ux.extu_frach<uy.extu_frach ||
ux.extu_fracl<uy.extu_fracl))) {
return x; /* |x|<|y| return x or x-y */
}
if(ux.extu_frach==uy.extu_frach && ux.extu_fracl==uy.extu_fracl) {
return Zero[sx]; /* |x|=|y| return x*0*/
}
}
/* determine ix = ilogb(x) */
if(ux.extu_exp == 0) { /* subnormal x */
ux.extu_ld *= 0x1.0p512;
ix = ux.extu_exp - (BIAS + 512);
} else {
ix = ux.extu_exp - BIAS;
}
/* determine iy = ilogb(y) */
if(uy.extu_exp == 0) { /* subnormal y */
uy.extu_ld *= 0x1.0p512;
iy = uy.extu_exp - (BIAS + 512);
} else {
iy = uy.extu_exp - BIAS;
}
/* set up {hx,lx}, {hy,ly} and align y to x */
hx = SET_NBIT(ux.extu_frach);
hy = SET_NBIT(uy.extu_frach);
lx = ux.extu_fracl;
ly = uy.extu_fracl;
/* fix point fmod */
n = ix - iy;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;}
else {
if ((hz|lz)==0) /* return sign(x)*0 */
return Zero[sx];
hx = hz+hz+(lz>>MANL_SHIFT); lx = lz+lz;
}
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) /* return sign(x)*0 */
return Zero[sx];
while(hx<(1LL<<HFRAC_BITS)) { /* normalize x */
hx = hx+hx+(lx>>MANL_SHIFT); lx = lx+lx;
iy -= 1;
}
ux.extu_frach = hx; /* The mantissa is truncated here if needed. */
ux.extu_fracl = lx;
if (iy < LDBL_MIN_EXP) {
ux.extu_exp = iy + (BIAS + 512);
ux.extu_ld *= 0x1p-512;
} else {
ux.extu_exp = iy + BIAS;
}
x = ux.extu_ld * one; /* create necessary signal */
return x; /* exact output */
}
#endif /* __HAVE_LONG_DOUBLE */

184
lib/libm/src/e_sqrtl.c Normal file
View File

@@ -0,0 +1,184 @@
/*-
* Copyright (c) 2007 Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/e_sqrtl.c 176720 2008-03-02 01:47:58Z das $");
#endif
__RCSID("$NetBSD: e_sqrtl.c,v 1.4 2013/11/22 20:15:06 martin Exp $");
#include <machine/ieee.h>
#include <float.h>
#include "math.h"
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
#ifdef HAVE_FENV_H
#include <fenv.h>
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#endif
#ifdef HAVE_FENV_H
/* Return (x + ulp) for normal positive x. Assumes no overflow. */
static inline long double
inc(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
if (++ux.extu_fracl == 0) {
if (++ux.extu_frach == 0) {
ux.extu_exp++;
ux.extu_frach |= LDBL_NBIT;
}
}
return (ux.extu_ld);
}
/* Return (x - ulp) for normal positive x. Assumes no underflow. */
static inline long double
dec(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
if (ux.extu_fracl-- == 0) {
if (ux.extu_frach-- == LDBL_NBIT) {
ux.extu_exp--;
ux.extu_frach |= LDBL_NBIT;
}
}
return (ux.extu_ld);
}
/*
* This is slow, but simple and portable. You should use hardware sqrt
* if possible.
*/
long double
__ieee754_sqrtl(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
int k, r;
long double lo, xn;
fenv_t env;
/* If x = NaN, then sqrt(x) = NaN. */
/* If x = Inf, then sqrt(x) = Inf. */
/* If x = -Inf, then sqrt(x) = NaN. */
if (ux.extu_exp == LDBL_MAX_EXP * 2 - 1)
return (x * x + x);
/* If x = +-0, then sqrt(x) = +-0. */
if ((ux.extu_frach | ux.extu_fracl | ux.extu_exp) == 0)
return (x);
/* If x < 0, then raise invalid and return NaN */
if (ux.extu_sign)
return ((x - x) / (x - x));
feholdexcept(&env);
if (ux.extu_exp == 0) {
/* Adjust subnormal numbers. */
ux.extu_ld *= 0x1.0p514;
k = -514;
} else {
k = 0;
}
/*
* ux.extu_ld is a normal number, so break it into ux.extu_ld = e*2^n where
* ux.extu_ld = (2*e)*2^2k for odd n and ux.extu_ld = (4*e)*2^2k for even n.
*/
if ((ux.extu_exp - EXT_EXP_BIAS) & 1) { /* n is even. */
k += ux.extu_exp - EXT_EXP_BIAS - 1; /* 2k = n - 2. */
ux.extu_exp = EXT_EXP_BIAS + 1; /* ux.extu_ld in [2,4). */
} else {
k += ux.extu_exp - EXT_EXP_BIAS; /* 2k = n - 1. */
ux.extu_exp = EXT_EXP_BIAS; /* ux.extu_ld in [1,2). */
}
/*
* Newton's iteration.
* Split ux.extu_ld into a high and low part to achieve additional precision.
*/
xn = sqrt(ux.extu_ld); /* 53-bit estimate of sqrtl(x). */
#if LDBL_MANT_DIG > 100
xn = (xn + (ux.extu_ld / xn)) * 0.5; /* 106-bit estimate. */
#endif
lo = ux.extu_ld;
ux.extu_fracl = 0; /* Zero out lower bits. */
lo = (lo - ux.extu_ld) / xn; /* Low bits divided by xn. */
xn = xn + (ux.extu_ld / xn); /* High portion of estimate. */
ux.extu_ld = xn + lo; /* Combine everything. */
ux.extu_exp += (k >> 1) - 1;
feclearexcept(FE_INEXACT);
r = fegetround();
fesetround(FE_TOWARDZERO); /* Set to round-toward-zero. */
xn = x / ux.extu_ld; /* Chopped quotient (inexact?). */
if (!fetestexcept(FE_INEXACT)) { /* Quotient is exact. */
if (xn == ux.extu_ld) {
fesetenv(&env);
return (ux.extu_ld);
}
/* Round correctly for inputs like x = y**2 - ulp. */
xn = dec(xn); /* xn = xn - ulp. */
}
if (r == FE_TONEAREST) {
xn = inc(xn); /* xn = xn + ulp. */
} else if (r == FE_UPWARD) {
ux.extu_ld = inc(ux.extu_ld); /* ux.extu_ld = ux.extu_ld + ulp. */
xn = inc(xn); /* xn = xn + ulp. */
}
ux.extu_ld = ux.extu_ld + xn; /* Chopped sum. */
feupdateenv(&env); /* Restore env and raise inexact */
ux.extu_exp--;
return (ux.extu_ld);
}
#else
/*
* No fenv support:
* poor man's version: just use double
*/
long double
__ieee754_sqrtl(long double x)
{
return __ieee754_sqrt((double)x);
}
#endif
#endif

24
lib/libm/src/k_standard.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: k_standard.c,v 1.16 2010/09/01 10:44:28 drochner Exp $");
__RCSID("$NetBSD: k_standard.c,v 1.19 2013/11/19 19:24:34 joerg Exp $");
#endif
#include "math.h"
@@ -370,7 +370,7 @@ __kernel_standard(double x, double y, int type)
if (_LIB_VERSION == _SVID_)
exc.retval = -HUGE;
else
exc.retval = -HUGE_VAL;
exc.retval = zero/zero;
if (_LIB_VERSION == _POSIX_)
errno = EDOM;
else if (!matherr(&exc)) {
@@ -406,7 +406,7 @@ __kernel_standard(double x, double y, int type)
if (_LIB_VERSION == _SVID_)
exc.retval = -HUGE;
else
exc.retval = -HUGE_VAL;
exc.retval = zero/zero;
if (_LIB_VERSION == _POSIX_)
errno = EDOM;
else if (!matherr(&exc)) {
@@ -514,9 +514,15 @@ __kernel_standard(double x, double y, int type)
break;
case 26:
case 126:
case 226:
/* sqrt(x<0) */
exc.type = DOMAIN;
exc.name = type < 100 ? "sqrt" : "sqrtf";
if (type == 26)
exc.name = "sqrt";
else if (type == 126)
exc.name = "sqrtf";
else
exc.name = "sqrtl";
if (_LIB_VERSION == _SVID_)
exc.retval = zero;
else
@@ -532,9 +538,15 @@ __kernel_standard(double x, double y, int type)
break;
case 27:
case 127:
case 227:
/* fmod(x,0) */
exc.type = DOMAIN;
exc.name = type < 100 ? "fmod" : "fmodf";
if (type == 27)
exc.name = "fmod";
else if (type == 127)
exc.name = "fmodf";
else
exc.name = "fmodl";
if (_LIB_VERSION == _SVID_)
exc.retval = x;
else
@@ -800,7 +812,7 @@ __kernel_standard(double x, double y, int type)
if (_LIB_VERSION == _SVID_)
exc.retval = -HUGE;
else
exc.retval = -HUGE_VAL;
exc.retval = zero/zero;
if (_LIB_VERSION == _POSIX_)
errno = EDOM;
else if (!matherr(&exc)) {

175
lib/libm/src/ldbl_dummy.c Normal file
View File

@@ -0,0 +1,175 @@
/* $NetBSD: ldbl_dummy.c,v 1.1 2013/11/12 17:36:14 joerg Exp $ */
/*-
* Copyright (c) 2013 The NetBSD Foundation, Inc.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
/*
* Simple long double -> double wrappers for various transcendental functions.
* They work neither on the additional range of long double nor do they use
* the additional precision. They exist as stop gap fix for various programs
* picking up long double, e.g. via the C++ run time.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: ldbl_dummy.c,v 1.1 2013/11/12 17:36:14 joerg Exp $");
#include "namespace.h"
#include <math.h>
__weak_alias(atan2l, _atan2l)
__weak_alias(hypotl, _hypotl)
__weak_alias(logl, _logl)
__weak_alias(log10l, _log10l)
__weak_alias(expl, _expl)
__weak_alias(exp2l, _exp2l)
__weak_alias(powl, _powl)
__weak_alias(cosl, _cosl)
__weak_alias(sinl, _sinl)
__weak_alias(tanl, _tanl)
__weak_alias(coshl, _coshl)
__weak_alias(sinhl, _sinhl)
__weak_alias(tanhl, _tanhl)
__weak_alias(acosl, _acosl)
__weak_alias(asinl, _asinl)
__weak_alias(atanl, _atanl)
__weak_alias(acoshl, _acoshl)
__weak_alias(asinhl, _asinhl)
__weak_alias(atanhl, _atanhl)
long double
atan2l(long double y, long double x)
{
return atan2(y, x);
}
long double
hypotl(long double x, long double y)
{
return hypot(x, y);
}
long double
logl(long double x)
{
return log(x);
}
long double
log10l(long double x)
{
return log10(x);
}
long double
expl(long double x)
{
return exp(x);
}
long double
exp2l(long double x)
{
return exp2(x);
}
long double
powl(long double x, long double y)
{
return pow(x, y);
}
long double
cosl(long double x)
{
return cos(x);
}
long double
sinl(long double x)
{
return sin(x);
}
long double
tanl(long double x)
{
return tan(x);
}
long double
sinhl(long double x)
{
return sinh(x);
}
long double
coshl(long double x)
{
return cosh(x);
}
long double
tanhl(long double x)
{
return tanh(x);
}
long double
acosl(long double x)
{
return acos(x);
}
long double
asinl(long double x)
{
return asin(x);
}
long double
atanl(long double x)
{
return atan(x);
}
long double
asinhl(long double x)
{
return asinh(x);
}
long double
acoshl(long double x)
{
return acosh(x);
}
long double
atanhl(long double x)
{
return atanh(x);
}

59
lib/libm/src/math_private.h
View File

@@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
* $NetBSD: math_private.h,v 1.16 2010/09/16 20:39:50 drochner Exp $
* $NetBSD: math_private.h,v 1.20 2013/11/19 19:24:34 joerg Exp $
*/
#ifndef _MATH_PRIVATE_H_
@@ -48,6 +48,9 @@ typedef union
u_int32_t msw;
u_int32_t lsw;
} parts;
struct {
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
@@ -63,6 +66,9 @@ typedef union
u_int32_t lsw;
u_int32_t msw;
} parts;
struct {
u_int64_t w;
} xparts;
} ieee_double_shape_type;
#endif
@@ -77,6 +83,15 @@ do { \
(ix1) = ew_u.parts.lsw; \
} while (/*CONSTCOND*/0)
/* Get a 64-bit int from a double. */
#define EXTRACT_WORD64(ix,d) \
do { \
ieee_double_shape_type ew_u; \
ew_u.value = (d); \
(ix) = ew_u.xparts.w; \
} while (/*CONSTCOND*/0)
/* Get the more significant 32 bit int from a double. */
#define GET_HIGH_WORD(i,d) \
@@ -105,6 +120,15 @@ do { \
(d) = iw_u.value; \
} while (/*CONSTCOND*/0)
/* Set a double from a 64-bit int. */
#define INSERT_WORD64(d,ix) \
do { \
ieee_double_shape_type iw_u; \
iw_u.xparts.w = (ix); \
(d) = iw_u.value; \
} while (/*CONSTCOND*/0)
/* Set the more significant 32 bits of a double from an int. */
#define SET_HIGH_WORD(d,v) \
@@ -275,4 +299,37 @@ extern float __kernel_cosf __P((float,float));
extern float __kernel_tanf __P((float,float,int));
extern int __kernel_rem_pio2f __P((float*,float*,int,int,int,const int*));
/* ieee style elementary long double functions */
extern long double __ieee754_fmodl(long double, long double);
extern long double __ieee754_sqrtl(long double);
/*
* TRUNC() is a macro that sets the trailing 27 bits in the mantissa of an
* IEEE double variable to zero. It must be expression-like for syntactic
* reasons, and we implement this expression using an inline function
* instead of a pure macro to avoid depending on the gcc feature of
* statement-expressions.
*/
#define TRUNC(d) (_b_trunc(&(d)))
static __inline void
_b_trunc(volatile double *_dp)
{
uint32_t _lw;
GET_LOW_WORD(_lw, *_dp);
SET_LOW_WORD(*_dp, _lw & 0xf8000000);
}
struct Double {
double a;
double b;
};
/*
* Functions internal to the math package, yet not static.
*/
double __exp__D(double, double);
struct Double __log__D(double);
#endif /* _MATH_PRIVATE_H_ */

38
lib/libm/src/namespace.h
View File

@@ -1,14 +1,18 @@
/* $NetBSD: namespace.h,v 1.3 2010/04/23 19:17:07 drochner Exp $ */
/* $NetBSD: namespace.h,v 1.10 2013/11/19 19:24:34 joerg Exp $ */
#define atan2 _atan2
#define atan2f _atan2f
#define atan2l _atan2l
#define hypot _hypot
#define hypotf _hypotf
#define hypotl _hypotl
#define exp _exp
#define expf _expf
#define expl _expl
#define log _log
#define logf _logf
#define logl _logl
#if 0 /* not yet - need to review use in machdep code first */
#define sin _sin
@@ -20,15 +24,47 @@
#endif /* notyet */
#define sinh _sinh
#define sinhf _sinhf
#define sinhl _sinhl
#define cosh _cosh
#define coshf _coshf
#define coshl _coshl
#define asin _asin
#define asinf _asinf
#define asinl _asinl
#define casin _casin
#define casinf _casinf
#define casinl _casinl
#define catan _catan
#define catanf _catanf
#define catanl _catanl
#define scalbn _scalbn
#define scalbnf _scalbnf
#define scalbnl _scalbnl
#define scalbln _scalbln
#define scalblnf _scalblnf
#define scalblnl _scalblnl
#define sqrtl _sqrtl
#define cbrtl _cbrtl
#define ceill _ceill
#define floorl _floorl
#define roundl _roundl
#define fmodl _fmodl
#define truncl _truncl
#define exp2l _exp2l
#define cosl _cosl
#define sinl _sinl
#define tanl _tanl
#define powl _powl
#define coshl _coshl
#define sinhl _sinhl
#define acosl _acosl
#define atanl _atanl
#define asinhl _asinhl
#define acoshl _acoshl
#define tanhl _tanhl
#define atanhl _atanhl
#define log10l _log10l

8
lib/libm/src/s_cbrt.c
View File

@@ -12,12 +12,18 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_cbrt.c,v 1.11 2002/05/26 22:01:54 wiz Exp $");
__RCSID("$NetBSD: s_cbrt.c,v 1.13 2013/11/20 12:29:13 joerg Exp $");
#endif
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_cbrtl, cbrt)
__weak_alias(cbrtl, _cbrtl)
#endif
/* cbrt(x)
* Return cube root of x
*/

145
lib/libm/src/s_cbrtl.c Normal file
View File

@@ -0,0 +1,145 @@
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
* Copyright (c) 2009-2011, Bruce D. Evans, Steven G. Kargl, David Schultz.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* The argument reduction and testing for exceptional cases was
* written by Steven G. Kargl with input from Bruce D. Evans
* and David A. Schultz.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_cbrtl.c,v 1.1 2013/11/19 19:24:34 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/s_cbrtl.c 238924 2012-07-30 21:58:28Z kargl $");
#endif
#include "namespace.h"
#include <machine/ieee.h>
#include <float.h>
#include "math.h"
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
__weak_alias(cbrtl, _cbrtl)
static const unsigned
B1 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
long double
cbrtl(long double x)
{
union ieee_ext_u ux, vx;
long double r, s, t, w;
double dr, dt, dx;
float ft, fx;
uint32_t hx;
int k;
ux.extu_ld = x;
/*
* If x = +-Inf, then cbrt(x) = +-Inf.
* If x = NaN, then cbrt(x) = NaN.
*/
if (ux.extu_exp == EXT_EXP_INFNAN)
return (x + x);
if ((ux.extu_frach | ux.extu_fracl | ux.extu_exp) == 0)
return (x);
vx.extu_ld = 1;
vx.extu_ext.ext_sign = ux.extu_ext.ext_sign;
ux.extu_ext.ext_sign = 0;
if (ux.extu_exp == 0) {
/* Adjust subnormal numbers. */
ux.extu_ld *= 0x1.0p514;
k = ux.extu_exp - EXT_EXP_BIAS - 514;
} else {
k = ux.extu_exp - EXT_EXP_BIAS;
}
ux.extu_exp = EXT_EXP_BIAS;
x = ux.extu_ld;
switch (k % 3) {
case 1:
case -2:
x = 2*x;
k--;
break;
case 2:
case -1:
x = 4*x;
k -= 2;
break;
}
vx.extu_exp = EXT_EXP_BIAS + k / 3;
/*
* The following is the guts of s_cbrtf, with the handling of
* special values removed and extra care for accuracy not taken,
* but with most of the extra accuracy not discarded.
*/
/* ~5-bit estimate: */
fx = x;
GET_FLOAT_WORD(hx, fx);
SET_FLOAT_WORD(ft, ((hx & 0x7fffffff) / 3 + B1));
/* ~16-bit estimate: */
dx = x;
dt = ft;
dr = dt * dt * dt;
dt = dt * (dx + dx + dr) / (dx + dr + dr);
/* ~47-bit estimate: */
dr = dt * dt * dt;
dt = dt * (dx + dx + dr) / (dx + dr + dr);
#if LDBL_MANT_DIG == 64
/*
* dt is cbrtl(x) to ~47 bits (after x has been reduced to 1 <= x < 8).
* Round it away from zero to 32 bits (32 so that t*t is exact, and
* away from zero for technical reasons).
*/
volatile double vd2 = 0x1.0p32;
volatile double vd1 = 0x1.0p-31;
#define vd ((long double)vd2 + vd1)
t = dt + vd - 0x1.0p32;
#elif LDBL_MANT_DIG == 113
/*
* Round dt away from zero to 47 bits. Since we don't trust the 47,
* add 2 47-bit ulps instead of 1 to round up. Rounding is slow and
* might be avoidable in this case, since on most machines dt will
* have been evaluated in 53-bit precision and the technical reasons
* for rounding up might not apply to either case in cbrtl() since
* dt is much more accurate than needed.
*/
t = dt + 0x2.0p-46 + 0x1.0p60L - 0x1.0p60;
#else
#error "Unsupported long double format"
#endif
/*
* Final step Newton iteration to 64 or 113 bits with
* error < 0.667 ulps
*/
s=t*t; /* t*t is exact */
r=x/s; /* error <= 0.5 ulps; |r| < |t| */
w=t+t; /* t+t is exact */
r=(r-t)/(w+r); /* r-t is exact; w+r ~= 3*t */
t=t+t*r; /* error <= 0.5 + 0.5/3 + epsilon */
t *= vx.extu_ld;
return t;
}
#endif /* __HAVE_LONG_DOUBLE */

7
lib/libm/src/s_ceil.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_ceil.c,v 1.13 2009/02/16 01:22:18 lukem Exp $");
__RCSID("$NetBSD: s_ceil.c,v 1.14 2013/11/11 23:57:34 joerg Exp $");
#endif
/*
@@ -29,6 +29,11 @@ __RCSID("$NetBSD: s_ceil.c,v 1.13 2009/02/16 01:22:18 lukem Exp $");
static const double huge = 1.0e300;
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_ceill, ceil)
__weak_alias(ceill, ceil)
#endif
double
ceil(double x)
{

112
lib/libm/src/s_ceill.c Normal file
View File

@@ -0,0 +1,112 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* From: @(#)s_ceil.c 5.1 93/09/24
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_ceill.c,v 1.1 2013/11/11 23:57:34 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/s_ceill.c 176280 2008-02-14 15:10:34Z bde $");
#endif
/*
* ceill(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to ceill(x).
*/
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef __weak_alias
__weak_alias(ceill, _ceill)
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define MANH_SIZE (EXT_FRACHBITS + 1)
#define INC_MANH(ux, c) do { \
uint64_t oi = ux.extu_frach; \
ux.extu_frach += (c); \
if (ux.extu_frach < oi) \
ux.extu_exp++; \
} while (0)
#else
#define MANH_SIZE EXT_FRACHBITS
#define INC_MANH(ux, c) do { \
uint64_t oi = ux.extu_frach; \
ux.extu_frach += (c); \
if (ux.extu_frach < oi) { \
ux.extu_exp++; \
ux.extu_frach |= 1llu << (EXT_FRACHBITS - 1); \
} \
} while (0)
#endif
static const long double huge = 1.0e300;
long double
ceill(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
int e = ux.extu_exp - LDBL_MAX_EXP + 1;
if (e < MANH_SIZE - 1) {
if (e < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0)
if (ux.extu_exp > 0 ||
(ux.extu_frach | ux.extu_fracl) != 0)
ux.extu_ld = ux.extu_sign ? -0.0 : 1.0;
} else {
uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
if (((ux.extu_frach & m) | ux.extu_fracl) == 0)
return (x); /* x is integral */
if (!ux.extu_sign) {
#ifdef LDBL_IMPLICIT_NBIT
if (e == 0)
ux.extu_exp++;
else
#endif
INC_MANH(ux, 1llu << (MANH_SIZE - e - 1));
}
if (huge + x > 0.0) { /* raise inexact flag */
ux.extu_frach &= ~m;
ux.extu_fracl = 0;
}
}
} else if (e < LDBL_MANT_DIG - 1) {
uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
if ((ux.extu_fracl & m) == 0)
return (x); /* x is integral */
if (!ux.extu_sign) {
if (e == MANH_SIZE - 1)
INC_MANH(ux, 1);
else {
uint64_t o = ux.extu_fracl;
ux.extu_fracl += 1llu << (LDBL_MANT_DIG - e - 1);
if (ux.extu_fracl < o) /* got a carry */
INC_MANH(ux, 1);
}
}
if (huge + x > 0.0) /* raise inexact flag */
ux.extu_fracl &= ~m;
}
return (ux.extu_ld);
}
#endif /* __HAVE_LONG_DOUBLE */

15
lib/libm/src/s_exp2.c
View File

@@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_exp2.c,v 1.2 2010/01/11 23:38:24 christos Exp $");
__RCSID("$NetBSD: s_exp2.c,v 1.3 2013/11/06 16:49:21 christos Exp $");
#ifdef __FBSDID
__FBSDID("$FreeBSD: src/lib/msun/src/s_exp2.c,v 1.7 2008/02/22 02:27:34 das Exp $");
#endif
@@ -343,9 +343,9 @@ static const double tbl[TBLSIZE * 2] = {
double
exp2(double x)
{
double r, t, twopk, twopkp1000, z;
double r, t, twopk, z;
uint32_t hx, ix, lx, i0;
int k;
int k, big;
/* Filter out exceptional cases. */
GET_HIGH_WORD(hx,x);
@@ -378,19 +378,20 @@ exp2(double x)
/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
t = tbl[i0]; /* exp2t[i0] */
z -= tbl[i0 + 1]; /* eps[i0] */
if (k >= -1021 << 20)
big = k >= -1021 << 20;
if (big)
INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
else
INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0);
INSERT_WORDS(twopk, 0x3ff00000 + k + (1000 << 20), 0);
r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
/* Scale by 2**(k>>20). */
if(k >= -1021 << 20) {
if (big) {
if (k == 1024 << 20)
return (r * 2.0 * 0x1p1023);
return (r * twopk);
} else {
return (r * twopkp1000 * twom1000);
return (r * twopk * twom1000);
}
}

6
lib/libm/src/s_fabs.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_fabs.c,v 1.10 2002/05/26 22:01:55 wiz Exp $");
__RCSID("$NetBSD: s_fabs.c,v 1.11 2013/11/29 22:16:10 joerg Exp $");
#endif
/*
@@ -22,6 +22,10 @@ __RCSID("$NetBSD: s_fabs.c,v 1.10 2002/05/26 22:01:55 wiz Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(fabsl, fabs)
#endif
double
fabs(double x)
{

6
lib/libm/src/s_fabsl.c
View File

@@ -1,4 +1,4 @@
/* $NetBSD: s_fabsl.c,v 1.2 2010/09/17 20:39:39 christos Exp $ */
/* $NetBSD: s_fabsl.c,v 1.3 2013/11/29 22:16:10 joerg Exp $ */
/*-
* Copyright (c) 2010 The NetBSD Foundation, Inc.
@@ -26,7 +26,7 @@
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_fabsl.c,v 1.2 2010/09/17 20:39:39 christos Exp $");
__RCSID("$NetBSD: s_fabsl.c,v 1.3 2013/11/29 22:16:10 joerg Exp $");
#include <math.h>
#include <machine/ieee.h>
@@ -35,7 +35,7 @@ __RCSID("$NetBSD: s_fabsl.c,v 1.2 2010/09/17 20:39:39 christos Exp $");
* fabsl(long double x)
* This function returns the absolute value of its argumetn x, |x|.
*/
#ifdef EXT_EXP_INFNAN
#ifdef __HAVE_LONG_DOUBLE
long double
fabsl(long double x)
{

4
lib/libm/src/s_finite.c
View File

@@ -12,11 +12,11 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_finite.c,v 1.11 2002/05/26 22:01:55 wiz Exp $");
__RCSID("$NetBSD: s_finite.c,v 1.12 2013/11/12 17:37:43 joerg Exp $");
#endif
/*
* finite(x) returns 1 is x is finite, else 0;
* finite(x) returns 1 if x is finite, else 0;
* no branching!
*/

4
lib/libm/src/s_finitef.c
View File

@@ -15,11 +15,11 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_finitef.c,v 1.7 2002/05/26 22:01:55 wiz Exp $");
__RCSID("$NetBSD: s_finitef.c,v 1.8 2013/11/12 17:37:43 joerg Exp $");
#endif
/*
* finitef(x) returns 1 is x is finite, else 0;
* finitef(x) returns 1 if x is finite, else 0;
* no branching!
*/

7
lib/libm/src/s_floor.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_floor.c,v 1.13 2009/02/16 01:27:36 lukem Exp $");
__RCSID("$NetBSD: s_floor.c,v 1.14 2013/11/11 23:57:34 joerg Exp $");
#endif
/*
@@ -27,6 +27,11 @@ __RCSID("$NetBSD: s_floor.c,v 1.13 2009/02/16 01:27:36 lukem Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_floorl, floor)
__weak_alias(floorl, floor)
#endif
static const double huge = 1.0e300;
double

112
lib/libm/src/s_floorl.c Normal file
View File

@@ -0,0 +1,112 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* From: @(#)s_floor.c 5.1 93/09/24
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_floorl.c,v 1.1 2013/11/11 23:57:34 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/s_floorl.c 176280 2008-02-14 15:10:34Z bde $");
#endif
/*
* floorl(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to floorl(x).
*/
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef __weak_alias
__weak_alias(floorl, _floorl)
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define MANH_SIZE (EXT_FRACHBITS + 1)
#define INC_MANH(ux, c) do { \
uint64_t oi = ux.extu_frach; \
ux.extu_frach += (c); \
if (ux.extu_frach < oi) \
ux.extu_exp++; \
} while (0)
#else
#define MANH_SIZE EXT_FRACHBITS
#define INC_MANH(ux, c) do { \
uint64_t oi = ux.extu_frach; \
ux.extu_frach += (c); \
if (ux.extu_frach < oi) { \
ux.extu_exp++; \
ux.extu_frach |= 1llu << (EXT_FRACHBITS - 1); \
} \
} while (0)
#endif
static const long double huge = 1.0e300;
long double
floorl(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
int e = ux.extu_exp - LDBL_MAX_EXP + 1;
if (e < MANH_SIZE - 1) {
if (e < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0)
if (ux.extu_exp > 0 ||
(ux.extu_frach | ux.extu_fracl) != 0)
ux.extu_ld = ux.extu_sign ? -1.0 : 0.0;
} else {
uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
if (((ux.extu_frach & m) | ux.extu_fracl) == 0)
return (x); /* x is integral */
if (ux.extu_sign) {
#ifdef LDBL_IMPLICIT_NBIT
if (e == 0)
ux.extu_exp++;
else
#endif
INC_MANH(ux, 1llu << (MANH_SIZE - e - 1));
}
if (huge + x > 0.0) { /* raise inexact flag */
ux.extu_frach &= ~m;
ux.extu_fracl = 0;
}
}
} else if (e < LDBL_MANT_DIG - 1) {
uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
if ((ux.extu_fracl & m) == 0)
return (x); /* x is integral */
if (ux.extu_sign) {
if (e == MANH_SIZE - 1)
INC_MANH(ux, 1);
else {
uint64_t o = ux.extu_fracl;
ux.extu_fracl += 1llu << (LDBL_MANT_DIG - e - 1);
if (ux.extu_fracl < o) /* got a carry */
INC_MANH(ux, 1);
}
}
if (huge + x > 0.0) /* raise inexact flag */
ux.extu_fracl &= ~m;
}
return (ux.extu_ld);
}
#endif /* __HAVE_LONG_DOUBLE */

293
lib/libm/src/s_fma.c Normal file
View File

@@ -0,0 +1,293 @@
/* $NetBSD: s_fma.c,v 1.6 2013/02/14 09:24:50 matt Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fma.c,v 1.8 2011/10/21 06:30:43 das Exp $");
#else
__RCSID("$NetBSD: s_fma.c,v 1.6 2013/02/14 09:24:50 matt Exp $");
#endif
#include <machine/ieee.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(fmal, fma)
#endif
/*
* A struct dd represents a floating-point number with twice the precision
* of a double. We maintain the invariant that "hi" stores the 53 high-order
* bits of the result.
*/
struct dd {
double hi;
double lo;
};
/*
* Compute a+b exactly, returning the exact result in a struct dd. We assume
* that both a and b are finite, but make no assumptions about their relative
* magnitudes.
*/
static inline struct dd
dd_add(double a, double b)
{
struct dd ret;
double s;
ret.hi = a + b;
s = ret.hi - a;
ret.lo = (a - (ret.hi - s)) + (b - s);
return (ret);
}
/*
* Compute a+b, with a small tweak: The least significant bit of the
* result is adjusted into a sticky bit summarizing all the bits that
* were lost to rounding. This adjustment negates the effects of double
* rounding when the result is added to another number with a higher
* exponent. For an explanation of round and sticky bits, see any reference
* on FPU design, e.g.,
*
* J. Coonen. An Implementation Guide to a Proposed Standard for
* Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
*/
static inline double
add_adjusted(double a, double b)
{
struct dd sum;
uint64_t hibits, lobits;
sum = dd_add(a, b);
if (sum.lo != 0) {
EXTRACT_WORD64(hibits, sum.hi);
if ((hibits & 1) == 0) {
/* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
EXTRACT_WORD64(lobits, sum.lo);
hibits += 1 - ((hibits ^ lobits) >> 62);
INSERT_WORD64(sum.hi, hibits);
}
}
return (sum.hi);
}
/*
* Compute ldexp(a+b, scale) with a single rounding error. It is assumed
* that the result will be subnormal, and care is taken to ensure that
* double rounding does not occur.
*/
static inline double
add_and_denormalize(double a, double b, int scale)
{
struct dd sum;
uint64_t hibits, lobits;
int bits_lost;
sum = dd_add(a, b);
/*
* If we are losing at least two bits of accuracy to denormalization,
* then the first lost bit becomes a round bit, and we adjust the
* lowest bit of sum.hi to make it a sticky bit summarizing all the
* bits in sum.lo. With the sticky bit adjusted, the hardware will
* break any ties in the correct direction.
*
* If we are losing only one bit to denormalization, however, we must
* break the ties manually.
*/
if (sum.lo != 0) {
EXTRACT_WORD64(hibits, sum.hi);
bits_lost = -((int)(hibits >> 52) & 0x7ff) - scale + 1;
if ((bits_lost != 1) ^ (int)(hibits & 1)) {
/* hibits += (int)copysign(1.0, sum.hi * sum.lo) */
EXTRACT_WORD64(lobits, sum.lo);
hibits += 1 - (((hibits ^ lobits) >> 62) & 2);
INSERT_WORD64(sum.hi, hibits);
}
}
return (ldexp(sum.hi, scale));
}
/*
* Compute a*b exactly, returning the exact result in a struct dd. We assume
* that both a and b are normalized, so no underflow or overflow will occur.
* The current rounding mode must be round-to-nearest.
*/
static inline struct dd
dd_mul(double a, double b)
{
static const double split = 0x1p27 + 1.0;
struct dd ret;
double ha, hb, la, lb, p, q;
p = a * split;
ha = a - p;
ha += p;
la = a - ha;
p = b * split;
hb = b - p;
hb += p;
lb = b - hb;
p = ha * hb;
q = ha * lb + la * hb;
ret.hi = p + q;
ret.lo = p - ret.hi + q + la * lb;
return (ret);
}
/*
* Fused multiply-add: Compute x * y + z with a single rounding error.
*
* We use scaling to avoid overflow/underflow, along with the
* canonical precision-doubling technique adapted from:
*
* Dekker, T. A Floating-Point Technique for Extending the
* Available Precision. Numer. Math. 18, 224-242 (1971).
*
* This algorithm is sensitive to the rounding precision. FPUs such
* as the i387 must be set in double-precision mode if variables are
* to be stored in FP registers in order to avoid incorrect results.
* This is the default on FreeBSD, but not on many other systems.
*
* Hardware instructions should be used on architectures that support it,
* since this implementation will likely be several times slower.
*/
double
fma(double x, double y, double z)
{
double xs, ys, zs, adj;
struct dd xy, r;
int oround;
int ex, ey, ez;
int spread;
/*
* Handle special cases. The order of operations and the particular
* return values here are crucial in handling special cases involving
* infinities, NaNs, overflows, and signed zeroes correctly.
*/
if (x == 0.0 || y == 0.0)
return (x * y + z);
if (z == 0.0)
return (x * y);
if (!isfinite(x) || !isfinite(y))
return (x * y + z);
if (!isfinite(z))
return (z);
xs = frexp(x, &ex);
ys = frexp(y, &ey);
zs = frexp(z, &ez);
oround = fegetround();
spread = ex + ey - ez;
/*
* If x * y and z are many orders of magnitude apart, the scaling
* will overflow, so we handle these cases specially. Rounding
* modes other than FE_TONEAREST are painful.
*/
if (spread < -DBL_MANT_DIG) {
feraiseexcept(FE_INEXACT);
if (!isnormal(z))
feraiseexcept(FE_UNDERFLOW);
switch (oround) {
case FE_TONEAREST:
return (z);
case FE_TOWARDZERO:
if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0))
return (z);
else
return (nextafter(z, 0));
case FE_DOWNWARD:
if ((x > 0.0) ^ (y < 0.0))
return (z);
else
return (nextafter(z, -INFINITY));
default: /* FE_UPWARD */
if ((x > 0.0) ^ (y < 0.0))
return (nextafter(z, INFINITY));
else
return (z);
}
}
if (spread <= DBL_MANT_DIG * 2)
zs = ldexp(zs, -spread);
else
zs = copysign(DBL_MIN, zs);
fesetround(FE_TONEAREST);
/*
* Basic approach for round-to-nearest:
*
* (xy.hi, xy.lo) = x * y (exact)
* (r.hi, r.lo) = xy.hi + z (exact)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
xy = dd_mul(xs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
if (r.hi == 0.0) {
/*
* When the addends cancel to 0, ensure that the result has
* the correct sign.
*/
fesetround(oround);
{
volatile double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexp(xy.lo, spread));
}
}
if (oround != FE_TONEAREST) {
/*
* There is no need to worry about double rounding in directed
* rounding modes.
*/
fesetround(oround);
adj = r.lo + xy.lo;
return (ldexp(r.hi + adj, spread));
}
adj = add_adjusted(r.lo, xy.lo);
if (spread + ilogb(r.hi) > -1023)
return (ldexp(r.hi + adj, spread));
else
return (add_and_denormalize(r.hi, adj, spread));
}

75
lib/libm/src/s_fmaf.c Normal file
View File

@@ -0,0 +1,75 @@
/* $NetBSD: s_fmaf.c,v 1.2 2013/02/11 01:29:58 christos Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmaf.c,v 1.3 2011/10/15 04:16:58 das Exp $");
#else
__RCSID("$NetBSD: s_fmaf.c,v 1.2 2013/02/11 01:29:58 christos Exp $");
#endif
#include <fenv.h>
#include "math.h"
#include "math_private.h"
/*
* Fused multiply-add: Compute x * y + z with a single rounding error.
*
* A double has more than twice as much precision than a float, so
* direct double-precision arithmetic suffices, except where double
* rounding occurs.
*/
float
fmaf(float x, float y, float z)
{
double xy, result;
uint32_t hr, lr;
xy = (double)x * y;
result = xy + z;
EXTRACT_WORDS(hr, lr, result);
/* Common case: The double precision result is fine. */
if ((lr & 0x1fffffff) != 0x10000000 || /* not a halfway case */
(hr & 0x7ff00000) == 0x7ff00000 || /* NaN */
result - xy == z || /* exact */
fegetround() != FE_TONEAREST) /* not round-to-nearest */
return (result);
/*
* If result is inexact, and exactly halfway between two float values,
* we need to adjust the low-order bit in the direction of the error.
*/
fesetround(FE_TOWARDZERO);
volatile double vxy = xy; /* XXX work around gcc CSE bug */
double adjusted_result = vxy + z;
fesetround(FE_TONEAREST);
if (result == adjusted_result)
SET_LOW_WORD(adjusted_result, lr + 1);
return (adjusted_result);
}

279
lib/libm/src/s_fmal.c Normal file
View File

@@ -0,0 +1,279 @@
/* $NetBSD: s_fmal.c,v 1.3 2013/02/12 21:40:19 martin Exp $ */
/*-
* Copyright (c) 2005-2011 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmal.c,v 1.7 2011/10/21 06:30:43 das Exp $");
#else
__RCSID("$NetBSD: s_fmal.c,v 1.3 2013/02/12 21:40:19 martin Exp $");
#endif
#include <machine/ieee.h>
#include <fenv.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
/*
* A struct dd represents a floating-point number with twice the precision
* of a long double. We maintain the invariant that "hi" stores the high-order
* bits of the result.
*/
struct dd {
long double hi;
long double lo;
};
/*
* Compute a+b exactly, returning the exact result in a struct dd. We assume
* that both a and b are finite, but make no assumptions about their relative
* magnitudes.
*/
static inline struct dd
dd_add(long double a, long double b)
{
struct dd ret;
long double s;
ret.hi = a + b;
s = ret.hi - a;
ret.lo = (a - (ret.hi - s)) + (b - s);
return (ret);
}
/*
* Compute a+b, with a small tweak: The least significant bit of the
* result is adjusted into a sticky bit summarizing all the bits that
* were lost to rounding. This adjustment negates the effects of double
* rounding when the result is added to another number with a higher
* exponent. For an explanation of round and sticky bits, see any reference
* on FPU design, e.g.,
*
* J. Coonen. An Implementation Guide to a Proposed Standard for
* Floating-Point Arithmetic. Computer, vol. 13, no. 1, Jan 1980.
*/
static inline long double
add_adjusted(long double a, long double b)
{
struct dd sum;
union ieee_ext_u u;
sum = dd_add(a, b);
if (sum.lo != 0) {
u.extu_ld = sum.hi;
if ((u.extu_ext.ext_fracl & 1) == 0)
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (sum.hi);
}
/*
* Compute ldexp(a+b, scale) with a single rounding error. It is assumed
* that the result will be subnormal, and care is taken to ensure that
* double rounding does not occur.
*/
static inline long double
add_and_denormalize(long double a, long double b, int scale)
{
struct dd sum;
int bits_lost;
union ieee_ext_u u;
sum = dd_add(a, b);
/*
* If we are losing at least two bits of accuracy to denormalization,
* then the first lost bit becomes a round bit, and we adjust the
* lowest bit of sum.hi to make it a sticky bit summarizing all the
* bits in sum.lo. With the sticky bit adjusted, the hardware will
* break any ties in the correct direction.
*
* If we are losing only one bit to denormalization, however, we must
* break the ties manually.
*/
if (sum.lo != 0) {
u.extu_ld = sum.hi;
bits_lost = -u.extu_ext.ext_exp - scale + 1;
if ((bits_lost != 1) ^ (int)(u.extu_ext.ext_fracl & 1))
sum.hi = nextafterl(sum.hi, INFINITY * sum.lo);
}
return (ldexp((double)sum.hi, scale));
}
/*
* Compute a*b exactly, returning the exact result in a struct dd. We assume
* that both a and b are normalized, so no underflow or overflow will occur.
* The current rounding mode must be round-to-nearest.
*/
static inline struct dd
dd_mul(long double a, long double b)
{
#if LDBL_MANT_DIG == 64
static const long double split = 0x1p32L + 1.0;
#elif LDBL_MANT_DIG == 113
static const long double split = 0x1p57L + 1.0;
#endif
struct dd ret;
long double ha, hb, la, lb, p, q;
p = a * split;
ha = a - p;
ha += p;
la = a - ha;
p = b * split;
hb = b - p;
hb += p;
lb = b - hb;
p = ha * hb;
q = ha * lb + la * hb;
ret.hi = p + q;
ret.lo = p - ret.hi + q + la * lb;
return (ret);
}
/*
* Fused multiply-add: Compute x * y + z with a single rounding error.
*
* We use scaling to avoid overflow/underflow, along with the
* canonical precision-doubling technique adapted from:
*
* Dekker, T. A Floating-Point Technique for Extending the
* Available Precision. Numer. Math. 18, 224-242 (1971).
*/
long double
fmal(long double x, long double y, long double z)
{
long double xs, ys, zs, adj;
struct dd xy, r;
int oround;
int ex, ey, ez;
int spread;
/*
* Handle special cases. The order of operations and the particular
* return values here are crucial in handling special cases involving
* infinities, NaNs, overflows, and signed zeroes correctly.
*/
if (x == 0.0 || y == 0.0)
return (x * y + z);
if (z == 0.0)
return (x * y);
if (!isfinite(x) || !isfinite(y))
return (x * y + z);
if (!isfinite(z))
return (z);
xs = frexpl(x, &ex);
ys = frexpl(y, &ey);
zs = frexpl(z, &ez);
oround = fegetround();
spread = ex + ey - ez;
/*
* If x * y and z are many orders of magnitude apart, the scaling
* will overflow, so we handle these cases specially. Rounding
* modes other than FE_TONEAREST are painful.
*/
if (spread < -LDBL_MANT_DIG) {
feraiseexcept(FE_INEXACT);
if (!isnormal(z))
feraiseexcept(FE_UNDERFLOW);
switch (oround) {
case FE_TONEAREST:
return (z);
case FE_TOWARDZERO:
if ((x > 0.0) ^ (y < 0.0) ^ (z < 0.0))
return (z);
else
return (nextafterl(z, 0));
case FE_DOWNWARD:
if ((x > 0.0) ^ (y < 0.0))
return (z);
else
return (nextafterl(z, (long double)-INFINITY));
default: /* FE_UPWARD */
if ((x > 0.0) ^ (y < 0.0))
return (nextafterl(z, (long double)INFINITY));
else
return (z);
}
}
if (spread <= LDBL_MANT_DIG * 2)
zs = ldexpl(zs, -spread);
else
zs = copysignl(LDBL_MIN, zs);
fesetround(FE_TONEAREST);
/*
* Basic approach for round-to-nearest:
*
* (xy.hi, xy.lo) = x * y (exact)
* (r.hi, r.lo) = xy.hi + z (exact)
* adj = xy.lo + r.lo (inexact; low bit is sticky)
* result = r.hi + adj (correctly rounded)
*/
xy = dd_mul(xs, ys);
r = dd_add(xy.hi, zs);
spread = ex + ey;
if (r.hi == 0.0) {
/*
* When the addends cancel to 0, ensure that the result has
* the correct sign.
*/
fesetround(oround);
{
volatile long double vzs = zs; /* XXX gcc CSE bug workaround */
return (xy.hi + vzs + ldexpl(xy.lo, spread));
}
}
if (oround != FE_TONEAREST) {
/*
* There is no need to worry about double rounding in directed
* rounding modes.
*/
fesetround(oround);
adj = r.lo + xy.lo;
return (ldexpl(r.hi + adj, spread));
}
adj = add_adjusted(r.lo, xy.lo);
if (spread + ilogbl(r.hi) > -16383)
return (ldexpl(r.hi + adj, spread));
else
return (add_and_denormalize(r.hi, adj, spread));
}
#endif

6
lib/libm/src/s_fmax.c
View File

@@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_fmax.c,v 1.2 2010/03/08 01:05:20 snj Exp $");
__RCSID("$NetBSD: s_fmax.c,v 1.3 2013/11/29 22:16:10 joerg Exp $");
#ifdef notdef
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmax.c,v 1.1 2004/06/30 07:04:01 das Exp $");
#endif
@@ -34,6 +34,10 @@ __FBSDID("$FreeBSD: src/lib/msun/src/s_fmax.c,v 1.1 2004/06/30 07:04:01 das Exp
#include <machine/ieee.h>
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(fmaxl, fmax)
#endif
double
fmax(double x, double y)
{

6
lib/libm/src/s_fmaxl.c
View File

@@ -25,20 +25,22 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_fmaxl.c,v 1.2 2010/03/08 01:05:20 snj Exp $");
__RCSID("$NetBSD: s_fmaxl.c,v 1.4 2013/11/29 22:16:10 joerg Exp $");
#ifdef notdef
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmaxl.c,v 1.1 2004/06/30 07:04:01 das Exp $");
#endif
#include <string.h>
#include <math.h>
#include <machine/ieee.h>
#ifdef EXT_EXP_INFNAN
#ifdef __HAVE_LONG_DOUBLE
long double
fmaxl(long double x, long double y)
{
union ieee_ext_u u[2];
memset(&u, 0, sizeof u);
u[0].extu_ld = x;
u[0].extu_ext.ext_frach &= ~0x80000000;
u[1].extu_ld = y;

6
lib/libm/src/s_fmin.c
View File

@@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_fmin.c,v 1.1 2009/10/04 22:04:30 christos Exp $");
__RCSID("$NetBSD: s_fmin.c,v 1.2 2013/11/29 22:16:10 joerg Exp $");
#ifdef notdef
__FBSDID("$FreeBSD: src/lib/msun/src/s_fmin.c,v 1.1 2004/06/30 07:04:01 das Exp $");
#endif
@@ -34,6 +34,10 @@ __FBSDID("$FreeBSD: src/lib/msun/src/s_fmin.c,v 1.1 2004/06/30 07:04:01 das Exp
#include <machine/ieee.h>
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(fminl, fmin)
#endif
double
fmin(double x, double y)
{

6
lib/libm/src/s_fminl.c
View File

@@ -25,20 +25,22 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_fminl.c,v 1.2 2010/03/08 01:05:20 snj Exp $");
__RCSID("$NetBSD: s_fminl.c,v 1.4 2013/11/29 22:16:10 joerg Exp $");
#ifdef notdef
__FBSDID("$FreeBSD: src/lib/msun/src/s_fminl.c,v 1.1 2004/06/30 07:04:01 das Exp $");
#endif
#include <math.h>
#include <string.h>
#include <machine/ieee.h>
#ifdef EXT_EXP_INFNAN
#ifdef __HAVE_LONG_DOUBLE
long double
fminl(long double x, long double y)
{
union ieee_ext_u u[2];
memset(&u, 0, sizeof u);
u[0].extu_ld = x;
u[0].extu_ext.ext_frach &= ~0x80000000;
u[1].extu_ld = y;

70
lib/libm/src/s_frexpl.c Normal file
View File

@@ -0,0 +1,70 @@
/* $NetBSD: s_frexpl.c,v 1.3 2013/02/12 21:40:19 martin Exp $ */
/*-
* Copyright (c) 2004-2005 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD: src/lib/msun/src/s_frexpl.c,v 1.1 2005/03/07 04:54:51 das Exp $
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_frexpl.c,v 1.3 2013/02/12 21:40:19 martin Exp $");
#include <machine/ieee.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
#if LDBL_MAX_EXP != 0x4000
#error "Unsupported long double format"
#endif
long double
frexpl(long double x, int *ex)
{
union ieee_ext_u u;
u.extu_ld = x;
switch (u.extu_ext.ext_exp) {
case 0: /* 0 or subnormal */
if ((u.extu_ext.ext_fracl | u.extu_ext.ext_frach) == 0) {
*ex = 0;
} else {
u.extu_ld *= 0x1.0p514;
*ex = u.extu_ext.ext_exp - 0x4200;
u.extu_ext.ext_exp = 0x3ffe;
}
break;
case 0x7fff: /* infinity or NaN; value of *ex is unspecified */
break;
default: /* normal */
*ex = u.extu_ext.ext_exp - 0x3ffe;
u.extu_ext.ext_exp = 0x3ffe;
break;
}
return (u.extu_ld);
}
#endif

10
lib/libm/src/s_ilogb.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_ilogb.c,v 1.12 2002/05/26 22:01:56 wiz Exp $");
__RCSID("$NetBSD: s_ilogb.c,v 1.14 2013/02/09 22:56:00 matt Exp $");
#endif
/* ilogb(double x)
@@ -24,6 +24,10 @@ __RCSID("$NetBSD: s_ilogb.c,v 1.12 2002/05/26 22:01:56 wiz Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(ilogbl,ilogb)
#endif
int
ilogb(double x)
{
@@ -34,7 +38,7 @@ ilogb(double x)
if(hx<0x00100000) {
GET_LOW_WORD(lx,x);
if((hx|lx)==0)
return 0x80000001; /* ilogb(0) = 0x80000001 */
return FP_ILOGB0; /* ilogb(0) = 0x80000001 */
else /* subnormal x */
if(hx==0) {
for (ix = -1043; lx>0; lx<<=1) ix -=1;
@@ -44,5 +48,5 @@ ilogb(double x)
return ix;
}
else if (hx<0x7ff00000) return (hx>>20)-1023;
else return 0x7fffffff;
else return FP_ILOGBNAN; /* inf too */
}

6
lib/libm/src/s_ilogbf.c
View File

@@ -15,7 +15,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_ilogbf.c,v 1.7 2002/05/26 22:01:56 wiz Exp $");
__RCSID("$NetBSD: s_ilogbf.c,v 1.8 2013/02/09 22:56:00 matt Exp $");
#endif
#include "math.h"
@@ -30,11 +30,11 @@ ilogbf(float x)
hx &= 0x7fffffff;
if(hx<0x00800000) {
if(hx==0)
return 0x80000001; /* ilogb(0) = 0x80000001 */
return FP_ILOGB0; /* ilogb(0) = 0x80000001 */
else /* subnormal x */
for (ix = -126,hx<<=8; hx>0; hx<<=1) ix -=1;
return ix;
}
else if (hx<0x7f800000) return (hx>>23)-127;
else return 0x7fffffff;
else return FP_ILOGBNAN; /* inf too */
}

76
lib/libm/src/s_ilogbl.c Normal file
View File

@@ -0,0 +1,76 @@
/* $NetBSD: s_ilogbl.c,v 1.2 2013/02/09 22:56:00 matt Exp $ */
/*-
* Copyright (c) 2011 The NetBSD Foundation, Inc.
* All rights reserved.
*
* This code is derived from software contributed to The NetBSD Foundation
* by Joerg Sonnenberger.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_ilogbl.c,v 1.2 2013/02/09 22:56:00 matt Exp $");
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#if LDBL_MANT_DIG == 64
#define FROM_UNDERFLOW 0x1p65L
#elif LDBL_MANT_DIG == 113
#define FROM_UNDERFLOW 0x1p114L
#else
#error Unsupported long double format
#endif
int
ilogbl(long double x)
{
union ieee_ext_u u;
if (x == 0.0L)
return FP_ILOGB0; /* ilogbl(0) = 0x80000001 */
u.extu_ld = x;
if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
return FP_ILOGBNAN; /* inf too */
if (u.extu_ext.ext_exp == 0) {
/*
* Scale denormalized numbers slightly,
* so that they are normal.
*/
u.extu_ld *= FROM_UNDERFLOW;
return u.extu_ext.ext_exp - EXT_EXP_BIAS - LDBL_MANT_DIG - 1;
}
return u.extu_ext.ext_exp - EXT_EXP_BIAS;
}
#endif

30
lib/libm/src/s_ldexp.c
View File

@@ -1,30 +0,0 @@
/* @(#)s_ldexp0.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_ldexp.c,v 1.11 2010/04/23 19:17:07 drochner Exp $");
#endif
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#include <errno.h>
double
ldexp(double value, int exp0)
{
if(!finite(value)||value==0.0) return value;
value = scalbn(value,exp0);
if(!finite(value)||value==0.0) errno = ERANGE;
return value;
}

33
lib/libm/src/s_ldexpf.c
View File

@@ -1,33 +0,0 @@
/* s_ldexp0f.c -- float version of s_ldexp0.c.
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
*/
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_ldexpf.c,v 1.8 2010/04/23 19:17:07 drochner Exp $");
#endif
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#include <errno.h>
float
ldexpf(float value, int exp0)
{
if(!finitef(value)||value==(float)0.0) return value;
value = scalbnf(value,exp0);
if(!finitef(value)||value==(float)0.0) errno = ERANGE;
return value;
}

6
lib/libm/src/s_logb.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_logb.c,v 1.11 2002/05/26 22:01:57 wiz Exp $");
__RCSID("$NetBSD: s_logb.c,v 1.12 2011/08/03 14:13:07 joerg Exp $");
#endif
/*
@@ -24,6 +24,10 @@ __RCSID("$NetBSD: s_logb.c,v 1.11 2002/05/26 22:01:57 wiz Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(logbl,logb)
#endif
double
logb(double x)
{

76
lib/libm/src/s_logbl.c Normal file
View File

@@ -0,0 +1,76 @@
/* $NetBSD: s_logbl.c,v 1.1 2011/08/03 14:13:07 joerg Exp $ */
/*-
* Copyright (c) 2011 The NetBSD Foundation, Inc.
* All rights reserved.
*
* This code is derived from software contributed to The NetBSD Foundation
* by Joerg Sonnenberger.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_logbl.c,v 1.1 2011/08/03 14:13:07 joerg Exp $");
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#if LDBL_MANT_DIG == 64
#define FROM_UNDERFLOW 0x1p65L
#elif LDBL_MANT_DIG == 113
#define FROM_UNDERFLOW 0x1p114L
#else
#error Unsupported long double format
#endif
long double
logbl(long double x)
{
union ieee_ext_u u;
if (x == 0.0L)
return -1.0L / fabsl(x); /* -HUGE_VALL + exception */
u.extu_ld = x;
if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
return fabsl(x); /* NaN or +Inf */
if (u.extu_ext.ext_exp == 0) {
/*
* Scale denormalized numbers slightly,
* so that they are normal.
*/
u.extu_ld *= FROM_UNDERFLOW;
return u.extu_ext.ext_exp - EXT_EXP_BIAS - LDBL_MANT_DIG - 1;
}
return u.extu_ext.ext_exp - EXT_EXP_BIAS;
}
#endif

124
lib/libm/src/s_nan.c Normal file
View File

@@ -0,0 +1,124 @@
/* $NetBSD: s_nan.c,v 1.2 2013/02/09 20:19:13 christos Exp $ */
/*-
* Copyright (c) 2007 David Schultz
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD: src/lib/msun/src/s_nan.c,v 1.2 2007/12/18 23:46:32 das Exp $
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_nan.c,v 1.2 2013/02/09 20:19:13 christos Exp $");
#include <sys/endian.h>
#include <ctype.h>
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <string.h>
#include "math_private.h"
/* XXX: Locale? not here? in <ctype.h>? */
static int
digittoint(char s) {
if (isdigit((unsigned char)s))
return s - '0';
if (islower((unsigned char)s))
return s - 'a' + 10;
return s - 'A' + 10;
}
/*
* Scan a string of hexadecimal digits (the format nan(3) expects) and
* make a bit array (using the local endianness). We stop when we
* encounter an invalid character, NUL, etc. If we overflow, we do
* the same as gcc's __builtin_nan(), namely, discard the high order bits.
*
* The format this routine accepts needs to be compatible with what is used
* in contrib/gdtoa/hexnan.c (for strtod/scanf) and what is used in
* __builtin_nan(). In fact, we're only 100% compatible for strings we
* consider valid, so we might be violating the C standard. But it's
* impossible to use nan(3) portably anyway, so this seems good enough.
*/
static void
_scan_nan(uint32_t *words, int num_words, const char *s)
{
int si; /* index into s */
int bitpos; /* index into words (in bits) */
memset(words, 0, num_words * sizeof(*words));
/* Allow a leading '0x'. (It's expected, but redundant.) */
if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
s += 2;
/* Scan forwards in the string, looking for the end of the sequence. */
for (si = 0; isxdigit((unsigned char)s[si]); si++)
continue;
/* Scan backwards, filling in the bits in words[] as we go. */
#if _BYTE_ORDER == _LITTLE_ENDIAN
for (bitpos = 0; bitpos < 32 * num_words; bitpos += 4) {
#else
for (bitpos = 32 * num_words - 4; bitpos >= 0; bitpos -= 4) {
#endif
if (--si < 0)
break;
words[bitpos / 32] |= digittoint(s[si]) << (bitpos % 32);
}
}
double
nan(const char *s)
{
union {
double d;
uint32_t bits[2];
} u;
_scan_nan(u.bits, 2, s);
#if _BYTE_ORDER == _LITTLE_ENDIAN
u.bits[1] |= 0x7ff80000;
#else
u.bits[0] |= 0x7ff80000;
#endif
return (u.d);
}
float
nanf(const char *s)
{
union {
float f;
uint32_t bits[1];
} u;
_scan_nan(u.bits, 1, s);
u.bits[0] |= 0x7fc00000;
return (u.f);
}
#if (LDBL_MANT_DIG == 53)
__weak_reference(nan, nanl);
#endif

64
lib/libm/src/s_nearbyint.c Normal file
View File

@@ -0,0 +1,64 @@
/* $NetBSD: s_nearbyint.c,v 1.3 2013/02/11 02:45:27 christos Exp $ */
/*-
* Copyright (c) 2004 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_nearbyint.c,v 1.2 2008/01/14 02:12:06 das Exp $");
#else
__RCSID("$NetBSD: s_nearbyint.c,v 1.3 2013/02/11 02:45:27 christos Exp $");
#endif
#include <machine/ieee.h>
#include <fenv.h>
#include <math.h>
/*
* We save and restore the floating-point environment to avoid raising
* an inexact exception. We can get away with using fesetenv()
* instead of feclearexcept()/feupdateenv() to restore the environment
* because the only exception defined for rint() is overflow, and
* rounding can't overflow as long as emax >= p.
*/
#define DECL(type, fn, rint) \
type \
fn(type x) \
{ \
type ret; \
fenv_t env; \
\
fegetenv(&env); \
ret = rint(x); \
fesetenv(&env); \
return (ret); \
}
DECL(double, nearbyint, rint)
DECL(float, nearbyintf, rintf)
#ifdef __HAVE_LONG_DOUBLE
DECL(long double, nearbyintl, rintl)
#endif

8
lib/libm/src/s_nextafter.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_nextafter.c,v 1.11 2002/05/26 22:01:57 wiz Exp $");
__RCSID("$NetBSD: s_nextafter.c,v 1.13 2013/07/18 22:32:53 matt Exp $");
#endif
/* IEEE functions
@@ -25,6 +25,10 @@ __RCSID("$NetBSD: s_nextafter.c,v 1.11 2002/05/26 22:01:57 wiz Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(nextafterl, nextafter)
#endif
double
nextafter(double x, double y)
{
@@ -39,7 +43,7 @@ nextafter(double x, double y)
if(((ix>=0x7ff00000)&&((ix-0x7ff00000)|lx)!=0) || /* x is nan */
((iy>=0x7ff00000)&&((iy-0x7ff00000)|ly)!=0)) /* y is nan */
return x+y;
if(x==y) return x; /* x=y, return x */
if(x==y) return y; /* x=y, return y */
if((ix|lx)==0) { /* x == 0 */
INSERT_WORDS(x,hy&0x80000000,1); /* return +-minsubnormal */
y = x*x;

4
lib/libm/src/s_nextafterf.c
View File

@@ -15,7 +15,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_nextafterf.c,v 1.7 2002/05/26 22:01:58 wiz Exp $");
__RCSID("$NetBSD: s_nextafterf.c,v 1.8 2011/04/18 15:59:09 drochner Exp $");
#endif
#include "math.h"
@@ -34,7 +34,7 @@ nextafterf(float x, float y)
if((ix>0x7f800000) || /* x is nan */
(iy>0x7f800000)) /* y is nan */
return x+y;
if(x==y) return x; /* x=y, return x */
if(x==y) return y; /* x=y, return y */
if(ix==0) { /* x == 0 */
SET_FLOAT_WORD(x,(hy&0x80000000)|1);/* return +-minsubnormal */
y = x*x;

14
lib/libm/src/s_nextafterl.c
View File

@@ -1,4 +1,4 @@
/* $NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $ */
/* $NetBSD: s_nextafterl.c,v 1.4 2013/07/18 22:31:13 matt Exp $ */
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
@@ -13,17 +13,23 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_nextafterl.c,v 1.2 2010/09/17 20:39:39 christos Exp $");
__RCSID("$NetBSD: s_nextafterl.c,v 1.4 2013/07/18 22:31:13 matt Exp $");
#include <float.h>
#include <math.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef EXT_EXP_INFNAN
#if LDBL_MAX_EXP != 0x4000
#error "Unsupported long double format"
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#endif
/*
* IEEE functions
* nextafterl(x,y)
@@ -83,7 +89,9 @@ nextafterl(long double x, long double y)
return x+x; /* overflow */
if (ux.extu_exp == 0) { /* underflow */
#ifndef LDBL_IMPLICIT_NBIT
mask_nbit_l(ux);
#endif
t = ux.extu_ld * ux.extu_ld;
if (t != ux.extu_ld) /* raise underflow flag */
return ux.extu_ld;
@@ -92,3 +100,5 @@ nextafterl(long double x, long double y)
return ux.extu_ld;
}
#endif
#endif /* __HAVE_LONG_DOUBLE */

8
lib/libm/src/s_nexttoward.c
View File

@@ -1,4 +1,4 @@
/* $NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos Exp $ */
/* $NetBSD: s_nexttoward.c,v 1.2 2013/08/21 13:03:56 martin Exp $ */
/* @(#)s_nextafter.c 5.1 93/09/24 */
/*
@@ -13,7 +13,7 @@
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos Exp $");
__RCSID("$NetBSD: s_nexttoward.c,v 1.2 2013/08/21 13:03:56 martin Exp $");
/*
* We assume that a long double has a 15-bit exponent. On systems
@@ -30,6 +30,10 @@ __RCSID("$NetBSD: s_nexttoward.c,v 1.1 2010/09/15 16:12:05 christos Exp $");
#error "Unsupported long double format"
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define LDBL_NBIT 0
#endif
/*
* The nexttoward() function is equivalent to nextafter() function,
* except that the second parameter shall have type long double and

69
lib/libm/src/s_nexttowardf.c Normal file
View File

@@ -0,0 +1,69 @@
/* $NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $ */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_nexttowardf.c,v 1.3 2011/02/10 07:38:38 das Exp $");
#else
__RCSID("$NetBSD: s_nexttowardf.c,v 1.3 2013/02/09 23:14:44 christos Exp $");
#endif
#include <string.h>
#include <float.h>
#include <machine/ieee.h>
#include "math.h"
#include "math_private.h"
#ifdef EXT_EXP_INFNAN
float
nexttowardf(float x, long double y)
{
volatile float t;
int32_t hx,ix;
union ieee_ext_u uy;
GET_FLOAT_WORD(hx,x);
ix = hx&0x7fffffff; /* |x| */
memset(&uy, 0, sizeof(uy));
uy.extu_ld = y;
uy.extu_ext.ext_frach &= ~0x80000000;
if((ix>0x7f800000) ||
(uy.extu_ext.ext_exp == EXT_EXP_INFNAN &&
(uy.extu_ext.ext_frach | uy.extu_ext.ext_fracl) != 0))
return x+y; /* x or y is nan */
if(x==y) return (float)y; /* x=y, return y */
if(ix==0) { /* x == 0 */
SET_FLOAT_WORD(x,(uy.extu_ext.ext_sign<<31)|1);/* return +-minsubnormal */
t = x*x;
if(t==x) return t; else return x; /* raise underflow flag */
}
if((hx >= 0) ^ (x < y)) /* x -= ulp */
hx -= 1;
else /* x += ulp */
hx += 1;
ix = hx&0x7f800000;
if(ix>=0x7f800000) return x+x; /* overflow */
if(ix<0x00800000) { /* underflow */
t = x*x;
if(t!=x) { /* raise underflow flag */
SET_FLOAT_WORD(x,hx);
return x;
}
}
SET_FLOAT_WORD(x,hx);
return x;
}
#endif

153
lib/libm/src/s_remquo.c Normal file
View File

@@ -0,0 +1,153 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
#include <float.h>
#include "math.h"
#include "math_private.h"
static const double Zero[] = {0.0, -0.0,};
/*
* Return the IEEE remainder and set *quo to the last n bits of the
* quotient, rounded to the nearest integer. We choose n=31 because
* we wind up computing all the integer bits of the quotient anyway as
* a side-effect of computing the remainder by the shift and subtract
* method. In practice, this is far more bits than are needed to use
* remquo in reduction algorithms.
*/
double
remquo(double x, double y, int *quo)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;
u_int32_t lx,ly,lz,q,sxy;
EXTRACT_WORDS(hx,lx,x);
EXTRACT_WORDS(hy,ly,y);
sxy = (hx ^ hy) & 0x80000000;
sx = hx&0x80000000; /* sign of x */
hx ^=sx; /* |x| */
hy &= 0x7fffffff; /* |y| */
/* purge off exception values */
if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */
((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */
return (x*y)/(x*y);
if(hx<=hy) {
if((hx<hy)||(lx<ly)) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */
}
if(lx==ly) {
*quo = 1;
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
}
}
/* determine ix = ilogb(x) */
if(hx<0x00100000) { /* subnormal x */
if(hx==0) {
for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
} else {
for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
}
} else ix = (hx>>20)-1023;
/* determine iy = ilogb(y) */
if(hy<0x00100000) { /* subnormal y */
if(hy==0) {
for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
} else {
for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
}
} else iy = (hy>>20)-1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -1022)
hx = 0x00100000|(0x000fffff&hx);
else { /* subnormal x, shift x to normal */
n = -1022-ix;
if(n<=31) {
hx = (hx<<n)|(lx>>(32-n));
lx <<= n;
} else {
hx = lx<<(n-32);
lx = 0;
}
}
if(iy >= -1022)
hy = 0x00100000|(0x000fffff&hy);
else { /* subnormal y, shift y to normal */
n = -1022-iy;
if(n<=31) {
hy = (hy<<n)|(ly>>(32-n));
ly <<= n;
} else {
hy = ly<<(n-32);
ly = 0;
}
}
/* fix point fmod */
n = ix - iy;
q = 0;
while(n--) {
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
else {hx = hz+hz+(lz>>31); lx = lz+lz; q++;}
q <<= 1;
}
hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
if(hz>=0) {hx=hz;lx=lz;q++;}
/* convert back to floating value and restore the sign */
if((hx|lx)==0) { /* return sign(x)*0 */
*quo = (sxy ? -q : q);
return Zero[(u_int32_t)sx>>31];
}
while(hx<0x00100000) { /* normalize x */
hx = hx+hx+(lx>>31); lx = lx+lx;
iy -= 1;
}
if(iy>= -1022) { /* normalize output */
hx = ((hx-0x00100000)|((iy+1023)<<20));
} else { /* subnormal output */
n = -1022 - iy;
if(n<=20) {
lx = (lx>>n)|((u_int32_t)hx<<(32-n));
hx >>= n;
} else if (n<=31) {
lx = (hx<<(32-n))|(lx>>n); hx = sx;
} else {
lx = hx>>(n-32); hx = sx;
}
}
fixup:
INSERT_WORDS(x,hx,lx);
y = fabs(y);
if (y < 0x1p-1021) {
if (x+x>y || (x+x==y && (q & 1))) {
q++;
x-=y;
}
} else if (x>0.5*y || (x==0.5*y && (q & 1))) {
q++;
x-=y;
}
GET_HIGH_WORD(hx,x);
SET_HIGH_WORD(x,hx^sx);
q &= 0x7fffffff;
*quo = (sxy ? -q : q);
return x;
}

120
lib/libm/src/s_remquof.c Normal file
View File

@@ -0,0 +1,120 @@
/* @(#)e_fmod.c 1.3 95/01/18 */
/*-
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
#include "math.h"
#include "math_private.h"
static const float Zero[] = {0.0, -0.0,};
/*
* Return the IEEE remainder and set *quo to the last n bits of the
* quotient, rounded to the nearest integer. We choose n=31 because
* we wind up computing all the integer bits of the quotient anyway as
* a side-effect of computing the remainder by the shift and subtract
* method. In practice, this is far more bits than are needed to use
* remquo in reduction algorithms.
*/
float
remquof(float x, float y, int *quo)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;
u_int32_t q,sxy;
GET_FLOAT_WORD(hx,x);
GET_FLOAT_WORD(hy,y);
sxy = (hx ^ hy) & 0x80000000;
sx = hx&0x80000000; /* sign of x */
hx ^=sx; /* |x| */
hy &= 0x7fffffff; /* |y| */
/* purge off exception values */
if(hy==0||hx>=0x7f800000||hy>0x7f800000) /* y=0,NaN;or x not finite */
return (x*y)/(x*y);
if(hx<hy) {
q = 0;
goto fixup; /* |x|<|y| return x or x-y */
} else if(hx==hy) {
*quo = 1;
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
}
/* determine ix = ilogb(x) */
if(hx<0x00800000) { /* subnormal x */
for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
} else ix = (hx>>23)-127;
/* determine iy = ilogb(y) */
if(hy<0x00800000) { /* subnormal y */
for (iy = -126,i=(hy<<8); i>0; i<<=1) iy -=1;
} else iy = (hy>>23)-127;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -126)
hx = 0x00800000|(0x007fffff&hx);
else { /* subnormal x, shift x to normal */
n = -126-ix;
hx <<= n;
}
if(iy >= -126)
hy = 0x00800000|(0x007fffff&hy);
else { /* subnormal y, shift y to normal */
n = -126-iy;
hy <<= n;
}
/* fix point fmod */
n = ix - iy;
q = 0;
while(n--) {
hz=hx-hy;
if(hz<0) hx = hx << 1;
else {hx = hz << 1; q++;}
q <<= 1;
}
hz=hx-hy;
if(hz>=0) {hx=hz;q++;}
/* convert back to floating value and restore the sign */
if(hx==0) { /* return sign(x)*0 */
*quo = (sxy ? -q : q);
return Zero[(u_int32_t)sx>>31];
}
while(hx<0x00800000) { /* normalize x */
hx <<= 1;
iy -= 1;
}
if(iy>= -126) { /* normalize output */
hx = ((hx-0x00800000)|((iy+127)<<23));
} else { /* subnormal output */
n = -126 - iy;
hx >>= n;
}
fixup:
SET_FLOAT_WORD(x,hx);
y = fabsf(y);
if (y < 0x1p-125f) {
if (x+x>y || (x+x==y && (q & 1))) {
q++;
x-=y;
}
} else if (x>0.5f*y || (x==0.5f*y && (q & 1))) {
q++;
x-=y;
}
GET_FLOAT_WORD(hx,x);
SET_FLOAT_WORD(x,hx^sx);
q &= 0x7fffffff;
*quo = (sxy ? -q : q);
return x;
}

95
lib/libm/src/s_rintl.c Normal file
View File

@@ -0,0 +1,95 @@
/* $NetBSD: s_rintl.c,v 1.5 2013/08/21 13:04:44 martin Exp $ */
/*-
* Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_rintl.c,v 1.5 2008/02/22 11:59:05 bde Exp $");
#else
__RCSID("$NetBSD: s_rintl.c,v 1.5 2013/08/21 13:04:44 martin Exp $");
#endif
#include <float.h>
#include <machine/ieee.h>
#include "math.h"
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
static const float
shift[2] = {
#if EXT_FRACBITS == 64
0x1.0p63, -0x1.0p63
#elif EXT_FRACBITS == 113
0x1.0p112, -0x1.0p112
#elif EXT_FRACBITS == 112
0x1.0p111, -0x1.0p111
#else
#error "Unsupported long double format"
#endif
};
static const float zero[2] = { 0.0, -0.0 };
long double
rintl(long double x)
{
union ieee_ext_u u;
uint32_t expsign;
int ex, sign;
u.extu_ld = x;
u.extu_ext.ext_frach &= ~0x80000000;
expsign = u.extu_ext.ext_sign;
ex = expsign & 0x7fff;
if (ex >= EXT_EXP_BIAS + EXT_FRACBITS - 1) {
if (ex == EXT_EXP_BIAS + EXT_FRACBITS)
return (x + x); /* Inf, NaN, or unsupported format */
return (x); /* finite and already an integer */
}
sign = expsign >> 15;
/*
* The following code assumes that intermediate results are
* evaluated in long double precision. If they are evaluated in
* greater precision, double rounding may occur, and if they are
* evaluated in less precision (as on i386), results will be
* wildly incorrect.
*/
x += shift[sign];
x -= shift[sign];
/*
* If the result is +-0, then it must have the same sign as x, but
* the above calculation doesn't always give this. Fix up the sign.
*/
if (ex < EXT_EXP_BIAS && x == 0.0L)
return (zero[sign]);
return (x);
}
#endif

7
lib/libm/src/s_round.c
View File

@@ -26,7 +26,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_round.c,v 1.2 2007/08/21 20:10:27 drochner Exp $");
__RCSID("$NetBSD: s_round.c,v 1.3 2013/11/11 23:57:34 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: src/lib/msun/src/s_round.c,v 1.1 2004/06/07 08:05:36 das Exp $");
#endif
@@ -34,6 +34,11 @@ __FBSDID("$FreeBSD: src/lib/msun/src/s_round.c,v 1.1 2004/06/07 08:05:36 das Exp
#include <math.h>
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_roundl, round)
__weak_alias(roundl, round)
#endif
double
round(double x)
{

63
lib/libm/src/s_roundl.c Normal file
View File

@@ -0,0 +1,63 @@
/*-
* Copyright (c) 2003, Steven G. Kargl
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice unmodified, this list of conditions, and the following
* disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
* IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
* OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
* INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
* NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
* DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
* THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
* THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_roundl.c,v 1.1 2013/11/11 23:57:34 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/s_roundl.c 153017 2005-12-02 13:45:06Z bde $");
#endif
#include "namespace.h"
#include <math.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef __weak_alias
__weak_alias(roundl, _roundl)
#endif
long double
roundl(long double x)
{
long double t;
if (!isfinite(x))
return (x);
if (x >= 0.0) {
t = floorl(x);
if (t - x <= -0.5)
t += 1.0;
return (t);
} else {
t = floorl(-x);
if (t + x <= -0.5)
t += 1.0;
return (-t);
}
}
#endif /* __HAVE_LONG_DOUBLE */

27
lib/libm/src/s_scalbn.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_scalbn.c,v 1.14 2010/04/23 19:17:07 drochner Exp $");
__RCSID("$NetBSD: s_scalbn.c,v 1.18 2013/05/20 19:40:09 joerg Exp $");
#endif
/*
@@ -26,18 +26,39 @@ __RCSID("$NetBSD: s_scalbn.c,v 1.14 2010/04/23 19:17:07 drochner Exp $");
#include "math.h"
#include "math_private.h"
#ifndef _LP64
__strong_alias(_scalbn, _scalbln)
#endif
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_scalbnl, _scalbn)
__strong_alias(_scalblnl, _scalbln)
__weak_alias(scalbnl, _scalbnl)
__weak_alias(scalblnl, _scalblnl)
#endif
#ifdef __weak_alias
__weak_alias(scalbn, _scalbn)
__weak_alias(scalbln, _scalbln)
__weak_alias(ldexp, _scalbn)
#endif
static const double
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
two54 = 0x1.0p54, /* 0x43500000, 0x00000000 */
twom54 = 0x1.0p-54, /* 0x3C900000, 0x00000000 */
huge = 1.0e+300,
tiny = 1.0e-300;
#ifdef _LP64
double
scalbn(double x, int n)
{
return scalbln(x, n);
}
#endif
double
scalbln(double x, long n)
{
int32_t k,hx,lx;
EXTRACT_WORDS(hx,lx,x);

21
lib/libm/src/s_scalbnf.c
View File

@@ -15,25 +15,36 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_scalbnf.c,v 1.9 2010/04/23 19:17:07 drochner Exp $");
__RCSID("$NetBSD: s_scalbnf.c,v 1.12 2013/05/20 19:40:09 joerg Exp $");
#endif
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#ifdef __weak_alias
__weak_alias(scalbnf, _scalbnf)
#ifndef _LP64
__strong_alias(_scalbnf, _scalblnf)
#endif
__weak_alias(scalbnf, _scalbnf)
__weak_alias(scalblnf, _scalblnf)
__weak_alias(ldexpf, _scalbnf)
static const float
two25 = 3.355443200e+07, /* 0x4c000000 */
twom25 = 2.9802322388e-08, /* 0x33000000 */
two25 = 0x1.0p25, /* 0x4c000000 */
twom25 = 0x1.0p-25, /* 0x33000000 */
huge = 1.0e+30,
tiny = 1.0e-30;
#ifdef _LP64
float
scalbnf(float x, int n)
{
return scalblnf(x, n);
}
#endif
float
scalblnf(float x, long n)
{
int32_t k,ix;
GET_FLOAT_WORD(ix,x);

116
lib/libm/src/s_scalbnl.c Normal file
View File

@@ -0,0 +1,116 @@
/* $NetBSD: s_scalbnl.c,v 1.9 2013/05/20 19:40:09 joerg Exp $ */
/*-
* Copyright (c) 2011 The NetBSD Foundation, Inc.
* All rights reserved.
*
* This code is derived from software contributed to The NetBSD Foundation
* by Joerg Sonnenberger.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
* ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
* TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
* PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
* BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_scalbnl.c,v 1.9 2013/05/20 19:40:09 joerg Exp $");
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef _LP64
long double
scalbnl(long double x, int n)
{
return scalblnl(x, n);
}
#else
__strong_alias(_scalbnl, _scalblnl)
#endif
__weak_alias(scalbnl, _scalbnl)
__weak_alias(scalblnl, _scalblnl)
__weak_alias(ldexpl, _scalbnl)
#if LDBL_MANT_DIG == 64
#define FROM_UNDERFLOW 0x1p65L
#define TO_UNDERFLOW 0x1p-65L
#elif LDBL_MANT_DIG == 113
#define FROM_UNDERFLOW 0x1p114L
#define TO_UNDERFLOW 0x1p-114L
#else
#error Unsupported long double format
#endif
long double
scalblnl(long double x, long n)
{
union ieee_ext_u u;
/* Trivial cases first */
if (n == 0 || x == 0.0L)
return x;
u.extu_ld = x;
/* NaN and infinite don't change either, but trigger exception */
if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
return x + x;
/* Protect against integer overflow in calculation of new exponent */
if (n > LDBL_MAX_EXP - LDBL_MIN_EXP + LDBL_MANT_DIG)
goto overflow;
if (n < LDBL_MIN_EXP - LDBL_MAX_EXP - LDBL_MANT_DIG)
goto underflow;
/* Scale denormalized numbers slightly, so that they are normal */
if (u.extu_ext.ext_exp == 0) {
u.extu_ld *= FROM_UNDERFLOW;
n -= LDBL_MANT_DIG + 1;
}
n += u.extu_ext.ext_exp;
if (n >= LDBL_MAX_EXP + EXT_EXP_BIAS)
goto overflow;
/* Positive exponent (incl. bias) means normal result */
if (n > 0) {
u.extu_ext.ext_exp = n;
return u.extu_ld;
}
/* Shift the exponent and let the multiply below handle subnormal */
n += LDBL_MANT_DIG + 1;
if (n <= 0)
goto underflow;
u.extu_ext.ext_exp = n;
return u.extu_ld * TO_UNDERFLOW;
underflow:
return LDBL_MIN * copysignl(LDBL_MIN, x);
overflow:
return LDBL_MAX * copysignl(LDBL_MAX, x);
}
#endif

47
lib/libm/src/s_tgammaf.c Normal file
View File

@@ -0,0 +1,47 @@
/*-
* Copyright (c) 2008 David Schultz <das@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
#if 0
__FBSDID("$FreeBSD: release/9.0.0/lib/msun/src/s_tgammaf.c 176388 2008-02-18 17:27:11Z das $");
#else
__RCSID("$NetBSD: s_tgammaf.c,v 1.1 2012/05/05 17:54:14 christos Exp $");
#endif
#include <math.h>
/*
* We simply call tgamma() rather than bloating the math library with
* a float-optimized version of it. The reason is that tgammaf() is
* essentially useless, since the function is superexponential and
* floats have very limited range.
*/
float
tgammaf(float x)
{
return (tgamma(x));
}

7
lib/libm/src/s_trunc.c
View File

@@ -15,7 +15,7 @@
__FBSDID("$FreeBSD: src/lib/msun/src/s_trunc.c,v 1.1 2004/06/20 09:25:43 das Exp $");
#endif
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: s_trunc.c,v 1.3 2008/04/25 22:21:53 christos Exp $");
__RCSID("$NetBSD: s_trunc.c,v 1.4 2013/11/13 12:58:11 joerg Exp $");
#endif
/*
@@ -30,6 +30,11 @@ __RCSID("$NetBSD: s_trunc.c,v 1.3 2008/04/25 22:21:53 christos Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_truncl, trunc)
__weak_alias(truncl, trunc)
#endif
static const double huge = 1.0e300;
double

79
lib/libm/src/s_truncl.c Normal file
View File

@@ -0,0 +1,79 @@
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
* From: @(#)s_floor.c 5.1 93/09/24
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: s_truncl.c,v 1.4 2013/11/13 12:58:11 joerg Exp $");
#if 0
__FBSDID("$FreeBSD: head/lib/msun/src/s_truncl.c 176280 2008-02-14 15:10:34Z bde $");
#endif
/*
* truncl(x)
* Return x rounded toward 0 to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to truncl(x).
*/
#include "namespace.h"
#include <float.h>
#include <math.h>
#include <stdint.h>
#include <machine/ieee.h>
#ifdef __HAVE_LONG_DOUBLE
#ifdef __weak_alias
__weak_alias(truncl, _truncl)
#endif
#ifdef LDBL_IMPLICIT_NBIT
#define MANH_SIZE (EXT_FRACHBITS + 1)
#else
#define MANH_SIZE EXT_FRACHBITS
#endif
static const long double huge = 1.0e300;
static const float zero[] = { 0.0, -0.0 };
long double
truncl(long double x)
{
union ieee_ext_u ux = { .extu_ld = x, };
int e = ux.extu_exp - LDBL_MAX_EXP + 1;
if (e < MANH_SIZE - 1) {
if (e < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0)
ux.extu_ld = zero[ux.extu_sign];
} else {
uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
if (((ux.extu_frach & m) | ux.extu_fracl) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
ux.extu_frach &= ~m;
ux.extu_fracl = 0;
}
}
} else if (e < LDBL_MANT_DIG - 1) {
uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
if ((ux.extu_fracl & m) == 0)
return (x); /* x is integral */
if (huge + x > 0.0) /* raise inexact flag */
ux.extu_fracl &= ~m;
}
return (ux.extu_ld);
}
#endif /* __HAVE_LONG_DOUBLE */

6
lib/libm/src/w_fmod.c
View File

@@ -12,7 +12,7 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: w_fmod.c,v 1.9 2002/05/26 22:02:00 wiz Exp $");
__RCSID("$NetBSD: w_fmod.c,v 1.10 2013/11/12 16:48:39 joerg Exp $");
#endif
/*
@@ -22,6 +22,10 @@ __RCSID("$NetBSD: w_fmod.c,v 1.9 2002/05/26 22:02:00 wiz Exp $");
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_fmodl, fmod)
__weak_alias(fmodl, fmod)
#endif
double
fmod(double x, double y) /* wrapper fmod */

46
lib/libm/src/w_fmodl.c Normal file
View File

@@ -0,0 +1,46 @@
/* @(#)w_fmod.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: w_fmodl.c,v 1.4 2013/11/21 13:41:10 martin Exp $");
/*
* wrapper fmodl(x,y)
*/
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
#ifdef __weak_alias
__weak_alias(fmodl, _fmodl)
#endif
long double
fmodl(long double x, long double y) /* wrapper fmod */
{
#ifdef _IEEE_LIBM
return __ieee754_fmodl(x,y);
#else
long double z;
z = __ieee754_fmodl(x,y);
if(_LIB_VERSION == _IEEE_ ||isnan(y)||isnan(x)) return z;
if(y==0.0) {
return __kernel_standard(x,y,227); /* fmod(x,0) */
} else
return z;
#endif
}
#endif

8
lib/libm/src/w_sqrt.c
View File

@@ -12,16 +12,22 @@
#include <sys/cdefs.h>
#if defined(LIBM_SCCS) && !defined(lint)
__RCSID("$NetBSD: w_sqrt.c,v 1.9 2002/05/26 22:02:03 wiz Exp $");
__RCSID("$NetBSD: w_sqrt.c,v 1.10 2013/11/19 19:24:34 joerg Exp $");
#endif
/*
* wrapper sqrt(x)
*/
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#ifndef __HAVE_LONG_DOUBLE
__strong_alias(_sqrtl, sqrt)
__weak_alias(sqrtl, _sqrtl)
#endif
double
sqrt(double x) /* wrapper sqrt */
{

44
lib/libm/src/w_sqrtl.c Normal file
View File

@@ -0,0 +1,44 @@
/* @(#)w_sqrt.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#include <sys/cdefs.h>
__RCSID("$NetBSD: w_sqrtl.c,v 1.2 2013/11/20 11:39:00 joerg Exp $");
/*
* wrapper sqrtl(x)
*/
#include "namespace.h"
#include "math.h"
#include "math_private.h"
#ifdef __HAVE_LONG_DOUBLE
__weak_alias(sqrtl, _sqrtl)
long double
sqrtl(long double x) /* wrapper sqrtl */
{
#ifdef _IEEE_LIBM
return __ieee754_sqrtl(x);
#else
long double z;
z = __ieee754_sqrtl(x);
if(_LIB_VERSION == _IEEE_ || isnan(x)) return z;
if(x<0.0) {
return __kernel_standard(x,x,226); /* sqrtl(negative) */
} else
return z;
#endif
}
#endif /* __HAVE_LONG_DOUBLE */