From 889d7d5c759cfec0b00d884a15976e94eb8bd43e Mon Sep 17 00:00:00 2001 From: Masashi Mizuno <29035331+unageek@users.noreply.github.com> Date: Sun, 4 Jul 2021 21:51:49 +0900 Subject: [PATCH] New example: compute exp(x) --- Cargo.toml | 4 ++ examples/exp.rs | 119 +++++++++++++++++++++++++++++++++ examples/gen_exp_constants.wls | 17 +++++ tools/gen_constants.wls | 18 +---- tools/interval.wl | 27 ++++++++ 5 files changed, 169 insertions(+), 16 deletions(-) create mode 100644 examples/exp.rs create mode 100755 examples/gen_exp_constants.wls create mode 100644 tools/interval.wl diff --git a/Cargo.toml b/Cargo.toml index 9e44113..ff9e83a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -63,6 +63,10 @@ optional = true name = "exact" required-features = ["gmp"] +[[example]] +name = "exp" +required-features = ["gmp"] + [[example]] name = "parse" required-features = ["gmp"] diff --git a/examples/exp.rs b/examples/exp.rs new file mode 100644 index 0000000..ccbde50 --- /dev/null +++ b/examples/exp.rs @@ -0,0 +1,119 @@ +// An interval extension of exp(x), exp2(x) and exp10(x). Not tightest, but simple. + +#![allow(clippy::approx_constant)] + +use inari::*; + +// To compute exp(x), change the base to 2 using the relation: +// +// exp(x) = 2^(x log_2 e). +// +// Now consider computing 2^x. Let D = [-1/2, 1/2]. +// Decompose x into an integer and a fractional parts: +// +// x = a + b, +// +// where a ∈ ℤ, b ∈ D. Therefore, 2^x = 2^a 2^b. The f64-interval extension of 2^a is trivial. +// From the Taylor series about b = 0, 2^b is enclosed by: +// +// N-1 (ln 2)^n (ln 2)^N +// 2^b ∈ { ∑ -------- b^n + -------- 2^β b^N | β ∈ D}, ∀b ∈ D. +// n=0 n! N! +// +// From the following relations, N = 14 would be sufficient: +// +// | (ln 2)^N | (ln 2)^N +// max | -------- 2^β b^N | ≤ -------- 2^(1/2) (1/2)^N, +// b,β∈D | N! | N! +// +// (ln 2)^N | +// -------- 2^(1/2) (1/2)^N | ≤ 2^-53. +// N! |N=14 +// +// References: +// +// - 大石 進一 (2018), 精度保証付き数値計算の基礎, コロナ社, ISBN 978-4-339-02887-4. +// - https://en.wikipedia.org/wiki/Taylor%27s_theorem#Explicit_formulas_for_the_remainder + +fn exp2_point_reduced(x: f64) -> Interval { + assert!(x.abs() <= 0.5); + + const N: usize = 14; + /// C[n] = (ln 2)^n / n! for n = 0, …, N - 1. + const C: [Interval; N] = [ + const_interval!(1.0, 1.0), + const_interval!(0.6931471805599453, 0.6931471805599454), + const_interval!(0.2402265069591007, 0.24022650695910072), + const_interval!(5.5504108664821576e-2, 5.550410866482158e-2), + const_interval!(9.618129107628477e-3, 9.618129107628479e-3), + const_interval!(1.3333558146428443e-3, 1.3333558146428445e-3), + const_interval!(1.540353039338161e-4, 1.5403530393381611e-4), + const_interval!(1.525273380405984e-5, 1.5252733804059841e-5), + const_interval!(1.3215486790144307e-6, 1.321548679014431e-6), + const_interval!(1.0178086009239699e-7, 1.01780860092397e-7), + const_interval!(7.0549116208011225e-9, 7.054911620801123e-9), + const_interval!(4.445538271870811e-10, 4.4455382718708116e-10), + const_interval!(2.5678435993488202e-11, 2.5678435993488206e-11), + const_interval!(1.3691488853904128e-12, 1.369148885390413e-12), + ]; + + let x = interval!(x, x).unwrap(); + // Initially, y = 2^[-1/2, 1/2] (ln 2)^N / N!. + let mut y = const_interval!(4.7932833733029885e-14, 9.586566746605978e-14); + for n in (0..N).rev() { + y = y.mul_add(x, C[n]); + } + y +} + +fn exp2_point(x: f64) -> Interval { + assert!(x.is_finite()); + + let mut a = x.floor(); + let mut b = x - a; + if b > 0.5 { + a += 1.0; + b -= 1.0; + } + + let exp2_a = if a < -1074.0 { + const_interval!(0.0, 5e-324) + } else if a > 1023.0 { + const_interval!(f64::MAX, f64::INFINITY) + } else { + interval!(a.exp2(), a.exp2()).unwrap() + }; + + let exp2_b = exp2_point_reduced(b); + + exp2_a * exp2_b +} + +fn exp2(x: Interval) -> Interval { + let inf = if x.inf() == f64::NEG_INFINITY { + 0.0 + } else { + exp2_point(x.inf()).inf() + }; + let sup = if x.sup() == f64::INFINITY { + f64::INFINITY + } else { + exp2_point(x.sup()).sup() + }; + interval!(inf, sup).unwrap() +} + +fn exp(x: Interval) -> Interval { + exp2(Interval::LOG2_E * x) +} + +fn exp10(x: Interval) -> Interval { + exp2(Interval::LOG2_10 * x) +} + +fn main() { + let x = interval!("[1.234567]").unwrap(); + println!("2^x ⊆ {:.15}", exp2(x)); + println!("e^x ⊆ {:.15}", exp(x)); + println!("10^x ⊆ {:.15}", exp10(x)); +} diff --git a/examples/gen_exp_constants.wls b/examples/gen_exp_constants.wls new file mode 100755 index 0000000..6278223 --- /dev/null +++ b/examples/gen_exp_constants.wls @@ -0,0 +1,17 @@ +#!/usr/bin/env wolframscript + +Import[FileNameJoin @ {DirectoryName[$InputFileName], "..", "tools", "interval.wl"}]; + +print[value_] := + Print[intervalExpression[value] <> ","]; + +n = 14; + +print[2^(-1/2) Log[2]^n / n!]; +print[2^(1/2) Log[2]^n / n!]; + +Print["----"]; + +For[i = 0, i < n, i++, + print[Log[2]^i / i!]; + ] diff --git a/tools/gen_constants.wls b/tools/gen_constants.wls index 8b68e83..85048e7 100755 --- a/tools/gen_constants.wls +++ b/tools/gen_constants.wls @@ -1,23 +1,9 @@ #!/usr/bin/env wolframscript -(* Format the source code in Wolfram Workbench before committing. *) - -toInterval[x_] := - Module[ {ulp, inf, sup}, - ulp = 2^(Floor[Log2[x]] - 52); - inf = Floor[x, ulp]; - sup = Ceiling[x, ulp]; - {inf, sup} - ]; - -toLiteral[x_] := - ExportString[x, "JSON"]; +Import[FileNameJoin @ {DirectoryName[$InputFileName], "interval.wl"}]; def[name_, value_] := - Module[ {inf, sup}, - {inf, sup} = toLiteral /@ toInterval[value]; - Print["pub const " <> name <> ": Self = const_interval!(" <> inf <> "," <> sup <> ");"]; - ]; + Print["pub const " <> name <> ": Self = " <> intervalExpression[value] <> ";"]; def["E", E]; def["FRAC_1_PI", 1/Pi]; diff --git a/tools/interval.wl b/tools/interval.wl new file mode 100644 index 0000000..b61ed87 --- /dev/null +++ b/tools/interval.wl @@ -0,0 +1,27 @@ +interval[x_] := + Module[ {max, ulp, inf, sup}, + max = SetPrecision[$MaxMachineNumber, Infinity]; + ulp = 2^(Max[Floor[Log2[Abs[x]]] - 52, -1074]); + inf = If[ x < -max, + -Infinity, + N @ Floor[x, ulp] + ]; + sup = If[ x > max, + Infinity, + N @ Ceiling[x, ulp] + ]; + {inf, sup} + ]; + +f64Expression[x_] := + Switch[x, + Infinity, "f64::INFINITY", + -Infinity, "f64::NEG_INFINITY", + _, ExportString[x, "JSON"] + ]; + +intervalExpression[x_] := + Module[ {inf, sup}, + {inf, sup} = f64Expression /@ interval[x]; + "const_interval!(" <> inf <> "," <> sup <> ")" + ];