From dc4f61ea8ba2bd57e848724802074f763df6e329 Mon Sep 17 00:00:00 2001 From: Alan Cox Date: Mon, 20 Jun 2016 20:02:04 +0100 Subject: [PATCH] libc: first batch of libm work --- Library/include/math.h | 88 +++++++++++++++++- Library/libs/API.list | 121 +++++++++++++++---------- Library/libs/Makefile | 13 ++- Library/libs/README.musl | 94 +++++++++++++++++++ Library/libs/acos.c | 102 +++++++++++++++++++++ Library/libs/acosf.c | 79 ++++++++++++++++ Library/libs/acosh.c | 55 ++++++++++++ Library/libs/acoshf.c | 46 ++++++++++ Library/libs/asin.c | 109 ++++++++++++++++++++++ Library/libs/asinf.c | 64 +++++++++++++ Library/libs/asinh.c | 56 ++++++++++++ Library/libs/asinhf.c | 49 ++++++++++ Library/libs/atan.c | 121 +++++++++++++++++++++++++ Library/libs/atan2.c | 120 +++++++++++++++++++++++++ Library/libs/atan2f.c | 93 +++++++++++++++++++ Library/libs/atanf.c | 96 ++++++++++++++++++++ Library/libs/atanh.c | 59 ++++++++++++ Library/libs/atanhf.c | 43 +++++++++ Library/libs/ceil.c | 82 +++++++++++++++++ Library/libs/ceilf.c | 54 +++++++++++ Library/libs/cos.c | 76 ++++++++++++++++ Library/libs/cosf.c | 74 +++++++++++++++ Library/libs/cosh.c | 75 ++++++++++++++++ Library/libs/coshf.c | 58 ++++++++++++ Library/libs/exp.c | 136 ++++++++++++++++++++++++++++ Library/libs/expf.c | 84 +++++++++++++++++ Library/libs/fabs.c | 13 +++ Library/libs/fabsf.c | 12 +++ Library/libs/floor.c | 82 +++++++++++++++++ Library/libs/floorf.c | 64 +++++++++++++ Library/libs/fmod.c | 146 ++++++++++++++++++++++++++++++ Library/libs/fmodf.c | 105 ++++++++++++++++++++++ Library/libs/frexp.c | 25 ++++++ Library/libs/frexpf.c | 25 ++++++ Library/libs/hypot.c | 124 +++++++++++++++++++++++++ Library/libs/hypotf.c | 87 ++++++++++++++++++ Library/libs/libm.h | 137 ++++++++++++++++++++++++++++ Library/libs/nan.c | 7 ++ Library/libs/nanf.c | 7 ++ Library/libs/scalbinf.c | 13 +++ Library/libs/scalbln.c | 13 +++ Library/libs/scalbn.c | 36 ++++++++ Library/libs/scalbnf.c | 35 ++++++++ Library/libs/sqrt.c | 190 +++++++++++++++++++++++++++++++++++++++ Library/libs/sqrtf.c | 88 ++++++++++++++++++ 45 files changed, 3200 insertions(+), 56 deletions(-) create mode 100644 Library/libs/README.musl create mode 100644 Library/libs/acos.c create mode 100644 Library/libs/acosf.c create mode 100644 Library/libs/acosh.c create mode 100644 Library/libs/acoshf.c create mode 100644 Library/libs/asin.c create mode 100644 Library/libs/asinf.c create mode 100644 Library/libs/asinh.c create mode 100644 Library/libs/asinhf.c create mode 100644 Library/libs/atan.c create mode 100644 Library/libs/atan2.c create mode 100644 Library/libs/atan2f.c create mode 100644 Library/libs/atanf.c create mode 100644 Library/libs/atanh.c create mode 100644 Library/libs/atanhf.c create mode 100644 Library/libs/ceil.c create mode 100644 Library/libs/ceilf.c create mode 100644 Library/libs/cos.c create mode 100644 Library/libs/cosf.c create mode 100644 Library/libs/cosh.c create mode 100644 Library/libs/coshf.c create mode 100644 Library/libs/exp.c create mode 100644 Library/libs/expf.c create mode 100644 Library/libs/fabs.c create mode 100644 Library/libs/fabsf.c create mode 100644 Library/libs/floor.c create mode 100644 Library/libs/floorf.c create mode 100644 Library/libs/fmod.c create mode 100644 Library/libs/fmodf.c create mode 100644 Library/libs/frexp.c create mode 100644 Library/libs/frexpf.c create mode 100644 Library/libs/hypot.c create mode 100644 Library/libs/hypotf.c create mode 100644 Library/libs/libm.h create mode 100644 Library/libs/nan.c create mode 100644 Library/libs/nanf.c create mode 100644 Library/libs/scalbinf.c create mode 100644 Library/libs/scalbln.c create mode 100644 Library/libs/scalbn.c create mode 100644 Library/libs/scalbnf.c create mode 100644 Library/libs/sqrt.c create mode 100644 Library/libs/sqrtf.c diff --git a/Library/include/math.h b/Library/include/math.h index af4a8c6e..669c9698 100644 --- a/Library/include/math.h +++ b/Library/include/math.h @@ -1,9 +1,10 @@ -/* - */ - #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); @@ -27,8 +28,87 @@ extern double tanh(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 diff --git a/Library/libs/API.list b/Library/libs/API.list index 372044ee..66f4a79a 100644 --- a/Library/libs/API.list +++ b/Library/libs/API.list @@ -1,23 +1,37 @@ 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 @@ -25,21 +39,23 @@ clearerr 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 @@ -48,8 +64,10 @@ execve 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 @@ -61,8 +79,10 @@ fgetgrent 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 @@ -71,12 +91,13 @@ fputs 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 @@ -105,7 +126,9 @@ gmtime CODE gsignal OBSOLETE hcreate hdestroy -hypot +hsearch +hypot CODE +hypotf CODE ioctl CODE isalnum CODE isalpha CODE @@ -119,22 +142,23 @@ ispunct 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 @@ -149,20 +173,21 @@ memcpy CODE 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 @@ -176,27 +201,27 @@ read 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 @@ -214,11 +239,11 @@ strspn CODE strtod CODE strtok CODE strtol CODE -swab *MISSING* +swab CODE sync CODE system CODE -tan -tanh +tan ULY +tanh ULY tdelete tempnam tfind @@ -248,6 +273,6 @@ vfprintf CODE vsprintf CODE wait CODE write CODE -y0 -y1 -yn +y0 ULY (in j0) +y1 ULY (in j1) +yn ULY (in jn) diff --git a/Library/libs/Makefile b/Library/libs/Makefile index d5e55e33..667632eb 100644 --- a/Library/libs/Makefile +++ b/Library/libs/Makefile @@ -21,11 +21,11 @@ SRC_CRT0NS = crt0nostdio$(PLATFORM).s 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 @@ -36,8 +36,8 @@ SRC_C += getpw.c __getpwent.c getpwnam.c getpwuid.c gets.c gettimeofday.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 @@ -68,6 +68,11 @@ SRC_HARD += regexp.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') diff --git a/Library/libs/README.musl b/Library/libs/README.musl new file mode 100644 index 00000000..c6c7008f --- /dev/null +++ b/Library/libs/README.musl @@ -0,0 +1,94 @@ +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/* diff --git a/Library/libs/acos.c b/Library/libs/acos.c new file mode 100644 index 00000000..c2db25c2 --- /dev/null +++ b/Library/libs/acos.c @@ -0,0 +1,102 @@ +/* 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 +#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); + } +} diff --git a/Library/libs/acosf.c b/Library/libs/acosf.c new file mode 100644 index 00000000..b2879e2a --- /dev/null +++ b/Library/libs/acosf.c @@ -0,0 +1,79 @@ +/* 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 +#include + +#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); + } +} diff --git a/Library/libs/acosh.c b/Library/libs/acosh.c new file mode 100644 index 00000000..3006da75 --- /dev/null +++ b/Library/libs/acosh.c @@ -0,0 +1,55 @@ +/* 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 +#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)); + } +} diff --git a/Library/libs/acoshf.c b/Library/libs/acoshf.c new file mode 100644 index 00000000..beae1dc1 --- /dev/null +++ b/Library/libs/acoshf.c @@ -0,0 +1,46 @@ +/* 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 +#include +#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)); + } +} diff --git a/Library/libs/asin.c b/Library/libs/asin.c new file mode 100644 index 00000000..d5a7f0fb --- /dev/null +++ b/Library/libs/asin.c @@ -0,0 +1,109 @@ +/* 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 +#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; +} diff --git a/Library/libs/asinf.c b/Library/libs/asinf.c new file mode 100644 index 00000000..3bd91f8a --- /dev/null +++ b/Library/libs/asinf.c @@ -0,0 +1,64 @@ +/* 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 +#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; +} diff --git a/Library/libs/asinh.c b/Library/libs/asinh.c new file mode 100644 index 00000000..705df8b6 --- /dev/null +++ b/Library/libs/asinh.c @@ -0,0 +1,56 @@ +/* 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 +#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; +} diff --git a/Library/libs/asinhf.c b/Library/libs/asinhf.c new file mode 100644 index 00000000..98b36c62 --- /dev/null +++ b/Library/libs/asinhf.c @@ -0,0 +1,49 @@ +/* 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 +#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; +} diff --git a/Library/libs/atan.c b/Library/libs/atan.c new file mode 100644 index 00000000..369df18a --- /dev/null +++ b/Library/libs/atan.c @@ -0,0 +1,121 @@ +/* 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 +#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; +} diff --git a/Library/libs/atan2.c b/Library/libs/atan2.c new file mode 100644 index 00000000..7d225ca4 --- /dev/null +++ b/Library/libs/atan2.c @@ -0,0 +1,120 @@ +/* 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 +#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(-,-) */ + } +} diff --git a/Library/libs/atan2f.c b/Library/libs/atan2f.c new file mode 100644 index 00000000..b8522c67 --- /dev/null +++ b/Library/libs/atan2f.c @@ -0,0 +1,93 @@ +/* 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 +#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(-,-) */ + } +} diff --git a/Library/libs/atanf.c b/Library/libs/atanf.c new file mode 100644 index 00000000..4bccd35f --- /dev/null +++ b/Library/libs/atanf.c @@ -0,0 +1,96 @@ +/* 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 +#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; +} diff --git a/Library/libs/atanh.c b/Library/libs/atanh.c new file mode 100644 index 00000000..0f7a3058 --- /dev/null +++ b/Library/libs/atanh.c @@ -0,0 +1,59 @@ +/* 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 +#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; +} diff --git a/Library/libs/atanhf.c b/Library/libs/atanhf.c new file mode 100644 index 00000000..a2b54fda --- /dev/null +++ b/Library/libs/atanhf.c @@ -0,0 +1,43 @@ +/* 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 +#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; +} diff --git a/Library/libs/ceil.c b/Library/libs/ceil.c new file mode 100644 index 00000000..19555180 --- /dev/null +++ b/Library/libs/ceil.c @@ -0,0 +1,82 @@ +/* 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; +} diff --git a/Library/libs/ceilf.c b/Library/libs/ceilf.c new file mode 100644 index 00000000..fec945b6 --- /dev/null +++ b/Library/libs/ceilf.c @@ -0,0 +1,54 @@ +/* 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; +} diff --git a/Library/libs/cos.c b/Library/libs/cos.c new file mode 100644 index 00000000..3dde803e --- /dev/null +++ b/Library/libs/cos.c @@ -0,0 +1,76 @@ +/* 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 +#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); + } +} diff --git a/Library/libs/cosf.c b/Library/libs/cosf.c new file mode 100644 index 00000000..a8792393 --- /dev/null +++ b/Library/libs/cosf.c @@ -0,0 +1,74 @@ +/* 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 +#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); + } +} diff --git a/Library/libs/cosh.c b/Library/libs/cosh.c new file mode 100644 index 00000000..9c3d7e98 --- /dev/null +++ b/Library/libs/cosh.c @@ -0,0 +1,75 @@ +/* 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 +#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; +} diff --git a/Library/libs/coshf.c b/Library/libs/coshf.c new file mode 100644 index 00000000..51e8042f --- /dev/null +++ b/Library/libs/coshf.c @@ -0,0 +1,58 @@ +/* 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 +#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; +} diff --git a/Library/libs/exp.c b/Library/libs/exp.c new file mode 100644 index 00000000..edcef465 --- /dev/null +++ b/Library/libs/exp.c @@ -0,0 +1,136 @@ +/* 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); +} diff --git a/Library/libs/expf.c b/Library/libs/expf.c new file mode 100644 index 00000000..e266248a --- /dev/null +++ b/Library/libs/expf.c @@ -0,0 +1,84 @@ +/* + * 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 +#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); +} diff --git a/Library/libs/fabs.c b/Library/libs/fabs.c new file mode 100644 index 00000000..18596239 --- /dev/null +++ b/Library/libs/fabs.c @@ -0,0 +1,13 @@ +/* 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; +} diff --git a/Library/libs/fabsf.c b/Library/libs/fabsf.c new file mode 100644 index 00000000..128da8f9 --- /dev/null +++ b/Library/libs/fabsf.c @@ -0,0 +1,12 @@ +/* From MUSL */ + +#include "libm.h" + +float fabsf(float x) +{ + union fshape u; + + u.value = x; + u.bits &= (uint32_t)-1 / 2; + return u.value; +} diff --git a/Library/libs/floor.c b/Library/libs/floor.c new file mode 100644 index 00000000..ecb9dde8 --- /dev/null +++ b/Library/libs/floor.c @@ -0,0 +1,82 @@ +/* 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; +} diff --git a/Library/libs/floorf.c b/Library/libs/floorf.c new file mode 100644 index 00000000..23fa5e7d --- /dev/null +++ b/Library/libs/floorf.c @@ -0,0 +1,64 @@ +/* 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; +} diff --git a/Library/libs/fmod.c b/Library/libs/fmod.c new file mode 100644 index 00000000..44d24dd6 --- /dev/null +++ b/Library/libs/fmod.c @@ -0,0 +1,146 @@ +/* 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 +#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<>(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<>(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 */ +} diff --git a/Library/libs/fmodf.c b/Library/libs/fmodf.c new file mode 100644 index 00000000..52b97790 --- /dev/null +++ b/Library/libs/fmodf.c @@ -0,0 +1,105 @@ +/* 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 +#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<= -126) + hy = 0x00800000|(0x007fffff&hy); + else { /* subnormal y, shift y to normal */ + n = -126-iy; + hy = hy<>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 */ +} diff --git a/Library/libs/frexp.c b/Library/libs/frexp.c new file mode 100644 index 00000000..509e92f9 --- /dev/null +++ b/Library/libs/frexp.c @@ -0,0 +1,25 @@ +/* From MUSL */ + +#include +#include + +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; +} diff --git a/Library/libs/frexpf.c b/Library/libs/frexpf.c new file mode 100644 index 00000000..35e97b71 --- /dev/null +++ b/Library/libs/frexpf.c @@ -0,0 +1,25 @@ +/* From MUSL */ + +#include +#include + +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; +} diff --git a/Library/libs/hypot.c b/Library/libs/hypot.c new file mode 100644 index 00000000..573ce223 --- /dev/null +++ b/Library/libs/hypot.c @@ -0,0 +1,124 @@ +/* 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 +#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; +} diff --git a/Library/libs/hypotf.c b/Library/libs/hypotf.c new file mode 100644 index 00000000..00ea7d4b --- /dev/null +++ b/Library/libs/hypotf.c @@ -0,0 +1,87 @@ +/* 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 +#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; +} diff --git a/Library/libs/libm.h b/Library/libs/libm.h new file mode 100644 index 00000000..b42eed52 --- /dev/null +++ b/Library/libs/libm.h @@ -0,0 +1,137 @@ +/* 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 +#include +#include + +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 diff --git a/Library/libs/nan.c b/Library/libs/nan.c new file mode 100644 index 00000000..525ff160 --- /dev/null +++ b/Library/libs/nan.c @@ -0,0 +1,7 @@ +#include + +double nan(const char *s) +{ + return __NaN; +} + \ No newline at end of file diff --git a/Library/libs/nanf.c b/Library/libs/nanf.c new file mode 100644 index 00000000..78536f8f --- /dev/null +++ b/Library/libs/nanf.c @@ -0,0 +1,7 @@ +#include + +float nanf(const char *s) +{ + return __NaN; +} + \ No newline at end of file diff --git a/Library/libs/scalbinf.c b/Library/libs/scalbinf.c new file mode 100644 index 00000000..8c717c21 --- /dev/null +++ b/Library/libs/scalbinf.c @@ -0,0 +1,13 @@ +/* From MUSL */ + +#include +#include + +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); +} diff --git a/Library/libs/scalbln.c b/Library/libs/scalbln.c new file mode 100644 index 00000000..350c91f3 --- /dev/null +++ b/Library/libs/scalbln.c @@ -0,0 +1,13 @@ +/* From MUSL */ + +#include +#include + +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); +} diff --git a/Library/libs/scalbn.c b/Library/libs/scalbn.c new file mode 100644 index 00000000..a20be3c7 --- /dev/null +++ b/Library/libs/scalbn.c @@ -0,0 +1,36 @@ +/* From MUSL */ + +#include +#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; +} diff --git a/Library/libs/scalbnf.c b/Library/libs/scalbnf.c new file mode 100644 index 00000000..70e10e21 --- /dev/null +++ b/Library/libs/scalbnf.c @@ -0,0 +1,35 @@ +/* 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; +} diff --git a/Library/libs/sqrt.c b/Library/libs/sqrt.c new file mode 100644 index 00000000..bf62cfc8 --- /dev/null +++ b/Library/libs/sqrt.c @@ -0,0 +1,190 @@ +/* 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 +#include +#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; +} diff --git a/Library/libs/sqrtf.c b/Library/libs/sqrtf.c new file mode 100644 index 00000000..3d485bca --- /dev/null +++ b/Library/libs/sqrtf.c @@ -0,0 +1,88 @@ +/* + * 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 +#include +#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; +} -- 2.34.1