Skip to content

refactor: sinh overflows in non-nearest rounding #1421

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Jul 30, 2020
Merged

refactor: sinh overflows in non-nearest rounding #1421

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 13 additions & 15 deletions std/assembly/math.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,12 +63,13 @@ function R(z: f64): f64 { // Rational approximation of (asin(x)-x)/x^3
/** @internal */
// @ts-ignore: decorator
@inline
function expo2(x: f64): f64 { // exp(x)/2 for x >= log(DBL_MAX)
function expo2(x: f64, sign: f64): f64 { // exp(x)/2 for x >= log(DBL_MAX)
const // see: musl/src/math/__expo2.c
k = <u32>2043,
kln2 = reinterpret<f64>(0x40962066151ADD8B); // 0x1.62066151add8bp+10
var scale = reinterpret<f64>(<u64>((<u32>0x3FF + k / 2) << 20) << 32);
return NativeMath.exp(x - kln2) * scale * scale;
// in directed rounding correct sign before rounding or overflow is important
return NativeMath.exp(x - kln2) * (sign * scale) * scale;
}

/** @internal */
Expand Down Expand Up @@ -759,7 +760,7 @@ export namespace NativeMath {
t = exp(x);
return 0.5 * (t + 1 / t);
}
t = expo2(x);
t = expo2(x, 1);
return t;
}

Expand Down Expand Up @@ -1470,18 +1471,16 @@ export namespace NativeMath {
var u = reinterpret<u64>(x) & 0x7FFFFFFFFFFFFFFF;
var absx = reinterpret<f64>(u);
var w = <u32>(u >> 32);
var t: f64;
var h = builtin_copysign(0.5, x);
if (w < 0x40862E42) {
t = expm1(absx);
let t = expm1(absx);
if (w < 0x3FF00000) {
if (w < 0x3FF00000 - (26 << 20)) return x;
return h * (2 * t - t * t / (t + 1));
}
return h * (t + t / (t + 1));
}
t = 2 * h * expo2(absx);
return t;
return expo2(absx, 2 * h);
}

// @ts-ignore: decorator
Expand Down Expand Up @@ -1760,12 +1759,13 @@ function Rf(z: f32): f32 { // Rational approximation of (asin(x)-x)/x^3

// @ts-ignore: decorator
@inline
function expo2f(x: f32): f32 { // exp(x)/2 for x >= log(DBL_MAX)
function expo2f(x: f32, sign: f32): f32 { // exp(x)/2 for x >= log(DBL_MAX)
const // see: musl/src/math/__expo2f.c
k = <u32>235,
kln2 = reinterpret<f32>(0x4322E3BC); // 0x1.45c778p+7f
var scale = reinterpret<f32>(<u32>(0x7F + (k >> 1)) << 23);
return NativeMathf.exp(x - kln2) * scale * scale;
// in directed rounding correct sign before rounding or overflow is important
return NativeMathf.exp(x - kln2) * (sign * scale) * scale;
}

// @ts-ignore: decorator
Expand Down Expand Up @@ -2253,7 +2253,7 @@ export namespace NativeMathf {
// return 0.5 * (t + 1 / t);
return 0.5 * t + 0.5 / t;
}
return expo2f(x);
return expo2f(x, 1);
}

// @ts-ignore: decorator
Expand Down Expand Up @@ -2758,18 +2758,16 @@ export namespace NativeMathf {
export function sinh(x: f32): f32 { // see: musl/src/math/sinhf.c
var u = reinterpret<u32>(x) & 0x7FFFFFFF;
var absx = reinterpret<f32>(u);
var t: f32;
var h = builtin_copysign<f32>(0.5, x);
if (u < 0x42B17217) {
t = expm1(absx);
let t = expm1(absx);
if (u < 0x3F800000) {
if (u < 0x3F800000 - (12 << 23)) return x;
return h * (2 * t - t * t / (t + 1));
}
return h * (t + t / (t + 1));
}
t = 2 * h * expo2f(absx);
return t;
return expo2f(absx, 2 * h);
}

// @ts-ignore: decorator
Expand Down Expand Up @@ -3225,4 +3223,4 @@ switch (e) {
...
case 63:
}
*/
*/
2 changes: 1 addition & 1 deletion tests/compiler/std/array.optimized.wat
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4236,7 +4236,7 @@
if
i32.const 0
i32.const 6464
i32.const 1398
i32.const 1399
i32.const 5
call $~lib/builtins/abort
unreachable
Expand Down
2 changes: 1 addition & 1 deletion tests/compiler/std/array.untouched.wat
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7415,7 +7415,7 @@
if
i32.const 0
i32.const 5456
i32.const 1398
i32.const 1399
i32.const 5
call $~lib/builtins/abort
unreachable
Expand Down
18 changes: 9 additions & 9 deletions tests/compiler/std/math.optimized.wat
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8432,7 +8432,7 @@
if
i32.const 0
i32.const 3616
i32.const 1398
i32.const 1399
i32.const 5
call $~lib/builtins/abort
unreachable
Expand Down Expand Up @@ -9672,18 +9672,18 @@
f64.mul
return
end
f64.const 2
local.get $2
f64.mul
local.get $1
f64.const 1416.0996898839683
f64.sub
call $~lib/math/NativeMath.exp
f64.const 2247116418577894884661631e283
f64.const 2
local.get $2
f64.mul
f64.const 2247116418577894884661631e283
f64.mul
f64.mul
f64.const 2247116418577894884661631e283
f64.mul
)
(func $std/math/test_sinh (param $0 f64) (param $1 f64) (param $2 f64) (result i32)
local.get $0
Expand Down Expand Up @@ -9760,18 +9760,18 @@
f32.mul
return
end
f32.const 2
local.get $3
f32.mul
local.get $1
f32.const 162.88958740234375
f32.sub
call $~lib/math/NativeMathf.exp
f32.const 1661534994731144841129758e11
f32.const 2
local.get $3
f32.mul
f32.const 1661534994731144841129758e11
f32.mul
f32.mul
f32.const 1661534994731144841129758e11
f32.mul
)
(func $std/math/test_sinhf (param $0 f32) (param $1 f32) (param $2 f32) (result i32)
local.get $0
Expand Down
Loading