feat(r1): math

r1/math.ts

@@ -0,0 +1,35 @@
+/** Computes the IEEE 754 floating-point remainder of x / y. */
+export const remainder = (x: number, y: number): number => {
+  if (isNaN(x) || isNaN(y) || !isFinite(x) || y === 0) return NaN
+
+  const quotient = x / y
+  let n = Math.round(quotient)
+
+  // When quotient is exactly halfway between two integers, round to the nearest even integer
+  if (Math.abs(quotient - n) === 0.5) n = 2 * Math.round(quotient / 2)
+
+  const rem = x - n * y
+  return !rem ? Math.sign(x) * 0 : rem
+}
+
+/** Returns the next representable floating-point value after x towards y. */
+export const nextAfter = (x: number, y: number): number => {
+  if (isNaN(x) || isNaN(y)) return NaN
+  if (x === y) return y
+  if (x === 0) return y > 0 ? Number.MIN_VALUE : -Number.MIN_VALUE
+
+  const buffer = new ArrayBuffer(8)
+  const view = new DataView(buffer)
+
+  view.setFloat64(0, x, true)
+  let bits = view.getBigUint64(0, true)
+
+  if (x > 0 === y > x) bits += 1n
+  else bits -= 1n
+
+  view.setBigUint64(0, bits, true)
+  return view.getFloat64(0, true)
+}
+
+/** Returns true IFF a is within epsilon distance of b. */
+export const float64Near = (a: number, b: number, epsilon: number) => Math.abs(a - b) <= epsilon;

r1/math_test.ts

@@ -0,0 +1,35 @@
+import test from 'node:test'
+import { equal, ok } from 'node:assert/strict'
+import { remainder, nextAfter, float64Near } from './math'
+
+test('remainder', t => {
+  equal(remainder(5.1, 2), -0.9000000000000004)
+  equal(remainder(5.5, 2), -0.5)
+  equal(remainder(-5.5, 2), 0.5)
+  equal(remainder(5, 2.5), 0)
+  equal(remainder(5.1, 0), NaN)
+  equal(remainder(Infinity, 2), NaN)
+  equal(remainder(5, NaN), NaN)
+  equal(remainder(0, 2), 0)
+  equal(remainder(5, 2), 1)
+  equal(remainder(-5, 2), -1)
+})
+
+test('nextAfter', t => {
+  equal(nextAfter(0, 1), 5e-324)
+  equal(nextAfter(0, -1), -5e-324)
+  equal(nextAfter(1, 2), 1.0000000000000002)
+  equal(nextAfter(1, 0), 0.9999999999999999)
+  equal(nextAfter(1, 1), 1)
+  equal(nextAfter(Number.MAX_VALUE, Infinity), Infinity)
+  equal(nextAfter(-Number.MAX_VALUE, -Infinity), -Infinity)
+  equal(nextAfter(Number.NaN, 1), NaN)
+  equal(nextAfter(1, Number.NaN), NaN)
+})
+
+test('float64Near', t => {
+  ok(float64Near(0, 0, 0))
+  ok(float64Near(1e-10, 1e-10*2, 1e-10))
+  ok(!float64Near(1e-10, 1e-9, 1e-10))
+  ok(!float64Near(1e-5, 1e-4, 1e-5/10))
+})

s1/angle.ts

@@ -1,5 +1,6 @@
 import type { Angle } from './_types'
 import { DEGREE } from './angle_constants'
+import { remainder } from '../r1/math'
 
 /**
  * Angle represents a 1D angle. The internal representation is a double precision value in radians, so conversion to and from radians is exact.
@@ -70,20 +71,6 @@ export const e7 = (a: Angle): number => round(degrees(a) * 1e7)
 /** Returns the absolute value of the angle. */
 export const abs = (a: Angle): Angle => Math.abs(a)
 
-/** Computes the IEEE 754 floating-point remainder of x / y. */
-function remainder(x: number, y: number): number {
-  if (isNaN(x) || isNaN(y) || !isFinite(x) || y === 0) return NaN
-
-  const quotient = x / y
-  let n = Math.round(quotient)
-
-  // When quotient is exactly halfway between two integers, round to the nearest even integer
-  if (Math.abs(quotient - n) === 0.5) n = 2 * Math.round(quotient / 2)
-
-  const rem = x - n * y
-  return !rem ? Math.sign(x) * 0 : rem
-}
-
 /** Returns an equivalent angle in (-π, π]. */
 export const normalized = (a: Angle): Angle => {
   let rad = remainder(a, 2 * Math.PI)