encounter/tp
0
0
Fork 0
mirror of https://github.com/encounter/tp.git synced 2026-07-11 06:18:46 -07:00
Code Issues Packages Projects Releases Wiki Activity

clean up dolphin files / work on some rels (#212)

Browse Source 
* d_a_alldie / d_a_tboxSw / d_a_tag_gstart / d_a_tag_hstop

* dolphin OS work / cleanup

* dolphin GX work / cleanup

* finish changing dolphin files to C

* more files into C

* match rest of MSL_C math functions

* more dolphin files converted to C

* remove asm

* d_bg_w work

* remove asm

* d_a_alink work / kytag14
This commit is contained in:
TakaRikka
2022-11-11 11:09:48 -07:00
committed by GitHub
parent f03b959831
commit 1114b13da8
665 changed files with 14862 additions and 30501 deletions
Download Patch File Download Diff File Expand all files Collapse all files
libs/MSL_C/Makefile
+2
View File
@@ -121,6 +121,8 @@ LIBMSL_C_PPCEABI_BARE_H_A_O_FILES := \
LIBMSL_C_PPCEABI_BARE_H_A_CFLAGS := \
-O4,p \
-lang=c \
-fp_contract on \
-use_lmw_stmw on \
LIBMSL_C_PPCEABI_BARE_H_A_LDFLAGS := \
-nodefaults \
libs/MSL_C/Math/Double_precision/e_acos.c
+101 -84
View File
@@ -1,89 +1,106 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_acos
//
#include "MSL_C/Math/Double_precision/e_acos.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)e_acos.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/* __ieee754_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
*/
void __ieee754_acos();
#include "fdlibm.h"
#include "MSL_C/math.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
one= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
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 */
void sqrt();
extern u32 __float_nan;
//
// Declarations:
//
/* ############################################################################################## */
/* 80456678-80456680 004C78 0008+00 1/1 0/0 0/0 .sdata2 @83 */
SECTION_SDATA2 static u8 lit_83[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
/* 80456680-80456688 004C80 0008+00 1/1 0/0 0/0 .sdata2 @84 */
SECTION_SDATA2 static f64 lit_84 = 3.141592653589793;
/* 80456688-80456690 004C88 0008+00 1/1 0/0 0/0 .sdata2 @85 */
SECTION_SDATA2 static f64 lit_85 = 1.5707963267948966;
/* 80456690-80456698 004C90 0008+00 1/1 0/0 0/0 .sdata2 @86 */
SECTION_SDATA2 static f64 lit_86 = 6.123233995736766e-17;
/* 80456698-804566A0 004C98 0008+00 1/1 0/0 0/0 .sdata2 @87 */
SECTION_SDATA2 static f64 lit_87 = 1.0 / 6.0;
/* 804566A0-804566A8 004CA0 0008+00 1/1 0/0 0/0 .sdata2 @88 */
SECTION_SDATA2 static f64 lit_88 = -0.3255658186224009;
/* 804566A8-804566B0 004CA8 0008+00 1/1 0/0 0/0 .sdata2 @89 */
SECTION_SDATA2 static f64 lit_89 = 0.20121253213486293;
/* 804566B0-804566B8 004CB0 0008+00 1/1 0/0 0/0 .sdata2 @90 */
SECTION_SDATA2 static f64 lit_90 = -0.04005553450067941;
/* 804566B8-804566C0 004CB8 0008+00 1/1 0/0 0/0 .sdata2 @91 */
SECTION_SDATA2 static f64 lit_91 = 0.0007915349942898145;
/* 804566C0-804566C8 004CC0 0008+00 1/1 0/0 0/0 .sdata2 @92 */
SECTION_SDATA2 static f64 lit_92 = 3.479331075960212e-05;
/* 804566C8-804566D0 004CC8 0008+00 1/1 0/0 0/0 .sdata2 @93 */
SECTION_SDATA2 static f64 lit_93 = 1.0;
/* 804566D0-804566D8 004CD0 0008+00 1/1 0/0 0/0 .sdata2 @94 */
SECTION_SDATA2 static f64 lit_94 = -2.403394911734414;
/* 804566D8-804566E0 004CD8 0008+00 1/1 0/0 0/0 .sdata2 @95 */
SECTION_SDATA2 static f64 lit_95 = 2.0209457602335057;
/* 804566E0-804566E8 004CE0 0008+00 1/1 0/0 0/0 .sdata2 @96 */
SECTION_SDATA2 static f64 lit_96 = -0.6882839716054533;
/* 804566E8-804566F0 004CE8 0008+00 1/1 0/0 0/0 .sdata2 @97 */
SECTION_SDATA2 static f64 lit_97 = 0.07703815055590194;
/* 804566F0-804566F8 004CF0 0008+00 1/1 0/0 0/0 .sdata2 @98 */
SECTION_SDATA2 static f64 lit_98 = 0.5;
/* 804566F8-80456700 004CF8 0008+00 1/1 0/0 0/0 .sdata2 @99 */
SECTION_SDATA2 static f64 lit_99 = 2.0;
/* 80369274-803694B0 363BB4 023C+00 0/0 1/1 0/0 .text __ieee754_acos */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_acos() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_acos/__ieee754_acos.s"
}
#pragma pop
#ifdef __STDC__
double __ieee754_acos(double x)
#else
double __ieee754_acos(x)
double x;
#endif
{
double z,p,q,r,w,s,c,df;
int hx,ix;
hx = __HI(x);
ix = hx&0x7fffffff;
if(ix>=0x3ff00000) { /* |x| >= 1 */
if(((ix-0x3ff00000)|__LO(x))==0) { /* |x|==1 */
if(hx>0) return 0.0; /* acos(1) = 0 */
else return pi+2.0*pio2_lo; /* acos(-1)= pi */
}
return *__float_nan; /* acos(|x|>1) is NaN */
}
if(ix<0x3fe00000) { /* |x| < 0.5 */
if(ix<=0x3c600000) return pio2_hi+pio2_lo;/*if|x|<2**-57*/
z = x*x;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+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 = (one+x)*0.5;
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+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 = (one-x)*0.5;
s = sqrt(z);
df = s;
__LO(df) = 0;
c = (z-df*df)/(s+df);
p = z*(pS0+z*(pS1+z*(pS2+z*(pS3+z*(pS4+z*pS5)))));
q = one+z*(qS1+z*(qS2+z*(qS3+z*qS4)));
r = p/q;
w = r*s+c;
return 2.0*(df+w);
}
}
libs/MSL_C/Math/Double_precision/e_asin.c
+112 -82
View File
@@ -1,87 +1,117 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_asin
//
#include "MSL_C/Math/Double_precision/e_asin.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)e_asin.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/* __ieee754_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.
*
*/
void __ieee754_asin();
#include "fdlibm.h"
#include "MSL_C/math.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
one
= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
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 */
/* coefficient 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 */
void sqrt();
extern u32 __float_nan;
//
// Declarations:
//
/* ############################################################################################## */
/* 80456700-80456708 004D00 0008+00 1/1 0/0 0/0 .sdata2 @94 */
SECTION_SDATA2 static f64 lit_94 = 1.5707963267948966;
/* 80456708-80456710 004D08 0008+00 1/1 0/0 0/0 .sdata2 @95 */
SECTION_SDATA2 static f64 lit_95 = 6.123233995736766e-17;
/* 80456710-80456718 004D10 0008+00 1/1 0/0 0/0 .sdata2 @96 */
SECTION_SDATA2 static f64 lit_96 = 1e+300;
/* 80456718-80456720 004D18 0008+00 1/1 0/0 0/0 .sdata2 @97 */
SECTION_SDATA2 static f64 lit_97 = 1.0;
/* 80456720-80456728 004D20 0008+00 1/1 0/0 0/0 .sdata2 @98 */
SECTION_SDATA2 static f64 lit_98 = 1.0 / 6.0;
/* 80456728-80456730 004D28 0008+00 1/1 0/0 0/0 .sdata2 @99 */
SECTION_SDATA2 static f64 lit_99 = -0.3255658186224009;
/* 80456730-80456738 004D30 0008+00 1/1 0/0 0/0 .sdata2 @100 */
SECTION_SDATA2 static f64 lit_100 = 0.20121253213486293;
/* 80456738-80456740 004D38 0008+00 1/1 0/0 0/0 .sdata2 @101 */
SECTION_SDATA2 static f64 lit_101 = -0.04005553450067941;
/* 80456740-80456748 004D40 0008+00 1/1 0/0 0/0 .sdata2 @102 */
SECTION_SDATA2 static f64 lit_102 = 0.0007915349942898145;
/* 80456748-80456750 004D48 0008+00 1/1 0/0 0/0 .sdata2 @103 */
SECTION_SDATA2 static f64 lit_103 = 3.479331075960212e-05;
/* 80456750-80456758 004D50 0008+00 1/1 0/0 0/0 .sdata2 @104 */
SECTION_SDATA2 static f64 lit_104 = -2.403394911734414;
/* 80456758-80456760 004D58 0008+00 1/1 0/0 0/0 .sdata2 @105 */
SECTION_SDATA2 static f64 lit_105 = 2.0209457602335057;
/* 80456760-80456768 004D60 0008+00 1/1 0/0 0/0 .sdata2 @106 */
SECTION_SDATA2 static f64 lit_106 = -0.6882839716054533;
/* 80456768-80456770 004D68 0008+00 1/1 0/0 0/0 .sdata2 @107 */
SECTION_SDATA2 static f64 lit_107 = 0.07703815055590194;
/* 80456770-80456778 004D70 0008+00 1/1 0/0 0/0 .sdata2 @108 */
SECTION_SDATA2 static f64 lit_108 = 0.5;
/* 80456778-80456780 004D78 0008+00 1/1 0/0 0/0 .sdata2 @109 */
SECTION_SDATA2 static f64 lit_109 = 2.0;
/* 80456780-80456788 004D80 0008+00 1/1 0/0 0/0 .sdata2 @110 */
SECTION_SDATA2 static f64 lit_110 = 0.7853981633974483;
/* 803694B0-803696E8 363DF0 0238+00 0/0 1/1 0/0 .text __ieee754_asin */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_asin() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_asin/__ieee754_asin.s"
}
#pragma pop
#ifdef __STDC__
double __ieee754_asin(double x)
#else
double __ieee754_asin(x) double x;
#endif
{
double t, w, p, q, c, r, s;
int hx, ix;
hx = __HI(x);
ix = hx & 0x7fffffff;
if (ix >= 0x3ff00000) { /* |x|>= 1 */
if (((ix - 0x3ff00000) | __LO(x)) == 0)
/* asin(1)=+-pi/2 with inexact */
return x * pio2_hi + x * pio2_lo;
return *__float_nan; /* asin(|x|>1) is NaN */
} else if (ix < 0x3fe00000) { /* |x|<0.5 */
if (ix < 0x3e400000) { /* if |x| < 2**-27 */
if (huge + x > one)
return x; /* return x with inexact if x!=0*/
} else
t = x * x;
p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
w = p / q;
return x + x * w;
}
/* 1> |x|>= 0.5 */
w = one - fabs(x);
t = w * 0.5;
p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
q = one + 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;
__LO(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;
else
return -t;
}
libs/MSL_C/Math/Double_precision/e_atan2.c
+134 -61
View File
@@ -1,70 +1,143 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_atan2
//
#include "MSL_C/Math/Double_precision/e_atan2.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)e_atan2.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*
*/
//
// Forward References:
//
/* __ieee754_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.
*/
void __ieee754_atan2();
#include "fdlibm.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
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 */
void atan();
#ifdef __STDC__
double __ieee754_atan2(double y, double x)
#else
double __ieee754_atan2(y, x) double y, x;
#endif
{
double z;
int k, m, hx, hy, ix, iy;
unsigned lx, ly;
//
// Declarations:
//
hx = __HI(x);
ix = hx & 0x7fffffff;
lx = __LO(x);
hy = __HI(y);
iy = hy & 0x7fffffff;
ly = __LO(y);
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) */
/* ############################################################################################## */
/* 80456788-80456790 004D88 0008+00 1/1 0/0 0/0 .sdata2 @145 */
SECTION_SDATA2 static f64 lit_145 = 3.141592653589793;
/* 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;
/* 80456790-80456798 004D90 0008+00 1/1 0/0 0/0 .sdata2 @146 */
SECTION_SDATA2 static f64 lit_146 = -3.141592653589793;
/* 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;
/* 80456798-804567A0 004D98 0008+00 1/1 0/0 0/0 .sdata2 @147 */
SECTION_SDATA2 static f64 lit_147 = -1.5707963267948966;
/* 804567A0-804567A8 004DA0 0008+00 1/1 0/0 0/0 .sdata2 @148 */
SECTION_SDATA2 static f64 lit_148 = 1.5707963267948966;
/* 804567A8-804567B0 004DA8 0008+00 1/1 0/0 0/0 .sdata2 @149 */
SECTION_SDATA2 static f64 lit_149 = 0.7853981633974483;
/* 804567B0-804567B8 004DB0 0008+00 1/1 0/0 0/0 .sdata2 @150 */
SECTION_SDATA2 static f64 lit_150 = -0.7853981633974483;
/* 804567B8-804567C0 004DB8 0008+00 1/1 0/0 0/0 .sdata2 @151 */
SECTION_SDATA2 static f64 lit_151 = 2.356194490192345;
/* 804567C0-804567C8 004DC0 0008+00 1/1 0/0 0/0 .sdata2 @152 */
SECTION_SDATA2 static f64 lit_152 = -2.356194490192345;
/* 804567C8-804567D0 004DC8 0008+00 1/1 0/0 0/0 .sdata2 @153 */
SECTION_SDATA2 static u8 lit_153[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
/* 804567D0-804567D8 004DD0 0008+00 1/1 0/0 0/0 .sdata2 @154 */
SECTION_SDATA2 static f64 lit_154 = -0.0;
/* 804567D8-804567E0 004DD8 0008+00 1/1 0/0 0/0 .sdata2 @155 */
SECTION_SDATA2 static f64 lit_155 = 1.2246467991473532e-16;
/* 803696E8-80369978 364028 0290+00 0/0 1/1 0/0 .text __ieee754_atan2 */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_atan2() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_atan2/__ieee754_atan2.s"
}
#pragma pop
/* 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:
__HI(z) ^= 0x80000000;
return z; /* atan(-,+) */
case 2:
return pi - (z - pi_lo); /* atan(+,-) */
default: /* case 3 */
return (z - pi_lo) - pi; /* atan(-,-) */
}
}
libs/MSL_C/Math/Double_precision/e_exp.c
+153 -95
View File
@@ -1,104 +1,162 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_exp
//
#include "MSL_C/Math/Double_precision/e_exp.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)e_exp.c 1.6 04/04/22 */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/* __ieee754_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 Remes 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*R1(r)
* = 1 + r + ----------- (for better accuracy)
* 2 - R1(r)
* where
* 2 4 10
* R1(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 > 7.09782712893383973096e+02 then exp(x) overflow
* if x < -7.45133219101941108420e+02 then exp(x) underflow
*
* 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.
*/
void __ieee754_exp();
#include "fdlibm.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
one = 1.0,
halF[2] = {0.5,-0.5,},
huge = 1.0e+300,
twom1000= 9.33263618503218878990e-302, /* 2**-1000=0x01700000,0*/
o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
-6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
-1.90821492927058770002e-10,},/* 0xbdea39ef, 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 */
//
// Declarations:
//
#ifdef __STDC__
double __ieee754_exp(double x) /* default IEEE double exp */
#else
double __ieee754_exp(x) /* default IEEE double exp */
double x;
#endif
{
double y, hi, lo, c, t;
int k, xsb;
unsigned hx;
/* ############################################################################################## */
/* 803A2340-803A2350 02E9A0 0010+00 1/1 0/0 0/0 .rodata halF */
SECTION_RODATA static u8 const halF[16] = {
0x3F, 0xE0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0xBF, 0xE0, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
COMPILER_STRIP_GATE(0x803A2340, &halF);
hx = __HI(x); /* high word of x */
xsb = (hx >> 31) & 1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* 803A2350-803A2360 02E9B0 0010+00 0/1 0/0 0/0 .rodata ln2HI */
#pragma push
#pragma force_active on
SECTION_RODATA static u8 const ln2HI[16] = {
0x3F, 0xE6, 0x2E, 0x42, 0xFE, 0xE0, 0x00, 0x00, 0xBF, 0xE6, 0x2E, 0x42, 0xFE, 0xE0, 0x00, 0x00,
};
COMPILER_STRIP_GATE(0x803A2350, &ln2HI);
#pragma pop
/* filter out non-finite argument */
if (hx >= 0x40862E42) { /* if |x|>=709.78... */
if (hx >= 0x7ff00000) {
if (((hx & 0xfffff) | __LO(x)) != 0)
return x + x; /* NaN */
else
return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */
}
if (x > o_threshold)
return huge * huge; /* overflow */
if (x < u_threshold)
return twom1000 * twom1000; /* underflow */
}
/* 803A2360-803A2370 02E9C0 0010+00 0/1 0/0 0/0 .rodata ln2LO */
#pragma push
#pragma force_active on
SECTION_RODATA static u8 const ln2LO[16] = {
0x3D, 0xEA, 0x39, 0xEF, 0x35, 0x79, 0x3C, 0x76, 0xBD, 0xEA, 0x39, 0xEF, 0x35, 0x79, 0x3C, 0x76,
};
COMPILER_STRIP_GATE(0x803A2360, &ln2LO);
#pragma pop
/* argument reduction */
if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */
hi = x - ln2HI[xsb];
lo = ln2LO[xsb];
k = 1 - xsb - xsb;
} else {
k = (int)(invln2 * x + halF[xsb]);
t = k;
hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
lo = t * ln2LO[0];
}
x = hi - lo;
} else if (hx < 0x3e300000) { /* when |x|<2**-28 */
if (huge + x > one)
return one + x; /* trigger inexact */
} else
k = 0;
/* 804567E0-804567E8 004DE0 0008+00 1/1 0/0 0/0 .sdata2 @115 */
SECTION_SDATA2 static u8 lit_115[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
/* 804567E8-804567F0 004DE8 0008+00 1/1 0/0 0/0 .sdata2 @116 */
SECTION_SDATA2 static f64 lit_116 = 709.782712893384;
/* 804567F0-804567F8 004DF0 0008+00 1/1 0/0 0/0 .sdata2 @117 */
SECTION_SDATA2 static f64 lit_117 = DOUBLE_INF;
/* 804567F8-80456800 004DF8 0008+00 1/1 0/0 0/0 .sdata2 @118 */
SECTION_SDATA2 static f64 lit_118 = -745.1332191019411;
/* 80456800-80456808 004E00 0008+00 1/1 0/0 0/0 .sdata2 @119 */
SECTION_SDATA2 static f64 lit_119 = 1.4426950408889634;
/* 80456808-80456810 004E08 0008+00 1/1 0/0 0/0 .sdata2 @120 */
SECTION_SDATA2 static f64 lit_120 = 1e+300;
/* 80456810-80456818 004E10 0008+00 1/1 0/0 0/0 .sdata2 @121 */
SECTION_SDATA2 static f64 lit_121 = 1.0;
/* 80456818-80456820 004E18 0008+00 1/1 0/0 0/0 .sdata2 @122 */
SECTION_SDATA2 static f64 lit_122 = 0.16666666666666602;
/* 80456820-80456828 004E20 0008+00 1/1 0/0 0/0 .sdata2 @123 */
SECTION_SDATA2 static f64 lit_123 = -0.0027777777777015593;
/* 80456828-80456830 004E28 0008+00 1/1 0/0 0/0 .sdata2 @124 */
SECTION_SDATA2 static f64 lit_124 = 6.613756321437934e-05;
/* 80456830-80456838 004E30 0008+00 1/1 0/0 0/0 .sdata2 @125 */
SECTION_SDATA2 static f64 lit_125 = -1.6533902205465252e-06;
/* 80456838-80456840 004E38 0008+00 1/1 0/0 0/0 .sdata2 @126 */
SECTION_SDATA2 static f64 lit_126 = 4.1381367970572385e-08;
/* 80456840-80456848 004E40 0008+00 1/1 0/0 0/0 .sdata2 @127 */
SECTION_SDATA2 static f64 lit_127 = 2.0;
/* 80456848-80456850 004E48 0008+00 1/1 0/0 0/0 .sdata2 @128 */
SECTION_SDATA2 static f64 lit_128 = 9.332636185032189e-302;
/* 80456850-80456858 004E50 0008+00 1/1 0/0 0/0 .sdata2 @131 */
SECTION_SDATA2 static f64 lit_131 = 4503601774854144.0 /* cast s32 to float */;
/* 80369978-80369B9C 3642B8 0224+00 0/0 1/1 0/0 .text __ieee754_exp */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_exp() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_exp/__ieee754_exp.s"
}
#pragma pop
/* x is now in primary range */
t = x * x;
c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
if (k == 0)
return one - ((x * c) / (c - 2.0) - x);
else
y = one - ((lo - (x * c) / (2.0 - c)) - hi);
if (k >= -1021) {
__HI(y) += (k << 20); /* add k to y's exponent */
return y;
} else {
__HI(y) += ((k + 1000) << 20); /* add k to y's exponent */
return y * twom1000;
}
}
libs/MSL_C/Math/Double_precision/e_pow.c
+380 -138
View File
File diff suppressed because it is too large Load Diff
libs/MSL_C/Math/Double_precision/e_rem_pio2.c
+171 -95
View File
@@ -1,105 +1,181 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_rem_pio2
//
#include "MSL_C/Math/Double_precision/e_rem_pio2.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)e_rem_pio2.c 1.4 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*
*/
//
// Forward References:
//
/* __ieee754_rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2()
*/
void __ieee754_rem_pio2();
#include "fdlibm.h"
#include "MSL_C/math.h"
//
// External References:
//
void __kernel_rem_pio2();
//
// Declarations:
//
/* ############################################################################################## */
/* 803A23B0-803A24B8 02EA10 0108+00 1/1 0/0 0/0 .rodata two_over_pi */
SECTION_RODATA static u8 const two_over_pi[264] = {
0x00, 0xA2, 0xF9, 0x83, 0x00, 0x6E, 0x4E, 0x44, 0x00, 0x15, 0x29, 0xFC, 0x00, 0x27, 0x57, 0xD1,
0x00, 0xF5, 0x34, 0xDD, 0x00, 0xC0, 0xDB, 0x62, 0x00, 0x95, 0x99, 0x3C, 0x00, 0x43, 0x90, 0x41,
0x00, 0xFE, 0x51, 0x63, 0x00, 0xAB, 0xDE, 0xBB, 0x00, 0xC5, 0x61, 0xB7, 0x00, 0x24, 0x6E, 0x3A,
0x00, 0x42, 0x4D, 0xD2, 0x00, 0xE0, 0x06, 0x49, 0x00, 0x2E, 0xEA, 0x09, 0x00, 0xD1, 0x92, 0x1C,
0x00, 0xFE, 0x1D, 0xEB, 0x00, 0x1C, 0xB1, 0x29, 0x00, 0xA7, 0x3E, 0xE8, 0x00, 0x82, 0x35, 0xF5,
0x00, 0x2E, 0xBB, 0x44, 0x00, 0x84, 0xE9, 0x9C, 0x00, 0x70, 0x26, 0xB4, 0x00, 0x5F, 0x7E, 0x41,
0x00, 0x39, 0x91, 0xD6, 0x00, 0x39, 0x83, 0x53, 0x00, 0x39, 0xF4, 0x9C, 0x00, 0x84, 0x5F, 0x8B,
0x00, 0xBD, 0xF9, 0x28, 0x00, 0x3B, 0x1F, 0xF8, 0x00, 0x97, 0xFF, 0xDE, 0x00, 0x05, 0x98, 0x0F,
0x00, 0xEF, 0x2F, 0x11, 0x00, 0x8B, 0x5A, 0x0A, 0x00, 0x6D, 0x1F, 0x6D, 0x00, 0x36, 0x7E, 0xCF,
0x00, 0x27, 0xCB, 0x09, 0x00, 0xB7, 0x4F, 0x46, 0x00, 0x3F, 0x66, 0x9E, 0x00, 0x5F, 0xEA, 0x2D,
0x00, 0x75, 0x27, 0xBA, 0x00, 0xC7, 0xEB, 0xE5, 0x00, 0xF1, 0x7B, 0x3D, 0x00, 0x07, 0x39, 0xF7,
0x00, 0x8A, 0x52, 0x92, 0x00, 0xEA, 0x6B, 0xFB, 0x00, 0x5F, 0xB1, 0x1F, 0x00, 0x8D, 0x5D, 0x08,
0x00, 0x56, 0x03, 0x30, 0x00, 0x46, 0xFC, 0x7B, 0x00, 0x6B, 0xAB, 0xF0, 0x00, 0xCF, 0xBC, 0x20,
0x00, 0x9A, 0xF4, 0x36, 0x00, 0x1D, 0xA9, 0xE3, 0x00, 0x91, 0x61, 0x5E, 0x00, 0xE6, 0x1B, 0x08,
0x00, 0x65, 0x99, 0x85, 0x00, 0x5F, 0x14, 0xA0, 0x00, 0x68, 0x40, 0x8D, 0x00, 0xFF, 0xD8, 0x80,
0x00, 0x4D, 0x73, 0x27, 0x00, 0x31, 0x06, 0x06, 0x00, 0x15, 0x56, 0xCA, 0x00, 0x73, 0xA8, 0xC9,
0x00, 0x60, 0xE2, 0x7B, 0x00, 0xC0, 0x8C, 0x6B,
};
COMPILER_STRIP_GATE(0x803A23B0, &two_over_pi);
/* 803A24B8-803A2538 02EB18 0080+00 1/1 0/0 0/0 .rodata npio2_hw */
SECTION_RODATA static u8 const npio2_hw[128] = {
0x3F, 0xF9, 0x21, 0xFB, 0x40, 0x09, 0x21, 0xFB, 0x40, 0x12, 0xD9, 0x7C, 0x40, 0x19, 0x21, 0xFB,
0x40, 0x1F, 0x6A, 0x7A, 0x40, 0x22, 0xD9, 0x7C, 0x40, 0x25, 0xFD, 0xBB, 0x40, 0x29, 0x21, 0xFB,
0x40, 0x2C, 0x46, 0x3A, 0x40, 0x2F, 0x6A, 0x7A, 0x40, 0x31, 0x47, 0x5C, 0x40, 0x32, 0xD9, 0x7C,
0x40, 0x34, 0x6B, 0x9C, 0x40, 0x35, 0xFD, 0xBB, 0x40, 0x37, 0x8F, 0xDB, 0x40, 0x39, 0x21, 0xFB,
0x40, 0x3A, 0xB4, 0x1B, 0x40, 0x3C, 0x46, 0x3A, 0x40, 0x3D, 0xD8, 0x5A, 0x40, 0x3F, 0x6A, 0x7A,
0x40, 0x40, 0x7E, 0x4C, 0x40, 0x41, 0x47, 0x5C, 0x40, 0x42, 0x10, 0x6C, 0x40, 0x42, 0xD9, 0x7C,
0x40, 0x43, 0xA2, 0x8C, 0x40, 0x44, 0x6B, 0x9C, 0x40, 0x45, 0x34, 0xAC, 0x40, 0x45, 0xFD, 0xBB,
0x40, 0x46, 0xC6, 0xCB, 0x40, 0x47, 0x8F, 0xDB, 0x40, 0x48, 0x58, 0xEB, 0x40, 0x49, 0x21, 0xFB,
};
COMPILER_STRIP_GATE(0x803A24B8, &npio2_hw);
/* 80456968-80456970 004F68 0008+00 1/1 0/0 0/0 .sdata2 @145 */
SECTION_SDATA2 static u8 lit_145[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
#ifdef __STDC__
static const int two_over_pi[] = {
#else
static int two_over_pi[] = {
#endif
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7,
0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C,
0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11,
0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7,
0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E,
0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
/* 80456970-80456978 004F70 0008+00 1/1 0/0 0/0 .sdata2 @146 */
SECTION_SDATA2 static f64 lit_146 = 1.5707963267341256;
#ifdef __STDC__
static const int npio2_hw[] = {
#else
static int npio2_hw[] = {
#endif
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, 0x402921FB, 0x402C463A, 0x402F6A7A, 0x4031475C,
0x4032D97C, 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB, 0x403AB41B, 0x403C463A, 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C,
0x4042106C, 0x4042D97C, 0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 0x404858EB, 0x404921FB,
};
/* 80456978-80456980 004F78 0008+00 1/1 0/0 0/0 .sdata2 @147 */
SECTION_SDATA2 static f64 lit_147 = 6.077100506506192e-11;
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
/* 80456980-80456988 004F80 0008+00 1/1 0/0 0/0 .sdata2 @148 */
SECTION_SDATA2 static f64 lit_148 = 6.077100506303966e-11;
#ifdef __STDC__
static const double
#else
static double
#endif
zero
= 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
/* 80456988-80456990 004F88 0008+00 1/1 0/0 0/0 .sdata2 @149 */
SECTION_SDATA2 static f64 lit_149 = 2.0222662487959506e-21;
#ifdef __STDC__
int __ieee754_rem_pio2(double x, double* y)
#else
int __ieee754_rem_pio2(x, y) double x, y[];
#endif
{
double z, w, t, r, fn;
double tx[3];
int e0, i, j, nx, n, ix, hx;
/* 80456990-80456998 004F90 0008+00 1/1 0/0 0/0 .sdata2 @150 */
SECTION_SDATA2 static f64 lit_150 = 0.5;
/* 80456998-804569A0 004F98 0008+00 1/1 0/0 0/0 .sdata2 @151 */
SECTION_SDATA2 static f64 lit_151 = 0.6366197723675814;
/* 804569A0-804569A8 004FA0 0008+00 1/1 0/0 0/0 .sdata2 @152 */
SECTION_SDATA2 static f64 lit_152 = 2.0222662487111665e-21;
/* 804569A8-804569B0 004FA8 0008+00 1/1 0/0 0/0 .sdata2 @153 */
SECTION_SDATA2 static f64 lit_153 = 8.4784276603689e-32;
/* 804569B0-804569B8 004FB0 0008+00 1/1 0/0 0/0 .sdata2 @154 */
SECTION_SDATA2 static f64 lit_154 = 16777216.0;
/* 804569B8-804569C0 004FB8 0008+00 1/1 0/0 0/0 .sdata2 @157 */
SECTION_SDATA2 static f64 lit_157 = 4503601774854144.0 /* cast s32 to float */;
/* 8036A708-8036AAA8 365048 03A0+00 0/0 3/3 0/0 .text __ieee754_rem_pio2 */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_rem_pio2() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_rem_pio2/__ieee754_rem_pio2.s"
}
#pragma pop
hx = __HI(x); /* high word of x */
ix = hx & 0x7fffffff;
if (ix <= 0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{
y[0] = x;
y[1] = 0;
return 0;
}
if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
if (hx > 0) {
z = x - pio2_1;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z - pio2_1t;
y[1] = (z - y[0]) - pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z - y[0]) - pio2_2t;
}
return 1;
} else { /* negative x */
z = x + pio2_1;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z + pio2_1t;
y[1] = (z - y[0]) + pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z - y[0]) + pio2_2t;
}
return -1;
}
}
if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
t = fabs(x);
n = (int)(t * invpio2 + half);
fn = (double)n;
r = t - fn * pio2_1;
w = fn * pio2_1t; /* 1st round good to 85 bit */
if (n < 32 && ix != npio2_hw[n - 1]) {
y[0] = r - w; /* quick check no cancellation */
} else {
j = ix >> 20;
y[0] = r - w;
i = j - (((__HI(y[0])) >> 20) & 0x7ff);
if (i > 16) { /* 2nd iteration needed, good to 118 */
t = r;
r = t - fn * pio2_2;
w = fn * pio2_2t - ((t - r) - fn * pio2_2);
y[0] = r - w;
i = j - (((__HI(y[0])) >> 20) & 0x7ff);
if (i > 49) { /* 3rd iteration need, 151 bits acc */
t = r; /* will cover all possible cases */
w = fn * pio2_3;
r = t - w;
w = fn * pio2_3t - ((t - r) - w);
y[0] = r - w;
}
}
}
y[1] = (r - y[0]) - w;
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
} else
return n;
}
/*
* all other (large) arguments
*/
if (ix >= 0x7ff00000) { /* x is inf or NaN */
y[0] = y[1] = x - x;
return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
__LO(z) = __LO(x);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
__HI(z) = ix - (e0 << 20);
for (i = 0; i < 2; i++) {
tx[i] = (double)((int)(z));
z = (z - tx[i]) * two24;
}
tx[2] = z;
nx = 3;
while (tx[nx - 1] == zero)
nx--; /* skip zero term */
n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
}
return n;
}
libs/MSL_C/Math/Double_precision/e_sqrt.c
+453 -30
View File
@@ -1,39 +1,462 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/e_sqrt
//
/* @(#)e_sqrt.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
#include "MSL_C/Math/Double_precision/e_sqrt.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* __ieee754_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
*
* Other methods : see the appended file at the end of the program below.
*---------------
*/
//
// Forward References:
//
#include "fdlibm.h"
#include "errno.h"
#include "MSL_C/math.h"
void __ieee754_sqrt();
#ifdef __STDC__
static const double one = 1.0, tiny = 1.0e-300;
#else
static double one = 1.0, tiny = 1.0e-300;
#endif
//
// External References:
//
#ifdef __STDC__
double __ieee754_sqrt(double x)
#else
double __ieee754_sqrt(x) double x;
#endif
{
double z;
int sign = (int)0x80000000;
unsigned r, t1, s1, ix1, q1;
int ix0, s0, q, m, t, i;
extern u32 __float_nan;
extern u8 errno[4 + 4 /* padding */];
ix0 = __HI(x); /* high word of x */
ix1 = __LO(x); /* low word of x */
//
// Declarations:
//
/* take care of Inf and NaN */
if ((ix0 & 0x7ff00000) == 0x7ff00000) {
errno = 33;
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 */
else if (ix0 < 0) {
errno = 33;
return *__float_nan;
} /* 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] */
/* ############################################################################################## */
/* 80456B48-80456B50 005148 0008+00 1/1 0/0 0/0 .sdata2 @164 */
SECTION_SDATA2 static f64 lit_164 = 1.0;
/* 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 */
/* 8036C7A0-8036C9C4 3670E0 0224+00 0/0 1/1 0/0 .text __ieee754_sqrt */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __ieee754_sqrt() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/e_sqrt/__ieee754_sqrt.s"
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 += 1;
ix0 -= t;
if (ix1 < t1)
ix0 -= 1;
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 = one - tiny; /* trigger inexact flag */
if (z >= one) {
z = one + tiny;
if (q1 == (unsigned)0xffffffff) {
q1 = 0;
q += 1;
} else if (z > one) {
if (q1 == (unsigned)0xfffffffe)
q += 1;
q1 += 2;
} else
q1 += (q1 & 1);
}
}
ix0 = (q >> 1) + 0x3fe00000;
ix1 = q1 >> 1;
if ((q & 1) == 1)
ix1 |= sign;
ix0 += (m << 20);
__HI(z) = ix0;
__LO(z) = ix1;
return z;
}
#pragma pop
/*
Other methods (use floating-point arithmetic)
-------------
(This is a copy of a drafted paper by Prof W. Kahan
and K.C. Ng, written in May, 1986)
Two algorithms are given here to implement sqrt(x)
(IEEE double precision arithmetic) in software.
Both supply sqrt(x) correctly rounded. The first algorithm (in
Section A) uses newton iterations and involves four divisions.
The second one uses reciproot iterations to avoid division, but
requires more multiplications. Both algorithms need the ability
to chop results of arithmetic operations instead of round them,
and the INEXACT flag to indicate when an arithmetic operation
is executed exactly with no roundoff error, all part of the
standard (IEEE 754-1985). The ability to perform shift, add,
subtract and logical AND operations upon 32-bit words is needed
too, though not part of the standard.
A. sqrt(x) by Newton Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
1 11 52 ...widths
------------------------------------------------------
x: |s| e | f |
------------------------------------------------------
msb lsb msb lsb ...order
------------------------ ------------------------
x0: |s| e | f1 | x1: | f2 |
------------------------ ------------------------
By performing shifts and subtracts on x0 and x1 (both regarded
as integers), we obtain an 8-bit approximation of sqrt(x) as
follows.
k := (x0>>1) + 0x1ff80000;
y0 := k - T1[31&(k>>15)]. ... y ~ sqrt(x) to 8 bits
Here k is a 32-bit integer and T1[] is an integer array containing
correction terms. Now magically the floating value of y (y's
leading 32-bit word is y0, the value of its trailing word is 0)
approximates sqrt(x) to almost 8-bit.
Value of T1:
static int T1[32]= {
0, 1024, 3062, 5746, 9193, 13348, 18162, 23592,
29598, 36145, 43202, 50740, 58733, 67158, 75992, 85215,
83599, 71378, 60428, 50647, 41945, 34246, 27478, 21581,
16499, 12183, 8588, 5674, 3403, 1742, 661, 130,};
(2) Iterative refinement
Apply Heron's rule three times to y, we have y approximates
sqrt(x) to within 1 ulp (Unit in the Last Place):
y := (y+x/y)/2 ... almost 17 sig. bits
y := (y+x/y)/2 ... almost 35 sig. bits
y := y-(y-x/y)/2 ... within 1 ulp
Remark 1.
Another way to improve y to within 1 ulp is:
y := (y+x/y) ... almost 17 sig. bits to 2*sqrt(x)
y := y - 0x00100006 ... almost 18 sig. bits to sqrt(x)
2
(x-y )*y
y := y + 2* ---------- ...within 1 ulp
2
3y + x
This formula has one division fewer than the one above; however,
it requires more multiplications and additions. Also x must be
scaled in advance to avoid spurious overflow in evaluating the
expression 3y*y+x. Hence it is not recommended uless division
is slow. If division is very slow, then one should use the
reciproot algorithm given in section B.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
I := FALSE; ... reset INEXACT flag I
R := RZ; ... set rounding mode to round-toward-zero
z := x/y; ... chopped quotient, possibly inexact
If(not I) then { ... if the quotient is exact
if(z=y) {
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
} else {
z := z - ulp; ... special rounding
}
}
i := TRUE; ... sqrt(x) is inexact
If (r=RN) then z=z+ulp ... rounded-to-nearest
If (r=RP) then { ... round-toward-+inf
y = y+ulp; z=z+ulp;
}
y := y+z; ... chopped sum
y0:=y0-0x00100000; ... y := y/2 is correctly rounded.
I := i; ... restore inexact flag
R := r; ... restore rounded mode
return sqrt(x):=y.
(4) Special cases
Square root of +inf, +-0, or NaN is itself;
Square root of a negative number is NaN with invalid signal.
B. sqrt(x) by Reciproot Iteration
(1) Initial approximation
Let x0 and x1 be the leading and the trailing 32-bit words of
a floating point number x (in IEEE double format) respectively
(see section A). By performing shifs and subtracts on x0 and y0,
we obtain a 7.8-bit approximation of 1/sqrt(x) as follows.
k := 0x5fe80000 - (x0>>1);
y0:= k - T2[63&(k>>14)]. ... y ~ 1/sqrt(x) to 7.8 bits
Here k is a 32-bit integer and T2[] is an integer array
containing correction terms. Now magically the floating
value of y (y's leading 32-bit word is y0, the value of
its trailing word y1 is set to zero) approximates 1/sqrt(x)
to almost 7.8-bit.
Value of T2:
static int T2[64]= {
0x1500, 0x2ef8, 0x4d67, 0x6b02, 0x87be, 0xa395, 0xbe7a, 0xd866,
0xf14a, 0x1091b,0x11fcd,0x13552,0x14999,0x15c98,0x16e34,0x17e5f,
0x18d03,0x19a01,0x1a545,0x1ae8a,0x1b5c4,0x1bb01,0x1bfde,0x1c28d,
0x1c2de,0x1c0db,0x1ba73,0x1b11c,0x1a4b5,0x1953d,0x18266,0x16be0,
0x1683e,0x179d8,0x18a4d,0x19992,0x1a789,0x1b445,0x1bf61,0x1c989,
0x1d16d,0x1d77b,0x1dddf,0x1e2ad,0x1e5bf,0x1e6e8,0x1e654,0x1e3cd,
0x1df2a,0x1d635,0x1cb16,0x1be2c,0x1ae4e,0x19bde,0x1868e,0x16e2e,
0x1527f,0x1334a,0x11051,0xe951, 0xbe01, 0x8e0d, 0x5924, 0x1edd,};
(2) Iterative refinement
Apply Reciproot iteration three times to y and multiply the
result by x to get an approximation z that matches sqrt(x)
to about 1 ulp. To be exact, we will have
-1ulp < sqrt(x)-z<1.0625ulp.
... set rounding mode to Round-to-nearest
y := y*(1.5-0.5*x*y*y) ... almost 15 sig. bits to 1/sqrt(x)
y := y*((1.5-2^-30)+0.5*x*y*y)... about 29 sig. bits to 1/sqrt(x)
... special arrangement for better accuracy
z := x*y ... 29 bits to sqrt(x), with z*y<1
z := z + 0.5*z*(1-z*y) ... about 1 ulp to sqrt(x)
Remark 2. The constant 1.5-2^-30 is chosen to bias the error so that
(a) the term z*y in the final iteration is always less than 1;
(b) the error in the final result is biased upward so that
-1 ulp < sqrt(x) - z < 1.0625 ulp
instead of |sqrt(x)-z|<1.03125ulp.
(3) Final adjustment
By twiddling y's last bit it is possible to force y to be
correctly rounded according to the prevailing rounding mode
as follows. Let r and i be copies of the rounding mode and
inexact flag before entering the square root program. Also we
use the expression y+-ulp for the next representable floating
numbers (up and down) of y. Note that y+-ulp = either fixed
point y+-1, or multiply y by nextafter(1,+-inf) in chopped
mode.
R := RZ; ... set rounding mode to round-toward-zero
switch(r) {
case RN: ... round-to-nearest
if(x<= z*(z-ulp)...chopped) z = z - ulp; else
if(x<= z*(z+ulp)...chopped) z = z; else z = z+ulp;
break;
case RZ:case RM: ... round-to-zero or round-to--inf
R:=RP; ... reset rounding mod to round-to-+inf
if(x<z*z ... rounded up) z = z - ulp; else
if(x>=(z+ulp)*(z+ulp) ...rounded up) z = z+ulp;
break;
case RP: ... round-to-+inf
if(x>(z+ulp)*(z+ulp)...chopped) z = z+2*ulp; else
if(x>z*z ...chopped) z = z+ulp;
break;
}
Remark 3. The above comparisons can be done in fixed point. For
example, to compare x and w=z*z chopped, it suffices to compare
x1 and w1 (the trailing parts of x and w), regarding them as
two's complement integers.
...Is z an exact square root?
To determine whether z is an exact square root of x, let z1 be the
trailing part of z, and also let x0 and x1 be the leading and
trailing parts of x.
If ((z1&0x03ffffff)!=0) ... not exact if trailing 26 bits of z!=0
I := 1; ... Raise Inexact flag: z is not exact
else {
j := 1 - [(x0>>20)&1] ... j = logb(x) mod 2
k := z1 >> 26; ... get z's 25-th and 26-th
fraction bits
I := i or (k&j) or ((k&(j+j+1))!=(x1&3));
}
R:= r ... restore rounded mode
return sqrt(x):=z.
If multiplication is cheaper then the foregoing red tape, the
Inexact flag can be evaluated by
I := i;
I := (z*z!=x) or I.
Note that z*z can overwrite I; this value must be sensed if it is
True.
Remark 4. If z*z = x exactly, then bit 25 to bit 0 of z1 must be
zero.
--------------------
z1: | f2 |
--------------------
bit 31 bit 0
Further more, bit 27 and 26 of z1, bit 0 and 1 of x1, and the odd
or even of logb(x) have the following relations:
-------------------------------------------------
bit 27,26 of z1 bit 1,0 of x1 logb(x)
-------------------------------------------------
00 00 odd and even
01 01 even
10 10 odd
10 00 even
11 01 even
-------------------------------------------------
(4) Special cases (see (4) of Section A).
*/
libs/MSL_C/Math/Double_precision/k_cos.c
+88 -55
View File
@@ -1,60 +1,93 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/k_cos
//
#include "MSL_C/Math/Double_precision/k_cos.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)k_cos.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/*
* __kernel_cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
* 3. cos(x) is approximated by a polynomial of degree 14 on
* [0,pi/4]
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) = 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy when x > 0.3, let qx = |x|/4 with
* the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
* Then
* cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
* Note that 1-qx and (x*x/2-qx) is EXACT here, and the
* magnitude of the latter is at least a quarter of x*x/2,
* thus, reducing the rounding error in the subtraction.
*/
void __kernel_cos();
#include "fdlibm.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
one
= 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
//
// Declarations:
//
/* ############################################################################################## */
/* 804569C0-804569C8 004FC0 0008+00 1/1 0/0 0/0 .sdata2 @65 */
SECTION_SDATA2 static f64 lit_65 = 1.0;
/* 804569C8-804569D0 004FC8 0008+00 1/1 0/0 0/0 .sdata2 @66 */
SECTION_SDATA2 static f64 lit_66 = 0.0416666666666666;
/* 804569D0-804569D8 004FD0 0008+00 1/1 0/0 0/0 .sdata2 @67 */
SECTION_SDATA2 static f64 lit_67 = -0.001388888888887411;
/* 804569D8-804569E0 004FD8 0008+00 1/1 0/0 0/0 .sdata2 @68 */
SECTION_SDATA2 static f64 lit_68 = 2.480158728947673e-05;
/* 804569E0-804569E8 004FE0 0008+00 1/1 0/0 0/0 .sdata2 @69 */
SECTION_SDATA2 static f64 lit_69 = -2.7557314351390663e-07;
/* 804569E8-804569F0 004FE8 0008+00 1/1 0/0 0/0 .sdata2 @70 */
SECTION_SDATA2 static f64 lit_70 = 2.087572321298175e-09;
/* 804569F0-804569F8 004FF0 0008+00 1/1 0/0 0/0 .sdata2 @71 */
SECTION_SDATA2 static f64 lit_71 = -1.1359647557788195e-11;
/* 804569F8-80456A00 004FF8 0008+00 1/1 0/0 0/0 .sdata2 @72 */
SECTION_SDATA2 static f64 lit_72 = 0.5;
/* 80456A00-80456A08 005000 0008+00 1/1 0/0 0/0 .sdata2 @73 */
SECTION_SDATA2 static f64 lit_73 = 0.28125;
/* 8036AAA8-8036AB9C 3653E8 00F4+00 0/0 2/2 0/0 .text __kernel_cos */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __kernel_cos() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/k_cos/__kernel_cos.s"
}
#pragma pop
#ifdef __STDC__
double __kernel_cos(double x, double y)
#else
double __kernel_cos(x, y) double x, y;
#endif
{
double a, hz, z, r, qx;
int ix;
ix = __HI(x) & 0x7fffffff; /* ix = |x|'s high word*/
if (ix < 0x3e400000) { /* if x < 2**27 */
if (((int)x) == 0)
return one; /* generate inexact */
}
z = x * x;
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
if (ix < 0x3FD33333) /* if |x| < 0.3 */
return one - (0.5 * z - (z * r - x * y));
else {
if (ix > 0x3fe90000) { /* x > 0.78125 */
qx = 0.28125;
} else {
__HI(qx) = ix - 0x00200000; /* x/4 */
__LO(qx) = 0;
}
hz = 0.5 * z - qx;
a = one - qx;
return a - (hz - (z * r - x * y));
}
}
libs/MSL_C/Math/Double_precision/k_rem_pio2.c
+338 -64
View File
@@ -1,79 +1,353 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/k_rem_pio2
//
#include "MSL_C/Math/Double_precision/k_rem_pio2.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)k_rem_pio2.c 1.2 95/01/04 */
/* $Id: k_rem_pio2.c,v 1.2.14.1 2002/01/31 15:24:13 ceciliar Exp $ */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/*
* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
* double x[],y[]; int e0,nx,prec; int ipio2[];
*
* __kernel_rem_pio2 return the last three digits of N with
* y = x - N*pi/2
* so that |y| < pi/2.
*
* The method is to compute the integer (mod 8) and fraction parts of
* (2/pi)*x without doing the full multiplication. In general we
* skip the part of the product that are known to be a huge integer (
* more accurately, = 0 mod 8 ). Thus the number of operations are
* independent of the exponent of the input.
*
* (2/pi) is represented by an array of 24-bit integers in ipio2[].
*
* Input parameters:
* x[] The input value (must be positive) is broken into nx
* pieces of 24-bit integers in double precision format.
* x[i] will be the i-th 24 bit of x. The scaled exponent
* of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
* match x's up to 24 bits.
*
* Example of breaking a double positive z into x[0]+x[1]+x[2]:
* e0 = ilogb(z)-23
* z = ldexp(z,-e0)
* for i = 0,1,2
* x[i] = floor(z)
* z = (z-x[i])*2**24
*
*
* y[] ouput result in an array of double precision numbers.
* The dimension of y[] is:
* 24-bit precision 1
* 53-bit precision 2
* 64-bit precision 2
* 113-bit precision 3
* The actual value is the sum of them. Thus for 113-bit
* precison, one may have to do something like:
*
* long double t,w,r_head, r_tail;
* t = (long double)y[2] + (long double)y[1];
* w = (long double)y[0];
* r_head = t+w;
* r_tail = w - (r_head - t);
*
* e0 The exponent of x[0]
*
* nx dimension of x[]
*
* prec an integer indicating the precision:
* 0 24 bits (single)
* 1 53 bits (double)
* 2 64 bits (extended)
* 3 113 bits (quad)
*
* ipio2[]
* integer array, contains the (24*i)-th to (24*i+23)-th
* bit of 2/pi after binary point. The corresponding
* floating value is
*
* ipio2[i] * 2^(-24(i+1)).
*
* External function:
* double ldexp(), floor();
*
*
* Here is the description of some local variables:
*
* jk jk+1 is the initial number of terms of ipio2[] needed
* in the computation. The recommended value is 2,3,4,
* 6 for single, double, extended,and quad.
*
* jz local integer variable indicating the number of
* terms of ipio2[] used.
*
* jx nx - 1
*
* jv index for pointing to the suitable ipio2[] for the
* computation. In general, we want
* ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
* is an integer. Thus
* e0-3-24*jv >= 0 or (e0-3)/24 >= jv
* Hence jv = max(0,(e0-3)/24).
*
* jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
*
* q[] double array with integral value, representing the
* 24-bits chunk of the product of x and 2/pi.
*
* q0 the corresponding exponent of q[0]. Note that the
* exponent for q[i] would be q0-24*i.
*
* PIo2[] double precision array, obtained by cutting pi/2
* into 24 bits chunks.
*
* f[] ipio2[] in floating point
*
* iq[] integer array by breaking up q[] in 24-bits chunk.
*
* fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
*
* ih integer. If >0 it indicates q[] is >= 0.5, hence
* it also indicates the *sign* of the result.
*
*/
void __kernel_rem_pio2();
/*
* 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.
*/
//
// External References:
//
#include "fdlibm.h"
void _savefpr_25();
void _restfpr_25();
void floor();
void ldexp();
#ifdef __STDC__
static const int init_jk[] = { 2, 3, 4, 6 };
/* initial value for jk */ /*- cc 020130 -*/
#else
static int init_jk[] = { 2, 3, 4, 6 }; /*- cc 020130 -*/
#endif
//
// Declarations:
//
/* ############################################################################################## */
/* 803A2538-803A2548 02EB98 0010+00 1/1 0/0 0/0 .rodata init_jk */
SECTION_RODATA static u8 const init_jk[16] = {
0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x03, 0x00, 0x00, 0x00, 0x04, 0x00, 0x00, 0x00, 0x06,
};
COMPILER_STRIP_GATE(0x803A2538, &init_jk);
/* 803A2548-803A2588 02EBA8 0040+00 1/1 0/0 0/0 .rodata PIo2 */
SECTION_RODATA static u8 const PIo2[64] = {
0x3F, 0xF9, 0x21, 0xFB, 0x40, 0x00, 0x00, 0x00, 0x3E, 0x74, 0x44, 0x2D, 0x00, 0x00, 0x00, 0x00,
0x3C, 0xF8, 0x46, 0x98, 0x80, 0x00, 0x00, 0x00, 0x3B, 0x78, 0xCC, 0x51, 0x60, 0x00, 0x00, 0x00,
0x39, 0xF0, 0x1B, 0x83, 0x80, 0x00, 0x00, 0x00, 0x38, 0x7A, 0x25, 0x20, 0x40, 0x00, 0x00, 0x00,
0x36, 0xE3, 0x82, 0x22, 0x80, 0x00, 0x00, 0x00, 0x35, 0x69, 0xF3, 0x1D, 0x00, 0x00, 0x00, 0x00,
};
COMPILER_STRIP_GATE(0x803A2548, &PIo2);
/* 80456A08-80456A10 005008 0008+00 1/1 0/0 0/0 .sdata2 @436 */
SECTION_SDATA2 static u8 lit_436[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
#ifdef __STDC__
static const double PIo2[] = {
#else
static double PIo2[] = {
#endif
1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
};
/* 80456A10-80456A18 005010 0008+00 1/1 0/0 0/0 .sdata2 @437 */
SECTION_SDATA2 static f64 lit_437 = 5.960464477539063e-08;
#ifdef __STDC__
static const double
#else
static double
#endif
zero
= 0.0,
one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
/* 80456A18-80456A20 005018 0008+00 1/1 0/0 0/0 .sdata2 @438 */
SECTION_SDATA2 static f64 lit_438 = 16777216.0;
#ifdef __STDC__
int __kernel_rem_pio2(double* x, double* y, int e0, int nx, int prec, const int* ipio2) /*- cc 020130 -*/
#else
int __kernel_rem_pio2(x, y, e0, nx, prec, ipio2) /*- cc 020130 -*/
double x[],
y[];
int e0, nx, prec;
int ipio2[]; /*- cc 020130 -*/
#endif
{
int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih; /*- cc 020130 -*/
double z, fw, f[20], fq[20], q[20];
/* 80456A20-80456A28 005020 0008+00 1/1 0/0 0/0 .sdata2 @439 */
SECTION_SDATA2 static f64 lit_439 = 8.0;
/* initialize jk*/
jk = init_jk[prec];
jp = jk;
/* 80456A28-80456A30 005028 0008+00 1/1 0/0 0/0 .sdata2 @440 */
SECTION_SDATA2 static f64 lit_440 = 0.125;
/* determine jx,jv,q0, note that 3>q0 */
jx = nx - 1;
jv = (e0 - 3) / 24;
if (jv < 0)
jv = 0;
q0 = e0 - 24 * (jv + 1);
/* 80456A30-80456A38 005030 0008+00 1/1 0/0 0/0 .sdata2 @441 */
SECTION_SDATA2 static f64 lit_441 = 0.5;
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv - jx;
m = jx + jk;
for (i = 0; i <= m; i++, j++)
f[i] = (j < 0) ? zero : (double)ipio2[j];
/* 80456A38-80456A40 005038 0008+00 1/1 0/0 0/0 .sdata2 @442 */
SECTION_SDATA2 static f64 lit_442 = 1.0;
/* compute q[0],q[1],...q[jk] */
for (i = 0; i <= jk; i++) {
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
/* 80456A40-80456A48 005040 0008+00 1/1 0/0 0/0 .sdata2 @445 */
SECTION_SDATA2 static f64 lit_445 = 4503601774854144.0 /* cast s32 to float */;
jz = jk;
recompute:
/* distill q[] into iq[] reversingly */
for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
fw = (double)((int)(twon24 * z)); /*- cc 020130 -*/
iq[i] = (int)(z - two24 * fw); /*- cc 020130 -*/
z = q[j - 1] + fw;
}
/* 8036AB9C-8036B9F0 3654DC 0E54+00 0/0 1/1 0/0 .text __kernel_rem_pio2 */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __kernel_rem_pio2() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/k_rem_pio2/__kernel_rem_pio2.s"
}
#pragma pop
/* compute n */
z = ldexp(z, q0); /* actual value of z */
z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
n = (int)z; /*- cc 020130 -*/
z -= (double)n;
ih = 0;
if (q0 > 0) { /* need iq[jz-1] to determine n */
i = (iq[jz - 1] >> (24 - q0));
n += i;
iq[jz - 1] -= i << (24 - q0);
ih = iq[jz - 1] >> (23 - q0);
} else if (q0 == 0)
ih = iq[jz - 1] >> 23;
else if (z >= 0.5)
ih = 2;
if (ih > 0) { /* q > 0.5 */
n += 1;
carry = 0;
for (i = 0; i < jz; i++) { /* compute 1-q */
j = iq[i];
if (carry == 0) {
if (j != 0) {
carry = 1;
iq[i] = 0x1000000 - j;
}
} else
iq[i] = 0xffffff - j;
}
if (q0 > 0) { /* rare case: chance is 1 in 12 */
switch (q0) {
case 1:
iq[jz - 1] &= 0x7fffff;
break;
case 2:
iq[jz - 1] &= 0x3fffff;
break;
}
}
if (ih == 2) {
z = one - z;
if (carry != 0)
z -= ldexp(one, q0);
}
}
/* check if recomputation is needed */
if (z == zero) {
j = 0;
for (i = jz - 1; i >= jk; i--)
j |= iq[i];
if (j == 0) { /* need recomputation */
for (k = 1; iq[jk - k] == 0; k++)
; /* k = no. of terms needed */
for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
f[jx + i] = (double)ipio2[jv + i];
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
jz += k;
goto recompute;
}
}
/* chop off zero terms */
if (z == 0.0) {
jz -= 1;
q0 -= 24;
while (iq[jz] == 0) {
jz--;
q0 -= 24;
}
} else { /* break z into 24-bit if necessary */
z = ldexp(z, -q0);
if (z >= two24) {
fw = (double)((int)(twon24 * z)); /*- cc 020130 -*/
iq[jz] = (int)(z - two24 * fw); /*- cc 020130 -*/
jz += 1;
q0 += 24;
iq[jz] = (int)fw; /*- cc 020130 -*/
} else
iq[jz] = (int)z; /*- cc 020130 -*/
}
/* convert integer "bit" chunk to floating-point value */
fw = ldexp(one, q0);
for (i = jz; i >= 0; i--) {
q[i] = fw * (double)iq[i];
fw *= twon24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
for (i = jz; i >= 0; i--) {
for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
fw += PIo2[k] * q[i + k];
fq[jz - i] = fw;
}
/* compress fq[] into y[] */
switch (prec) {
case 0:
fw = 0.0;
for (i = jz; i >= 0; i--)
fw += fq[i];
y[0] = (ih == 0) ? fw : -fw;
break;
case 1:
case 2:
fw = 0.0;
for (i = jz; i >= 0; i--)
fw += fq[i];
y[0] = (ih == 0) ? fw : -fw;
fw = fq[0] - fw;
for (i = 1; i <= jz; i++)
fw += fq[i];
y[1] = (ih == 0) ? fw : -fw;
break;
case 3: /* painful */
for (i = jz; i > 0; i--) {
fw = fq[i - 1] + fq[i];
fq[i] += fq[i - 1] - fw;
fq[i - 1] = fw;
}
for (i = jz; i > 1; i--) {
fw = fq[i - 1] + fq[i];
fq[i] += fq[i - 1] - fw;
fq[i - 1] = fw;
}
for (fw = 0.0, i = jz; i >= 2; i--)
fw += fq[i];
if (ih == 0) {
y[0] = fq[0];
y[1] = fq[1];
y[2] = fw;
} else {
y[0] = -fq[0];
y[1] = -fq[1];
y[2] = -fw;
}
}
return n & 7;
}
libs/MSL_C/Math/Double_precision/k_sin.c
+75 -49
View File
@@ -1,54 +1,80 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/k_sin
//
#include "MSL_C/Math/Double_precision/k_sin.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)k_sin.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
//
// Forward References:
//
/* __kernel_sin( x, y, iy)
* kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
*
* Algorithm
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
* 3. sin(x) is approximated by a polynomial of degree 13 on
* [0,pi/4]
* 3 13
* sin(x) ~ x + S1*x + ... + S6*x
* where
*
* |sin(x) 2 4 6 8 10 12 | -58
* |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
* | x |
*
* 4. sin(x+y) = sin(x) + sin'(x')*y
* ~ sin(x) + (1-x*x/2)*y
* For better accuracy, let
* 3 2 2 2 2
* r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
* then 3 2
* sin(x) = x + (S1*x + (x *(r-y/2)+y))
*/
void __kernel_sin();
#include "fdlibm.h"
//
// External References:
//
#ifdef __STDC__
static const double
#else
static double
#endif
half
= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
//
// Declarations:
//
/* ############################################################################################## */
/* 80456A48-80456A50 005048 0008+00 1/1 0/0 0/0 .sdata2 @60 */
SECTION_SDATA2 static f64 lit_60 = 0.00833333333332249;
/* 80456A50-80456A58 005050 0008+00 1/1 0/0 0/0 .sdata2 @61 */
SECTION_SDATA2 static f64 lit_61 = -0.0001984126982985795;
/* 80456A58-80456A60 005058 0008+00 1/1 0/0 0/0 .sdata2 @62 */
SECTION_SDATA2 static f64 lit_62 = 2.7557313707070068e-06;
/* 80456A60-80456A68 005060 0008+00 1/1 0/0 0/0 .sdata2 @63 */
SECTION_SDATA2 static f64 lit_63 = -2.5050760253406863e-08;
/* 80456A68-80456A70 005068 0008+00 1/1 0/0 0/0 .sdata2 @64 */
SECTION_SDATA2 static f64 lit_64 = 1.58969099521155e-10;
/* 80456A70-80456A78 005070 0008+00 1/1 0/0 0/0 .sdata2 @65 */
SECTION_SDATA2 static f64 lit_65 = -0.16666666666666632;
/* 80456A78-80456A80 005078 0008+00 1/1 0/0 0/0 .sdata2 @66 */
SECTION_SDATA2 static f64 lit_66 = 0.5;
/* 8036B9F0-8036BA90 366330 00A0+00 0/0 2/2 0/0 .text __kernel_sin */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __kernel_sin() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/k_sin/__kernel_sin.s"
}
#pragma pop
#ifdef __STDC__
double __kernel_sin(double x, double y, int iy)
#else
double __kernel_sin(x, y, iy) double x, y;
int iy; /* iy=0 if y is zero */
#endif
{
double z, r, v;
int ix;
ix = __HI(x) & 0x7fffffff; /* high word of x */
if (ix < 0x3e400000) /* |x| < 2**-27 */
{
if ((int)x == 0)
return x;
} /* generate inexact */
z = x * x;
v = z * x;
r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
if (iy == 0)
return x + v * (S1 + z * r);
else
return x - ((z * (half * y - v * r) - y) - v * S1);
}
libs/MSL_C/Math/Double_precision/k_tan.c
+167 -54
View File
@@ -1,68 +1,181 @@
//===========================================================================
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/k_tan
// k_tan.c
//
#include "MSL_C/Math/Double_precision/k_tan.h"
#include "dol2asm.h"
#include "dolphin/types.h"
// Part of the standard mathematical function library
//
// Forward References:
//===========================================================================
//####ECOSGPLCOPYRIGHTBEGIN####
// -------------------------------------------
// This file is part of eCos, the Embedded Configurable Operating System.
// Copyright (C) 1998, 1999, 2000, 2001, 2002 Red Hat, Inc.
//
void __kernel_tan();
// eCos is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the Free
// Software Foundation; either version 2 or (at your option) any later version.
//
// External References:
// eCos is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License along
// with eCos; if not, write to the Free Software Foundation, Inc.,
// 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA.
//
// Declarations:
// As a special exception, if other files instantiate templates or use macros
// or inline functions from this file, or you compile this file and link it
// with other works to produce a work based on this file, this file does not
// by itself cause the resulting work to be covered by the GNU General Public
// License. However the source code for this file must still be made available
// in accordance with section (3) of the GNU General Public License.
//
// This exception does not invalidate any other reasons why a work based on
// this file might be covered by the GNU General Public License.
//
// Alternative licenses for eCos may be arranged by contacting Red Hat, Inc.
// at http://sources.redhat.com/ecos/ecos-license/
// -------------------------------------------
//####ECOSGPLCOPYRIGHTEND####
//===========================================================================
//#####DESCRIPTIONBEGIN####
//
// Author(s): jlarmour
// Contributors: jlarmour
// Date: 1998-02-13
// Purpose:
// Description:
// Usage:
//
//####DESCRIPTIONEND####
//
//===========================================================================
/* ############################################################################################## */
/* 803A2588-803A25F0 02EBE8 0068+00 1/1 0/0 0/0 .rodata T */
SECTION_RODATA static u8 const T[104] = {
0x3F, 0xD5, 0x55, 0x55, 0x55, 0x55, 0x55, 0x63, 0x3F, 0xC1, 0x11, 0x11, 0x11, 0x10, 0xFE,
0x7A, 0x3F, 0xAB, 0xA1, 0xBA, 0x1B, 0xB3, 0x41, 0xFE, 0x3F, 0x96, 0x64, 0xF4, 0x84, 0x06,
0xD6, 0x37, 0x3F, 0x82, 0x26, 0xE3, 0xE9, 0x6E, 0x84, 0x93, 0x3F, 0x6D, 0x6D, 0x22, 0xC9,
0x56, 0x03, 0x28, 0x3F, 0x57, 0xDB, 0xC8, 0xFE, 0xE0, 0x83, 0x15, 0x3F, 0x43, 0x44, 0xD8,
0xF2, 0xF2, 0x65, 0x01, 0x3F, 0x30, 0x26, 0xF7, 0x1A, 0x8D, 0x10, 0x68, 0x3F, 0x14, 0x7E,
0x88, 0xA0, 0x37, 0x92, 0xA6, 0x3F, 0x12, 0xB8, 0x0F, 0x32, 0xF0, 0xA7, 0xE9, 0xBE, 0xF3,
0x75, 0xCB, 0xDB, 0x60, 0x53, 0x73, 0x3E, 0xFB, 0x2A, 0x70, 0x74, 0xBF, 0x7A, 0xD4,
};
COMPILER_STRIP_GATE(0x803A2588, &T);
// Derived from code with the following copyright
/* 80456A80-80456A88 005080 0008+00 1/1 0/0 0/0 .sdata2 @94 */
SECTION_SDATA2 static f64 lit_94 = 1.0;
/* @(#)k_tan.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*/
/* 80456A88-80456A90 005088 0008+00 1/1 0/0 0/0 .sdata2 @95 */
SECTION_SDATA2 static f64 lit_95 = -1.0;
/* __kernel_tan( x, y, k )
* kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input k indicates whether tan (if k=1) or
* -1/tan (if k= -1) is returned.
*
* Algorithm
* 1. Since tan(-x) = -tan(x), we need only to consider positive x.
* 2. if x < 2^-28 (hx<0x3e300000 0), return x with inexact if x!=0.
* 3. tan(x) is approximated by a odd polynomial of degree 27 on
* [0,0.67434]
* 3 27
* tan(x) ~ x + T1*x + ... + T13*x
* where
*
* |tan(x) 2 4 26 | -59.2
* |----- - (1+T1*x +T2*x +.... +T13*x )| <= 2
* | x |
*
* Note: tan(x+y) = tan(x) + tan'(x)*y
* ~ tan(x) + (1+x*x)*y
* Therefore, for better accuracy in computing tan(x+y), let
* 3 2 2 2 2
* r = x *(T2+x *(T3+x *(...+x *(T12+x *T13))))
* then
* 3 2
* tan(x+y) = x + (T1*x + (x *(r+y)+y))
*
* 4. For x in [0.67434,pi/4], let y = pi/4 - x, then
* tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
* = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
*/
/* 80456A90-80456A98 005090 0008+00 1/1 0/0 0/0 .sdata2 @96 */
SECTION_SDATA2 static f64 lit_96 = 0.7853981633974483;
#include "fdlibm.h"
#include "MSL_C/math.h"
static const double one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pio4 = 7.85398163397448278999e-01, /* 0x3FE921FB, 0x54442D18 */
pio4lo = 3.06161699786838301793e-17, /* 0x3C81A626, 0x33145C07 */
T[] = {
3.33333333333334091986e-01, /* 0x3FD55555, 0x55555563 */
1.33333333333201242699e-01, /* 0x3FC11111, 0x1110FE7A */
5.39682539762260521377e-02, /* 0x3FABA1BA, 0x1BB341FE */
2.18694882948595424599e-02, /* 0x3F9664F4, 0x8406D637 */
8.86323982359930005737e-03, /* 0x3F8226E3, 0xE96E8493 */
3.59207910759131235356e-03, /* 0x3F6D6D22, 0xC9560328 */
1.45620945432529025516e-03, /* 0x3F57DBC8, 0xFEE08315 */
5.88041240820264096874e-04, /* 0x3F4344D8, 0xF2F26501 */
2.46463134818469906812e-04, /* 0x3F3026F7, 0x1A8D1068 */
7.81794442939557092300e-05, /* 0x3F147E88, 0xA03792A6 */
7.14072491382608190305e-05, /* 0x3F12B80F, 0x32F0A7E9 */
-1.85586374855275456654e-05, /* 0xBEF375CB, 0xDB605373 */
2.59073051863633712884e-05, /* 0x3EFB2A70, 0x74BF7AD4 */
};
/* 80456A98-80456AA0 005098 0008+00 1/1 0/0 0/0 .sdata2 @97 */
SECTION_SDATA2 static f64 lit_97 = 3.061616997868383e-17;
/* 80456AA0-80456AA8 0050A0 0008+00 1/1 0/0 0/0 .sdata2 @98 */
SECTION_SDATA2 static u8 lit_98[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
/* 80456AA8-80456AB0 0050A8 0008+00 1/1 0/0 0/0 .sdata2 @99 */
SECTION_SDATA2 static f64 lit_99 = 2.0;
/* 80456AB0-80456AB8 0050B0 0008+00 1/1 0/0 0/0 .sdata2 @101 */
SECTION_SDATA2 static f64 lit_101 = 4503601774854144.0 /* cast s32 to float */;
/* 8036BA90-8036BCA4 3663D0 0214+00 0/0 1/1 0/0 .text __kernel_tan */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void __kernel_tan() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/k_tan/__kernel_tan.s"
double __kernel_tan(double x, double y, int iy)
{
double z, r, v, w, s;
int ix, hx;
hx = __HI(x); /* high word of x */
ix = hx & 0x7fffffff; /* high word of |x| */
if (ix < 0x3e300000) /* x < 2**-28 */
{
if ((int)x == 0) { /* generate inexact */
if (((ix | __LO(x)) | (iy + 1)) == 0) {
double ret = fabs(x);
return one / ret;
} else
return (iy == 1) ? x : -one / x;
}
}
if (ix >= 0x3FE59428) { /* |x|>=0.6744 */
if (hx < 0) {
x = -x;
y = -y;
}
z = pio4 - x;
w = pio4lo - y;
x = z + w;
y = 0.0;
}
z = x * x;
w = z * z;
/* Break x^5*(T[1]+x^2*T[2]+...) into
* x^5(T[1]+x^4*T[3]+...+x^20*T[11]) +
* x^5(x^2*(T[2]+x^4*T[4]+...+x^22*[T12]))
*/
r = T[1] + w * (T[3] + w * (T[5] + w * (T[7] + w * (T[9] + w * T[11]))));
v = z * (T[2] + w * (T[4] + w * (T[6] + w * (T[8] + w * (T[10] + w * T[12])))));
s = z * x;
r = y + z * (s * (r + v) + y);
r += T[0] * s;
w = x + r;
if (ix >= 0x3FE59428) {
v = (double)iy;
return (double)(1 - ((hx >> 30) & 2)) * (v - 2.0 * (x - (w * w / (w + v) - r)));
}
if (iy == 1)
return w;
else { /* if allow error up to 2 ulp,
simply return -1.0/(x+r) here */
/* compute -1.0/(x+r) accurately */
double a, t;
z = w;
__LO(z) = 0;
v = r - (z - x); /* z+v = r+x */
t = a = -1.0 / w; /* a = -1.0/w */
__LO(t) = 0;
s = 1.0 + t * z;
return t + a * (s + t * v);
}
}
#pragma pop
// EOF k_tan.c
libs/MSL_C/Math/Double_precision/s_atan.c
+131 -67
View File
@@ -1,79 +1,143 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/s_atan
//
#include "MSL_C/Math/Double_precision/s_atan.h"
#include "dol2asm.h"
#include "dolphin/types.h"
/* @(#)s_atan.c 1.3 95/01/18 */
/*
* ====================================================
* 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.
* ====================================================
*
*/
//
// Forward References:
//
/* 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.
*/
void atan();
#include "fdlibm.h"
//
// External References:
//
//
// Declarations:
//
/* ############################################################################################## */
/* 803A25F0-803A2610 02EC50 0020+00 1/1 0/0 0/0 .rodata atanhi */
SECTION_RODATA static u8 const atanhi[32] = {
0x3F, 0xDD, 0xAC, 0x67, 0x05, 0x61, 0xBB, 0x4F, 0x3F, 0xE9, 0x21, 0xFB, 0x54, 0x44, 0x2D, 0x18,
0x3F, 0xEF, 0x73, 0x0B, 0xD2, 0x81, 0xF6, 0x9B, 0x3F, 0xF9, 0x21, 0xFB, 0x54, 0x44, 0x2D, 0x18,
#ifdef __STDC__
static const double atanhi[] = {
#else
static double atanhi[] = {
#endif
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 */
};
COMPILER_STRIP_GATE(0x803A25F0, &atanhi);
/* 803A2610-803A2630 02EC70 0020+00 0/1 0/0 0/0 .rodata atanlo */
#pragma push
#pragma force_active on
SECTION_RODATA static u8 const atanlo[32] = {
0x3C, 0x7A, 0x2B, 0x7F, 0x22, 0x2F, 0x65, 0xE2, 0x3C, 0x81, 0xA6, 0x26, 0x33, 0x14, 0x5C, 0x07,
0x3C, 0x70, 0x07, 0x88, 0x7A, 0xF0, 0xCB, 0xBD, 0x3C, 0x91, 0xA6, 0x26, 0x33, 0x14, 0x5C, 0x07,
#ifdef __STDC__
static const double atanlo[] = {
#else
static double atanlo[] = {
#endif
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 */
};
COMPILER_STRIP_GATE(0x803A2610, &atanlo);
#pragma pop
/* 803A2630-803A2688 02EC90 0058+00 0/1 0/0 0/0 .rodata aT */
#pragma push
#pragma force_active on
SECTION_RODATA static u8 const aT[88] = {
0x3F, 0xD5, 0x55, 0x55, 0x55, 0x55, 0x55, 0x0D, 0xBF, 0xC9, 0x99, 0x99, 0x99, 0x98, 0xEB,
0xC4, 0x3F, 0xC2, 0x49, 0x24, 0x92, 0x00, 0x83, 0xFF, 0xBF, 0xBC, 0x71, 0xC6, 0xFE, 0x23,
0x16, 0x71, 0x3F, 0xB7, 0x45, 0xCD, 0xC5, 0x4C, 0x20, 0x6E, 0xBF, 0xB3, 0xB0, 0xF2, 0xAF,
0x74, 0x9A, 0x6D, 0x3F, 0xB1, 0x0D, 0x66, 0xA0, 0xD0, 0x3D, 0x51, 0xBF, 0xAD, 0xDE, 0x2D,
0x52, 0xDE, 0xFD, 0x9A, 0x3F, 0xA9, 0x7B, 0x4B, 0x24, 0x76, 0x0D, 0xEB, 0xBF, 0xA2, 0xB4,
0x44, 0x2C, 0x6A, 0x6C, 0x2F, 0x3F, 0x90, 0xAD, 0x3A, 0xE3, 0x22, 0xDA, 0x11,
#ifdef __STDC__
static const double aT[] = {
#else
static double aT[] = {
#endif
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 */
};
COMPILER_STRIP_GATE(0x803A2630, &aT);
#pragma pop
/* 80456AB8-80456AC0 0050B8 0008+00 1/1 0/0 0/0 .sdata2 @115 */
SECTION_SDATA2 static f64 lit_115 = 1e+300;
#ifdef __STDC__
static const double
#else
static double
#endif
one
= 1.0,
huge = 1.0e300;
/* 80456AC0-80456AC8 0050C0 0008+00 1/1 0/0 0/0 .sdata2 @116 */
SECTION_SDATA2 static f64 lit_116 = 1.0;
#ifdef __STDC__
double atan(double x)
#else
double atan(x) double x;
#endif
{
double w, s1, s2, z;
int ix, hx, id;
/* 80456AC8-80456AD0 0050C8 0008+00 1/1 0/0 0/0 .sdata2 @117 */
SECTION_SDATA2 static f64 lit_117 = 2.0;
/* 80456AD0-80456AD8 0050D0 0008+00 1/1 0/0 0/0 .sdata2 @118 */
SECTION_SDATA2 static f64 lit_118 = 1.5;
/* 80456AD8-80456AE0 0050D8 0008+00 1/1 0/0 0/0 .sdata2 @119 */
SECTION_SDATA2 static f64 lit_119 = -1.0;
/* 8036BCA4-8036BEBC 3665E4 0218+00 0/0 2/2 0/0 .text atan */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void atan() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/s_atan/atan.s"
}
#pragma pop
hx = __HI(x);
ix = hx & 0x7fffffff;
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (__LO(x) != 0)))
return x + x; /* NaN */
if (hx > 0)
return atanhi[3] + atanlo[3];
else
return -atanhi[3] - atanlo[3];
}
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e200000) { /* |x| < 2^-29 */
if (huge + x > one)
return x; /* raise inexact */
}
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 - one) / (2.0 + x);
} else { /* 11/16<=|x|< 19/16 */
id = 1;
x = (x - one) / (x + one);
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2;
x = (x - 1.5) / (one + 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);
else {
z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
return (hx < 0) ? -z : z;
}
}
libs/MSL_C/Math/Double_precision/s_ceil.c
+53 -53
View File
@@ -34,57 +34,57 @@ double ceil(double x)
double ceil(x) double x;
#endif
{
int i0, i1, j0;
unsigned i, j;
i0 = __HI(x);
i1 = __LO(x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
if (j0 < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0) { /* return 0*sign(x) if |x|<1 */
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)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
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 = ((unsigned)(0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
if (i0 > 0) {
if (j0 == 20)
i0 += 1;
else {
j = i1 + (1 << (52 - j0));
if (j < i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
}
__HI(x) = i0;
__LO(x) = i1;
return x;
int i0, i1, j0;
unsigned i, j;
i0 = __HI(x);
i1 = __LO(x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
if (j0 < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0) { /* return 0*sign(x) if |x|<1 */
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)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
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 = ((unsigned)(0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
if (i0 > 0) {
if (j0 == 20)
i0 += 1;
else {
j = i1 + (1 << (52 - j0));
if (j < i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
}
__HI(x) = i0;
__LO(x) = i1;
return x;
}
libs/MSL_C/Math/Double_precision/s_ldexp.c
+55 -48
View File
@@ -1,52 +1,59 @@
//
// Generated By: dol2asm
// Translation Unit: Math/Double_precision/s_ldexp
//
/* @(#)s_ldexp.c 1.2 95/01/04 */
/* $Id: s_ldexp.c,v 1.3.14.1 2002/01/31 15:24:14 ceciliar Exp $ */
/*
* ====================================================
* 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 "MSL_C/Math/Double_precision/s_ldexp.h"
#include "dol2asm.h"
#include "dolphin/types.h"
#include "fdlibm.h"
#include "MSL_C/math.h" /* for isfinite macro */
static const double
//
// Forward References:
//
two54
= 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
big = 1.0e+300, tiny = 1.0e-300;
void ldexp();
double ldexp(double x, int n)
{
s32 k, hx, lx; /*- cc 020130 -*/
if (!isfinite(x) || x == 0.0)
return x;
//
// External References:
//
void copysign();
//
// Declarations:
//
/* ############################################################################################## */
/* 80456B10-80456B18 005110 0008+00 1/1 0/0 0/0 .sdata2 @91 */
SECTION_SDATA2 static u8 lit_91[8] = {
0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
};
/* 80456B18-80456B20 005118 0008+00 1/1 0/0 0/0 .sdata2 @92 */
SECTION_SDATA2 static f64 lit_92 = 1.8014398509481984e+16;
/* 80456B20-80456B28 005120 0008+00 1/1 0/0 0/0 .sdata2 @93 */
SECTION_SDATA2 static f64 lit_93 = 1e-300;
/* 80456B28-80456B30 005128 0008+00 1/1 0/0 0/0 .sdata2 @94 */
SECTION_SDATA2 static f64 lit_94 = 1e+300;
/* 80456B30-80456B38 005130 0008+00 1/1 0/0 0/0 .sdata2 @95 */
SECTION_SDATA2 static f64 lit_95 = 5.551115123125783e-17;
/* 8036C2D0-8036C494 366C10 01C4+00 0/0 3/3 0/0 .text ldexp */
#pragma push
#pragma optimization_level 0
#pragma optimizewithasm off
asm void ldexp() {
nofralloc
#include "asm/MSL_C/Math/Double_precision/s_ldexp/ldexp.s"
}
#pragma pop
hx = __HI(x);
lx = __LO(x);
k = (hx & 0x7ff00000) >> 20; /* extract exponent */
if (k == 0) { /* 0 or subnormal x */
if ((lx | (hx & 0x7fffffff)) == 0)
return x; /* +-0 */
x *= two54;
hx = __HI(x);
k = ((hx & 0x7ff00000) >> 20) - 54;
if (n < -50000)
return tiny * x; /*underflow*/
}
if (k == 0x7ff)
return x + x; /* NaN or Inf */
k = k + n;
if (k > 0x7fe)
return big * copysign(big, x); /* overflow */
if (k > 0) /* normal result */
{
__HI(x) = (hx & 0x800fffff) | (k << 20);
return x;
}
if (k <= -54)
if (n > 50000) /* in case integer overflow in n+k */
return big * copysign(big, x); /*overflow*/
else
return tiny * copysign(tiny, x); /*underflow*/
k += 54; /* subnormal result */
__HI(x) = (hx & 0x800fffff) | (k << 20);
return x * twom54;
}