-/*
- */
-
#ifndef _MATH_H
#define _MATH_H
+#ifndef double
+
+/* Compiler with full float/double support */
+
extern double fabs(double);
extern double floor(double);
extern double ceil(double);
extern double log(double);
extern double log10(double);
+extern double log1p(double);
+extern double log2(double);
+extern double logb(double);
+
+extern double nan(const char *__tagp);
+
extern double pow(double, double);
extern double exp(double);
-extern double eval_poly(double, double *, int);
+extern double scalbln(double, long);
+extern double scalbn(double, int);
+
+extern double hypot(double, double);
+
+#else
+
+/* We have double defined as float .. so fix up the support routines */
+#define fabs(a) fabsf(a)
+#define floor(a) floorf(a)
+#define ceil(a) ceilf(a)
+#define modf(a,b) modff(a,b)
+#define frexp(a,b) frexpf(a,b)
+#define ldexp(a,b) ldexpf(a,b)
+/* FIXME atof equivalence */
+
+#define sqrt(a) sqrtf(a)
+
+#define sin(a) sinf(a)
+#define cos(a) cosf(a)
+#define tan(a) tanf(a)
+#define asin(a) asinf(a)
+#define acos(a) acosf(a)
+#define atan(a) atanf(a)
+#define atan2(a,b) atan2f(a,b)
+#define sinh(a) sinhf(a)
+#define cosh(a) coshf(a)
+#define tanh(a) tanhf(a)
+
+#define log(a) logf(a)
+#define log10(a) log10f(a)
+#define log1p(a) log1pf(a)
+#define log2(a) log2f(a)
+#define logb(a) logbf(a)
+
+#define nan(a) nanf(a)
+
+#define pow(a,b) powf(a,b)
+#define exp(a) expf(a)
+#define hypot(a,b) hypotf(a,b)
+
+#define scalbln(a,b) scalblnf(a,b)
+#define scalbn(a,b) scalblf(a,b)
+
+#endif
+
+extern float acosf(float);
+extern float acoshf(float);
+extern float asinf(float);
+extern float asinhf(float);
+extern float atanf(float);
+extern float atan2f(float);
+extern float atanhf(float);
+extern float ceilf(float);
+extern float fabsf(float);
+extern float floorf(float);
+extern float frexpf(float, int *);
+extern float hypotf(float, float);
+extern float logf(float);
+extern float log10f(float);
+extern float log1pf(float);
+extern float log2f(float);
+extern float logbf(float);
+extern float nanf(const char *__tagp);
+extern float scalblf(float, long);
+extern float scalbnf(float, int);
+extern float sqrtf(float);
+
+/* FIXME: sort out the right NaN's */
+#define __sNaN 0x1.fffff0p128
+#define __NaN 0x1.fffff0p128
+
+#define __FINFINITY 1e40f
+
#endif
Initial target API (SYS5.2 equivalent - ish)
-a64l
+Z80 = implementation for float form in the SDCC library and we should eventually
+use that not our own if possible
+
+a64l CODE
abort FIXME
abs CODE
access CODE
-acos Z80
+acos CODE
+acosf Z80 | CODE
+acosh CODE
+acoshf CODE (check Z80)
advance
alarm CODE
asctime CODE
-asin Z80
+asin CODE
+asinf Z80 | CODE
+asinh CODE
+asinhf CODE (check Z80)
assert CODE
-atan Z80
-atan2 Z80
+atan CODE
+atanf Z80 | CODE
+atan2 CODE
+atan2f Z80 | CODE
+atanh CODE
+atanhf CODE (check Z80)
atof CODE
atoi CODE
atol CODE
-bsearch
+bsearch CODE
calloc CODE
-ceil Z80
+ceil CODE
+ceilf Z80 | CODE
chdir CODE
chmod CODE
chown CODE
clock CODE
close CODE
compile CODE (but no to spec)
-cos Z80
-cosh Z80
+cos CODE
+cosf Z80 | CODE
+cosh CODE
+coshf Z80 | CODE
creat CODE
crypt CODE (but tea not DES)
ctermid CODE
cuserid CODE
-drand48
+drand48 CODE
dup CODE
-encrypt
+encrypt NO (obsolete single DES)
endgrent CODE
endpwent CODE
endutent CODE
-erand48
-erf
-erfc
+erand48 CODE
+erf ULY
+erfc ULY
execl CODE
execle *MISSING ?*
execlp CODE
execvp CODE
exit CODE
_exit CODE
-exp
-fabs Z80
+exp CODE
+expf CODE | Z80
+fabs CODE
+fabsf CODE | Z80
fclose CODE
fcntl CODE
fdopen CODE
fgetpwent CODE
fgets CODE
fileno CODE
-floor Z80
-fmod
+floor CODE
+floorf Z80 | CODE
+fmod CODE
+fmodf CODE
fopen CODE
fork CODE
fprintf CODE
fread CODE
free CODE
freopen CODE
-frexp
+frexp CODE
+frexpf CODE
fscanf CODE
fseek CODE
fstat CODE
ftell CODE
-ftw
+ftw https://opensource.apple.com/source/Libc/Libc-825.26/gen/nftw.c ?
fwrite CODE
gamma
getc CODE
gsignal OBSOLETE
hcreate
hdestroy
-hypot
+hsearch
+hypot CODE
+hypotf CODE
ioctl CODE
isalnum CODE
isalpha CODE
isspace CODE
isupper CODE
isxdigit CODE
-j0
-j1
-jn
-jrand48
+j0 ULY
+j1 ULY
+jn ULY
+jrand48 CODE
kill CODE
-l64a
-lcong48
-ldexp
+l64a CODE
+lcong48 CODE
+ldexp ULY
lfind CODE
link CODE
localtime CODE
lockf NO (flock)
-log Z80
-log10
+log CODE | Z80
+logf ULY
+log10 ULY
longjmp CODE
-lrand48
+lrand48 CODE
lsearch CODE
lseek CODE
mallinfo
memset CODE
mknod CODE
mktemp
-modf Z80
+modf ULY
+modff Z80 | ULY
monitor
mount CODE
-mrand48
+mrand48 CODE
nlist
-nrand48
+nrand48 CODE
open CODE
pause CODE
pclose CODE
perror CODE
pipe CODE
popen CODE
-pow Z80
-printf CODE (check all fmts included) - do libfl version
+pow Z80 | ULY
+printf CODE (check all fmts included)
putc CODE
putchar CODE
putenv CODE
realloc CODE
rewind CODE
scanf CODE
-seed48
+seed48 CODE
setbuf CODE
setgid CODE
setgrent CODE
setjmp CODE
-setkey
+setkey NO (obsolete, single DES)
setpgrp CODE
setuid CODE
setutent CODE
setvbuf CODE
sgetl
signal CODE (but kernel bits to do)
-sin
-sinh
+sin ULY
+sinh ULY
sleep CODE
sprintf CODE
sputl
sqrt Z80
srand CODE
-srand48
-ssignal
+srand48 CODE
+ssignal OBSOLETE
stat CODE
step
stime CODE
strtod CODE
strtok CODE
strtol CODE
-swab *MISSING*
+swab CODE
sync CODE
system CODE
-tan
-tanh
+tan ULY
+tanh ULY
tdelete
tempnam
tfind
vsprintf CODE
wait CODE
write CODE
-y0
-y1
-yn
+y0 ULY (in j0)
+y1 ULY (in j1)
+yn ULY (in jn)
OBJ_CRT0NS = $(SRC_CRT0NS:.s=.rel)
SRC_ASM = enter.s htonl-z80.s htons-z80.s
OBJ_ASM = $(SRC_ASM:.s=.rel)
-SRC_C = __argv.c abort.c asctime.c assert.c atexit.c atoi_small.c
+SRC_C = __argv.c a64l.c abort.c asctime.c assert.c atexit.c atoi_small.c
SRC_C += bcmp.c bcopy.c bsearch.c bzero.c calloc.c cfree.c clock.c closedir.c
SRC_C += closedir_r.c clock_gettime.c clock_getres.c clock_settime.c
SRC_C += creat.c crypt.c ctermid.c ctime.c cuserid.c
-SRC_C += difftime.c err.c errno.c error.c
+SRC_C += difftime.c drand48.c err.c errno.c error.c
SRC_C += execl.c execv.c execvp.c exit.c
SRC_C += fclose.c fflush.c fgetc.c fgetgrent.c fgetpwent.c
SRC_C += fgetpos.c fgets.c fopen.c fork.c fprintf.c fputc.c fputs.c fread.c
SRC_C += getw.c gmtime.c gmtime_r.c grent.c
SRC_C += inet_addr.c inet_aton.c inet_network.c inet_ntoa.c inet_ntop.c inet_pton.c
SRC_C += index.c initgroups.c isatty.c killpg.c
-SRC_C += libintl.c
-SRC_C += localtim.c localtim_r.c lseek.c lsearch.c lstat.c ltoa.c ltostr.c
+SRC_C += l64a.c libintl.c localtim.c localtim_r.c
+SRC_C += lrand48.c lseek.c lsearch.c lstat.c ltoa.c ltostr.c
SRC_C += malloc.c mkfifo.c mkstemps.c nanosleep.c opendir.c opendir_r.c
SRC_C += pause.c perror.c
SRC_C += popen.c printf.c putenv.c putgetch.c putpwent.c putw.c pwent.c qsort.c
#SRC_C += strncpy.c strpbrk.c strrchr.c strspn.c strstr.c strtok.c
#SRC_C += memchr.c memcmp.c memcpy.c memset.c
+SRC_C += acosf.c acoshf.c asinf.c asinhf.c atan2f.c atanf.c atanhf.c
+SRC_C += ceilf.c expf.c
+SRC_C += fabsf.c floorf.c fmodf.c frexpf.c hypotf.c
+SRC_C += scalbnf.c scalbinf.c sqrtf.c
+
SRC_CT += termcap.c
SRC_CURS = $(shell find curses -name '*.c')
--- /dev/null
+Elements of this C library as marked are taken from MUSL
+
+musl as a whole is licensed under the following standard MIT license:
+
+Copyright © 2005-2012 Rich Felker
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+
+
+Portions of this software are contributed or derived from software
+authored by third parties. Complete details on the copyright status of
+all code included in musl follows below:
+
+
+The TRE regular expression implementation (src/regex/reg* and
+src/regex/tre*) is Copyright © 2001-2008 Ville Laurikari and licensed
+under a 2-clause BSD license (license text in the source files). The
+included version has been heavily modified by Rich Felker in 2012, in
+the interests of size, simplicity, and namespace cleanliness.
+
+Most of the math library code (src/math/* and src/complex/*) is
+Copyright © 1993,2004 Sun Microsystems or
+Copyright © 2003-2011 David Schultz or
+Copyright © 2003-2009 Steven G. Kargl or
+Copyright © 2003-2009 Bruce D. Evans or
+Copyright © 2008 Stephen L. Moshier
+and labelled as such. All have been licensed under extremely
+permissive terms. See the comments in the individual files for
+details.
+
+The implementation of DES for crypt (src/misc/crypt_des.c) is
+Copyright © 1994 David Burren. It is licensed under a BSD license.
+
+The implementation of blowfish crypt (src/misc/crypt_blowfish.c) was
+originally written by Solar Designer and placed into the public
+domain. The code also comes with a fallback permissive license for use
+in jurisdictions that may not recognize the public domain.
+
+The smoothsort implementation (src/stdlib/qsort.c) is Copyright © 2011
+Valentin Ochs and is licensed under an MIT-style license.
+
+The BSD PRNG implementation (src/prng/random.c) and XSI search API
+(src/search/*.c) functions are Copyright © 2011 Szabolcs Nagy and
+licensed under following terms: "Permission to use, copy, modify,
+and/or distribute this code for any purpose with or without fee is
+hereby granted. There is no warranty."
+
+The x86_64 port was written by Nicholas J. Kain. Several files (crt)
+were released into the public domain; others are licensed under the
+standard MIT license terms at the top of this file. See individual
+files for their copyright status.
+
+The mips and microblaze ports were originally written by Richard
+Pennington for use in the ellcc project. The original code was adapted
+by Rich Felker for build system and code conventions during upstream
+integration. It is licensed under the standard MIT terms.
+
+The powerpc port was also originally written by Richard Pennington,
+and later supplemented and integrated by John Spencer. It is licensed
+under the standard MIT terms.
+
+All other files which have no copyright comments are original works
+Copyright © 2005-2012 Rich Felker, the main author of this library.
+The decision to exclude such comments is intentional, as it should be
+possible to carry around the complete source code on tiny storage
+media. All public header files (include/*) should be treated as Public
+Domain as they intentionally contain no content which can be covered
+by copyright. Some source modules may fall in this category as well.
+If you believe that a file is so trivial that it should be in the
+Public Domain, please contact me and, if I agree, I will explicitly
+release it from copyright.
+
+The following files are trivial, in my opinion not copyrightable in
+the first place, and hereby explicitly released to the Public Domain:
+
+All public headers: include/*
+Startup files: crt/*
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_acos.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* acos(x)
+ * Method :
+ * acos(x) = pi/2 - asin(x)
+ * acos(-x) = pi/2 + asin(x)
+ * For |x|<=0.5
+ * acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
+ * For x>0.5
+ * acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
+ * = 2asin(sqrt((1-x)/2))
+ * = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
+ * = 2f + (2c + 2s*z*R(z))
+ * where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
+ * for f so that f+c ~ sqrt(z).
+ * For x<-0.5
+ * acos(x) = pi - 2asin(sqrt((1-|x|)/2))
+ * = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
+ *
+ * Special cases:
+ * if x is NaN, return x itself;
+ * if |x|>1, return NaN with invalid signal.
+ *
+ * Function needed: sqrt
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double
+pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
+pio2_hi = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+// FIXME
+static const volatile double
+pio2_lo = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
+static const double
+pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double acos(double x)
+{
+ double z,p,q,r,w,s,c,df;
+ int32_t hx,ix;
+
+ GET_HIGH_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x3ff00000) { /* |x| >= 1 */
+ uint32_t lx;
+
+ GET_LOW_WORD(lx,x);
+ if ((ix-0x3ff00000 | lx) == 0) { /* |x|==1 */
+ if (hx > 0) return 0.0; /* acos(1) = 0 */
+ return pi + 2.0*pio2_lo; /* acos(-1)= pi */
+ }
+ return (x-x)/(x-x); /* acos(|x|>1) is NaN */
+ }
+ if (ix < 0x3fe00000) { /* |x| < 0.5 */
+ if (ix <= 0x3c600000) /* |x| < 2**-57 */
+ return pio2_hi + pio2_lo;
+ z = x*x;
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ r = p/q;
+ return pio2_hi - (x - (pio2_lo-x*r));
+ } else if (hx < 0) { /* x < -0.5 */
+ z = (1.0+x)*0.5;
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ s = sqrt(z);
+ r = p/q;
+ w = r*s-pio2_lo;
+ return pi - 2.0*(s+w);
+ } else { /* x > 0.5 */
+ z = (1.0-x)*0.5;
+ s = sqrt(z);
+ df = s;
+ SET_LOW_WORD(df,0);
+ c = (z-df*df)/(s+df);
+ p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
+ q = 1.0+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
+ r = p/q;
+ w = r*s+c;
+ return 2.0*(df+w);
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.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 <math.h>
+#include <signal.h>
+
+#include "libm.h"
+
+static const float
+pi = 3.1415925026e+00, /* 0x40490fda */
+pio2_hi = 1.5707962513e+00; /* 0x3fc90fda */
+static const volatile float
+pio2_lo = 7.5497894159e-08; /* 0x33a22168 */
+static const float
+pS0 = 1.6666586697e-01,
+pS1 = -4.2743422091e-02,
+pS2 = -8.6563630030e-03,
+qS1 = -7.0662963390e-01;
+
+float acosf(float x)
+{
+ float z,p,q,r,w,s,c,df;
+ int32_t hx,ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x3f800000) { /* |x| >= 1 */
+ if (ix == 0x3f800000) { /* |x| == 1 */
+ if (hx > 0) return 0.0f; /* acos(1) = 0 */
+ return pi + 2.0f*pio2_lo; /* acos(-1)= pi */
+ }
+ raise(SIGFPE);
+ return __NaN;
+/* return (x-x)/(x-x);*/ /* acos(|x|>1) is NaN */
+ }
+ if (ix < 0x3f000000) { /* |x| < 0.5 */
+ if (ix <= 0x32800000) /* |x| < 2**-26 */
+ return pio2_hi + pio2_lo;
+ z = x*x;
+ p = z*(pS0+z*(pS1+z*pS2));
+ q = 1.0f+z*qS1;
+ r = p/q;
+ return pio2_hi - (x - (pio2_lo-x*r));
+ } else if (hx < 0) { /* x < -0.5 */
+ z = (1.0f+x)*0.5f;
+ p = z*(pS0+z*(pS1+z*pS2));
+ q = 1.0f+z*qS1;
+ s = sqrtf(z);
+ r = p/q;
+ w = r*s-pio2_lo;
+ return pi - 2.0f*(s+w);
+ } else { /* x > 0.5 */
+ int32_t idf;
+
+ z = (1.0f-x)*0.5f;
+ s = sqrtf(z);
+ df = s;
+ GET_FLOAT_WORD(idf,df);
+ SET_FLOAT_WORD(df,idf&0xfffff000);
+ c = (z-df*df)/(s+df);
+ p = z*(pS0+z*(pS1+z*pS2));
+ q = 1.0f+z*qS1;
+ r = p/q;
+ w = r*s+c;
+ return 2.0f*(df+w);
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_acosh.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+/* acosh(x)
+ * Method :
+ * Based on
+ * acosh(x) = log [ x + sqrt(x*x-1) ]
+ * we have
+ * acosh(x) := log(x)+ln2, if x is large; else
+ * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else
+ * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1.
+ *
+ * Special cases:
+ * acosh(x) is NaN with signal if x<1.
+ * acosh(NaN) is NaN without signal.
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double
+ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
+
+double acosh(double x)
+{
+ double t;
+ int32_t hx;
+ uint32_t lx;
+
+ EXTRACT_WORDS(hx, lx, x);
+ if (hx < 0x3ff00000) { /* x < 1 */
+ return (x-x)/(x-x);
+ } else if (hx >= 0x41b00000) { /* x > 2**28 */
+ if (hx >= 0x7ff00000) /* x is inf of NaN */
+ return x+x;
+ return log(x) + ln2; /* acosh(huge) = log(2x) */
+ } else if ((hx-0x3ff00000 | lx) == 0) {
+ return 0.0; /* acosh(1) = 0 */
+ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
+ t = x*x;
+ return log(2.0*x - 1.0/(x+sqrt(t-1.0)));
+ } else { /* 1 < x < 2 */
+ t = x-1.0;
+ return log1p(t + sqrt(2.0*t+t*t));
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_acoshf.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 <math.h>
+#include <signal.h>
+#include "libm.h"
+
+static const float
+ln2 = 6.9314718246e-01; /* 0x3f317218 */
+
+float acoshf(float x)
+{
+ float t;
+ int32_t hx;
+
+ GET_FLOAT_WORD(hx, x);
+ if (hx < 0x3f800000) { /* x < 1 */
+ raise(SIGFPE);
+ return __NaN;
+/* return (x-x)/(x-x); */
+ } else if (hx >= 0x4d800000) { /* x > 2**28 */
+ if (hx >= 0x7f800000) /* x is inf of NaN */
+ return x + x;
+ return logf(x) + ln2; /* acosh(huge)=log(2x) */
+ } else if (hx == 0x3f800000) {
+ return 0.0f; /* acosh(1) = 0 */
+ } else if (hx > 0x40000000) { /* 2**28 > x > 2 */
+ t = x*x;
+ return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f)));
+ } else { /* 1 < x < 2 */
+ t = x-1.0f;
+ return log1pf(t + sqrtf(2.0f*t+t*t));
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_asin.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* asin(x)
+ * Method :
+ * Since asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ * we approximate asin(x) on [0,0.5] by
+ * asin(x) = x + x*x^2*R(x^2)
+ * where
+ * R(x^2) is a rational approximation of (asin(x)-x)/x^3
+ * and its remez error is bounded by
+ * |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
+ *
+ * For x in [0.5,1]
+ * asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ * Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
+ * then for x>0.98
+ * asin(x) = pi/2 - 2*(s+s*z*R(z))
+ * = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
+ * For x<=0.98, let pio4_hi = pio2_hi/2, then
+ * f = hi part of s;
+ * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
+ * and
+ * asin(x) = pi/2 - 2*(s+s*z*R(z))
+ * = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
+ * = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
+ *
+ * Special cases:
+ * if x is NaN, return x itself;
+ * if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double
+huge = 1.000e+300,
+pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
+pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
+pio4_hi = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
+/* coefficients for R(x^2) */
+pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
+pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
+pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
+pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
+pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
+pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
+qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
+qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
+qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
+qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+double asin(double x)
+{
+ double t=0.0,w,p,q,c,r,s;
+ int32_t hx,ix;
+
+ GET_HIGH_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x3ff00000) { /* |x|>= 1 */
+ uint32_t lx;
+
+ GET_LOW_WORD(lx, x);
+ if ((ix-0x3ff00000 | lx) == 0)
+ /* asin(1) = +-pi/2 with inexact */
+ return x*pio2_hi + x*pio2_lo;
+ return (x-x)/(x-x); /* asin(|x|>1) is NaN */
+ } else if (ix < 0x3fe00000) { /* |x|<0.5 */
+ if (ix < 0x3e500000) { /* if |x| < 2**-26 */
+ if (huge+x > 1.0)
+ return x; /* return x with inexact if x!=0*/
+ }
+ t = x*x;
+ p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+ q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+ w = p/q;
+ return x + x*w;
+ }
+ /* 1 > |x| >= 0.5 */
+ w = 1.0 - fabs(x);
+ t = w*0.5;
+ p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5)))));
+ q = 1.0+t*(qS1+t*(qS2+t*(qS3+t*qS4)));
+ s = sqrt(t);
+ if (ix >= 0x3FEF3333) { /* if |x| > 0.975 */
+ w = p/q;
+ t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
+ } else {
+ w = s;
+ SET_LOW_WORD(w,0);
+ c = (t-w*w)/(s+w);
+ r = p/q;
+ p = 2.0*s*r-(pio2_lo-2.0*c);
+ q = pio4_hi - 2.0*w;
+ t = pio4_hi - (p-q);
+ }
+ if (hx > 0)
+ return t;
+ return -t;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.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 <math.h>
+#include "libm.h"
+
+static const float
+huge = 1.000e+30,
+/* coefficients for R(x^2) */
+pS0 = 1.6666586697e-01,
+pS1 = -4.2743422091e-02,
+pS2 = -8.6563630030e-03,
+qS1 = -7.0662963390e-01;
+
+static const double
+pio2 = 1.570796326794896558e+00;
+
+float asinf(float x)
+{
+ double s;
+ float t,w,p,q;
+ int32_t hx,ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x3f800000) { /* |x| >= 1 */
+ if (ix == 0x3f800000) /* |x| == 1 */
+ return x*pio2; /* asin(+-1) = +-pi/2 with inexact */
+ return __sNaN; /* asin(|x|>1) is NaN */
+ } else if (ix < 0x3f000000) { /* |x|<0.5 */
+ if (ix < 0x39800000) { /* |x| < 2**-12 */
+ if (huge+x > 1.0f)
+ return x; /* return x with inexact if x!=0 */
+ }
+ t = x*x;
+ p = t*(pS0+t*(pS1+t*pS2));
+ q = 1.0f+t*qS1;
+ w = p/q;
+ return x + x*w;
+ }
+ /* 1 > |x| >= 0.5 */
+ w = 1.0f - fabsf(x);
+ t = w*0.5f;
+ p = t*(pS0+t*(pS1+t*pS2));
+ q = 1.0f+t*qS1;
+ s = sqrt(t);
+ w = p/q;
+ t = pio2-2.0*(s+s*w);
+ if (hx > 0)
+ return t;
+ return -t;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_asinh.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* asinh(x)
+ * Method :
+ * Based on
+ * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ]
+ * we have
+ * asinh(x) := x if 1+x*x=1,
+ * := sign(x)*(log(x)+ln2)) for large |x|, else
+ * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else
+ * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2)))
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double
+ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
+huge= 1.00000000000000000000e+300;
+
+double asinh(double x)
+{
+ double t,w;
+ int32_t hx,ix;
+
+ GET_HIGH_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x7ff00000) /* x is inf or NaN */
+ return x+x;
+ if (ix < 0x3e300000) { /* |x| < 2**-28 */
+ /* return x inexact except 0 */
+ if (huge+x > 1.0)
+ return x;
+ }
+ if (ix > 0x41b00000) { /* |x| > 2**28 */
+ w = log(fabs(x)) + ln2;
+ } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
+ t = fabs(x);
+ w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t));
+ } else { /* 2.0 > |x| > 2**-28 */
+ t = x*x;
+ w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t)));
+ }
+ if (hx > 0)
+ return w;
+ return -w;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_asinhf.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 <math.h>
+#include "libm.h"
+
+static const float
+ln2 = 6.9314718246e-01, /* 0x3f317218 */
+huge= 1.0000000000e+30;
+
+float asinhf(float x)
+{
+ float t,w;
+ int32_t hx,ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x7f800000) /* x is inf or NaN */
+ return x+x;
+ if (ix < 0x31800000) { /* |x| < 2**-28 */
+ /* return x inexact except 0 */
+ if (huge+x > 1.0f)
+ return x;
+ }
+ if (ix > 0x4d800000) { /* |x| > 2**28 */
+ w = logf(fabsf(x)) + ln2;
+ } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */
+ t = fabsf(x);
+ w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t));
+ } else { /* 2.0 > |x| > 2**-28 */
+ t = x*x;
+ w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t)));
+ }
+ if (hx > 0)
+ return w;
+ return -w;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_atan.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* atan(x)
+ * Method
+ * 1. Reduce x to positive by atan(x) = -atan(-x).
+ * 2. According to the integer k=4t+0.25 chopped, t=x, the argument
+ * is further reduced to one of the following intervals and the
+ * arctangent of t is evaluated by the corresponding formula:
+ *
+ * [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
+ * [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
+ * [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
+ * [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
+ * [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double atanhi[] = {
+ 4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
+ 7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
+ 9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
+ 1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
+};
+
+static const double atanlo[] = {
+ 2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
+ 3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
+ 1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
+ 6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
+};
+
+static const double aT[] = {
+ 3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
+ -1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
+ 1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
+ -1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
+ 9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
+ -7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
+ 6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
+ -5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
+ 4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
+ -3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
+ 1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
+};
+
+static const double huge = 1.0e300;
+
+double atan(double x)
+{
+ double w,s1,s2,z;
+ int32_t ix,hx,id;
+
+ GET_HIGH_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x44100000) { /* if |x| >= 2^66 */
+ uint32_t low;
+
+ GET_LOW_WORD(low, x);
+ if (ix > 0x7ff00000 ||
+ (ix == 0x7ff00000 && low != 0)) /* NaN */
+ return x+x;
+ if (hx > 0)
+ return atanhi[3] + *(volatile double *)&atanlo[3];
+ else
+ return -atanhi[3] - *(volatile double *)&atanlo[3];
+ }
+ if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
+ if (ix < 0x3e400000) { /* |x| < 2^-27 */
+ /* raise inexact */
+ if (huge+x > 1.0)
+ return x;
+ }
+ id = -1;
+ } else {
+ x = fabs(x);
+ if (ix < 0x3ff30000) { /* |x| < 1.1875 */
+ if (ix < 0x3fe60000) { /* 7/16 <= |x| < 11/16 */
+ id = 0;
+ x = (2.0*x-1.0)/(2.0+x);
+ } else { /* 11/16 <= |x| < 19/16 */
+ id = 1;
+ x = (x-1.0)/(x+1.0);
+ }
+ } else {
+ if (ix < 0x40038000) { /* |x| < 2.4375 */
+ id = 2;
+ x = (x-1.5)/(1.0+1.5*x);
+ } else { /* 2.4375 <= |x| < 2^66 */
+ id = 3;
+ x = -1.0/x;
+ }
+ }
+ }
+ /* end of argument reduction */
+ z = x*x;
+ w = z*z;
+ /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+ s1 = z*(aT[0]+w*(aT[2]+w*(aT[4]+w*(aT[6]+w*(aT[8]+w*aT[10])))));
+ s2 = w*(aT[1]+w*(aT[3]+w*(aT[5]+w*(aT[7]+w*aT[9]))));
+ if (id < 0)
+ return x - x*(s1+s2);
+ z = atanhi[id] - (x*(s1+s2) - atanlo[id] - x);
+ return hx < 0 ? -z : z;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+/* atan2(y,x)
+ * Method :
+ * 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
+ * 2. Reduce x to positive by (if x and y are unexceptional):
+ * ARG (x+iy) = arctan(y/x) ... if x > 0,
+ * ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
+ *
+ * Special cases:
+ *
+ * ATAN2((anything), NaN ) is NaN;
+ * ATAN2(NAN , (anything) ) is NaN;
+ * ATAN2(+-0, +(anything but NaN)) is +-0 ;
+ * ATAN2(+-0, -(anything but NaN)) is +-pi ;
+ * ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
+ * ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
+ * ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
+ * ATAN2(+-INF,+INF ) is +-pi/4 ;
+ * ATAN2(+-INF,-INF ) is +-3pi/4;
+ * ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+#include <math.h>
+#include "libm.h"
+
+// FIXME
+static const volatile double
+tiny = 1.0e-300;
+static const double
+pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
+pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
+pi = 3.1415926535897931160E+00; /* 0x400921FB, 0x54442D18 */
+static const volatile double
+pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
+
+double atan2(double y, double x)
+{
+ double z;
+ int32_t k,m,hx,hy,ix,iy;
+ uint32_t lx,ly;
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff;
+ EXTRACT_WORDS(hy, ly, y);
+ iy = hy & 0x7fffffff;
+ if ((ix|((lx|-lx)>>31)) > 0x7ff00000 ||
+ (iy|((ly|-ly)>>31)) > 0x7ff00000) /* x or y is NaN */
+ return x+y;
+ if ((hx-0x3ff00000 | lx) == 0) /* x = 1.0 */
+ return atan(y);
+ m = ((hy>>31)&1) | ((hx>>30)&2); /* 2*sign(x)+sign(y) */
+
+ /* when y = 0 */
+ if ((iy|ly) == 0) {
+ switch(m) {
+ case 0:
+ case 1: return y; /* atan(+-0,+anything)=+-0 */
+ case 2: return pi+tiny; /* atan(+0,-anything) = pi */
+ case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
+ }
+ }
+ /* when x = 0 */
+ if ((ix|lx) == 0)
+ return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
+ /* when x is INF */
+ if (ix == 0x7ff00000) {
+ if (iy == 0x7ff00000) {
+ switch(m) {
+ case 0: return pi_o_4+tiny; /* atan(+INF,+INF) */
+ case 1: return -pi_o_4-tiny; /* atan(-INF,+INF) */
+ case 2: return 3.0*pi_o_4+tiny; /* atan(+INF,-INF) */
+ case 3: return -3.0*pi_o_4-tiny; /* atan(-INF,-INF) */
+ }
+ } else {
+ switch(m) {
+ case 0: return 0.0; /* atan(+...,+INF) */
+ case 1: return -0.0; /* atan(-...,+INF) */
+ case 2: return pi+tiny; /* atan(+...,-INF) */
+ case 3: return -pi-tiny; /* atan(-...,-INF) */
+ }
+ }
+ }
+ /* when y is INF */
+ if (iy == 0x7ff00000)
+ return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
+
+ /* compute y/x */
+ k = (iy-ix)>>20;
+ if (k > 60) { /* |y/x| > 2**60 */
+ z = pi_o_2+0.5*pi_lo;
+ m &= 1;
+ } else if (hx < 0 && k < -60) /* 0 > |y|/x > -2**-60 */
+ z = 0.0;
+ else /* safe to do y/x */
+ z = atan(fabs(y/x));
+ switch (m) {
+ case 0: return z; /* atan(+,+) */
+ case 1: return -z; /* atan(-,+) */
+ case 2: return pi - (z-pi_lo); /* atan(+,-) */
+ default: /* case 3 */
+ return (z-pi_lo) - pi; /* atan(-,-) */
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_atan2f.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 <math.h>
+#include "libm.h"
+
+static const volatile float
+tiny = 1.0e-30;
+static const float
+pi_o_4 = 7.8539818525e-01, /* 0x3f490fdb */
+pi_o_2 = 1.5707963705e+00, /* 0x3fc90fdb */
+pi = 3.1415927410e+00; /* 0x40490fdb */
+static const volatile float
+pi_lo = -8.7422776573e-08; /* 0xb3bbbd2e */
+
+float atan2f(float y, float x)
+{
+ float z;
+ int32_t k,m,hx,hy,ix,iy;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ GET_FLOAT_WORD(hy, y);
+ iy = hy & 0x7fffffff;
+ if (ix > 0x7f800000 || iy > 0x7f800000) /* x or y is NaN */
+ return x+y;
+ if (hx == 0x3f800000) /* x=1.0 */
+ return atanf(y);
+ m = ((hy>>31)&1) | ((hx>>30)&2); /* 2*sign(x)+sign(y) */
+
+ /* when y = 0 */
+ if (iy == 0) {
+ switch (m) {
+ case 0:
+ case 1: return y; /* atan(+-0,+anything)=+-0 */
+ case 2: return pi+tiny; /* atan(+0,-anything) = pi */
+ case 3: return -pi-tiny; /* atan(-0,-anything) =-pi */
+ }
+ }
+ /* when x = 0 */
+ if (ix == 0)
+ return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
+ /* when x is INF */
+ if (ix == 0x7f800000) {
+ if (iy == 0x7f800000) {
+ switch (m) {
+ case 0: return pi_o_4+tiny; /* atan(+INF,+INF) */
+ case 1: return -pi_o_4-tiny; /* atan(-INF,+INF) */
+ case 2: return 3.0f*pi_o_4+tiny; /*atan(+INF,-INF)*/
+ case 3: return -3.0f*pi_o_4-tiny; /*atan(-INF,-INF)*/
+ }
+ } else {
+ switch (m) {
+ case 0: return 0.0f; /* atan(+...,+INF) */
+ case 1: return -0.0f; /* atan(-...,+INF) */
+ case 2: return pi+tiny; /* atan(+...,-INF) */
+ case 3: return -pi-tiny; /* atan(-...,-INF) */
+ }
+ }
+ }
+ /* when y is INF */
+ if (iy == 0x7f800000)
+ return hy < 0 ? -pi_o_2-tiny : pi_o_2+tiny;
+
+ /* compute y/x */
+ k = (iy-ix)>>23;
+ if (k > 26) { /* |y/x| > 2**26 */
+ z = pi_o_2 + 0.5f*pi_lo;
+ m &= 1;
+ } else if (k < -26 && hx < 0) /* 0 > |y|/x > -2**-26 */
+ z = 0.0;
+ else /* safe to do y/x */
+ z = atanf(fabsf(y/x));
+ switch (m) {
+ case 0: return z; /* atan(+,+) */
+ case 1: return -z; /* atan(-,+) */
+ case 2: return pi - (z-pi_lo); /* atan(+,-) */
+ default: /* case 3 */
+ return (z-pi_lo) - pi; /* atan(-,-) */
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.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 <math.h>
+#include "libm.h"
+
+static const float atanhi[] = {
+ 4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
+ 7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
+ 9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
+ 1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
+};
+
+static const float atanlo[] = {
+ 5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
+ 3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
+ 3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
+ 7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
+};
+
+static const float aT[] = {
+ 3.3333328366e-01,
+ -1.9999158382e-01,
+ 1.4253635705e-01,
+ -1.0648017377e-01,
+ 6.1687607318e-02,
+};
+
+static const float huge = 1.0e30;
+
+float atanf(float x)
+{
+ float w,s1,s2,z;
+ int32_t ix,hx,id;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix >= 0x4c800000) { /* if |x| >= 2**26 */
+ if (ix > 0x7f800000) /* NaN */
+ return x+x;
+ if (hx > 0)
+ return atanhi[3] + *(volatile float *)&atanlo[3];
+ else
+ return -atanhi[3] - *(volatile float *)&atanlo[3];
+ }
+ if (ix < 0x3ee00000) { /* |x| < 0.4375 */
+ if (ix < 0x39800000) { /* |x| < 2**-12 */
+ /* raise inexact */
+ if(huge+x>1.0f)
+ return x;
+ }
+ id = -1;
+ } else {
+ x = fabsf(x);
+ if (ix < 0x3f980000) { /* |x| < 1.1875 */
+ if (ix < 0x3f300000) { /* 7/16 <= |x| < 11/16 */
+ id = 0;
+ x = (2.0f*x - 1.0f)/(2.0f + x);
+ } else { /* 11/16 <= |x| < 19/16 */
+ id = 1;
+ x = (x - 1.0f)/(x + 1.0f);
+ }
+ } else {
+ if (ix < 0x401c0000) { /* |x| < 2.4375 */
+ id = 2;
+ x = (x - 1.5f)/(1.0f + 1.5f*x);
+ } else { /* 2.4375 <= |x| < 2**26 */
+ id = 3;
+ x = -1.0f/x;
+ }
+ }
+ }
+ /* end of argument reduction */
+ z = x*x;
+ w = z*z;
+ /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+ s1 = z*(aT[0]+w*(aT[2]+w*aT[4]));
+ s2 = w*(aT[1]+w*aT[3]);
+ if (id < 0)
+ return x - x*(s1+s2);
+ z = atanhi[id] - ((x*(s1+s2) - atanlo[id]) - x);
+ return hx < 0 ? -z : z;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_atanh.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ */
+/* atanh(x)
+ * Method :
+ * 1.Reduced x to positive by atanh(-x) = -atanh(x)
+ * 2.For x>=0.5
+ * 1 2x x
+ * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------)
+ * 2 1 - x 1 - x
+ *
+ * For x<0.5
+ * atanh(x) = 0.5*log1p(2x+2x*x/(1-x))
+ *
+ * Special cases:
+ * atanh(x) is NaN if |x| > 1 with signal;
+ * atanh(NaN) is that NaN with no signal;
+ * atanh(+-1) is +-INF with signal.
+ *
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double huge = 1e300;
+
+double atanh(double x)
+{
+ double t;
+ int32_t hx,ix;
+ uint32_t lx;
+
+ EXTRACT_WORDS(hx, lx, x);
+ ix = hx & 0x7fffffff;
+ if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */
+ return __sNaN;
+ if (ix == 0x3ff00000)
+ return __Nan;
+ if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */
+ return x;
+ SET_HIGH_WORD(x, ix);
+ if (ix < 0x3fe00000) { /* x < 0.5 */
+ t = x+x;
+ t = 0.5*log1p(t + t*x/(1.0-x));
+ } else
+ t = 0.5*log1p((x+x)/(1.0-x));
+ if (hx >= 0)
+ return t;
+ return -t;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_atanhf.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 <math.h>
+#include "libm.h"
+
+static const float huge = 1e30;
+
+float atanhf(float x)
+{
+ float t;
+ int32_t hx,ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix > 0x3f800000) /* |x| > 1 */
+ return __sNaN;
+ if (ix == 0x3f800000)
+ return __NaN;
+ if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */
+ return x;
+ SET_FLOAT_WORD(x, ix);
+ if (ix < 0x3f000000) { /* x < 0.5 */
+ t = x+x;
+ t = 0.5f*log1pf(t + t*x/(1.0f-x));
+ } else
+ t = 0.5f*log1pf((x+x)/(1.0f-x));
+ if (hx >= 0)
+ return t;
+ return -t;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ceil.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/*
+ * ceil(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to ceil(x).
+ */
+
+#include "libm.h"
+
+static const double huge = 1.0e300;
+
+double ceil(double x)
+{
+ int32_t i0,i1,j0;
+ uint32_t i,j;
+
+ EXTRACT_WORDS(i0, i1, x);
+ // FIXME signed shift
+ j0 = ((i0>>20)&0x7ff) - 0x3ff;
+ if (j0 < 20) {
+ if (j0 < 0) {
+ /* raise inexact if x != 0 */
+ if (huge+x > 0.0) {
+ if (i0 < 0) {
+ i0 = 0x80000000;
+ i1=0;
+ } else if ((i0|i1) != 0) {
+ i0=0x3ff00000;
+ i1=0;
+ }
+ }
+ } else {
+ i = 0x000fffff>>j0;
+ if (((i0&i)|i1) == 0) /* x is integral */
+ return x;
+ /* raise inexact flag */
+ if (huge+x > 0.0) {
+ if (i0 > 0)
+ i0 += 0x00100000>>j0;
+ i0 &= ~i;
+ i1 = 0;
+ }
+ }
+ } else if (j0 > 51) {
+ if (j0 == 0x400) /* inf or NaN */
+ return x+x;
+ return x; /* x is integral */
+ } else {
+ i = (uint32_t)0xffffffff>>(j0-20);
+ if ((i1&i) == 0)
+ return x; /* x is integral */
+ /* raise inexact flag */
+ if (huge+x > 0.0) {
+ if (i0 > 0) {
+ if (j0 == 20)
+ i0 += 1;
+ else {
+ j = i1 + (1<<(52-j0));
+ if (j < i1) /* got a carry */
+ i0 += 1;
+ i1 = j;
+ }
+ }
+ i1 &= ~i;
+ }
+ }
+ INSERT_WORDS(x, i0, i1);
+ return x;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_ceilf.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 "libm.h"
+
+static const float huge = 1.0e30;
+
+float ceilf(float x)
+{
+ int32_t i0,j0;
+ uint32_t i;
+
+ GET_FLOAT_WORD(i0, x);
+ j0 = ((i0>>23)&0xff) - 0x7f;
+ if (j0 < 23) {
+ if (j0 < 0) {
+ /* raise inexact if x != 0 */
+ if (huge+x > 0.0f) {
+ if (i0 < 0)
+ i0 = 0x80000000;
+ else if(i0 != 0)
+ i0 = 0x3f800000;
+ }
+ } else {
+ i = 0x007fffff>>j0;
+ if ((i0&i) == 0)
+ return x; /* x is integral */
+ /* raise inexact flag */
+ if (huge+x > 0.0f) {
+ if (i0 > 0)
+ i0 += 0x00800000>>j0;
+ i0 &= ~i;
+ }
+ }
+ } else {
+ if (j0 == 0x80) /* inf or NaN */
+ return x+x;
+ return x; /* x is integral */
+ }
+ SET_FLOAT_WORD(x, i0);
+ return x;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cos.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* cos(x)
+ * Return cosine function of x.
+ *
+ * kernel function:
+ * __sin ... sine function on [-pi/4,pi/4]
+ * __cos ... cosine function on [-pi/4,pi/4]
+ * __rem_pio2 ... argument reduction routine
+ *
+ * Method.
+ * Let S,C and T denote the sin, cos and tan respectively on
+ * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
+ * in [-pi/4 , +pi/4], and let n = k mod 4.
+ * We have
+ *
+ * n sin(x) cos(x) tan(x)
+ * ----------------------------------------------------------
+ * 0 S C T
+ * 1 C -S -1/T
+ * 2 -S -C T
+ * 3 -C S -1/T
+ * ----------------------------------------------------------
+ *
+ * Special cases:
+ * Let trig be any of sin, cos, or tan.
+ * trig(+-INF) is NaN, with signals;
+ * trig(NaN) is that NaN;
+ *
+ * Accuracy:
+ * TRIG(x) returns trig(x) nearly rounded
+ */
+
+#include <math.h>
+#include "libm.h"
+
+double cos(double x)
+{
+ double y[2],z=0.0;
+ int32_t n, ix;
+
+ GET_HIGH_WORD(ix, x);
+
+ /* |x| ~< pi/4 */
+ ix &= 0x7fffffff;
+ if (ix <= 0x3fe921fb) {
+ if (ix < 0x3e46a09e) /* if x < 2**-27 * sqrt(2) */
+ /* raise inexact if x != 0 */
+ if ((int)x == 0)
+ return 1.0;
+ return __cos(x, z);
+ }
+
+ /* cos(Inf or NaN) is NaN */
+ if (ix >= 0x7ff00000)
+ return x-x;
+
+ /* argument reduction needed */
+ n = __rem_pio2(x, y);
+ switch (n&3) {
+ case 0: return __cos(y[0], y[1]);
+ case 1: return -__sin(y[0], y[1], 1);
+ case 2: return -__cos(y[0], y[1]);
+ default:
+ return __sin(y[0], y[1], 1);
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_cosf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * Optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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 <math.h>
+#include "libm.h"
+
+/* Small multiples of pi/2 rounded to double precision. */
+static const double
+c1pio2 = 1*M_PI_2, /* 0x3FF921FB, 0x54442D18 */
+c2pio2 = 2*M_PI_2, /* 0x400921FB, 0x54442D18 */
+c3pio2 = 3*M_PI_2, /* 0x4012D97C, 0x7F3321D2 */
+c4pio2 = 4*M_PI_2; /* 0x401921FB, 0x54442D18 */
+
+float cosf(float x)
+{
+ double y;
+ int32_t n, hx, ix;
+
+ GET_FLOAT_WORD(hx, x);
+ ix = hx & 0x7fffffff;
+ if (ix <= 0x3f490fda) { /* |x| ~<= pi/4 */
+ if (ix < 0x39800000) /* |x| < 2**-12 */
+ if ((int)x == 0) /* raise inexact if x != 0 */
+ return 1.0;
+ return __cosdf(x);
+ }
+ if (ix <= 0x407b53d1) { /* |x| ~<= 5*pi/4 */
+ if (ix > 0x4016cbe3) /* |x| ~> 3*pi/4 */
+ return -__cosdf(hx > 0 ? x-c2pio2 : x+c2pio2);
+ else {
+ if (hx > 0)
+ return __sindf(c1pio2 - x);
+ else
+ return __sindf(x + c1pio2);
+ }
+ }
+ if (ix <= 0x40e231d5) { /* |x| ~<= 9*pi/4 */
+ if (ix > 0x40afeddf) /* |x| ~> 7*pi/4 */
+ return __cosdf(hx > 0 ? x-c4pio2 : x+c4pio2);
+ else {
+ if (hx > 0)
+ return __sindf(x - c3pio2);
+ else
+ return __sindf(-c3pio2 - x);
+ }
+ }
+
+ /* cos(Inf or NaN) is NaN */
+ if (ix >= 0x7f800000)
+ return x-x;
+
+ /* general argument reduction needed */
+ n = __rem_pio2f(x,&y);
+ switch (n&3) {
+ case 0: return __cosdf(y);
+ case 1: return __sindf(-y);
+ case 2: return -__cosdf(y);
+ default:
+ return __sindf(y);
+ }
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_cosh.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* cosh(x)
+ * Method :
+ * mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
+ * 1. Replace x by |x| (cosh(x) = cosh(-x)).
+ * 2.
+ * [ exp(x) - 1 ]^2
+ * 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
+ * 2*exp(x)
+ *
+ * exp(x) + 1/exp(x)
+ * ln2/2 <= x <= 22 : cosh(x) := -------------------
+ * 2
+ * 22 <= x <= lnovft : cosh(x) := exp(x)/2
+ * lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
+ * ln2ovft < x : cosh(x) := huge*huge (overflow)
+ *
+ * Special cases:
+ * cosh(x) is |x| if x is +INF, -INF, or NaN.
+ * only cosh(0)=1 is exact for finite x.
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double huge = 1.0e300;
+
+double cosh(double x)
+{
+ double t, w;
+ int32_t ix;
+
+ GET_HIGH_WORD(ix, x);
+ ix &= 0x7fffffff;
+
+ /* x is INF or NaN */
+ if (ix >= 0x7ff00000)
+ return x*x;
+
+ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+ if (ix < 0x3fd62e43) {
+ t = expm1(fabs(x));
+ w = 1.0+t;
+ if (ix < 0x3c800000)
+ return w; /* cosh(tiny) = 1 */
+ return 1.0 + (t*t)/(w+w);
+ }
+
+ /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
+ if (ix < 0x40360000) {
+ t = exp(fabs(x));
+ return 0.5*t + 0.5/t;
+ }
+
+ /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
+ if (ix < 0x40862E42)
+ return 0.5*exp(fabs(x));
+
+ /* |x| in [log(maxdouble), overflowthresold] */
+ if (ix <= 0x408633CE)
+ return __expo2(fabs(x));
+
+ /* |x| > overflowthresold, cosh(x) overflow */
+ return huge*huge;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_coshf.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 <math.h>
+#include "libm.h"
+
+static const float huge = 1.0e30;
+
+float coshf(float x)
+{
+ float t, w;
+ int32_t ix;
+
+ GET_FLOAT_WORD(ix, x);
+ ix &= 0x7fffffff;
+
+ /* x is INF or NaN */
+ if (ix >= 0x7f800000)
+ return x*x;
+
+ /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+ if (ix < 0x3eb17218) {
+ t = expm1f(fabsf(x));
+ w = 1.0f+t;
+ if (ix<0x39800000)
+ return 1.0f; /* cosh(tiny) = 1 */
+ return 1.0f + (t*t)/(w+w);
+ }
+
+ /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
+ if (ix < 0x41100000) {
+ t = expf(fabsf(x));
+ return 0.5f*t + 0.5f/t;
+ }
+
+ /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */
+ if (ix < 0x42b17217)
+ return 0.5f*expf(fabsf(x));
+
+ /* |x| in [log(maxfloat), overflowthresold] */
+ if (ix <= 0x42b2d4fc)
+ return __expo2f(fabsf(x));
+
+ /* |x| > overflowthresold, cosh(x) overflow */
+ return huge*huge;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ * 1. Argument reduction:
+ * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ * Given x, find r and integer k such that
+ *
+ * x = k*ln2 + r, |r| <= 0.5*ln2.
+ *
+ * Here r will be represented as r = hi-lo for better
+ * accuracy.
+ *
+ * 2. Approximation of exp(r) by a special rational function on
+ * the interval [0,0.34658]:
+ * Write
+ * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ * We use a special Remez algorithm on [0,0.34658] to generate
+ * a polynomial of degree 5 to approximate R. The maximum error
+ * of this polynomial approximation is bounded by 2**-59. In
+ * other words,
+ * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ * (where z=r*r, and the values of P1 to P5 are listed below)
+ * and
+ * | 5 | -59
+ * | 2.0+P1*z+...+P5*z - R(z) | <= 2
+ * | |
+ * The computation of exp(r) thus becomes
+ * 2*r
+ * exp(r) = 1 + ----------
+ * R(r) - r
+ * r*c(r)
+ * = 1 + r + ----------- (for better accuracy)
+ * 2 - c(r)
+ * where
+ * 2 4 10
+ * c(r) = r - (P1*r + P2*r + ... + P5*r ).
+ *
+ * 3. Scale back to obtain exp(x):
+ * From step 1, we have
+ * exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ * exp(INF) is INF, exp(NaN) is NaN;
+ * exp(-INF) is 0, and
+ * for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ * according to an error analysis, the error is always less than
+ * 1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ * For IEEE double
+ * if x > 709.782712893383973096 then exp(x) overflows
+ * if x < -745.133219101941108420 then exp(x) underflows
+ */
+
+#include "libm.h"
+
+static const double
+half[2] = {0.5,-0.5},
+ln2hi = 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
+ln2lo = 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
+invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
+P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
+P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
+P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
+P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
+P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+double exp(double x)
+{
+ double hi, lo, c, xx;
+ int k, sign;
+ uint32_t hx;
+
+ GET_HIGH_WORD(hx, x);
+ sign = hx>>31;
+ hx &= 0x7fffffff; /* high word of |x| */
+
+ /* special cases */
+ if (hx >= 0x40862e42) { /* if |x| >= 709.78... */
+ if (isnan(x))
+ return x;
+ if (hx == 0x7ff00000 && sign) /* -inf */
+ return 0;
+ if (x > 709.782712893383973096) {
+ /* overflow if x!=inf */
+ x = (double)(0x1p1023 * x);
+ return x;
+ }
+ if (x < -745.13321910194110842) {
+ /* underflow */
+ x = (double)(0x1p-1000 * 0x1p-1000);
+ return x;
+ }
+ }
+
+ /* argument reduction */
+ if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
+ if (hx >= 0x3ff0a2b2) /* if |x| >= 1.5 ln2 */
+ k = (int)(invln2*x + half[sign]);
+ else
+ k = 1 - sign - sign;
+ hi = x - k*ln2hi; /* k*ln2hi is exact here */
+ lo = k*ln2lo;
+ x = (double)(hi - lo);
+ } else if (hx > 0x3e300000) { /* if |x| > 2**-28 */
+ k = 0;
+ hi = x;
+ lo = 0;
+ } else {
+ volatile float dummy;
+ /* inexact if x!=0 */
+ dummy = 0x1p1023 + x; /* ensure we evaluate and maybe fault */
+ return 1 + x;
+ }
+
+ /* x is now in primary range */
+ xx = x*x;
+ c = x - xx*(P1+xx*(P2+xx*(P3+xx*(P4+xx*P5))));
+ x = 1 + (x*c/(2-c) - lo + hi);
+ if (k == 0)
+ return x;
+ return scalbn(x, k);
+}
--- /dev/null
+/*
+ * 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 <math.h>
+#include "libm.h"
+
+static const float
+half[2] = {0.5,-0.5},
+ln2hi = 6.9314575195e-1f, /* 0x3f317200 */
+ln2lo = 1.4286067653e-6f, /* 0x35bfbe8e */
+invln2 = 1.4426950216e+0f, /* 0x3fb8aa3b */
+/*
+ * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
+ * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74
+ */
+P1 = 1.6666625440e-1f, /* 0xaaaa8f.0p-26 */
+P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */
+
+float expf(float x)
+{
+ float hi, lo, c, xx;
+ int k, sign;
+ uint32_t hx;
+
+ GET_FLOAT_WORD(hx, x);
+ sign = hx >> 31; /* sign bit of x */
+ hx &= 0x7fffffff; /* high word of |x| */
+
+ /* special cases */
+ if (hx >= 0x42b17218) { /* if |x| >= 88.722839f or NaN */
+ if (hx > 0x7f800000) /* NaN */
+ return x;
+ if (!sign) {
+ /* overflow if x!=inf */
+ x = (float)(x * 0x1p127f);
+ return x;
+ }
+ if (hx == 0x7f800000) /* -inf */
+ return 0;
+ if (hx >= 0x42cff1b5) { /* x <= -103.972084f */
+ /* underflow */
+ x = (float)(0x1p-100f*0x1p-100f);
+ return x;
+ }
+ }
+
+ /* argument reduction */
+ if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
+ if (hx > 0x3f851592) /* if |x| > 1.5 ln2 */
+ k = invln2*x + half[sign];
+ else
+ k = 1 - sign - sign;
+ hi = x - k*ln2hi; /* k*ln2hi is exact here */
+ lo = k*ln2lo;
+ x = (float)(hi - lo);
+ } else if (hx > 0x39000000) { /* |x| > 2**-14 */
+ k = 0;
+ hi = x;
+ lo = 0;
+ } else {
+ /* raise inexact */
+ volatile float dummy = 0x1p127f + x;
+ return 1 + x;
+ }
+
+ /* x is now in primary range */
+ xx = x*x;
+ c = x - xx*(P1+xx*P2);
+ x = 1 + (x*c/(2-c) - lo + hi);
+ if (k == 0)
+ return x;
+ return scalbnf(x, k);
+}
--- /dev/null
+/* From MUSL */
+
+#include "math.h"
+#include "libm.h"
+
+double fabs(double x)
+{
+ union dshape u;
+
+ u.value = x;
+ u.bits &= (uint64_t)-1 / 2;
+ return u.value;
+}
--- /dev/null
+/* From MUSL */
+
+#include "libm.h"
+
+float fabsf(float x)
+{
+ union fshape u;
+
+ u.value = x;
+ u.bits &= (uint32_t)-1 / 2;
+ return u.value;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_floor.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/*
+ * floor(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to floor(x).
+ */
+
+#include "libm.h"
+
+static const double huge = 1.0e300;
+
+double floor(double x)
+{
+ int32_t i0,i1,j0;
+ uint32_t i,j;
+
+ EXTRACT_WORDS(i0, i1, x);
+ // FIXME: signed shift
+ j0 = ((i0>>20)&0x7ff) - 0x3ff;
+ if (j0 < 20) {
+ if (j0 < 0) { /* |x| < 1 */
+ /* raise inexact if x != 0 */
+ if (huge+x > 0.0) {
+ if (i0 >= 0) { /* x >= 0 */
+ i0 = i1 = 0;
+ } else if (((i0&0x7fffffff)|i1) != 0) {
+ i0 = 0xbff00000;
+ i1 = 0;
+ }
+ }
+ } else {
+ i = 0x000fffff>>j0;
+ if (((i0&i)|i1) == 0)
+ return x; /* x is integral */
+ /* raise inexact flag */
+ if (huge+x > 0.0) {
+ if (i0 < 0)
+ i0 += 0x00100000>>j0;
+ i0 &= ~i;
+ i1 = 0;
+ }
+ }
+ } else if (j0 > 51) {
+ if (j0 == 0x400)
+ return x+x; /* inf or NaN */
+ else
+ return x; /* x is integral */
+ } else {
+ i = (uint32_t)0xffffffff>>(j0-20);
+ if ((i1&i) == 0)
+ return x; /* x is integral */
+ /* raise inexact flag */
+ if (huge+x > 0.0) {
+ if (i0 < 0) {
+ if (j0 == 20)
+ i0++;
+ else {
+ j = i1+(1<<(52-j0));
+ if (j < i1)
+ i0++; /* got a carry */
+ i1 = j;
+ }
+ }
+ i1 &= ~i;
+ }
+ }
+ INSERT_WORDS(x, i0, i1);
+ return x;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/s_floorf.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.
+ * ====================================================
+ */
+/*
+ * floorf(x)
+ * Return x rounded toward -inf to integral value
+ * Method:
+ * Bit twiddling.
+ * Exception:
+ * Inexact flag raised if x not equal to floorf(x).
+ */
+
+#include "libm.h"
+
+static const float huge = 1.0e30;
+
+float floorf(float x)
+{
+ int32_t i0,j0;
+ uint32_t i;
+
+ GET_FLOAT_WORD(i0, x);
+ // FIXME: signed shift
+ j0 = ((i0>>23)&0xff) - 0x7f;
+ if (j0 < 23) {
+ if (j0 < 0) { /* |x| < 1 */
+ /* raise inexact if x != 0 */
+ if (huge+x > 0.0f) {
+ if (i0 >= 0) /* x >= 0 */
+ i0 = 0;
+ else if ((i0&0x7fffffff) != 0)
+ i0 = 0xbf800000;
+ }
+ } else {
+ i = 0x007fffff>>j0;
+ if ((i0&i) == 0)
+ return x; /* x is integral */
+ /* raise inexact flag */
+ if (huge+x > 0.0f) {
+ if (i0 < 0)
+ i0 += 0x00800000>>j0;
+ i0 &= ~i;
+ }
+ }
+ } else {
+ if (j0 == 0x80) /* inf or NaN */
+ return x+x;
+ else
+ return x; /* x is integral */
+ }
+ SET_FLOAT_WORD(x, i0);
+ return x;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_fmod.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/*
+ * fmod(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const double Zero[] = {0.0, -0.0,};
+
+double fmod(double x, double y)
+{
+ int32_t n,hx,hy,hz,ix,iy,sx,i;
+ uint32_t lx,ly,lz;
+
+ EXTRACT_WORDS(hx, lx, x);
+ EXTRACT_WORDS(hy, ly, y);
+ 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) /* |x| < |y| */
+ return x;
+ if (lx == ly) /* |x| = |y|, return x*0 */
+ return Zero[(uint32_t)sx>>31];
+ }
+
+ /* 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;
+ 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 {
+ if ((hz|lz) == 0) /* return sign(x)*0 */
+ return Zero[(uint32_t)sx>>31];
+ hx = hz+hz+(lz>>31);
+ 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[(uint32_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));
+ INSERT_WORDS(x, hx|sx, lx);
+ } else { /* subnormal output */
+ n = -1022 - iy;
+ if (n <= 20) {
+ lx = (lx>>n)|((uint32_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;
+ }
+ INSERT_WORDS(x, hx|sx, lx);
+ }
+ return x; /* exact output */
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_fmodf.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.
+ * ====================================================
+ */
+/*
+ * fmodf(x,y)
+ * Return x mod y in exact arithmetic
+ * Method: shift and subtract
+ */
+
+#include <math.h>
+#include "libm.h"
+
+static const float Zero[] = {0.0, -0.0,};
+
+float fmodf(float x, float y)
+{
+ int32_t n,hx,hy,hz,ix,iy,sx,i;
+
+ GET_FLOAT_WORD(hx, x);
+ GET_FLOAT_WORD(hy, y);
+ sx = hx & 0x80000000; /* sign of x */
+ hx ^= sx; /* |x| */
+ hy &= 0x7fffffff; /* |y| */
+
+ /* purge off exception values */
+ if (hy == 0 || hx >= 0x7f800000 || /* y=0,or x not finite */
+ hy > 0x7f800000) /* or y is NaN */
+ return (x*y)/(x*y);
+ if (hx < hy) /* |x| < |y| */
+ return x;
+ if (hx == hy) /* |x| = |y|, return x*0 */
+ return Zero[(uint32_t)sx>>31];
+
+ /* 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 = hx<<n;
+ }
+ if (iy >= -126)
+ hy = 0x00800000|(0x007fffff&hy);
+ else { /* subnormal y, shift y to normal */
+ n = -126-iy;
+ hy = hy<<n;
+ }
+
+ /* fix point fmod */
+ n = ix - iy;
+ while (n--) {
+ hz = hx-hy;
+ if (hz<0)
+ hx = hx+hx;
+ else {
+ if(hz == 0) /* return sign(x)*0 */
+ return Zero[(uint32_t)sx>>31];
+ hx = hz+hz;
+ }
+ }
+ hz = hx-hy;
+ if (hz >= 0)
+ hx = hz;
+
+ /* convert back to floating value and restore the sign */
+ if (hx == 0) /* return sign(x)*0 */
+ return Zero[(uint32_t)sx>>31];
+ while (hx < 0x00800000) { /* normalize x */
+ hx = hx+hx;
+ iy -= 1;
+ }
+ if (iy >= -126) { /* normalize output */
+ hx = ((hx-0x00800000)|((iy+127)<<23));
+ SET_FLOAT_WORD(x, hx|sx);
+ } else { /* subnormal output */
+ n = -126 - iy;
+ hx >>= n;
+ SET_FLOAT_WORD(x, hx|sx);
+ }
+ return x; /* exact output */
+}
--- /dev/null
+/* From MUSL */
+
+#include <math.h>
+#include <stdint.h>
+
+double frexp(double x, int *e)
+{
+ union { double d; uint64_t i; } y = { x };
+ int ee = y.i>>52 & 0x7ff;
+
+ if (!ee) {
+ if (x) {
+ x = frexp(x*0x1p64, e);
+ *e -= 64;
+ } else *e = 0;
+ return x;
+ } else if (ee == 0x7ff) {
+ return x;
+ }
+
+ *e = ee - 0x3fe;
+ y.i &= 0x800fffffffffffffull;
+ y.i |= 0x3fe0000000000000ull;
+ return y.d;
+}
--- /dev/null
+/* From MUSL */
+
+#include <math.h>
+#include <stdint.h>
+
+float frexpf(float x, int *e)
+{
+ union { float f; uint32_t i; } y = { x };
+ int ee = y.i>>23 & 0xff;
+
+ if (!ee) {
+ if (x) {
+ x = frexpf(x*0x1p64, e);
+ *e -= 64;
+ } else *e = 0;
+ return x;
+ } else if (ee == 0xff) {
+ return x;
+ }
+
+ *e = ee - 0x7e;
+ y.i &= 0x807ffffful;
+ y.i |= 0x3f000000ul;
+ return y.f;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_hypot.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* hypot(x,y)
+ *
+ * Method :
+ * If (assume round-to-nearest) z=x*x+y*y
+ * has error less than sqrt(2)/2 ulp, then
+ * sqrt(z) has error less than 1 ulp (exercise).
+ *
+ * So, compute sqrt(x*x+y*y) with some care as
+ * follows to get the error below 1 ulp:
+ *
+ * Assume x>y>0;
+ * (if possible, set rounding to round-to-nearest)
+ * 1. if x > 2y use
+ * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
+ * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
+ * 2. if x <= 2y use
+ * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
+ * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
+ * y1= y with lower 32 bits chopped, y2 = y-y1.
+ *
+ * NOTE: scaling may be necessary if some argument is too
+ * large or too tiny
+ *
+ * Special cases:
+ * hypot(x,y) is INF if x or y is +INF or -INF; else
+ * hypot(x,y) is NAN if x or y is NAN.
+ *
+ * Accuracy:
+ * hypot(x,y) returns sqrt(x^2+y^2) with error less
+ * than 1 ulps (units in the last place)
+ */
+
+#include <math.h>
+#include "libm.h"
+
+double hypot(double x, double y)
+{
+ double a,b,t1,t2,y1,y2,w;
+ int32_t j,k,ha,hb;
+
+ GET_HIGH_WORD(ha, x);
+ ha &= 0x7fffffff;
+ GET_HIGH_WORD(hb, y);
+ hb &= 0x7fffffff;
+ if (hb > ha) {
+ a = y;
+ b = x;
+ j=ha; ha=hb; hb=j;
+ } else {
+ a = x;
+ b = y;
+ }
+ a = fabs(a);
+ b = fabs(b);
+ if (ha - hb > 0x3c00000) /* x/y > 2**60 */
+ return a+b;
+ k = 0;
+ if (ha > 0x5f300000) { /* a > 2**500 */
+ if(ha >= 0x7ff00000) { /* Inf or NaN */
+ uint32_t low;
+ /* Use original arg order iff result is NaN; quieten sNaNs. */
+ w = fabs(x+0.0) - fabs(y+0.0);
+ GET_LOW_WORD(low, a);
+ if (((ha&0xfffff)|low) == 0) w = a;
+ GET_LOW_WORD(low, b);
+ if (((hb^0x7ff00000)|low) == 0) w = b;
+ return w;
+ }
+ /* scale a and b by 2**-600 */
+ ha -= 0x25800000; hb -= 0x25800000; k += 600;
+ SET_HIGH_WORD(a, ha);
+ SET_HIGH_WORD(b, hb);
+ }
+ if (hb < 0x20b00000) { /* b < 2**-500 */
+ if (hb <= 0x000fffff) { /* subnormal b or 0 */
+ uint32_t low;
+ GET_LOW_WORD(low, b);
+ if ((hb|low) == 0)
+ return a;
+ t1 = 0;
+ SET_HIGH_WORD(t1, 0x7fd00000); /* t1 = 2^1022 */
+ b *= t1;
+ a *= t1;
+ k -= 1022;
+ } else { /* scale a and b by 2^600 */
+ ha += 0x25800000; /* a *= 2^600 */
+ hb += 0x25800000; /* b *= 2^600 */
+ k -= 600;
+ SET_HIGH_WORD(a, ha);
+ SET_HIGH_WORD(b, hb);
+ }
+ }
+ /* medium size a and b */
+ w = a - b;
+ if (w > b) {
+ t1 = 0;
+ SET_HIGH_WORD(t1, ha);
+ t2 = a-t1;
+ w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
+ } else {
+ a = a + a;
+ y1 = 0;
+ SET_HIGH_WORD(y1, hb);
+ y2 = b - y1;
+ t1 = 0;
+ SET_HIGH_WORD(t1, ha+0x00100000);
+ t2 = a - t1;
+ w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+ }
+ if (k)
+ w = scalbn(w, k);
+ return w;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_hypotf.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 <math.h>
+#include "libm.h"
+
+float hypotf(float x, float y)
+{
+ float a,b,t1,t2,y1,y2,w;
+ int32_t j,k,ha,hb;
+
+ GET_FLOAT_WORD(ha,x);
+ ha &= 0x7fffffff;
+ GET_FLOAT_WORD(hb,y);
+ hb &= 0x7fffffff;
+ if (hb > ha) {
+ a = y;
+ b = x;
+ j=ha; ha=hb; hb=j;
+ } else {
+ a = x;
+ b = y;
+ }
+ a = fabsf(a);
+ b = fabsf(b);
+ if (ha - hb > 0xf000000) /* x/y > 2**30 */
+ return a+b;
+ k = 0;
+ if (ha > 0x58800000) { /* a > 2**50 */
+ if(ha >= 0x7f800000) { /* Inf or NaN */
+ /* Use original arg order iff result is NaN; quieten sNaNs. */
+ w = fabsf(x+0.0f) - fabsf(y+0.0f);
+ if (ha == 0x7f800000) w = a;
+ if (hb == 0x7f800000) w = b;
+ return w;
+ }
+ /* scale a and b by 2**-68 */
+ ha -= 0x22000000; hb -= 0x22000000; k += 68;
+ SET_FLOAT_WORD(a, ha);
+ SET_FLOAT_WORD(b, hb);
+ }
+ if (hb < 0x26800000) { /* b < 2**-50 */
+ if (hb <= 0x007fffff) { /* subnormal b or 0 */
+ if (hb == 0)
+ return a;
+ SET_FLOAT_WORD(t1, 0x7e800000); /* t1 = 2^126 */
+ b *= t1;
+ a *= t1;
+ k -= 126;
+ } else { /* scale a and b by 2^68 */
+ ha += 0x22000000; /* a *= 2^68 */
+ hb += 0x22000000; /* b *= 2^68 */
+ k -= 68;
+ SET_FLOAT_WORD(a, ha);
+ SET_FLOAT_WORD(b, hb);
+ }
+ }
+ /* medium size a and b */
+ w = a - b;
+ if (w > b) {
+ SET_FLOAT_WORD(t1, ha&0xfffff000);
+ t2 = a - t1;
+ w = sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
+ } else {
+ a = a + a;
+ SET_FLOAT_WORD(y1, hb&0xfffff000);
+ y2 = b - y1;
+ SET_FLOAT_WORD(t1,(ha+0x00800000)&0xfffff000);
+ t2 = a - t1;
+ w = sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
+ }
+ if (k)
+ w = scalbnf(w, k);
+ return w;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/math_private.h */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+#ifndef _LIBM_H
+#define _LIBM_H
+
+#include <stdint.h>
+#include <float.h>
+#include <math.h>
+
+union fshape {
+ float value;
+ uint32_t bits;
+};
+
+
+/* Check platform actually has a real double type */
+#ifndef double
+
+union dshape {
+ double value;
+ uint64_t bits;
+};
+
+/* Get two 32 bit ints from a double. */
+#define EXTRACT_WORDS(hi,lo,d) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ (hi) = __u.bits >> 32; \
+ (lo) = (uint32_t)__u.bits; \
+} while (0)
+
+/* Get a 64 bit int from a double. */
+#define EXTRACT_WORD64(i,d) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ (i) = __u.bits; \
+} while (0)
+
+/* Get the more significant 32 bit int from a double. */
+#define GET_HIGH_WORD(i,d) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ (i) = __u.bits >> 32; \
+} while (0)
+
+/* Get the less significant 32 bit int from a double. */
+#define GET_LOW_WORD(i,d) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ (i) = (uint32_t)__u.bits; \
+} while (0)
+
+/* Set a double from two 32 bit ints. */
+#define INSERT_WORDS(d,hi,lo) \
+do { \
+ union dshape __u; \
+ __u.bits = ((uint64_t)(hi) << 32) | (uint32_t)(lo); \
+ (d) = __u.value; \
+} while (0)
+
+/* Set a double from a 64 bit int. */
+#define INSERT_WORD64(d,i) \
+do { \
+ union dshape __u; \
+ __u.bits = (i); \
+ (d) = __u.value; \
+} while (0)
+
+/* Set the more significant 32 bits of a double from an int. */
+#define SET_HIGH_WORD(d,hi) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ __u.bits &= 0xffffffff; \
+ __u.bits |= (uint64_t)(hi) << 32; \
+ (d) = __u.value; \
+} while (0)
+
+/* Set the less significant 32 bits of a double from an int. */
+#define SET_LOW_WORD(d,lo) \
+do { \
+ union dshape __u; \
+ __u.value = (d); \
+ __u.bits &= 0xffffffff00000000ull; \
+ __u.bits |= (uint32_t)(lo); \
+ (d) = __u.value; \
+} while (0)
+
+#endif
+
+/* Get a 32 bit int from a float. */
+#define GET_FLOAT_WORD(i,d) \
+do { \
+ union fshape __u; \
+ __u.value = (d); \
+ (i) = __u.bits; \
+} while (0)
+
+/* Set a float from a 32 bit int. */
+#define SET_FLOAT_WORD(d,i) \
+do { \
+ union fshape __u; \
+ __u.bits = (i); \
+ (d) = __u.value; \
+} while (0)
+
+/* fdlibm kernel functions */
+
+int __rem_pio2_large(double*,double*,int,int,int);
+
+int __rem_pio2(double,double*);
+double __sin(double,double,int);
+double __cos(double,double);
+double __tan(double,double,int);
+double __expo2(double);
+
+int __rem_pio2f(float,double*);
+float __sindf(double);
+float __cosdf(double);
+float __tandf(double,int);
+float __expo2f(float);
+
+#endif
--- /dev/null
+#include <math.h>
+
+double nan(const char *s)
+{
+ return __NaN;
+}
+
\ No newline at end of file
--- /dev/null
+#include <math.h>
+
+float nanf(const char *s)
+{
+ return __NaN;
+}
+
\ No newline at end of file
--- /dev/null
+/* From MUSL */
+
+#include <limits.h>
+#include <math.h>
+
+float scalblnf(float x, long n)
+{
+ if (n > INT_MAX)
+ n = INT_MAX;
+ else if (n < INT_MIN)
+ n = INT_MIN;
+ return scalbnf(x, n);
+}
--- /dev/null
+/* From MUSL */
+
+#include <limits.h>
+#include <math.h>
+
+double scalbln(double x, long n)
+{
+ if (n > INT_MAX)
+ n = INT_MAX;
+ else if (n < INT_MIN)
+ n = INT_MIN;
+ return scalbn(x, n);
+}
--- /dev/null
+/* From MUSL */
+
+#include <math.h>
+#include "libm.h"
+
+double scalbn(double x, int n)
+{
+ double scale;
+
+ if (n > 1023) {
+ x *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023) {
+ x *= 0x1p1023;
+ n -= 1023;
+ if (n > 1023) {
+ x = (double)(x * 0x1p1023);
+ return x;
+ }
+ }
+ } else if (n < -1022) {
+ x *= 0x1p-1022;
+ n += 1022;
+ if (n < -1022) {
+ x *= 0x1p-1022;
+ n += 1022;
+ if (n < -1022) {
+ x =(double)(x * 0x1p-1022);
+ return x;
+ }
+ }
+ }
+ INSERT_WORDS(scale, (uint32_t)(0x3ff+n)<<20, 0);
+ x = (double)(x * scale);
+ return x;
+}
--- /dev/null
+/* From MUSL */
+
+#include "libm.h"
+
+float scalbnf(float x, int n)
+{
+ float scale;
+
+ if (n > 127) {
+ x *= 0x1p127f;
+ n -= 127;
+ if (n > 127) {
+ x *= 0x1p127f;
+ n -= 127;
+ if (n > 127) {
+ x = (float)(x * 0x1p127f);
+ return x;
+ }
+ }
+ } else if (n < -126) {
+ x *= 0x1p-126f;
+ n += 126;
+ if (n < -126) {
+ x *= 0x1p-126f;
+ n += 126;
+ if (n < -126) {
+ x = (float)(x * 0x1p-126f);
+ return x;
+ }
+ }
+ }
+ SET_FLOAT_WORD(scale, (uint32_t)(0x7f+n)<<23);
+ x = (float)(x * scale);
+ return x;
+}
--- /dev/null
+/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* sqrt(x)
+ * Return correctly rounded sqrt.
+ * ------------------------------------------
+ * | Use the hardware sqrt if you have one |
+ * ------------------------------------------
+ * Method:
+ * Bit by bit method using integer arithmetic. (Slow, but portable)
+ * 1. Normalization
+ * Scale x to y in [1,4) with even powers of 2:
+ * find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
+ * sqrt(x) = 2^k * sqrt(y)
+ * 2. Bit by bit computation
+ * Let q = sqrt(y) truncated to i bit after binary point (q = 1),
+ * i 0
+ * i+1 2
+ * s = 2*q , and y = 2 * ( y - q ). (1)
+ * i i i i
+ *
+ * To compute q from q , one checks whether
+ * i+1 i
+ *
+ * -(i+1) 2
+ * (q + 2 ) <= y. (2)
+ * i
+ * -(i+1)
+ * If (2) is false, then q = q ; otherwise q = q + 2 .
+ * i+1 i i+1 i
+ *
+ * With some algebric manipulation, it is not difficult to see
+ * that (2) is equivalent to
+ * -(i+1)
+ * s + 2 <= y (3)
+ * i i
+ *
+ * The advantage of (3) is that s and y can be computed by
+ * i i
+ * the following recurrence formula:
+ * if (3) is false
+ *
+ * s = s , y = y ; (4)
+ * i+1 i i+1 i
+ *
+ * otherwise,
+ * -i -(i+1)
+ * s = s + 2 , y = y - s - 2 (5)
+ * i+1 i i+1 i i
+ *
+ * One may easily use induction to prove (4) and (5).
+ * Note. Since the left hand side of (3) contain only i+2 bits,
+ * it does not necessary to do a full (53-bit) comparison
+ * in (3).
+ * 3. Final rounding
+ * After generating the 53 bits result, we compute one more bit.
+ * Together with the remainder, we can decide whether the
+ * result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ * (it will never equal to 1/2ulp).
+ * The rounding mode can be detected by checking whether
+ * huge + tiny is equal to huge, and whether huge - tiny is
+ * equal to huge for some floating point number "huge" and "tiny".
+ *
+ * Special cases:
+ * sqrt(+-0) = +-0 ... exact
+ * sqrt(inf) = inf
+ * sqrt(-ve) = NaN ... with invalid signal
+ * sqrt(NaN) = NaN ... with invalid signal for signaling NaN
+ */
+
+#include <math.h>
+#include <signal.h>
+#include "libm.h"
+
+static const double tiny = 1.0e-300;
+
+double sqrt(double x)
+{
+ double z;
+ int32_t sign = (int)0x80000000;
+ int32_t ix0,s0,q,m,t,i;
+ uint32_t r,t1,s1,ix1,q1;
+
+ EXTRACT_WORDS(ix0, ix1, x);
+
+ /* take care of Inf and NaN */
+ if ((ix0&0x7ff00000) == 0x7ff00000) {
+ return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
+ }
+ /* take care of zero */
+ if (ix0 <= 0) {
+ if (((ix0&~sign)|ix1) == 0)
+ return x; /* sqrt(+-0) = +-0 */
+ if (ix0 < 0) {
+ raise(SIGFPE);
+ return __sNAN;
+ /* return (x-x)/(x-x); */ /* sqrt(-ve) = sNaN */
+ }
+ }
+ /* normalize x */
+ m = ix0>>20;
+ if (m == 0) { /* subnormal x */
+ while (ix0 == 0) {
+ m -= 21;
+ ix0 |= (ix1>>11);
+ ix1 <<= 21;
+ }
+ for (i=0; (ix0&0x00100000) == 0; i++)
+ ix0<<=1;
+ m -= i - 1;
+ ix0 |= ix1>>(32-i);
+ ix1 <<= i;
+ }
+ m -= 1023; /* unbias exponent */
+ ix0 = (ix0&0x000fffff)|0x00100000;
+ if (m & 1) { /* odd m, double x to make it even */
+ ix0 += ix0 + ((ix1&sign)>>31);
+ ix1 += ix1;
+ }
+ m >>= 1; /* m = [m/2] */
+
+ /* generate sqrt(x) bit by bit */
+ ix0 += ix0 + ((ix1&sign)>>31);
+ ix1 += ix1;
+ q = q1 = s0 = s1 = 0; /* [q,q1] = sqrt(x) */
+ r = 0x00200000; /* r = moving bit from right to left */
+
+ while (r != 0) {
+ t = s0 + r;
+ if (t <= ix0) {
+ s0 = t + r;
+ ix0 -= t;
+ q += r;
+ }
+ ix0 += ix0 + ((ix1&sign)>>31);
+ ix1 += ix1;
+ r >>= 1;
+ }
+
+ r = sign;
+ while (r != 0) {
+ t1 = s1 + r;
+ t = s0;
+ if (t < ix0 || (t == ix0 && t1 <= ix1)) {
+ s1 = t1 + r;
+ if ((t1&sign) == sign && (s1&sign) == 0)
+ s0++;
+ ix0 -= t;
+ if (ix1 < t1)
+ ix0--;
+ ix1 -= t1;
+ q1 += r;
+ }
+ ix0 += ix0 + ((ix1&sign)>>31);
+ ix1 += ix1;
+ r >>= 1;
+ }
+
+ /* use floating add to find out rounding direction */
+ if ((ix0|ix1) != 0) {
+ z = 1.0 - tiny; /* raise inexact flag */
+ if (z >= 1.0) {
+ z = 1.0 + tiny;
+ if (q1 == (uint32_t)0xffffffff) {
+ q1 = 0;
+ q++;
+ } else if (z > 1.0) {
+ if (q1 == (uint32_t)0xfffffffe)
+ q++;
+ q1 += 2;
+ } else
+ q1 += q1 & 1;
+ }
+ }
+ ix0 = (q>>1) + 0x3fe00000;
+ ix1 = q1>>1;
+ if (q&1)
+ ix1 |= sign;
+ ix0 += m << 20;
+ INSERT_WORDS(z, ix0, ix1);
+ return z;
+}
--- /dev/null
+/*
+ * 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 <math.h>
+#include <signal.h>
+#include "libm.h"
+
+static const float tiny = 1.0e-30;
+
+float sqrtf(float x)
+{
+ float z;
+ int32_t sign = (int)0x80000000;
+ int32_t ix,s,q,m,t,i;
+ uint32_t r;
+
+ GET_FLOAT_WORD(ix, x);
+
+ /* take care of Inf and NaN */
+ if ((ix&0x7f800000) == 0x7f800000)
+ return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
+
+ /* take care of zero */
+ if (ix <= 0) {
+ if ((ix&~sign) == 0)
+ return x; /* sqrt(+-0) = +-0 */
+ if (ix < 0) {
+ raise(SIGFPE);
+ return __sNaN;
+ }
+/* return (x-x)/(x-x); *//* sqrt(-ve) = sNaN */
+ }
+ /* normalize x */
+ m = ix>>23;
+ if (m == 0) { /* subnormal x */
+ for (i = 0; (ix&0x00800000) == 0; i++)
+ ix<<=1;
+ m -= i - 1;
+ }
+ m -= 127; /* unbias exponent */
+ ix = (ix&0x007fffff)|0x00800000;
+ if (m&1) /* odd m, double x to make it even */
+ ix += ix;
+ m >>= 1; /* m = [m/2] */
+
+ /* generate sqrt(x) bit by bit */
+ ix += ix;
+ q = s = 0; /* q = sqrt(x) */
+ r = 0x01000000; /* r = moving bit from right to left */
+
+ while (r != 0) {
+ t = s + r;
+ if (t <= ix) {
+ s = t+r;
+ ix -= t;
+ q += r;
+ }
+ ix += ix;
+ r >>= 1;
+ }
+
+ /* use floating add to find out rounding direction */
+ if (ix != 0) {
+ z = 1.0f - tiny; /* raise inexact flag */
+ if (z >= 1.0f) {
+ z = 1.0f + tiny;
+ if (z > 1.0f)
+ q += 2;
+ else
+ q += q & 1;
+ }
+ }
+ ix = (q>>1) + 0x3f000000;
+ ix += m << 23;
+ SET_FLOAT_WORD(z, ix);
+ return z;
+}
public / FUZIX.git / commitdiff