UpsilonNumworks/Upsilon
1
0
Fork 0
mirror of https://github.com/UpsilonNumworks/Upsilon.git synced 2026-03-18 21:30:38 +01:00
Code Issues Packages Projects Releases Wiki Activity

[liba] Add extra libm symbols (ldexp, fmod, trunc, atan2)

Browse Source
This commit is contained in:
Romain Goyet
2017-11-08 17:46:46 +01:00
parent 3adc5d9fac
commit 0ca52e201f
14 changed files with 558 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
liba/Makefile
View File

@@ -12,11 +12,13 @@ objs += $(addprefix liba/src/, \
fpclassify.o \
fpclassifyf.o \
ieee754.o \
ldexp.o \
malloc.o \
memcmp.o \
memcpy.o \
memmove.o \
memset.o \
nearbyint.o \
nearbyintf.o \
strcmp.o \
strchr.o \
@@ -30,9 +32,11 @@ objs += $(addprefix liba/src/external/openbsd/, \
e_acoshf.o \
e_asinf.o \
e_atanhf.o \
e_atan2.o \
e_atan2f.o \
e_coshf.o \
e_expf.o \
e_fmod.o \
e_fmodf.o \
e_lgammaf_r.o \
e_log10f.o \
@@ -54,7 +58,9 @@ objs += $(addprefix liba/src/external/openbsd/, \
s_fabsf.o \
s_floorf.o \
s_frexpf.o \
s_frexp.o \
s_log1pf.o \
s_modf.o \
s_modff.o \
s_roundf.o \
s_scalbnf.o \
@@ -62,6 +68,7 @@ objs += $(addprefix liba/src/external/openbsd/, \
s_sinf.o \
s_tanf.o \
s_tanhf.o \
s_trunc.o \
s_truncf.o \
w_lgammaf.o \
)
@@ -172,4 +179,5 @@ objs += $(addprefix liba/src/external/softfloat/src/, \
s_shortShiftRightJam64.o \
s_subMagsF64.o \
softfloat_state.o \
ui32_to_f64.o \
)

16
liba/include/math.h
View File

@@ -36,6 +36,7 @@ typedef double double_t;
#define FP_ZERO 0x10
#define fpclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)
#define signbit(x) __builtin_signbit(x)
#define isfinite(x) __builtin_isfinite(x)
#define isnormal(x) __builtin_isnormal(x)
#define isnan(x) __builtin_isnan(x)
@@ -65,7 +66,6 @@ float log1pf(float x);
float log10f(float x);
float logf(float x);
float modff(float value, float *iptr);
float nanf(const char *s);
float nearbyintf(float x);
float powf(float x, float y);
float roundf(float x);
@@ -82,6 +82,7 @@ double acosh(double x);
double asin(double x);
double asinh(double x);
double atan(double x);
double atan2(double y, double x);
double atanh(double x);
double ceil(double x);
double copysign(double x, double y);
@@ -91,12 +92,17 @@ double exp(double x);
double expm1(double x);
double fabs(double x);
double floor(double x);
double fmod(double x, double y);
double frexp(double x, int *eptr);
double lgamma(double x);
double lgamma_r(double x, int *signgamp);
double log1p(double x);
double log10(double x);
double log(double x);
double modf(double value, double *iptr);
double nearbyint(double x);
double pow(double x, double y);
double nearbyint(double x);
double round(double x);
double scalbn(double x, int n);
double sin(double x);
@@ -104,6 +110,7 @@ double sinh(double x);
double sqrt(double x);
double tan(double x);
double tanh(double x);
double trunc(double x);
/* The C99 standard says that any libc function can be re-declared as a macro.
* (See N1124 paragraph 7.1.4). This means that C files willing to actually
@@ -125,8 +132,8 @@ double tanh(double x);
#define expm1f(x) __builtin_expm1f(x)
#define fabsf(x) __builtin_fabsf(x)
#define floorf(x) __builtin_floorf(x)
#define frexpf(x, y) __builtin_frexpf(x, y)
#define fmodf(x, y) __builtin_fmodf(x, y)
#define frexpf(x, y) __builtin_frexpf(x, y)
#define ldexpf(x, n) __builtin_ldexpf(x, n)
#define lgammaf(x) __builtin_lgammaf(x)
#define lgammaf_r(x, signgamp) __builtin_lgammaf_r(x, signgamp)
@@ -150,6 +157,7 @@ double tanh(double x);
#define asin(x) __builtin_asin(x)
#define asinh(x) __builtin_asinh(x)
#define atan(x) __builtin_atan(x)
#define atan2(y, x) __builtin_atan2(y, x)
#define atanh(x) __builtin_atanh(x)
#define ceil(x) __builtin_ceil(x)
#define copysign(x, y) __builtin_copysign(x, y)
@@ -159,11 +167,14 @@ double tanh(double x);
#define expm1(x) __builtin_expm1(x)
#define fabs(x) __builtin_fabs(x)
#define floor(x) __builtin_floor(x)
#define fmod(x, y) __builtin_fmod(x, y)
#define ldexp(x, n) __builtin_ldexp(x, n)
#define lgamma(x) __builtin_lgamma(x)
#define lgamma_r(x, signgamp) __builtin_lgamma_r(x, signgamp)
#define log1p(x) __builtin_log1p(x)
#define log10(x) __builtin_log10(x)
#define log(x) __builtin_log(x)
#define nan(s) __builtin_nan(s)
#define pow(x, y) __builtin_pow(x, y)
#define round(x) __builtin_round(x)
#define scalbn(x, n) __builtin_scalbn(x, n)
@@ -172,6 +183,7 @@ double tanh(double x);
#define sqrt(x) __builtin_sqrt(x)
#define tan(x) __builtin_tan(x)
#define tanh(x) __builtin_tanh(x)
#define trunc(x) __builtin_trunc(x)
extern int signgam;

4
liba/src/aeabi-rt/double.c
View File

@@ -12,6 +12,10 @@ aeabi_double_t __aeabi_i2d(int i) {
return d(i32_to_f64(i));
}
aeabi_double_t __aeabi_ui2d(unsigned int i) {
return d(ui32_to_f64(i));
}
aeabi_float_t __aeabi_d2f(aeabi_double_t d) {
return f(f64_to_f32(f64(d)));
}

127
liba/src/external/openbsd/e_atan2.c vendored Normal file
View File

@@ -0,0 +1,127 @@
/* @(#)e_atan2.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* 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 <sys/cdefs.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
static const double
tiny = 1.0e-300,
zero = 0.0,
pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
double
atan2(double y, double x)
{
double z;
int32_t k,m,hx,hy,ix,iy;
u_int32_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) return atan(y); /* x=1.0 */
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 zero ; /* atan(+...,+INF) */
case 1: return -zero ; /* 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) z=pi_o_2+0.5*pi_lo; /* |y/x| > 2**60 */
else if(hx<0&&k<-60) z=0.0; /* |y|/x < -2**60 */
else z=atan(fabs(y/x)); /* safe to do y/x */
switch (m) {
case 0: return z ; /* atan(+,+) */
case 1: {
u_int32_t zh;
GET_HIGH_WORD(zh,z);
SET_HIGH_WORD(z,zh ^ 0x80000000);
}
return z ; /* atan(-,+) */
case 2: return pi-(z-pi_lo);/* atan(+,-) */
default: /* case 3 */
return (z-pi_lo)-pi;/* atan(-,-) */
}
}
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(atan2l, atan2);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

128
liba/src/external/openbsd/e_fmod.c vendored Normal file
View File

@@ -0,0 +1,128 @@
/* @(#)e_fmod.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* fmod(x,y)
* Return x mod y in exact arithmetic
* Method: shift and subtract
*/
#include "math.h"
#include "math_private.h"
static const double one = 1.0, Zero[] = {0.0, -0.0,};
double
fmod(double x, double y)
{
int32_t n,hx,hy,hz,ix,iy,sx,i;
u_int32_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)) return x; /* |x|<|y| return x */
if(lx==ly)
return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/
}
/* determine ix = ilogb(x) */
if(hx<0x00100000) { /* subnormal x */
if(hx==0) {
for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
} else {
for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
}
} else ix = (hx>>20)-1023;
/* determine iy = ilogb(y) */
if(hy<0x00100000) { /* subnormal y */
if(hy==0) {
for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
} else {
for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
}
} else iy = (hy>>20)-1023;
/* set up {hx,lx}, {hy,ly} and align y to x */
if(ix >= -1022)
hx = 0x00100000|(0x000fffff&hx);
else { /* subnormal x, shift x to normal */
n = -1022-ix;
if(n<=31) {
hx = (hx<<n)|(lx>>(32-n));
lx <<= n;
} else {
hx = lx<<(n-32);
lx = 0;
}
}
if(iy >= -1022)
hy = 0x00100000|(0x000fffff&hy);
else { /* subnormal y, shift y to normal */
n = -1022-iy;
if(n<=31) {
hy = (hy<<n)|(ly>>(32-n));
ly <<= n;
} else {
hy = ly<<(n-32);
ly = 0;
}
}
/* fix point fmod */
n = ix - iy;
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[(u_int32_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[(u_int32_t)sx>>31];
while(hx<0x00100000) { /* normalize x */
hx = hx+hx+(lx>>31); lx = lx+lx;
iy -= 1;
}
if(iy>= -1022) { /* normalize output */
hx = ((hx-0x00100000)|((iy+1023)<<20));
INSERT_WORDS(x,hx|sx,lx);
} else { /* subnormal output */
n = -1022 - iy;
if(n<=20) {
lx = (lx>>n)|((u_int32_t)hx<<(32-n));
hx >>= n;
} else if (n<=31) {
lx = (hx<<(32-n))|(lx>>n); hx = sx;
} else {
lx = hx>>(n-32); hx = sx;
}
INSERT_WORDS(x,hx|sx,lx);
x *= one; /* create necessary signal */
}
return x; /* exact output */
}

7
liba/src/external/openbsd/include/math.h vendored
View File

@@ -26,12 +26,12 @@
#undef fabsf
#undef floorf
#undef fmodf
#undef ldexpf
#undef lgammaf
#undef lgammaf_r
#undef log1pf
#undef log10f
#undef logf
#undef nanf
#undef nearbyintf
#undef powf
#undef roundf
@@ -41,12 +41,14 @@
#undef sqrtf
#undef tanf
#undef tanhf
#undef truncf
#undef acos
#undef acosh
#undef asin
#undef asinh
#undef atan
#undef atan2
#undef atanh
#undef ceil
#undef copysign
@@ -56,6 +58,8 @@
#undef expm1
#undef fabs
#undef floor
#undef fmod
#undef ldexp
#undef lgamma
#undef lgamma_r
#undef log1p
@@ -69,3 +73,4 @@
#undef sqrt
#undef tan
#undef tanh
#undef trunc

56
liba/src/external/openbsd/s_frexp.c vendored Normal file
View File

@@ -0,0 +1,56 @@
/* @(#)s_frexp.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* for non-zero x
* x = frexp(arg,&exp);
* return a double fp quantity x such that 0.5 <= |x| <1.0
* and the corresponding binary exponent "exp". That is
* arg = x*2^exp.
* If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
* with *exp=0.
*/
#include <sys/cdefs.h>
#include <float.h>
#include <math.h>
#include "math_private.h"
static const double
two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
double
frexp(double x, int *eptr)
{
int32_t hx, ix, lx;
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
*eptr = 0;
if(ix>=0x7ff00000||((ix|lx)==0)) return x; /* 0,inf,nan */
if (ix<0x00100000) { /* subnormal */
x *= two54;
GET_HIGH_WORD(hx,x);
ix = hx&0x7fffffff;
*eptr = -54;
}
*eptr += (ix>>20)-1022;
hx = (hx&0x800fffff)|0x3fe00000;
SET_HIGH_WORD(x,hx);
return x;
}
#if LDBL_MANT_DIG == 53
#ifdef __weak_alias
__weak_alias(frexpl, frexp);
#endif /* __weak_alias */
#endif /* LDBL_MANT_DIG == 53 */

71
liba/src/external/openbsd/s_modf.c vendored Normal file
View File

@@ -0,0 +1,71 @@
/* @(#)s_modf.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* modf(double x, double *iptr)
* return fraction part of x, and return x's integral part in *iptr.
* Method:
* Bit twiddling.
*
* Exception:
* No exception.
*/
#include "math.h"
#include "math_private.h"
static const double one = 1.0;
double
modf(double x, double *iptr)
{
int32_t i0,i1,j0;
u_int32_t i;
EXTRACT_WORDS(i0,i1,x);
j0 = ((i0>>20)&0x7ff)-0x3ff; /* exponent of x */
if(j0<20) { /* integer part in high x */
if(j0<0) { /* |x|<1 */
INSERT_WORDS(*iptr,i0&0x80000000,0); /* *iptr = +-0 */
return x;
} else {
i = (0x000fffff)>>j0;
if(((i0&i)|i1)==0) { /* x is integral */
u_int32_t high;
*iptr = x;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0&(~i),0);
return x - *iptr;
}
}
} else if (j0>51) { /* no fraction part */
u_int32_t high;
*iptr = x*one;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else { /* fraction part in low x */
i = ((u_int32_t)(0xffffffff))>>(j0-20);
if((i1&i)==0) { /* x is integral */
u_int32_t high;
*iptr = x;
GET_HIGH_WORD(high,x);
INSERT_WORDS(x,high&0x80000000,0); /* return +-0 */
return x;
} else {
INSERT_WORDS(*iptr,i0,i1&(~i));
return x - *iptr;
}
}
}

63
liba/src/external/openbsd/s_trunc.c vendored Normal file
View File

@@ -0,0 +1,63 @@
/* @(#)s_floor.c 5.1 93/09/24 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
#if 0
#include <sys/cdefs.h>
__FBSDID("$FreeBSD: src/lib/msun/src/s_trunc.c,v 1.1 2004/06/20 09:25:43 das Exp $");
#endif
/*
* trunc(x)
* Return x rounded toward 0 to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to trunc(x).
*/
#include "math.h"
#include "math_private.h"
static const double huge = 1.0e300;
double
trunc(double x)
{
int32_t i0,i1,jj0;
u_int32_t i;
EXTRACT_WORDS(i0,i1,x);
jj0 = ((i0>>20)&0x7ff)-0x3ff;
if(jj0<20) {
if(jj0<0) { /* raise inexact if x != 0 */
if(huge+x>0.0) {/* |x|<1, so return 0*sign(x) */
i0 &= 0x80000000U;
i1 = 0;
}
} else {
i = (0x000fffff)>>jj0;
if(((i0&i)|i1)==0) return x; /* x is integral */
if(huge+x>0.0) { /* raise inexact flag */
i0 &= (~i); i1=0;
}
}
} else if (jj0>51) {
if(jj0==0x400) return x+x; /* inf or NaN */
else return x; /* x is integral */
} else {
i = ((u_int32_t)(0xffffffff))>>(jj0-20);
if((i1&i)==0) return x; /* x is integral */
if(huge+x>0.0) /* raise inexact flag */
i1 &= (~i);
}
INSERT_WORDS(x,i0,i1);
return x;
}

1
liba/src/external/softfloat/port/softfloat.h vendored
View File

@@ -13,6 +13,7 @@ typedef struct { uint64_t v; } float64_t;
float32_t f64_to_f32(float64_t x);
float64_t f32_to_f64(float32_t x);
float64_t i32_to_f64(int32_t i);
float64_t ui32_to_f64(uint32_t i);
int_fast32_t f64_to_i32_r_minMag(float64_t x, bool exact);
bool f64_eq(float64_t x, float64_t y);
bool f64_le(float64_t x, float64_t y);

59
liba/src/external/softfloat/src/ui32_to_f64.c vendored Normal file
View File

@@ -0,0 +1,59 @@
/*============================================================================
This C source file is part of the SoftFloat IEEE Floating-Point Arithmetic
Package, Release 3c, by John R. Hauser.
Copyright 2011, 2012, 2013, 2014, 2015, 2016 The Regents of the University of
California. All Rights Reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
1. Redistributions of source code must retain the above copyright notice,
this list of conditions, and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright notice,
this list of conditions, and the following disclaimer in the documentation
and/or other materials provided with the distribution.
3. Neither the name of the University nor the names of its contributors may
be used to endorse or promote products derived from this software without
specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS "AS IS", AND ANY
EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, ARE
DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY
DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=============================================================================*/
#include <stdint.h>
#include "platform.h"
#include "internals.h"
#include "softfloat.h"
float64_t ui32_to_f64( uint32_t a )
{
uint_fast64_t uiZ;
int_fast8_t shiftDist;
union ui64_f64 uZ;
if ( ! a ) {
uiZ = 0;
} else {
shiftDist = softfloat_countLeadingZeros32( a ) + 21;
uiZ =
packToF64UI( 0, 0x432 - shiftDist, (uint_fast64_t) a<<shiftDist );
}
uZ.ui = uiZ;
return uZ.f;
}

9
liba/src/ldexp.c Normal file
View File

@@ -0,0 +1,9 @@
#include <math.h>
#undef ldexp
// For some reason, OpenBSD's libm defines ldexpf but doesn't define ldexp
double ldexp(double x, int n) {
return scalbn(x, n);
}

9
liba/src/nearbyint.c Normal file
View File

@@ -0,0 +1,9 @@
#include <math.h>
// See nearbyintf.c for comments
#undef nearbyint
double nearbyint(double x) {
return round(x);
}

6
liba/src/nearbyintf.c
View File

@@ -2,16 +2,16 @@
#undef nearbyintf
/* nearbyintf is not supposed to be the same as roundf according to openBSD
/* nearbyintf is not supposed to be the same as roundf according to OpenBSD's
* documentation. Indeed:
* - they round halfway cases to the nearest integer instead of away from zero
* (nearbyint(4.5) = 4 while round(4.5) = 5 for example).
* - nearbyint is guaranteed never to set floating point flags or trigger
* floating point exceptions (which are disabled by default)
* accroding to openBSD open manual page.
* accroding to OpenBSD open manual page.
*
* However, as nearbyintf is required by micropython only, in order to save
* space and time (nearbyintf openBSD 6.0 implementation required some other
* space and time (nearbyintf OpenBSD 6.0 implementation required some other
* function implementation), we used roundf. */
float nearbyintf(float x) {