diff --git a/runtime/src/main/cpp/math/COPYRIGHT b/runtime/src/main/cpp/math/COPYRIGHT new file mode 100644 index 00000000000..a8a38a56305 --- /dev/null +++ b/runtime/src/main/cpp/math/COPYRIGHT @@ -0,0 +1,128 @@ +musl as a whole is licensed under the following standard MIT license: + +---------------------------------------------------------------------- +Copyright © 2005-2014 Rich Felker, et al. + +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. +---------------------------------------------------------------------- + +Authors/contributors include: + +Anthony G. Basile +Arvid Picciani +Bobby Bingham +Boris Brezillon +Chris Spiegel +Emil Renner Berthing +Hiltjo Posthuma +Isaac Dunham +Jens Gustedt +Jeremy Huntwork +John Spencer +Justin Cormack +Luca Barbato +Luka Perkov +Michael Forney +Nicholas J. Kain +orc +Pascal Cuoq +Pierre Carrier +Rich Felker +Richard Pennington +Solar Designer +Strake +Szabolcs Nagy +Timo Teräs +Valentin Ochs +William Haddon + +Portions of this software are derived from third-party works licensed +under terms compatible with the above MIT license: + +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. + +Much 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 in comments in the individual source files. All +have been licensed under extremely permissive terms. + +The ARM memcpy code (src/string/armel/memcpy.s) is Copyright © 2008 +The Android Open Source Project and is licensed under a two-clause BSD +license. It was taken from Bionic libc, used on Android. + +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 +produced specifically for use as part of this library, written either +by Rich Felker, the main author of the library, or by one or more +contibutors listed above. Details on authorship of individual files +can be found in the git version control history of the project. The +omission of copyright and license comments in each file is in the +interest of source tree size. + +All public header files (include/* and arch/*/bits/*) 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 the authors and +request an explicit statement releasing it from copyright. + +The following files are trivial, believed not to be copyrightable in +the first place, and hereby explicitly released to the Public Domain: + +All public headers: include/*, arch/*/bits/* +Startup files: crt/* diff --git a/runtime/src/main/cpp/math/endian.h b/runtime/src/main/cpp/math/endian.h new file mode 100644 index 00000000000..1bd444518ab --- /dev/null +++ b/runtime/src/main/cpp/math/endian.h @@ -0,0 +1,82 @@ +#ifndef _ENDIAN_H +#define _ENDIAN_H + +#include + +#define __LITTLE_ENDIAN 1234 +#define __BIG_ENDIAN 4321 +#define __PDP_ENDIAN 3412 + +#if defined(__GNUC__) && defined(__BYTE_ORDER__) +#define __BYTE_ORDER __BYTE_ORDER__ +#else +#include +#endif + +#if defined(_GNU_SOURCE) || defined(_BSD_SOURCE) + +#define BIG_ENDIAN __BIG_ENDIAN +#define LITTLE_ENDIAN __LITTLE_ENDIAN +#define PDP_ENDIAN __PDP_ENDIAN +#define BYTE_ORDER __BYTE_ORDER + +#include + +static __inline uint16_t __bswap16(uint16_t __x) +{ + return __x<<8 | __x>>8; +} + +static __inline uint32_t __bswap32(uint32_t __x) +{ + return __x>>24 | __x>>8&0xff00 | __x<<8&0xff0000 | __x<<24; +} + +static __inline uint64_t __bswap64(uint64_t __x) +{ + return __bswap32(__x)+0ULL<<32 | __bswap32(__x>>32); +} + +#if __BYTE_ORDER == __LITTLE_ENDIAN +#define htobe16(x) __bswap16(x) +#define be16toh(x) __bswap16(x) +#define betoh16(x) __bswap16(x) +#define htobe32(x) __bswap32(x) +#define be32toh(x) __bswap32(x) +#define betoh32(x) __bswap32(x) +#define htobe64(x) __bswap64(x) +#define be64toh(x) __bswap64(x) +#define betoh64(x) __bswap64(x) +#define htole16(x) (uint16_t)(x) +#define le16toh(x) (uint16_t)(x) +#define letoh16(x) (uint16_t)(x) +#define htole32(x) (uint32_t)(x) +#define le32toh(x) (uint32_t)(x) +#define letoh32(x) (uint32_t)(x) +#define htole64(x) (uint64_t)(x) +#define le64toh(x) (uint64_t)(x) +#define letoh64(x) (uint64_t)(x) +#else +#define htobe16(x) (uint16_t)(x) +#define be16toh(x) (uint16_t)(x) +#define betoh16(x) (uint16_t)(x) +#define htobe32(x) (uint32_t)(x) +#define be32toh(x) (uint32_t)(x) +#define betoh32(x) (uint32_t)(x) +#define htobe64(x) (uint64_t)(x) +#define be64toh(x) (uint64_t)(x) +#define betoh64(x) (uint64_t)(x) +#define htole16(x) __bswap16(x) +#define le16toh(x) __bswap16(x) +#define letoh16(x) __bswap16(x) +#define htole32(x) __bswap32(x) +#define le32toh(x) __bswap32(x) +#define letoh32(x) __bswap32(x) +#define htole64(x) __bswap64(x) +#define le64toh(x) __bswap64(x) +#define letoh64(x) __bswap64(x) +#endif + +#endif + +#endif diff --git a/runtime/src/main/cpp/math/fmod.cpp b/runtime/src/main/cpp/math/fmod.cpp new file mode 100644 index 00000000000..6849722bac5 --- /dev/null +++ b/runtime/src/main/cpp/math/fmod.cpp @@ -0,0 +1,68 @@ +#include +#include + +double fmod(double x, double y) +{ + union {double f; uint64_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>52 & 0x7ff; + int ey = uy.i>>52 & 0x7ff; + int sx = ux.i>>63; + uint64_t i; + + /* in the followings uxi should be ux.i, but then gcc wrongly adds */ + /* float load/store to inner loops ruining performance and code size */ + uint64_t uxi = ux.i; + + if (uy.i<<1 == 0 || isnan(y) || ex == 0x7ff) + return (x*y)/(x*y); + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; + return x; + } + + /* normalize x and y */ + if (!ex) { + for (i = uxi<<12; i>>63 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1ULL >> 12; + uxi |= 1ULL << 52; + } + if (!ey) { + for (i = uy.i<<12; i>>63 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1ULL >> 12; + uy.i |= 1ULL << 52; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 63 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + uxi <<= 1; + } + i = uxi - uy.i; + if (i >> 63 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + for (; uxi>>52 == 0; uxi <<= 1, ex--); + + /* scale result */ + if (ex > 0) { + uxi -= 1ULL << 52; + uxi |= (uint64_t)ex << 52; + } else { + uxi >>= -ex + 1; + } + uxi |= (uint64_t)sx << 63; + ux.i = uxi; + return ux.f; +} diff --git a/runtime/src/main/cpp/math/fmodf.cpp b/runtime/src/main/cpp/math/fmodf.cpp new file mode 100644 index 00000000000..ff58f93365b --- /dev/null +++ b/runtime/src/main/cpp/math/fmodf.cpp @@ -0,0 +1,65 @@ +#include +#include + +float fmodf(float x, float y) +{ + union {float f; uint32_t i;} ux = {x}, uy = {y}; + int ex = ux.i>>23 & 0xff; + int ey = uy.i>>23 & 0xff; + uint32_t sx = ux.i & 0x80000000; + uint32_t i; + uint32_t uxi = ux.i; + + if (uy.i<<1 == 0 || isnan(y) || ex == 0xff) + return (x*y)/(x*y); + if (uxi<<1 <= uy.i<<1) { + if (uxi<<1 == uy.i<<1) + return 0*x; + return x; + } + + /* normalize x and y */ + if (!ex) { + for (i = uxi<<9; i>>31 == 0; ex--, i <<= 1); + uxi <<= -ex + 1; + } else { + uxi &= -1U >> 9; + uxi |= 1U << 23; + } + if (!ey) { + for (i = uy.i<<9; i>>31 == 0; ey--, i <<= 1); + uy.i <<= -ey + 1; + } else { + uy.i &= -1U >> 9; + uy.i |= 1U << 23; + } + + /* x mod y */ + for (; ex > ey; ex--) { + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + uxi <<= 1; + } + i = uxi - uy.i; + if (i >> 31 == 0) { + if (i == 0) + return 0*x; + uxi = i; + } + for (; uxi>>23 == 0; uxi <<= 1, ex--); + + /* scale result up */ + if (ex > 0) { + uxi -= 1U << 23; + uxi |= (uint32_t)ex << 23; + } else { + uxi >>= -ex + 1; + } + uxi |= sx; + ux.i = uxi; + return ux.f; +} diff --git a/runtime/src/main/cpp/math/libm.h b/runtime/src/main/cpp/math/libm.h new file mode 100644 index 00000000000..bef69d8cf52 --- /dev/null +++ b/runtime/src/main/cpp/math/libm.h @@ -0,0 +1,163 @@ +/* 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 +//#include "endian.h" + +#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 +#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t m; + uint16_t se; + } i; +}; +#elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 && __BYTE_ORDER == __LITTLE_ENDIAN +union ldshape { + long double f; + struct { + uint64_t lo; + uint32_t mid; + uint16_t top; + uint16_t se; + } i; + struct { + uint64_t lo; + uint64_t hi; + } i2; +}; +#else +#error Unsupported long double representation +#endif + +#ifdef __EMSCRIPTEN__ +/* + * asm.js doesn't have user-accessible floating point exceptions, so there's + * no point in trying to force expression evaluations to produce them. + */ +#define FORCE_EVAL(x) +#else +#define FORCE_EVAL(x) do { \ + if (sizeof(x) == sizeof(float)) { \ + volatile float __x; \ + __x = (x); \ + } else if (sizeof(x) == sizeof(double)) { \ + volatile double __x; \ + __x = (x); \ + } else { \ + volatile long double __x; \ + __x = (x); \ + } \ +} while(0) +#endif + +/* Get two 32 bit ints from a double. */ +#define EXTRACT_WORDS(hi,lo,d) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.f = (d); \ + (hi) = __u.i >> 32; \ + (lo) = (uint32_t)__u.i; \ +} while (0) + +/* Get the more significant 32 bit int from a double. */ +#define GET_HIGH_WORD(hi,d) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.f = (d); \ + (hi) = __u.i >> 32; \ +} while (0) + +/* Get the less significant 32 bit int from a double. */ +#define GET_LOW_WORD(lo,d) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.f = (d); \ + (lo) = (uint32_t)__u.i; \ +} while (0) + +/* Set a double from two 32 bit ints. */ +#define INSERT_WORDS(d,hi,lo) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.i = ((uint64_t)(hi)<<32) | (uint32_t)(lo); \ + (d) = __u.f; \ +} while (0) + +/* Set the more significant 32 bits of a double from an int. */ +#define SET_HIGH_WORD(d,hi) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.f = (d); \ + __u.i &= 0xffffffff; \ + __u.i |= (uint64_t)(hi) << 32; \ + (d) = __u.f; \ +} while (0) + +/* Set the less significant 32 bits of a double from an int. */ +#define SET_LOW_WORD(d,lo) \ +do { \ + union {double f; uint64_t i;} __u; \ + __u.f = (d); \ + __u.i &= 0xffffffff00000000ull; \ + __u.i |= (uint32_t)(lo); \ + (d) = __u.f; \ +} while (0) + +/* Get a 32 bit int from a float. */ +#define GET_FLOAT_WORD(w,d) \ +do { \ + union {float f; uint32_t i;} __u; \ + __u.f = (d); \ + (w) = __u.i; \ +} while (0) + +/* Set a float from a 32 bit int. */ +#define SET_FLOAT_WORD(d,w) \ +do { \ + union {float f; uint32_t i;} __u; \ + __u.i = (w); \ + (d) = __u.f; \ +} 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); + +int __rem_pio2l(long double, long double *); +long double __sinl(long double, long double, int); +long double __cosl(long double, long double); +long double __tanl(long double, long double, int); + +/* polynomial evaluation */ +long double __polevll(long double, const long double *, int); +long double __p1evll(long double, const long double *, int); + +#endif diff --git a/runtime/src/main/cpp/math/scalbn.cpp b/runtime/src/main/cpp/math/scalbn.cpp new file mode 100644 index 00000000000..530e07c79f0 --- /dev/null +++ b/runtime/src/main/cpp/math/scalbn.cpp @@ -0,0 +1,31 @@ +#include +#include + +double scalbn(double x, int n) +{ + union {double f; uint64_t i;} u; + double_t y = x; + + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) { + y *= 0x1p1023; + n -= 1023; + if (n > 1023) + n = 1023; + } + } else if (n < -1022) { + y *= 0x1p-1022; + n += 1022; + if (n < -1022) { + y *= 0x1p-1022; + n += 1022; + if (n < -1022) + n = -1022; + } + } + u.i = (uint64_t)(0x3ff+n)<<52; + x = y * u.f; + return x; +}