root: add F32::exp and F32::ln

cspell.json

@@ -14,6 +14,7 @@
     "eisel",
     "elementwise",
     "elif",
+    "expf",
     "extrinsics",
     "frameless",
     "freelist",
@@ -27,6 +28,7 @@
     "indexmap",
     "itertools",
     "lemire",
+    "logf",
     "miri",
     "msun",
     "muls",

root/numeric/F32/F32.vi

@@ -5,7 +5,12 @@ use util::{duplicate, erase};
 #[builtin = "F32"]
 pub type F32;
 
+mod common;
+
 pub mod F32 {
+  pub mod exp;
+  pub mod ln;
+
   pub const nan: F32 = 0.0 / 0.0;
   pub const inf: F32 = 1.0 / 0.0;
   pub const neg_inf: F32 = -inf;

root/numeric/F32/common.vi

@@ -0,0 +1,51 @@
+
+const num_mantissa_bits: N32 = Float::num_mantissa_bits[F32];
+const max_exponent: I32 = Float::max_exponent[F32];
+const exp_mask: N32 = (1[N32] << Float::num_exponent_bits[F32]) - 1;
+
+/// normalize returns a normal number y and exponent exp
+/// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
+pub fn .normalize(x: F32) -> (F32, I32) {
+  if x.to_bits() & (exp_mask << num_mantissa_bits) == 0 {
+    return (x * (1 << num_mantissa_bits) as F32, -num_mantissa_bits);
+  }
+  return (x, +0);
+}
+
+/// Calculates frac * 2 ^ exp.
+pub fn apply_exponent(frac: F32, exp: I32) -> F32 {
+  // special cases
+  if frac == 0.0 or frac == F32::inf or frac == F32::neg_inf or frac.is_nan() {
+    return frac;
+  }
+  let (frac, e) = frac.normalize();
+  exp += e;
+  let x = frac.to_bits();
+  exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - max_exponent;
+  if exp < -149 {
+    // underflow
+    return if frac < 0.0 {
+      -0.0
+    } else {
+      0.0
+    };
+  }
+  if exp > +127 {
+    // overflow
+    return if frac < 0.0 {
+      F32::neg_inf
+    } else {
+      F32::inf
+    };
+  }
+  let m = 1.0;
+  if exp < -126 {
+    // denormalize
+    exp += +24
+    // 2**-24
+    m = 1.0 / (1 << 24) as F32
+  }
+  x &= !(exp_mask << num_mantissa_bits);
+  x |= (exp + max_exponent) as N32 << num_mantissa_bits;
+  return m * F32::from_bits(x);
+}

root/numeric/F32/exp.vi

@@ -0,0 +1,82 @@
+
+// Adapted from https://github.com/rust-lang/compiler-builtins/blob/libm-v0.2.16/libm/src/math/exp.rs
+//
+// origin: FreeBSD /usr/src/lib/msun/src/e_expf.c
+//
+// Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+const ln2_hi: F32 = 6.9314575195e-01;
+const ln2_lo: F32 = 1.4286067653e-06;
+const inv_ln2: F32 = 1.4426950216e+00;
+
+const P1: F32 = 1.6666625440e-1;
+const P2: F32 = -2.7667332906e-3;
+
+const overflow: F32 = 88.7228394[F32];
+const underflow: F32 = -103.972084[F32];
+
+/// Exponential, base *e* (F32)
+///
+/// Calculate the exponential of `f`, that is, *e* raised to the power `f`
+/// (where *e* is the base of the natural system of logarithms, approximately 2.71828).
+pub fn .exp(f: F32) -> F32 {
+  when {
+    f.is_nan() { f }
+    f >= overflow { F32::inf }
+    f <= underflow { 0.0 }
+    _ {
+      let hf = f.to_bits();
+      let is_negative = (hf >> 31) == 1;
+      // high word of |f|
+      hf &= 0x7fffffff;
+
+      // argument reduction
+      when {
+        hf > 0x3eb17218 {
+          // if |f| > 0.5 ln2
+          let k = when {
+            hf > 0x3f851592 {
+              if is_negative {
+                // if f < 1.5 ln2
+                (inv_ln2 * f - 0.5[F32]) as I32
+              } else {
+                // if f > 1.5 ln2
+                (inv_ln2 * f + 0.5[F32]) as I32
+              }
+            }
+            is_negative { -1 }
+            _ { +1 }
+          };
+          let kf = k as F32;
+          let hi = f - kf * ln2_hi;
+          // k*ln2hi is exact here
+          let lo = kf * ln2_lo;
+          f = hi - lo;
+          exp_poly_approximation(f, hi, lo, k)
+        }
+        hf > 0x39000000 { exp_poly_approximation(f, f, 0.0, +0) }
+        _ { 1.0 + f }
+      }
+    }
+  }
+}
+
+fn exp_poly_approximation(f: F32, hi: F32, lo: F32, k: I32) -> F32 {
+  let ff = f * f;
+  let c = f - ff * (P1 + ff * P2);
+  let y = 1.0 + (f * c / (2.0 - c) - lo + hi);
+  if k == +0 {
+    y
+  } else {
+    common::apply_exponent(y, k)
+  }
+}

root/numeric/F32/ln.vi

@@ -0,0 +1,73 @@
+
+// Adapted from https://github.com/rust-lang/compiler-builtins/blob/libm-v0.2.16/libm/src/math/logf.rs
+//
+// origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
+//
+// Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+
+const ln2_hi: F32 = 6.9313812256e-01;
+const ln2_lo: F32 = 9.0580006145e-06;
+// |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]).
+const lg1: F32 = 0.66666662693;
+const lg2: F32 = 0.40000972152;
+const lg3: F32 = 0.28498786688;
+const lg4: F32 = 0.24279078841;
+
+/// The natural logarithm of `f` (F32).
+pub fn .ln(f: F32) -> F32 {
+  // 0x1p25f === 2 ^ 25
+  let x1p25 = F32::from_bits(0x4c000000);
+  let bits = f.to_bits();
+  let k = +0;
+
+  when {
+    bits < 0x00800000 or (bits >> 31) != 0 {
+      // f < 2**-126 (f small or negative)
+      if bits << 1 == 0 {
+        // log(+-0)=-inf
+        return -1.0 / (f * f);
+      }
+      if bits as I32 < +0 {
+        return F32::nan;
+      }
+      // subnormal number, scale up f
+      k -= +25;
+      f *= x1p25;
+      bits = f.to_bits();
+    }
+    bits >= 0x7f800000 {
+      // f is infinite or NaN
+      return f;
+    }
+    bits == 0x3f800000 {
+      // f == 1.0
+      return 0.0;
+    }
+  }
+
+  // reduce f into [sqrt(2)/2, sqrt(2)]
+  bits += (0x3f800000 - 0x3f3504f3);
+  k += ((bits >> 23) as I32) - +0x7f;
+  bits = (bits & 0x007fffff) + 0x3f3504f3;
+  f = F32::from_bits(bits);
+
+  f -= 1.0;
+  let s = f / (2.0 + f);
+  let z = s * s;
+  let w = z * z;
+  let t1 = w * (lg2 + w * lg4);
+  let t2 = z * (lg1 + w * lg3);
+  let r = t2 + t1;
+  let hfsq = 0.5 * f * f;
+  let dk = k as F32;
+  s * (hfsq + r) + dk * ln2_lo - hfsq + f + dk * ln2_hi
+}

root/numeric/F64/common.vi

@@ -5,15 +5,12 @@ pub const log2e: F64 = 1.44269504088896338700e+00[F64];
 
 const num_mantissa_bits: N32 = Float::num_mantissa_bits[F64];
 const max_exponent: I32 = Float::max_exponent[F64];
-
 const exp_mask: N32 = (1[N32] << Float::num_exponent_bits[F64]) - 1;
 
-// normalize returns a normal number y and exponent exp
-// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
-pub fn normalize(x: F64) -> (F64, I32) {
-  // 2**-1022
-  const SmallestNormal: F64 = 2.2250738585072014e-308[F64];
-  if x.abs() < SmallestNormal {
+/// normalize returns a normal number y and exponent exp
+/// satisfying x == y × 2**exp. It assumes x is finite and non-zero.
+pub fn .normalize(x: F64) -> (F64, I32) {
+  if (x.to_bits() >> num_mantissa_bits) as N32 & exp_mask == 0 {
     return (x * (1[N64] << 52) as F64, -52);
   }
   return (x, +0);
@@ -33,7 +30,7 @@ pub fn fraction_exponent(f: F64) -> { frac: F64, exp: I32 } {
   if f == 0.0[F64] or f.is_nan() or f == F64::inf or f == F64::neg_inf {
     return { frac: f, exp: +0 };
   }
-  let (frac, exp) = normalize(f);
+  let (frac, exp) = f.normalize();
   let x = frac.to_bits();
   exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - (max_exponent - +1);
   x &= !(exp_mask as N64 << num_mantissa_bits);
@@ -48,7 +45,7 @@ pub fn apply_exponent(frac: F64, exp: I32) -> F64 {
   if frac == 0.0[F64] or frac == F64::inf or frac == F64::neg_inf or frac.is_nan() {
     return frac;
   }
-  let (frac, e) = normalize(frac);
+  let (frac, e) = frac.normalize();
   exp += e;
   let x = frac.to_bits();
   exp += ((x >> num_mantissa_bits) as N32 & exp_mask) as I32 - max_exponent;

tests/programs/f32_ops.vi

@@ -1,4 +1,9 @@
 
+use #root::rng::Pcg32;
+
+#[configurable]
+const max: N32 = 10_000;
+
 pub fn main(&io: &IO) {
   let interesting_floats = [
     +0.0,
@@ -30,10 +35,15 @@ pub fn main(&io: &IO) {
     1.0e10,
   ];
 
+  let rng = Pcg32::seeded("how many logs could a log calculator calculate if a log calculator could calculate logs");
+  for _ in 0..max {
+    interesting_floats.push_back(F32::from_bits(rng.gen_n32()));
+  }
+
   for f in interesting_floats {
-    let sqrt = f.sqrt();
-    io.println("{f} = {f.to_bits().to_hex()}");
-    io.println("sqrt({f}) = {sqrt} ({sqrt.to_bits().to_hex()})");
+    io.println("sqrt {f.to_bits()} {f.sqrt().to_bits()}");
+    io.println("exp {f.to_bits()} {f.exp().to_bits()}");
+    io.println("ln {f.to_bits()} {f.ln().to_bits()}");
     io.println("");
   }
 }

tests/programs/f64_ops.vi

@@ -1,4 +1,9 @@
 
+use #root::rng::Pcg32;
+
+#[configurable]
+const max: N32 = 10_000;
+
 pub fn main(&io: &IO) {
   let interesting_floats = [
     +0.0[F64],
@@ -33,14 +38,15 @@ pub fn main(&io: &IO) {
     1.0e10[F64],
   ];
 
+  let rng = Pcg32::seeded("how many logs could a log calculator calculate if a log calculator could calculate logs");
+  for _ in 0..max {
+    interesting_floats.push_back(F64::from_bits(N64(rng.gen_n32(), rng.gen_n32())));
+  }
+
   for f in interesting_floats {
-    io.println("{f} = {f.to_bits().to_hex()}");
-    let sqrt = f.sqrt();
-    io.println("sqrt({f}) = {sqrt} ({sqrt.to_bits().to_hex()})");
-    let exp = f.exp();
-    io.println("exp({f}) = {exp} ({exp.to_bits().to_hex()})");
-    let ln = f.ln();
-    io.println("ln({f}) = {ln} ({ln.to_bits().to_hex()})");
+    io.println("sqrt {f.to_bits()} {f.sqrt().to_bits()}");
+    io.println("exp {f.to_bits()} {f.exp().to_bits()}");
+    io.println("ln {f.to_bits()} {f.ln().to_bits()}");
     io.println("");
   }
 }