Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
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
791 changes: 457 additions & 334 deletions core/math/math.odin

Large diffs are not rendered by default.

125 changes: 88 additions & 37 deletions core/math/math_basic.odin
Original file line number Diff line number Diff line change
Expand Up @@ -42,88 +42,139 @@ foreign _ {
}

@(require_results)
sin_f16 :: proc "contextless" (θ: f16) -> f16 {
sin_f16 :: proc "contextless" (θ: f16, loc := #caller_location) -> f16 {
validate_finite(θ, loc)
return _sin_f16(θ)
}
@(require_results)
sin_f32 :: proc "contextless" (θ: f32) -> f32 {
sin_f32 :: proc "contextless" (θ: f32, loc := #caller_location) -> f32 {
validate_finite(θ, loc)
return _sin_f32(θ)
}
@(require_results)
sin_f64 :: proc "contextless" (θ: f64) -> f64 {
sin_f64 :: proc "contextless" (θ: f64, loc := #caller_location) -> f64 {
validate_finite(θ, loc)
return _sin_f64(θ)
}

@(require_results)
cos_f16 :: proc "contextless" (θ: f16) -> f16 {
cos_f16 :: proc "contextless" (θ: f16, loc := #caller_location) -> f16 {
validate_finite(θ, loc)
return _cos_f16(θ)
}
@(require_results)
cos_f32 :: proc "contextless" (θ: f32) -> f32 {
cos_f32 :: proc "contextless" (θ: f32, loc := #caller_location) -> f32 {
validate_finite(θ, loc)
return _cos_f32(θ)
}
@(require_results)
cos_f64 :: proc "contextless" (θ: f64) -> f64 {
cos_f64 :: proc "contextless" (θ: f64, loc := #caller_location) -> f64 {
validate_finite(θ, loc)
return _cos_f64(θ)
}

@(require_results)
pow_f16 :: proc "contextless" (x, power: f16) -> f16 {
return _pow_f16(x, power)
pow_f16 :: proc "contextless" (x, power: f16, loc := #caller_location) -> f16 {
validate_finite(x, loc)
validate_finite(power, loc)
result := _pow_f16(x, power)
validate_finite(result, loc)
return result
}
@(require_results)
pow_f32 :: proc "contextless" (x, power: f32) -> f32 {
return _pow_f32(x, power)
pow_f32 :: proc "contextless" (x, power: f32, loc := #caller_location) -> f32 {
validate_finite(x, loc)
validate_finite(power, loc)
result := _pow_f32(x, power)
validate_finite(result, loc)
return result
}
@(require_results)
pow_f64 :: proc "contextless" (x, power: f64) -> f64 {
return _pow_f64(x, power)
pow_f64 :: proc "contextless" (x, power: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)
validate_finite(power, loc)
result := _pow_f64(x, power)
validate_finite(result, loc)
return result
}

@(require_results)
fmuladd_f16 :: proc "contextless" (a, b, c: f16) -> f16 {
return _fmuladd_f16(a, b, c)
fmuladd_f16 :: proc "contextless" (a, b, c: f16, loc := #caller_location) -> f16 {
validate_finite(a, loc)
validate_finite(b, loc)
validate_finite(c, loc)
result := _fmuladd_f16(a, b, c)
validate_finite(result, loc)
return result
}
@(require_results)
fmuladd_f32 :: proc "contextless" (a, b, c: f32) -> f32 {
return _fmuladd_f32(a, b, c)
fmuladd_f32 :: proc "contextless" (a, b, c: f32, loc := #caller_location) -> f32 {
validate_finite(a, loc)
validate_finite(b, loc)
validate_finite(c, loc)
result := _fmuladd_f32(a, b, c)
validate_finite(result, loc)
return result
}
@(require_results)
fmuladd_f64 :: proc "contextless" (a, b, c: f64) -> f64 {
return _fmuladd_f64(a, b, c)
fmuladd_f64 :: proc "contextless" (a, b, c: f64, loc := #caller_location) -> f64 {
validate_finite(a, loc)
validate_finite(b, loc)
validate_finite(c, loc)
result := _fmuladd_f64(a, b, c)
validate_finite(result, loc)
return result
}

@(require_results)
exp_f16 :: proc "contextless" (x: f16) -> f16 {
return _exp_f16(x)
exp_f16 :: proc "contextless" (x: f16, loc := #caller_location) -> f16 {
validate_finite(x, loc)
result := _exp_f16(x)
validate_finite(result, loc)
return result
}
@(require_results)
exp_f32 :: proc "contextless" (x: f32) -> f32 {
return _exp_f32(x)
exp_f32 :: proc "contextless" (x: f32, loc := #caller_location) -> f32 {
validate_finite(x, loc)
result := _exp_f32(x)
validate_finite(result, loc)
return result
}
@(require_results)
exp_f64 :: proc "contextless" (x: f64) -> f64 {
return _exp_f64(x)
exp_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)
result := _exp_f64(x)
validate_finite(result, loc)
return result
}


@(require_results)
sqrt_f16 :: proc "contextless" (x: f16) -> f16 {
sqrt_f16 :: proc "contextless" (x: f16, loc := #caller_location) -> f16 {
validate_finite(x, loc)
validation_assert(x >= 0, loc)
return intrinsics.sqrt(x)
}
@(require_results)
sqrt_f32 :: proc "contextless" (x: f32) -> f32 {
sqrt_f32 :: proc "contextless" (x: f32, loc := #caller_location) -> f32 {
validate_finite(x, loc)
validation_assert(x >= 0, loc)
return intrinsics.sqrt(x)
}
@(require_results)
sqrt_f64 :: proc "contextless" (x: f64) -> f64 {
sqrt_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)
validation_assert(x >= 0, loc)
return intrinsics.sqrt(x)
}



@(require_results)
ln_f64 :: proc "contextless" (x: f64) -> f64 {
ln_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)
validation_assert(x > 0, loc)

// The original C code, the long comment, and the constants
// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c
// and came with this notice.
Expand Down Expand Up @@ -205,7 +256,7 @@ ln_f64 :: proc "contextless" (x: f64) -> f64 {
}

// reduce
f1, ki := frexp(x)
f1, ki := frexp(x, loc)
if f1 < SQRT_TWO/2 {
f1 *= 2
ki -= 1
Expand All @@ -224,14 +275,14 @@ ln_f64 :: proc "contextless" (x: f64) -> f64 {
return k*LN2_HI - ((hfsq - (s*(hfsq+R) + k*LN2_LO)) - f)
}

@(require_results) ln_f16 :: proc "contextless" (x: f16) -> f16 { return #force_inline f16(ln_f64(f64(x))) }
@(require_results) ln_f32 :: proc "contextless" (x: f32) -> f32 { return #force_inline f32(ln_f64(f64(x))) }
@(require_results) ln_f16le :: proc "contextless" (x: f16le) -> f16le { return #force_inline f16le(ln_f64(f64(x))) }
@(require_results) ln_f16be :: proc "contextless" (x: f16be) -> f16be { return #force_inline f16be(ln_f64(f64(x))) }
@(require_results) ln_f32le :: proc "contextless" (x: f32le) -> f32le { return #force_inline f32le(ln_f64(f64(x))) }
@(require_results) ln_f32be :: proc "contextless" (x: f32be) -> f32be { return #force_inline f32be(ln_f64(f64(x))) }
@(require_results) ln_f64le :: proc "contextless" (x: f64le) -> f64le { return #force_inline f64le(ln_f64(f64(x))) }
@(require_results) ln_f64be :: proc "contextless" (x: f64be) -> f64be { return #force_inline f64be(ln_f64(f64(x))) }
@(require_results) ln_f16 :: proc "contextless" (x: f16, loc := #caller_location) -> f16 { return #force_inline f16(ln_f64(f64(x), loc)) }
@(require_results) ln_f32 :: proc "contextless" (x: f32, loc := #caller_location) -> f32 { return #force_inline f32(ln_f64(f64(x), loc)) }
@(require_results) ln_f16le :: proc "contextless" (x: f16le, loc := #caller_location) -> f16le { return #force_inline f16le(ln_f64(f64(x), loc)) }
@(require_results) ln_f16be :: proc "contextless" (x: f16be, loc := #caller_location) -> f16be { return #force_inline f16be(ln_f64(f64(x), loc)) }
@(require_results) ln_f32le :: proc "contextless" (x: f32le, loc := #caller_location) -> f32le { return #force_inline f32le(ln_f64(f64(x), loc)) }
@(require_results) ln_f32be :: proc "contextless" (x: f32be, loc := #caller_location) -> f32be { return #force_inline f32be(ln_f64(f64(x), loc)) }
@(require_results) ln_f64le :: proc "contextless" (x: f64le, loc := #caller_location) -> f64le { return #force_inline f64le(ln_f64(f64(x), loc)) }
@(require_results) ln_f64be :: proc "contextless" (x: f64be, loc := #caller_location) -> f64be { return #force_inline f64be(ln_f64(f64(x), loc)) }
ln :: proc{
ln_f16, ln_f16le, ln_f16be,
ln_f32, ln_f32le, ln_f32be,
Expand Down
81 changes: 53 additions & 28 deletions core/math/math_basic_js.odin
Original file line number Diff line number Diff line change
Expand Up @@ -8,46 +8,71 @@ foreign import "odin_env"
@(default_calling_convention="c")
foreign odin_env {
@(link_name="sin", require_results)
sin_f64 :: proc(θ: f64) -> f64 ---
_sin_f64 :: proc(θ: f64) -> f64 ---
@(link_name="cos", require_results)
cos_f64 :: proc(θ: f64) -> f64 ---
_cos_f64 :: proc(θ: f64) -> f64 ---
@(link_name="pow", require_results)
pow_f64 :: proc(x, power: f64) -> f64 ---
_pow_f64 :: proc(x, power: f64) -> f64 ---
@(link_name="fmuladd", require_results)
fmuladd_f64 :: proc(a, b, c: f64) -> f64 ---
_fmuladd_f64 :: proc(a, b, c: f64) -> f64 ---
@(link_name="ln", require_results)
ln_f64 :: proc(x: f64) -> f64 ---
_ln_f64 :: proc(x: f64) -> f64 ---
@(link_name="exp", require_results)
exp_f64 :: proc(x: f64) -> f64 ---
_exp_f64 :: proc(x: f64) -> f64 ---
}

@(require_results)
sin_f64 :: proc "contextless" (θ: f64, loc := #caller_location) -> f64 {
return _sin_f64(θ)
}
@(require_results)
cos_f64 :: proc "contextless" (θ: f64, loc := #caller_location) -> f64 {
return _cos_f64(θ)
}
@(require_results)
pow_f64 :: proc "contextless" (x, power: f64, loc := #caller_location) -> f64 {
return _pow_f64(x, power)
}
@(require_results)
fmuladd_f64 :: proc "contextless" (a, b, c: f64, loc := #caller_location) -> f64 {
return _fmuladd_f64(a, b, c)
}
@(require_results)
ln_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
return _ln_f64(x)
}
@(require_results)
exp_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
return _exp_f64(x)
}

@(require_results)
sqrt_f64 :: proc "contextless" (x: f64) -> f64 {
return intrinsics.sqrt(x)
}

@(require_results) sqrt_f16 :: proc "c" (x: f16) -> f16 { return f16(sqrt_f64(f64(x))) }
@(require_results) sin_f16 :: proc "c" (θ: f16) -> f16 { return f16(sin_f64(f64(θ))) }
@(require_results) cos_f16 :: proc "c" (θ: f16) -> f16 { return f16(cos_f64(f64(θ))) }
@(require_results) pow_f16 :: proc "c" (x, power: f16) -> f16 { return f16(pow_f64(f64(x), f64(power))) }
@(require_results) fmuladd_f16 :: proc "c" (a, b, c: f16) -> f16 { return f16(fmuladd_f64(f64(a), f64(a), f64(c))) }
@(require_results) ln_f16 :: proc "c" (x: f16) -> f16 { return f16(ln_f64(f64(x))) }
@(require_results) exp_f16 :: proc "c" (x: f16) -> f16 { return f16(exp_f64(f64(x))) }

@(require_results) sqrt_f32 :: proc "c" (x: f32) -> f32 { return f32(sqrt_f64(f64(x))) }
@(require_results) sin_f32 :: proc "c" (θ: f32) -> f32 { return f32(sin_f64(f64(θ))) }
@(require_results) cos_f32 :: proc "c" (θ: f32) -> f32 { return f32(cos_f64(f64(θ))) }
@(require_results) pow_f32 :: proc "c" (x, power: f32) -> f32 { return f32(pow_f64(f64(x), f64(power))) }
@(require_results) fmuladd_f32 :: proc "c" (a, b, c: f32) -> f32 { return f32(fmuladd_f64(f64(a), f64(a), f64(c))) }
@(require_results) ln_f32 :: proc "c" (x: f32) -> f32 { return f32(ln_f64(f64(x))) }
@(require_results) exp_f32 :: proc "c" (x: f32) -> f32 { return f32(exp_f64(f64(x))) }

@(require_results) ln_f16le :: proc "contextless" (x: f16le) -> f16le { return #force_inline f16le(ln_f64(f64(x))) }
@(require_results) ln_f16be :: proc "contextless" (x: f16be) -> f16be { return #force_inline f16be(ln_f64(f64(x))) }
@(require_results) ln_f32le :: proc "contextless" (x: f32le) -> f32le { return #force_inline f32le(ln_f64(f64(x))) }
@(require_results) ln_f32be :: proc "contextless" (x: f32be) -> f32be { return #force_inline f32be(ln_f64(f64(x))) }
@(require_results) ln_f64le :: proc "contextless" (x: f64le) -> f64le { return #force_inline f64le(ln_f64(f64(x))) }
@(require_results) ln_f64be :: proc "contextless" (x: f64be) -> f64be { return #force_inline f64be(ln_f64(f64(x))) }
@(require_results) sqrt_f16 :: proc "c" (x: f16, loc := #caller_location) -> f16 { return f16(sqrt_f64(f64(x), loc)) }
@(require_results) sin_f16 :: proc "c" (θ: f16, loc := #caller_location) -> f16 { return f16(sin_f64(f64(θ), loc)) }
@(require_results) cos_f16 :: proc "c" (θ: f16, loc := #caller_location) -> f16 { return f16(cos_f64(f64(θ), loc)) }
@(require_results) pow_f16 :: proc "c" (x, power: f16, loc := #caller_location) -> f16 { return f16(pow_f64(f64(x), f64(power), loc)) }
@(require_results) fmuladd_f16 :: proc "c" (a, b, c: f16, loc := #caller_location) -> f16 { return f16(fmuladd_f64(f64(a), f64(a), f64(c), loc)) }
@(require_results) ln_f16 :: proc "c" (x: f16, loc := #caller_location) -> f16 { return f16(ln_f64(f64(x), loc)) }
@(require_results) exp_f16 :: proc "c" (x: f16, loc := #caller_location) -> f16 { return f16(exp_f64(f64(x), loc)) }

@(require_results) sqrt_f32 :: proc "c" (x: f32, loc := #caller_location) -> f32 { return f32(sqrt_f64(f64(x), loc)) }
@(require_results) sin_f32 :: proc "c" (θ: f32, loc := #caller_location) -> f32 { return f32(sin_f64(f64(θ), loc)) }
@(require_results) cos_f32 :: proc "c" (θ: f32, loc := #caller_location) -> f32 { return f32(cos_f64(f64(θ), loc)) }
@(require_results) pow_f32 :: proc "c" (x, power: f32, loc := #caller_location) -> f32 { return f32(pow_f64(f64(x), f64(power), loc)) }
@(require_results) fmuladd_f32 :: proc "c" (a, b, c: f32, loc := #caller_location) -> f32 { return f32(fmuladd_f64(f64(a), f64(a), f64(c), loc)) }
@(require_results) ln_f32 :: proc "c" (x: f32, loc := #caller_location) -> f32 { return f32(ln_f64(f64(x), loc)) }
@(require_results) exp_f32 :: proc "c" (x: f32, loc := #caller_location) -> f32 { return f32(exp_f64(f64(x), loc)) }

@(require_results) ln_f16le :: proc "contextless" (x: f16le, loc := #caller_location) -> f16le { return #force_inline f16le(ln_f64(f64(x), loc)) }
@(require_results) ln_f16be :: proc "contextless" (x: f16be, loc := #caller_location) -> f16be { return #force_inline f16be(ln_f64(f64(x), loc)) }
@(require_results) ln_f32le :: proc "contextless" (x: f32le, loc := #caller_location) -> f32le { return #force_inline f32le(ln_f64(f64(x), loc)) }
@(require_results) ln_f32be :: proc "contextless" (x: f32be, loc := #caller_location) -> f32be { return #force_inline f32be(ln_f64(f64(x), loc)) }
@(require_results) ln_f64le :: proc "contextless" (x: f64le, loc := #caller_location) -> f64le { return #force_inline f64le(ln_f64(f64(x), loc)) }
@(require_results) ln_f64be :: proc "contextless" (x: f64be, loc := #caller_location) -> f64be { return #force_inline f64be(ln_f64(f64(x), loc)) }
ln :: proc{
ln_f16, ln_f16le, ln_f16be,
ln_f32, ln_f32le, ln_f32be,
Expand Down
36 changes: 20 additions & 16 deletions core/math/math_erf.odin
Original file line number Diff line number Diff line change
Expand Up @@ -117,15 +117,17 @@ erf :: proc{
erf_f64,
}

@(require_results) erf_f16 :: proc "contextless" (x: f16) -> f16 { return f16(erf_f64(f64(x))) }
@(require_results) erf_f16le :: proc "contextless" (x: f16le) -> f16le { return f16le(erf_f64(f64(x))) }
@(require_results) erf_f16be :: proc "contextless" (x: f16be) -> f16be { return f16be(erf_f64(f64(x))) }
@(require_results) erf_f32 :: proc "contextless" (x: f32) -> f32 { return f32(erf_f64(f64(x))) }
@(require_results) erf_f32le :: proc "contextless" (x: f32le) -> f32le { return f32le(erf_f64(f64(x))) }
@(require_results) erf_f32be :: proc "contextless" (x: f32be) -> f32be { return f32be(erf_f64(f64(x))) }
@(require_results) erf_f16 :: proc "contextless" (x: f16, loc := #caller_location) -> f16 { return f16(erf_f64(f64(x), loc)) }
@(require_results) erf_f16le :: proc "contextless" (x: f16le, loc := #caller_location) -> f16le { return f16le(erf_f64(f64(x), loc)) }
@(require_results) erf_f16be :: proc "contextless" (x: f16be, loc := #caller_location) -> f16be { return f16be(erf_f64(f64(x), loc)) }
@(require_results) erf_f32 :: proc "contextless" (x: f32, loc := #caller_location) -> f32 { return f32(erf_f64(f64(x), loc)) }
@(require_results) erf_f32le :: proc "contextless" (x: f32le, loc := #caller_location) -> f32le { return f32le(erf_f64(f64(x), loc)) }
@(require_results) erf_f32be :: proc "contextless" (x: f32be, loc := #caller_location) -> f32be { return f32be(erf_f64(f64(x), loc)) }

@(require_results)
erf_f64 :: proc "contextless" (x: f64) -> f64 {
erf_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)

erx :: 0h3FEB0AC160000000
// Coefficients for approximation to erf in [0, 0.84375]
efx :: 0h3FC06EBA8214DB69
Expand Down Expand Up @@ -251,7 +253,7 @@ erf_f64 :: proc "contextless" (x: f64) -> f64 {
S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
}
z := transmute(f64)(0xffffffff00000000 & transmute(u64)x) // pseudo-single (20-bit) precision x
r := exp(-z*z-0.5625) * exp((z-x)*(z+x)+R/S)
r := exp(-z*z-0.5625, loc) * exp((z-x)*(z+x)+R/S, loc)
if sign {
return r/x - 1
}
Expand All @@ -269,15 +271,17 @@ erfc :: proc{
erfc_f64,
}

@(require_results) erfc_f16 :: proc "contextless" (x: f16) -> f16 { return f16(erfc_f64(f64(x))) }
@(require_results) erfc_f16le :: proc "contextless" (x: f16le) -> f16le { return f16le(erfc_f64(f64(x))) }
@(require_results) erfc_f16be :: proc "contextless" (x: f16be) -> f16be { return f16be(erfc_f64(f64(x))) }
@(require_results) erfc_f32 :: proc "contextless" (x: f32) -> f32 { return f32(erfc_f64(f64(x))) }
@(require_results) erfc_f32le :: proc "contextless" (x: f32le) -> f32le { return f32le(erfc_f64(f64(x))) }
@(require_results) erfc_f32be :: proc "contextless" (x: f32be) -> f32be { return f32be(erfc_f64(f64(x))) }
@(require_results) erfc_f16 :: proc "contextless" (x: f16, loc := #caller_location) -> f16 { return f16(erfc_f64(f64(x), loc)) }
@(require_results) erfc_f16le :: proc "contextless" (x: f16le, loc := #caller_location) -> f16le { return f16le(erfc_f64(f64(x), loc)) }
@(require_results) erfc_f16be :: proc "contextless" (x: f16be, loc := #caller_location) -> f16be { return f16be(erfc_f64(f64(x), loc)) }
@(require_results) erfc_f32 :: proc "contextless" (x: f32, loc := #caller_location) -> f32 { return f32(erfc_f64(f64(x), loc)) }
@(require_results) erfc_f32le :: proc "contextless" (x: f32le, loc := #caller_location) -> f32le { return f32le(erfc_f64(f64(x), loc)) }
@(require_results) erfc_f32be :: proc "contextless" (x: f32be, loc := #caller_location) -> f32be { return f32be(erfc_f64(f64(x), loc)) }

@(require_results)
erfc_f64 :: proc "contextless" (x: f64) -> f64 {
erfc_f64 :: proc "contextless" (x: f64, loc := #caller_location) -> f64 {
validate_finite(x, loc)

erx :: 0h3FEB0AC160000000
// Coefficients for approximation to erf in [0, 0.84375]
efx :: 0h3FC06EBA8214DB69
Expand Down Expand Up @@ -399,7 +403,7 @@ erfc_f64 :: proc "contextless" (x: f64) -> f64 {
S = 1 + s*(sb1+s*(sb2+s*(sb3+s*(sb4+s*(sb5+s*(sb6+s*sb7))))))
}
z := transmute(f64)(0xffffffff00000000 & transmute(u64)x) // pseudo-single (20-bit) precision x
r := exp(-z*z-0.5625) * exp((z-x)*(z+x)+R/S)
r := exp(-z*z-0.5625, loc) * exp((z-x)*(z+x)+R/S, loc)
if sign {
return 2 - r/x
}
Expand Down
Loading
Loading