Brought some musl math functions for wasm32.

This commit is contained in:
Alexander Gorshenev
2017-08-03 17:32:56 +03:00
committed by alexander-gorshenev
parent c5d78b5b54
commit 9ec5b24832
6 changed files with 537 additions and 0 deletions
+128
View File
@@ -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/*
+82
View File
@@ -0,0 +1,82 @@
#ifndef _ENDIAN_H
#define _ENDIAN_H
#include <features.h>
#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 <bits/endian.h>
#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 <stdint.h>
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
+68
View File
@@ -0,0 +1,68 @@
#include <math.h>
#include <stdint.h>
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;
}
+65
View File
@@ -0,0 +1,65 @@
#include <math.h>
#include <stdint.h>
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;
}
+163
View File
@@ -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 <stdint.h>
#include <float.h>
#include <math.h>
//#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
+31
View File
@@ -0,0 +1,31 @@
#include <math.h>
#include <stdint.h>
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;
}