maelstrom/milsko
1
0
Fork 0
forked from pyrite-dev/milsko
Code Pull requests Activity

format

Browse source 
git-svn-id: http://svn2.nishi.boats/svn/milsko/trunk@550 b9cfdab3-6d41-4d17-bbe4-086880011989
This commit is contained in:
Nishi 2025-11-01 08:05:32 +00:00
commit 22e75cc25d
7 changed files with 290 additions and 277 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/math/nbsd_math.h
View file

@ -20,8 +20,7 @@
#include <Mw/BaseTypes.h> #include <Mw/BaseTypes.h>
typedef union _ieee_double_shape_type typedef union _ieee_double_shape_type {
{
double value; double value;
struct struct
{ {
@ -82,6 +81,4 @@ do { \
(d) = sl_u.value; \ (d) = sl_u.value; \
} while(0) } while(0)
#endif #endif

52
src/math/nbsd_pow.c
View file

@ -60,15 +60,20 @@
*/ */
static const double static const double
bp[] = {1.0, 1.5,}, bp[] = {
dp_h[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */ 1.0,
dp_l[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */ 1.5,
zero = 0.0, },
one = 1.0, dp_h[] = {
two = 2.0, 0.0,
two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */ 5.84962487220764160156e-01,
hugev = 1.0e300, }, /* 0x3FE2B803, 0x40000000 */
tinyv = 1.0e-300, dp_l[] = {
0.0,
1.35003920212974897128e-08,
}, /* 0x3E4CFDEB, 0x43CFD006 */
zero = 0.0, one = 1.0, two = 2.0, two53 = 9007199254740992.0, /* 0x43400000, 0x00000000 */
hugev = 1.0e300, tinyv = 1.0e-300,
/* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */ /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */ L1 = 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */ L2 = 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
@ -93,8 +98,7 @@ ivln2_h = 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/ ivln2_l = 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
double double
nbsd_pow(double x, double y) nbsd_pow(double x, double y) {
{
double z, ax, z_h, z_l, p_h, p_l; double z, ax, z_h, z_l, p_h, p_l;
double yy1, t1, t2, r, s, t, u, v, w; double yy1, t1, t2, r, s, t, u, v, w;
MwI32 i, j, k, yisint, n; MwI32 i, j, k, yisint, n;
@ -103,7 +107,8 @@ nbsd_pow(double x, double y)
EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hx, lx, x);
EXTRACT_WORDS(hy, ly, y); EXTRACT_WORDS(hy, ly, y);
ix = hx&0x7fffffff; iy = hy&0x7fffffff; ix = hx & 0x7fffffff;
iy = hy & 0x7fffffff;
/* y==zero: x**0 = 1 */ /* y==zero: x**0 = 1 */
if((iy | ly) == 0) return one; if((iy | ly) == 0) return one;
@ -147,7 +152,9 @@ nbsd_pow(double x, double y)
return (hy < 0) ? -y : zero; return (hy < 0) ? -y : zero;
} }
if(iy == 0x3ff00000) { /* y is +-1 */ if(iy == 0x3ff00000) { /* y is +-1 */
if(hy<0) return one/x; else return x; if(hy < 0) return one / x;
else
return x;
} }
if(hy == 0x40000000) return x * x; /* y is 2 */ if(hy == 0x40000000) return x * x; /* y is 2 */
if(hy == 0x3fe00000) { /* y is 0.5 */ if(hy == 0x3fe00000) { /* y is 0.5 */
@ -206,15 +213,23 @@ nbsd_pow(double x, double y)
double ss, s2, s_h, s_l, t_h, t_l; double ss, s2, s_h, s_l, t_h, t_l;
n = 0; n = 0;
/* take care subnormal number */ /* take care subnormal number */
if(ix<0x00100000) if(ix < 0x00100000) {
{ax *= two53; n -= 53; GET_HIGH_WORD(ix,ax); } ax *= two53;
n -= 53;
GET_HIGH_WORD(ix, ax);
}
n += ((ix) >> 20) - 0x3ff; n += ((ix) >> 20) - 0x3ff;
j = ix & 0x000fffff; j = ix & 0x000fffff;
/* determine interval */ /* determine interval */
ix = j | 0x3ff00000; /* normalize ix */ ix = j | 0x3ff00000; /* normalize ix */
if(j <= 0x3988E) k = 0; /* |x|<sqrt(3/2) */ if(j <= 0x3988E) k = 0; /* |x|<sqrt(3/2) */
else if(j<0xBB67A) k=1; /* |x|<sqrt(3) */ else if(j < 0xBB67A)
else {k=0;n+=1;ix -= 0x00100000;} k = 1; /* |x|<sqrt(3) */
else {
k = 0;
n += 1;
ix -= 0x00100000;
}
SET_HIGH_WORD(ax, ix); SET_HIGH_WORD(ax, ix);
/* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */ /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
@ -300,6 +315,7 @@ nbsd_pow(double x, double y)
GET_HIGH_WORD(j, z); GET_HIGH_WORD(j, z);
j += (n << 20); j += (n << 20);
if((j >> 20) <= 0) z = scalbn(z, n); /* subnormal output */ if((j >> 20) <= 0) z = scalbn(z, n); /* subnormal output */
else SET_HIGH_WORD(z,j); else
SET_HIGH_WORD(z, j);
return s * z; return s * z;
} }