📄 math32lb.h
字号:
c = arg1f32; // save x
arg1f32 *= 3.32192809489; // 1/log10(2) = 3.32192809489
arg1f32 += 0.5; // k = [ x / log10(2) + .5 ]
arg1f32 = floor32( arg1f32);
float32 e = arg1f32; // save float k
FpFlags &= ~0x40;
// floating point to integer conversion
xexp = arg1f32; // k = [ x / ln(2) + .5 ]
FpFlags |= 0x40;
arg1f32 = e * -0.301025390625; // c1
d = arg1f32 + c;
arg1f32 = e * -4.6050389811952113E-6; // c2
arg1f32 += d;
d = arg1f32; // save f
if (!(d.midH8 & 0x80)) {
// minimax approximation on [0,log10(2)/2]
arg1f32 *= 6.388992868121E-1;/* EXP1032H5 */
arg1f32 += 1.154596329197E0; /* EXP1032H4 */
arg1f32 *= d;
arg1f32 += 2.035920309947E0; /* EXP1032H3 */
arg1f32 *= d;
arg1f32 += 2.650909138708E0; /* EXP1032H2 */
arg1f32 *= d;
arg1f32 += 2.302585504840E0; /* EXP1032H1 */
}
else {
// minimax approximation on [-log10(2)/2,0]
arg1f32 *= 4.544952589676E-1;/* EXP1032L5 */
arg1f32 += 1.157459289066E0; /* EXP1032L4 */
arg1f32 *= d;
arg1f32 += 2.033640565225E0; /* EXP1032L3 */
arg1f32 *= d;
arg1f32 += 2.650914554552E0; /* EXP1032L2 */
arg1f32 *= d;
arg1f32 += 2.302584716116E0; /* EXP1032L1 */
}
arg1f32 *= d;
if (!(savedFlags & 0x40))
FpFlags &= ~0x40;
arg1f32 += 1.0; /* EXP1032H0/EXP1032L0 */
arg1f32.high8 += xexp;
goto RETURN;
EXP1:
arg1f32 = 1.0; // return 10**x = 1.0
goto RETURN;
DOMERR:
FpDomainError = 1; // domain error
RETURN:
return arg1f32;
}
float32 exp( sharedM float32 arg1f32)
// Maximum argument : 88.7228391117 = log(2**128)
// Minimum argument : -87.3365447506 = log(2**-126)
{
sharedM float32 arg2f32;
float32 c, d;
char xexp;
if (( arg1f32.high8 - 94) & 0x80)
goto EXP1; // return e**x = 1
W = 133 - arg1f32.high8;
if (!Carry)
goto DOMERR;
if (!Zero_)
goto ARGOK;
if (!(arg1f32.midH8 & 0x80)) {
W = 49 - arg1f32.midH8;
if (!Carry)
goto DOMERR;
if (!Zero_)
goto ARGOK;
W = 114 - arg1f32.midL8;
if (!Carry)
goto DOMERR;
if (!Zero_)
goto ARGOK;
W = 24 - arg1f32.low8;
if (!Carry)
goto DOMERR;
}
else {
W = 174 - arg1f32.midH8;
if (!Carry)
goto DOMERR;
if (!Zero_)
goto ARGOK;
W = 172 - arg1f32.midL8;
if (!Carry)
goto DOMERR;
if (!Zero_)
goto ARGOK;
W = 80 - arg1f32.low8;
if (!Carry)
goto DOMERR;
}
ARGOK:
char savedFlags = FpFlags;
FpFlags |= 0x40; // enable rounding
// Range reduction
c = arg1f32; // save x
arg1f32 *= 1.44269504089; // 1/ln(2) = 1.44269504089
arg1f32 += 0.5; // k = [ x / ln(2) + .5 ]
arg1f32 = floor32( arg1f32);
float32 e = arg1f32; // save float k
// floating point to integer conversion
xexp = arg1f32; // k = [ x / ln(2) + .5 ]
arg1f32 = e * -0.693359375; // c1
d = arg1f32 + c;
arg1f32 = e * 0.00021219444005; // c2
arg1f32 += d;
d = arg1f32; // save f
if (!(d.midH8 & 0x80)) {
arg1f32 *= 0.989943653774E-2; /* EXP32H5 */
arg1f32 += 0.410473706887E-1; /* EXP32H4 */
arg1f32 *= d;
arg1f32 += 0.166777360103; /* EXP32H3 */
arg1f32 *= d;
arg1f32 += 0.499991163105; /* EXP32H2 */
arg1f32 *= d;
arg1f32 += 1.00000025499; /* EXP32H1 */
}
else {
arg1f32 *= 0.699995870637E-2; /* EXP32L5 */
arg1f32 += 0.411548782678E-1; /* EXP32L4 */
arg1f32 *= d;
arg1f32 += 0.166574299807; /* EXP32L3 */
arg1f32 *= d;
arg1f32 += 0.499992371926; /* EXP32L2 */
arg1f32 *= d;
arg1f32 += 0.999999766814; /* EXP32L1 */
}
arg1f32 *= d;
arg1f32 += 1.0; /* EXP32H0 or EXP32L0 */
arg1f32.high8 += xexp;
if (!(savedFlags & 0x40))
FpFlags &= ~0x40; // restore rounding flag
goto RETURN;
EXP1:
arg1f32 = 1.0; // return 10**x = 1.0
goto RETURN;
DOMERR:
FpDomainError = 1; // domain error
RETURN:
return arg1f32;
}
float32 cosin32( sharedM float32 arg1f32, sharedM float32 arg2f32, char cosinus)
{
float32 c, d;
char j;
char savedFlags = FpFlags; // save rounding flag
FpFlags |= 0x40; // enable rounding
char csign = 0; // initialize sign
if (!cosinus && (arg1f32.midH8 & 0x80))
csign |= 0x80;
arg1f32.midH8 &= ~0x80; // use |x|
// loss threshold check: arg1f32 <= +512.0, arg1f32 >= -512.0
if (arg1f32.high8 >= 0x88)
FpDomainError = 1; // domain error
c = arg1f32; // save |x|
arg1f32 *= 1.27323954474; // 4/pi
FpFlags &= ~0x40;
// y = [ |x| * (4/pi) ]
arg1f32.low24 = arg1f32; // floating point to integer conversion
FpFlags |= 0x40;
if (arg1f32.low8 & 0x1)
arg1f32.low24 += 1; // make arg1f32 even
// j = y mod 8
j = arg1f32.low8 & 7; // 0,2,4,6
if (j >= 4) {
csign ^= 128;
j -= 4;
}
arg1f32 = arg1f32.low24; // integer to floating point conversion
// save y in DARG
W = arg1f32.high8;
d.high8 = W;
#if __CoreSet__ / 100 == 17
if (W == 0) // NB17
goto ZEQX;
#else
if (Zero_)
goto ZEQX;
#endif
d.low24 = arg1f32.low24;
arg1f32 *= -7.853851318359375e-01;
c = arg1f32 + c; // z1 = |x| - y * (p1)
arg1f32 = d * -1.303153112530708e-05;
c = arg1f32 + c; // z2 = z1 - y * (p2)
arg1f32 = d * -3.038469476363188e-11;
arg1f32 += c; // z3 = z2 - y * (p3)
// save z in c
W = arg1f32.high8; c.high8 = W; arg2f32.high8 = W;
W = arg1f32.midH8; c.midH8 = W; arg2f32.midH8 = W;
W = arg1f32.midL8; c.midL8 = W; arg2f32.midL8 = W;
W = arg1f32.low8; c.low8 = W; arg2f32.low8 = W;
goto POLYNOM;
ZEQX:
W = c.high8; arg1f32.high8 = W; arg2f32.high8 = W;
W = c.midH8; arg1f32.midH8 = W; arg2f32.midH8 = W;
W = c.midL8; arg1f32.midL8 = W; arg2f32.midL8 = W;
W = c.low8; arg1f32.low8 = W; arg2f32.low8 = W;
POLYNOM:
arg1f32 *= arg2f32; // z * z
d = arg1f32; // save z * z
if ((( rr( j) ^ j) & 0x1) ^ cosinus) {
// minimax polynomial coefficients for
// cos(z) = 1 -.5*z**2 + z**4*q(z**2) on [-pi/4,pi/4]
arg1f32 *= 2.443315706066392E-5; /* COS32D2 */
arg1f32 += -1.388731625438419E-3; /* COS32D1 */
arg1f32 *= d;
arg1f32 += 4.166664568297614E-2; /* COS32D0 */
arg1f32 *= d;
arg1f32 *= d;
arg2f32 = d;
if (arg2f32.high8)
arg2f32.high8 -= 1; // arg2f32 *= 0.5
arg1f32 -= arg2f32;
FpFlags &= ~0x40;
arg1f32 += 1.0;
}
else {
// minimax polynomial coefficients for
// sin(z) = z+z*(z**2)*p(z**2) on [-pi/4,pi/4]
arg1f32 *= 2.718121647219611E-6; /* SIN32D3 */
arg1f32 += -1.983931227180460E-4; /* SIN32D2 */
arg1f32 *= d;
arg1f32 += 8.333329304850749E-3; /* SIN32D1 */
arg1f32 *= d;
arg1f32 += -1.666666664079712E-1; /* SIN32D0 */
arg1f32 *= d;
arg1f32 *= c;
FpFlags &= ~0x40;
arg1f32 += c;
}
W = 128; // LSB+1
if (cosinus && (j & 0x2)) // NB17
csign ^= W;
if (csign & 0x80)
arg1f32.midH8 ^= W;
if (savedFlags & 0x40)
FpFlags |= 0x40; // restore rounding flag
return arg1f32;
}
float32 cos( sharedM float32 arg1f32)
{
sharedM float32 arg2f32; // allocation 'trick'
arg1f32 = cosin32( arg1f32, arg2f32, 1);
return arg1f32;
}
float32 sin( sharedM float32 arg1f32)
{
sharedM float32 arg2f32; // allocation 'trick'
arg1f32 = cosin32( arg1f32, arg2f32, 0);
return arg1f32;
}
float32 sqrt( sharedM float32 arg1f32)
{
sharedM float32 arg2f32; // allocation 'trick'
if (arg1f32.midH8 & 0x80) // test for negative argument
goto DOMERR;
if (!arg1f32.high8) // return if argument zero
goto RETURN;
char cexp = arg1f32.high8; // save exponent in cexp
char savedFLAGS = FpFlags;
FpFlags |= 0x40; // enable rounding
arg1f32.high8 = FpBIAS; // compute z
arg1f32 -= 1.0; // z = f - 1.0
float32 d = arg1f32;
// Q(z)
arg1f32 += 3.294831307E1; /* SQRT32Q2 */
arg1f32 *= d;
arg1f32 += 1.333544312E2; /* SQRT32Q1 */
arg1f32 *= d;
arg1f32 += 1.210947497E2; /* SQRT32Q0 */
float32 e = arg1f32; // save Q(z)
arg1f32 = d; // restore z
// P(z)
arg1f32 *= 7.370062896E0; /* SQRT32P2 */
arg1f32 += 5.154073142E1; /* SQRT32P1 */
arg1f32 *= d;
arg1f32 += 6.054736157E1; /* SQRT32P0 */
arg1f32 /= e; // P(z)/Q(z)
arg1f32 *= d; // z*P(z)/Q(z)
arg1f32 += 1.0; // sqrt(1+z)=1+z*P(z)/Q(z)
if (!(cexp & 0x1)) // is cexp even or odd?
arg1f32 *= 1.41421356237; // sqrt(2)
// divide exponent by two
cexp += 127;
arg1f32.high8 = rr( cexp);
if (!(savedFLAGS & 0x40))
FpFlags &= ~0x40;
goto RETURN;
DOMERR:
FpDomainError = 1; // domain error
RETURN:
return arg1f32;
}
#pragma library 0
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -