Okay I need to commit this now before I forget.

benches/benches.rs

@@ -21,7 +21,7 @@ extern crate test;
 
 use noise::Seed;
 use noise::{perlin2, perlin3, perlin4};
-use noise::{open_simplex2, open_simplex3};
+use noise::{open_simplex2, open_simplex3, open_simplex4};
 use noise::{value2, value3, value4};
 use noise::{cell2_range, cell3_range, cell4_range};
 use noise::{cell2_range_inv, cell3_range_inv, cell4_range_inv};
@@ -61,6 +61,12 @@ fn bench_open_simplex3(bencher: &mut Bencher) {
     bencher.iter(|| open_simplex3(black_box(&seed), black_box(&[42.0f64, 37.0, 26.0])));
 }
 
+#[bench]
+fn bench_open_simplex4(bencher: &mut Bencher) {
+    let seed = Seed::new(0);
+    bencher.iter(|| open_simplex4(black_box(&seed), black_box(&[42.0f64, 37.0, 26.0, 128.0])));
+}
+
 #[bench]
 fn bench_value2(bencher: &mut Bencher) {
     let seed = Seed::new(0);
@@ -247,6 +253,18 @@ fn bench_open_simplex3_64x64(bencher: &mut Bencher) {
     });
 }
 
+#[bench]
+fn bench_open_simplex4_64x64(bencher: &mut Bencher) {
+    let seed = Seed::new(0);
+    bencher.iter(|| {
+        for y in 0..64 {
+            for x in 0..64 {
+                black_box(open_simplex4(black_box(&seed), &[x as f64, y as f64, x as f64, y as f64]));
+            }
+        }
+    });
+}
+
 #[bench]
 fn bench_value2_64x64(bencher: &mut Bencher) {
     let seed = Seed::new(0);

src/lib.rs

@@ -33,7 +33,7 @@ pub use seed::Seed;
 pub use math::{Point2, Point3, Point4};
 pub use perlin::{perlin2, perlin3, perlin4};
 pub use value::{value2, value3, value4};
-pub use open_simplex::{open_simplex2, open_simplex3};
+pub use open_simplex::{open_simplex2, open_simplex3, open_simplex4};
 pub use brownian::{Brownian2, Brownian3, Brownian4};
 
 pub use cell::{range_sqr_euclidian2, range_sqr_euclidian3, range_sqr_euclidian4};

src/open_simplex.rs

@@ -17,13 +17,16 @@
 //! http://uniblock.tumblr.com/post/97868843242/noise
 
 use num_traits::Float;
+use std::ops::{Add};
 
 use {gradient, math, Seed};
 
 const STRETCH_CONSTANT_2D: f64 = -0.211324865405187; //(1/sqrt(2+1)-1)/2;
 const SQUISH_CONSTANT_2D: f64 = 0.366025403784439; //(sqrt(2+1)-1)/2;
 const STRETCH_CONSTANT_3D: f64 = -1.0 / 6.0; //(1/Math.sqrt(3+1)-1)/3;
 const SQUISH_CONSTANT_3D: f64 = 1.0 / 3.0; //(Math.sqrt(3+1)-1)/3;
+const STRETCH_CONSTANT_4D: f64 = -0.138196601125011; //(Math.sqrt(4+1)-1)/4;
+const SQUISH_CONSTANT_4D: f64 = 0.309016994374947; //(Math.sqrt(4+1)-1)/4;
 
 const NORM_CONSTANT_2D: f32 = 1.0 / 14.0;
 const NORM_CONSTANT_3D: f32 = 1.0 / 14.0;
@@ -258,3 +261,221 @@ pub fn open_simplex3<T: Float>(seed: &Seed, point: &::Point3<T>) -> T {
 
     return value * math::cast(NORM_CONSTANT_3D);
 }
+
+/// 4-dimensional [OpenSimplex Noise](http://uniblock.tumblr.com/post/97868843242/noise)
+///
+/// This is a slower but higher quality form of gradient noise than
+/// `noise::perlin4`.
+pub fn open_simplex4<T: Float>(seed: &Seed, point: &math::Point4<T>) -> T {
+    fn gradient<T: Float>(seed: &Seed, stretched_floor: &math::Point4<T>, pos: &math::Point4<T>) -> T {
+        let zero: T = math::cast(0.0);
+        let attn = math::cast::<_, T>(2.0_f64) - math::dot4(*pos, *pos);
+        if attn > zero {
+            let index = seed.get4::<isize>(math::cast4::<_, isize>(*stretched_floor));
+            let vec = gradient::get4::<T>(index);
+            math::pow4(attn) * math::dot4(*pos, vec)
+        } else {
+            zero
+        }
+    }
+
+    // Constants.
+    let stretch_constant: T = math::cast(STRETCH_CONSTANT_4D);
+    let squish_constant: T = math::cast(SQUISH_CONSTANT_4D);
+    let zero: T = math::cast(0.0);
+    let one: T = math::cast(1.0);
+    let two: T = math::cast(2.0);
+    let three: T = math::cast(3.0);
+
+    // Place input coordinates on simplectic honeycomb.
+    let stretch_offset = math::fold4(*point, Add::add) * stretch_constant;
+    let stretched = math::map4(*point, |v| v + stretch_offset);
+
+    // Floor to get simplectic honeycomb coordinates of rhombo-hypercube
+    // super-cell origin.
+    let stretched_floor = math::map4(stretched, Float::floor);
+
+    // Skew out to get actual coordinates of stretched rhombo-hypercube origin.
+    // We'll need these later.
+    let squish_offset = math::fold4(stretched_floor, Add::add) * squish_constant;
+    let skewed_floor = math::map4(stretched_floor, |v| v + squish_offset);
+
+    // Compute simplectic honeycomb coordinates relative to rhombo-hypercube
+    // origin.
+    let rel_coords = math::sub4(stretched, stretched_floor);
+
+    // Sum those together to get a value that determines which region
+    // we're in.
+    let region_sum = math::fold4(rel_coords, Add::add);
+
+    // Position relative to origin point.
+    let pos0 = math::sub4(*point, skewed_floor);
+
+    let mut value = zero;
+    if region_sum <= one {
+        // We're inside the pentachoron (4-Simplex) at (0, 0, 0, 0)
+
+        // Contribution at (0, 0, 0, 0)
+        value = value + gradient(seed, &stretched_floor, &pos0);
+
+        // Contribution at (1, 0, 0, 0)
+        let pos1;
+        {
+            let vertex = math::add4(stretched_floor, [one, zero, zero, zero]);
+            pos1 = math::sub4(pos0, [
+                one + squish_constant,
+                squish_constant,
+                squish_constant,
+                squish_constant
+            ]);
+            value = value + gradient(seed, &vertex, &pos1);
+        }
+
+        // Contribution at (0, 1, 0, 0)
+        let pos2;
+        {
+            let vertex = math::add4(stretched_floor, [zero, one, zero, zero]);
+            pos2 = [
+                pos1[0] - one,
+                pos1[1] + one,
+                pos1[2],
+                pos1[3]
+            ];
+            value = value + gradient(seed, &vertex, &pos2);
+        }
+
+        // TODO
+		// Contribution at (0, 0, 1, 0)
+		// Contribution at (0, 0, 0, 1)
+    } else if region_sum >= three {
+        // We're inside the pentachoron (4-Simplex) at (1, 1, 1, 1)
+        let squish_constant_3 = three * squish_constant;
+
+        // Contribution at (1, 1, 1, 0)
+        let pos4;
+        {
+            let vertex = math::add4(stretched_floor, [one, one, one, zero]);
+            pos4 = math::sub4(pos0, [
+                one + squish_constant_3,
+                one + squish_constant_3,
+                one + squish_constant_3,
+                squish_constant_3
+            ]);
+            value = value + gradient(seed, &vertex, &pos4);
+        }
+
+        // Contribution at (1, 1, 0, 1)
+        let pos3;
+        {
+            let vertex = math::add4(stretched_floor, [one, one, zero, one]);
+            pos3 = [
+                pos4[0],
+                pos4[1],
+                pos4[2] + one,
+                pos4[3] - one
+            ];
+            value = value + gradient(seed, &vertex, &pos3);
+        }
+
+        // TODO
+		// Contribution at (1, 0, 1, 1)
+		// Contribution at (0, 1, 1, 1)
+    } else if region_sum <= two {
+        // We're inside the first dispentachoron (Rectified 4-Simplex)
+
+        // Contribution at (1, 0, 0, 0)
+        let pos1;
+        {
+            let vertex = math::add4(stretched_floor, [one, zero, zero, zero]);
+            pos1 = math::sub4(pos0, [
+                one + squish_constant,
+                squish_constant,
+                squish_constant,
+                squish_constant
+            ]);
+            value = value + gradient(seed, &vertex, &pos1);
+        }
+
+        // Contribution at (0, 1, 0, 0)
+        let pos2;
+        {
+            let vertex = math::add4(stretched_floor, [zero, one, zero, zero]);
+            pos2 = [
+                pos1[0] + one,
+                pos1[1] - one,
+                pos1[2],
+                pos1[3]
+            ];
+            value = value + gradient(seed, &vertex, &pos2);
+        }
+
+        // TODO
+        // Contribution at (0, 0, 1, 0)
+        // Contribution at (0, 0, 0, 1)
+        // Contribution at (1, 1, 0, 0)
+        // Contribution at (1, 0, 1, 0)
+        // Contribution at (1, 0, 0, 1)
+        // Contribution at (0, 1, 1, 0)
+        // Contribution at (0, 1, 0, 1)
+        // Contribution at (0, 0, 1, 1)
+    } else {
+        // We're inside the second dispentachoron (Rectified 4-Simplex)
+        let squish_constant_3 = three * squish_constant;
+
+        // Contribution at (1, 1, 1, 0)
+        let pos4;
+        {
+            let vertex = math::add4(stretched_floor, [one, one, one, zero]);
+            pos4 = math::sub4(pos0, [
+                one + squish_constant_3,
+                one + squish_constant_3,
+                one + squish_constant_3,
+                squish_constant_3
+            ]);
+            value = value + gradient(seed, &vertex, &pos4);
+        }
+
+        // Contribution at (1, 1, 0, 1)
+        let pos3;
+        {
+            let vertex = math::add4(stretched_floor, [zero, one, zero, zero]);
+            pos3 = [
+                pos4[0],
+                pos4[1],
+                pos4[2] + one,
+                pos4[3] - one
+            ];
+            value = value + gradient(seed, &vertex, &pos3);
+        }
+
+        // oh my god this is tiring.
+
+        // Contribution at (1, 0, 1, 1)
+        let pos2;
+        {
+            let vertex = math::add4(stretched_floor, [one, zero, one, one]);
+            pos2 = [
+                pos4[0],
+                pos4[1] + one,
+                pos4[2],
+                pos3[3]
+            ];
+            value = value + gradient(seed, &vertex, &pos2);
+        }
+
+        // Contribution at (0, 1, 1, 1)
+        let pos1;
+        {
+            let vertex = math::add4(stretched_floor, [zero, one, one, one]);
+            pos1 = [
+                pos4[0] + one,
+                pos4[1],
+                pos4[2],
+                pos3[3]
+            ];
+            value = value + gradient(seed, &vertex, &pos1);
+        }
+    }
+
+    value
+}