Fix density function inaccuracies (#44)

* Add clamp to PerlinNoise.wrap

* Add some Math.fround to spline calculations

* Fix tests

* Clamp more Math.floor calls

* Add more Math.fround

* Cleanup SimplexNoise sample2D

src/math/CubicSpline.ts

@@ -1,5 +1,5 @@
 import { Json } from '../util/index.js'
-import { binarySearch, lerp } from './Util.js'
+import { binarySearch, floatLerp } from './Util.js'
 
 export interface NumberFunction<C> {
 	compute(c: C): number,
@@ -80,24 +80,25 @@ export namespace CubicSpline {
 			const coordinate = this.coordinate.compute(c)
 			const i = binarySearch(0, this.locations.length, n => coordinate < this.locations[n]) - 1
 			const n = this.locations.length - 1
+			// TODO: use linear extend for this 
 			if (i < 0) {
-				return this.values[0].compute(c) + this.derivatives[0] * (coordinate - this.locations[0])  //TODO: us linear extend for this 
+				return Math.fround(this.values[0].compute(c) + Math.fround(this.derivatives[0] * Math.fround(coordinate - this.locations[0])))
 			}
 			if (i === n) {
-				return this.values[n].compute(c) + this.derivatives[n] * (coordinate - this.locations[n])  //TODO: us linear extend for this 
+				return Math.fround(this.values[n].compute(c) + Math.fround(this.derivatives[n] * Math.fround(coordinate - this.locations[n])))
 			}
 			const loc0 = this.locations[i]
 			const loc1 = this.locations[i + 1]
 			const der0 = this.derivatives[i]
 			const der1 = this.derivatives[i + 1]
-			const f = (coordinate - loc0) / (loc1 - loc0)
-	
+			const f = Math.fround(Math.fround(coordinate - loc0) / Math.fround(loc1 - loc0))
+
 			const val0 = this.values[i].compute(c)
 			const val1 = this.values[i + 1].compute(c)
-			
-			const f8 = der0 * (loc1 - loc0) - (val1 - val0)
-			const f9 = -der1 * (loc1 - loc0) + (val1 - val0)
-			const f10 = lerp(f, val0, val1) + f * (1.0 - f) * lerp(f, f8, f9)
+
+			const f8 = Math.fround(Math.fround(der0 * Math.fround(loc1 - loc0)) - Math.fround(val1 - val0))
+			const f9 = Math.fround(Math.fround(-der1 * Math.fround(loc1 - loc0)) + Math.fround(val1 - val0))
+			const f10 = Math.fround(floatLerp(f, val0, val1) + Math.fround(Math.fround(f * Math.fround(1.0 - f)) * floatLerp(f, f8, f9)))
 			return f10
 		}
 
@@ -114,11 +115,11 @@ export namespace CubicSpline {
 		}
 
 		public addPoint(location: number, value: number | CubicSpline<C>, derivative = 0) {
-			this.locations.push(location)
+			this.locations.push(Math.fround(location))
 			this.values.push(typeof value === 'number'
-				? new CubicSpline.Constant(value)
+				? new CubicSpline.Constant(Math.fround(value))
 				: value)
-			this.derivatives.push(derivative)
+			this.derivatives.push(Math.fround(derivative))
 			return this
 		}
 
@@ -159,7 +160,7 @@ export namespace CubicSpline {
 			for (var i = 0; i < lastIdx; ++i) {
 				const locationLeft = this.locations[i]
 				const locationRight = this.locations[i + 1]
-				const locationDelta = locationRight - locationLeft
+				const locationDelta = Math.fround(locationRight - locationLeft)
 				const splineLeft = this.values[i]
 				const splineRight = this.values[i + 1]
 				const minLeft = splineLeft.min()
@@ -169,18 +170,18 @@ export namespace CubicSpline {
 				const derivativeLeft = this.derivatives[i]
 				const derivativeRight = this.derivatives[i + 1]
 				if (derivativeLeft !== 0.0 || derivativeRight !== 0.0) {
-					const maxValueDeltaLeft = derivativeLeft * locationDelta
-					const maxValueDeltaRight = derivativeRight * locationDelta
+					const maxValueDeltaLeft = Math.fround(derivativeLeft * locationDelta)
+					const maxValueDeltaRight = Math.fround(derivativeRight * locationDelta)
 					const minValue = Math.min(minLeft, minRight)
 					const maxValue = Math.max(maxLeft, maxRight)
-					const minDeltaLeft = maxValueDeltaLeft - maxRight + minLeft
-					const maxDeltaLeft = maxValueDeltaLeft - minRight + maxLeft
-					const minDeltaRight = -maxValueDeltaRight + minRight - maxLeft
-					const maxDeltaRight = -maxValueDeltaRight + maxRight - minLeft
+					const minDeltaLeft = Math.fround(Math.fround(maxValueDeltaLeft - maxRight) + minLeft)
+					const maxDeltaLeft = Math.fround(Math.fround(maxValueDeltaLeft - minRight) + maxLeft)
+					const minDeltaRight = Math.fround(Math.fround(-maxValueDeltaRight + minRight) - maxLeft)
+					const maxDeltaRight = Math.fround(Math.fround(-maxValueDeltaRight + maxRight) - minLeft)
 					const minDelta = Math.min(minDeltaLeft, minDeltaRight)
 					const maxDelta = Math.max(maxDeltaLeft, maxDeltaRight)
-					splineMin = Math.min(splineMin, minValue + 0.25 * minDelta)
-					splineMax = Math.max(splineMax, maxValue + 0.25 * maxDelta)
+					splineMin = Math.min(splineMin, Math.fround(minValue + Math.fround(0.25 * minDelta)))
+					splineMax = Math.max(splineMax, Math.fround(maxValue + Math.fround(0.25 * maxDelta)))
 				}
 			}
 
@@ -191,7 +192,10 @@ export namespace CubicSpline {
 
 		private static linearExtend(location: number, locations: number[], value: number, derivatives: number[], useIndex: number) {
 			const derivative = derivatives[useIndex]
-			return derivative == 0.0 ? value : value + derivative * (location - locations[useIndex])
+			if (derivative == 0) {
+				return value
+			}
+			return Math.fround(value + Math.fround(derivative * Math.fround(location - locations[useIndex])))
 		}
 	}
 }

src/math/Util.ts

@@ -1,5 +1,10 @@
 import type { Random } from './random/index.js'
 
+const MIN_INT = -2147483648
+const MAX_INT = 2147483647
+const MIN_LONG = -9223372036854776000
+const MAX_LONG = 9223372036854776000
+
 export function square(x: number) {
 	return x * x
 }
@@ -12,6 +17,10 @@ export function lerp(a: number, b: number, c: number): number {
 	return b + a * (c - b)
 }
 
+export function floatLerp(a: number, b: number, c: number): number {
+	return Math.fround(b + Math.fround(a * Math.fround(c - b)))
+}
+
 export function lerp2(a: number, b: number, c: number, d: number, e: number, f: number): number {
 	return lerp(b, lerp(a, c, d), lerp(a, e, f))
 }
@@ -60,6 +69,14 @@ export function clampedMap(a: number, b: number, c: number, d: number, e: number
 	return clampedLerp(d, e, inverseLerp(a, b, c))
 }
 
+export function intFloor(a: number) {
+	return clamp(Math.floor(a), MIN_INT, MAX_INT)
+}
+
+export function longFloor(a: number) {
+	return clamp(Math.floor(a), MIN_LONG, MAX_LONG)
+}
+
 export function binarySearch(n: number, n2: number, predicate: (value: number) => boolean) {
 	let n3 = n2 - n
 	while (n3 > 0) {

src/math/noise/ImprovedNoise.ts

@@ -1,5 +1,5 @@
 import type { Random } from '../random/index.js'
-import { lerp3, smoothstep } from '../Util.js'
+import { intFloor, lerp3, smoothstep } from '../Util.js'
 import { SimplexNoise } from './SimplexNoise.js'
 
 export class ImprovedNoise {
@@ -29,17 +29,17 @@ export class ImprovedNoise {
 		const x2 = x + this.xo
 		const y2 = y + this.yo
 		const z2 = z + this.zo
-		const x3 = Math.floor(x2)
-		const y3 = Math.floor(y2)
-		const z3 = Math.floor(z2)
+		const x3 = intFloor(x2)
+		const y3 = intFloor(y2)
+		const z3 = intFloor(z2)
 		const x4 = x2 - x3
 		const y4 = y2 - y3
 		const z4 = z2 - z3
 
 		let y6 = 0
 		if (yScale !== 0) {
 			const t = yLimit >= 0 && yLimit < y4 ? yLimit : y4
-			y6 = Math.floor(t / yScale + 1e-7) * yScale
+			y6 = intFloor(t / yScale + 1e-7) * yScale
 		}
 
 		return this.sampleAndLerp(x3, y3, z3, x4, y4 - y6, z4, y4)

src/math/noise/PerlinNoise.ts

@@ -1,5 +1,6 @@
 import type { Random } from '../random/index.js'
 import { XoroshiroRandom } from '../random/index.js'
+import { longFloor } from '../Util.js'
 import { ImprovedNoise } from './ImprovedNoise.js'
 
 export class PerlinNoise {
@@ -79,6 +80,6 @@ export class PerlinNoise {
 	}
 
 	public static wrap(value: number) {
-		return value - Math.floor(value / 3.3554432E7 + 0.5) * 3.3554432E7
+		return value - longFloor(value / 3.3554432E7 + 0.5) * 3.3554432E7
 	}
 }

src/math/noise/SimplexNoise.ts

@@ -1,4 +1,5 @@
 import type { Random } from '../random/index.js'
+import { intFloor } from '../Util.js'
 
 export class SimplexNoise {
 	private static readonly GRADIENT = [[1, 1, 0], [-1, 1, 0], [1, -1, 0], [-1, -1, 0], [1, 0, 1], [-1, 0, 1], [1, 0, -1], [-1, 0, -1], [0, 1, 1], [0, -1, 1], [0, 1, -1], [0, -1, -1], [1, 1, 0], [0, -1, 1], [-1, 1, 0], [0, -1, -1]]
@@ -28,16 +29,16 @@ export class SimplexNoise {
 	}
 
 	public sample2D(d: number, d2: number) {
-		let d3
-		let n3
-		let d4
 		const d6 = (d + d2) * SimplexNoise.F2
-		const n4 = Math.floor(d + d6)
-		const d7 = n4 - (d3 = (n4 + (n3 = Math.floor(d2 + d6))) * SimplexNoise.G2)
+		const n4 = intFloor(d + d6)
+		const n3 = intFloor(d2 + d6)
+		const d3 = (n4 + n3) * SimplexNoise.G2
+		const d7 = n4 - d3
 		const d8 = d - d7
 		let a
 		let b
-		if (d8 > (d4 = d2 - (n3 - d3))) {
+		const d4 = d2 - (n3 - d3)
+		if (d8 > d4) {
 			a = 1
 			b = 0
 		} else {
@@ -61,9 +62,9 @@ export class SimplexNoise {
 
 	public sample(x: number, y: number, z: number) {
 		const d5 = (x + y + z) * 0.3333333333333333
-		const x2 = Math.floor(x + d5)
-		const y2 = Math.floor(y + d5)
-		const z2 = Math.floor(z + d5)
+		const x2 = intFloor(x + d5)
+		const y2 = intFloor(y + d5)
+		const z2 = intFloor(z + d5)
 		const d7 = (x2 + y2 + z2) * 0.16666666666666666
 		const x3 = x - (x2 - d7)
 		const y3 = y - (y2 - d7)

test/math/Spline.test.ts

@@ -11,10 +11,10 @@ describe('Spline', () => {
 			.addPoint(-0.44, -0.12)
 			.addPoint(-0.18, -0.12)
 
-		expect(spline.compute(-1.6)).toEqual(0.044)
-		expect(spline.compute(-0.7)).toEqual(-0.2222)
+		expect(spline.compute(-1.6)).toBeCloseTo(0.044, DELTA)
+		expect(spline.compute(-0.7)).toBeCloseTo(-0.2222, DELTA)
 		expect(spline.compute(-0.5)).toBeCloseTo(-0.21653879, DELTA)
-		expect(spline.compute(-0.2)).toEqual(-0.12)
+		expect(spline.compute(-0.2)).toBeCloseTo(-0.12, DELTA)
 	})
 
 	it('derivatives', () => {
@@ -25,7 +25,7 @@ describe('Spline', () => {
 			.addPoint(0.6, 0.4, 0.0)
 
 		expect(spline.compute(-0.1)).toBeCloseTo(-0.0022000019, DELTA)
-		expect(spline.compute(0)).toEqual(0.0178)
+		expect(spline.compute(0)).toBeCloseTo(0.0178, DELTA)
 		expect(spline.compute(0.31)).toBeCloseTo(0.24358201, DELTA)
 		expect(spline.compute(0.4)).toBeCloseTo(0.69171876, DELTA)
 	})

test/worldgen/DensityFunction.test.ts

@@ -6,6 +6,8 @@ import { DensityFunction as DF, NoiseGeneratorSettings, NoiseRouter, WorldgenReg
 import { RandomState } from '../../src/worldgen/RandomState.js'
 
 describe('DensityFunction', () => {
+	const DELTA = 1e-7
+
 	const ContextA = DF.context(1, 2, 3)
 	const ContextB = DF.context(2, 3, 4)
 	const ContextC = DF.context(12, -30, 1)
@@ -151,12 +153,12 @@ describe('DensityFunction', () => {
 			.addPoint(5, 0.2)
 			.addPoint(20, 0.7)
 		const fn2 = wrap(new DF.Spline(spline))
-		expect(fn2.compute(DF.context(0, 0, 0))).toEqual(1)
-		expect(fn2.compute(DF.context(0, 3.2, 0))).toEqual(0.4363904)
-		expect(fn2.compute(DF.context(0, 5, 0))).toEqual(0.2)
-		expect(fn2.compute(DF.context(0, 11, 0))).toEqual(0.376)
-		expect(fn2.compute(DF.context(0, 20, 0))).toEqual(0.7)
-		expect(fn2.compute(DF.context(0, 25, 0))).toEqual(0.7)
+		expect(fn2.compute(DF.context(0, 0, 0))).toBeCloseTo(1, DELTA)
+		expect(fn2.compute(DF.context(0, 3.2, 0))).toBeCloseTo(0.4363904, DELTA)
+		expect(fn2.compute(DF.context(0, 5, 0))).toBeCloseTo(0.2, DELTA)
+		expect(fn2.compute(DF.context(0, 11, 0))).toBeCloseTo(0.376, DELTA)
+		expect(fn2.compute(DF.context(0, 20, 0))).toBeCloseTo(0.7, DELTA)
+		expect(fn2.compute(DF.context(0, 25, 0))).toBeCloseTo(0.7, DELTA)
 	})
 
 	it('YClampedGradient', () => {