--wip-- [skip ci]

src/interpreter.rs

@@ -5,19 +5,44 @@ use std::simd::{
 };
 
 use crate::{
-    CSG,
+    BYTECODE_IMITATES_GLSL, CSG,
     ssa::{SSAInput, SSAInstruction, SSAOpcode},
 };
 
 pub type Value = std::simd::f32x8;
+pub type Mask = std::simd::mask32x8;
+
+pub const VALUE_NAN: Value = Value::splat(core::f32::NAN);
+pub const VALUE_1: Value = Value::splat(1.0);
+pub const VALUE_0: Value = Value::splat(0.0);
+pub const VALUE_05: Value = Value::splat(0.5);
+pub const VALUE_M1: Value = Value::splat(-1.0);
+pub const VALUE_2: Value = Value::splat(2.0);
+pub const VALUE_PI: Value = Value::splat(core::f32::consts::PI);
+
+fn glsign(f: Value) -> Value {
+    if BYTECODE_IMITATES_GLSL {
+        f.simd_eq(VALUE_0).select(f, f.signum())
+    } else {
+        f.signum()
+    }
+}
+
+fn glfract(f: Value) -> Value {
+    if BYTECODE_IMITATES_GLSL {
+        f - f.floor()
+    } else {
+        f.fract()
+    }
+}
 
 #[derive(Clone, Debug)]
-pub(crate) struct Interpreter<'csg> {
+pub(crate) struct PointInterpreter<'csg> {
     value_map: Vec<Value>,
     csg:       &'csg CSG,
 }
 
-impl Interpreter<'_> {
+impl PointInterpreter<'_> {
     fn load(&self, consider: SSAInput) -> Value {
         match consider {
             SSAInput::Register(r) => self.value_map[r as usize],
@@ -30,25 +55,16 @@ impl Interpreter<'_> {
     }
 }
 
-//imitate glsl, not rust
+impl<'csg> PointInterpreter<'csg> {
-fn sign(f: Value) -> Value {
-    f.simd_eq(Value::splat(0.0)).select(f, f.signum())
-}
-
-fn fract(f: Value) -> Value {
-    f - f.floor()
-}
-
-impl<'csg> Interpreter<'csg> {
     pub(crate) fn new(csg: &'csg CSG) -> Self {
-        Interpreter {
+        PointInterpreter {
-            value_map: vec![Value::splat(0.); csg.parts.last_output as usize],
+            value_map: vec![VALUE_0; csg.parts.last_output as usize],
             csg,
         }
     }
 
     fn clear_stacks(&mut self) {
-        self.value_map = vec![Value::splat(0.); self.csg.parts.last_output as usize];
+        self.value_map = vec![VALUE_0; self.csg.parts.last_output as usize];
     }
 
     fn param_one(&mut self, instruction: &SSAInstruction, func: impl Fn(Value) -> Value) {
@@ -93,7 +109,8 @@ impl<'csg> Interpreter<'csg> {
         }
     }
 
-    // cargo asm "erroccfisuvsc::interpreter::Interpreter::scene" --no-color --rust > scene.asm
+    // cargo asm "erroccfisuvsc::interpreter::Interpreter::scene" --no-color --rust
+    // > scene.asm
 
     pub(crate) fn scene(&mut self, px: Value, py: Value, pz: Value) -> Value {
         self.clear_stacks();
@@ -236,7 +253,7 @@ impl<'csg> Interpreter<'csg> {
 
                     match instruction.opcode.size {
                         1 => {
-                            self.store(instruction.outputs[0], Value::splat(1.0));
+                            self.store(instruction.outputs[0], VALUE_1);
                         },
                         2 => {
                             let max = val_a[0].simd_max(val_a[1]);
@@ -278,7 +295,7 @@ impl<'csg> Interpreter<'csg> {
                     self.param_one(instruction, |val_a| val_a.ceil());
                 },
                 SSAFract => {
-                    self.param_one(instruction, |val_a| fract(val_a));
+                    self.param_one(instruction, |val_a| glfract(val_a));
                 },
                 SSASin => {
                     self.param_one(instruction, |val_a| val_a.sin());
@@ -323,18 +340,16 @@ impl<'csg> Interpreter<'csg> {
                 },
                 SSASmoothMin => {
                     self.param_three(instruction, |d1, d2, k| {
-                        let h = (Value::splat(0.5) + (Value::splat(0.5) * (d2 - d1) / k))
+                        let h =
-                            .simd_clamp(Value::splat(0.0), Value::splat(1.0));
+                            (VALUE_05 + (VALUE_05 * (d2 - d1) / k)).simd_clamp(VALUE_0, VALUE_1);
-                        return ((d2 * (Value::splat(1.0) - h)) + (d1 * h))
+                        return ((d2 * (VALUE_1 - h)) + (d1 * h)) - k * h * (VALUE_1 - h);
-                            - k * h * (Value::splat(1.0) - h);
                     });
                 },
                 SSASmoothMax => {
                     self.param_three(instruction, |d1, d2, k| {
-                        let h = (Value::splat(0.5) - (Value::splat(0.5) * (d2 + d1) / k))
+                        let h =
-                            .simd_clamp(Value::splat(0.0), Value::splat(1.0));
+                            (VALUE_05 - (VALUE_05 * (d2 + d1) / k)).simd_clamp(VALUE_0, VALUE_1);
-                        return ((d2 * (Value::splat(1.0) - h)) + (-d1 * h))
+                        return ((d2 * (VALUE_1 - h)) + (-d1 * h)) + k * h * (VALUE_1 - h);
-                            + k * h * (Value::splat(1.0) - h);
                     });
                 },
                 SSASmoothMinMaterial => todo!(),
@@ -346,7 +361,7 @@ impl<'csg> Interpreter<'csg> {
                 },
                 SSAMix => {
                     self.param_three(instruction, |val_a, val_b, val_c| {
-                        (val_a * (Value::splat(1.0) - val_c)) + (val_b * val_c)
+                        (val_a * (VALUE_1 - val_c)) + (val_b * val_c)
                     });
                 },
                 SSAFMA => {
@@ -374,13 +389,13 @@ impl<'csg> Interpreter<'csg> {
                     let qx = pos_x.abs() - rad_x;
                     let qy = pos_y.abs() - rad_y;
                     let qz = pos_z.abs() - rad_z;
-                    let qxmax = qx.simd_max(Value::splat(0.0));
+                    let qxmax = qx.simd_max(VALUE_0);
-                    let qymax = qy.simd_max(Value::splat(0.0));
+                    let qymax = qy.simd_max(VALUE_0);
-                    let qzmax = qz.simd_max(Value::splat(0.0));
+                    let qzmax = qz.simd_max(VALUE_0);
                     self.store(
                         instruction.outputs[0],
                         ((qxmax * qxmax) + (qymax * qymax) + (qzmax * qzmax)).sqrt()
-                            + qx.simd_max(qy.simd_max(qz)).simd_min(Value::splat(0.0)),
+                            + qx.simd_max(qy.simd_max(qz)).simd_min(VALUE_0),
                     );
                 },
                 SSASDFTorus => {
@@ -397,22 +412,19 @@ impl<'csg> Interpreter<'csg> {
                 },
                 SSACompare => {
                     self.param_two(instruction, |val_a, val_b| {
-                        val_a.simd_gt(val_b).select(
-                            Value::splat(1.0),
                         val_a
-                                .simd_lt(val_b)
+                            .simd_gt(val_b)
-                                .select(Value::splat(-1.0), Value::splat(0.0)),
+                            .select(VALUE_1, val_a.simd_lt(val_b).select(VALUE_M1, VALUE_0))
-                        )
                     });
                 },
                 SSAAnd => {
                     self.param_two(instruction, |val_a, val_b| {
-                        val_a.simd_eq(Value::splat(0.0)).select(val_a, val_b)
+                        val_a.simd_eq(VALUE_0).select(val_a, val_b)
                     });
                 },
                 SSAOr => {
                     self.param_two(instruction, |val_a, val_b| {
-                        val_a.simd_eq(Value::splat(0.0)).select(val_b, val_a)
+                        val_a.simd_eq(VALUE_0).select(val_b, val_a)
                     });
                 },
                 SSARecip => {
@@ -420,14 +432,12 @@ impl<'csg> Interpreter<'csg> {
                 },
                 SSANot => {
                     self.param_one(instruction, |val_a| {
-                        val_a
+                        val_a.simd_eq(VALUE_0).select(VALUE_1, VALUE_0)
-                            .simd_eq(Value::splat(0.0))
-                            .select(Value::splat(1.0), Value::splat(0.0))
                     });
                 },
-                SSAStop => return Value::splat(0.0),
+                SSAStop => return VALUE_0,
             }
         }
-        return Value::splat(f32::NAN);
+        return VALUE_NAN;
     }
 }

src/interval.rs

@@ -0,0 +1,621 @@
+// Based on https://github.com/mkeeter/fidget/blob/main/fidget/src/core/types/interval.rs
+
+use std::simd::{
+    StdFloat,
+    cmp::{SimdPartialEq, SimdPartialOrd},
+    num::SimdFloat,
+};
+
+use crate::interpreter::{
+    Mask, VALUE_0, VALUE_1, VALUE_2, VALUE_05, VALUE_M1, VALUE_NAN, VALUE_PI, Value,
+};
+
+/// Stores a range, with conservative calculations to guarantee that it always
+/// contains the actual value.
+///
+/// # Warning
+/// This implementation does not set rounding modes, so it may not be _perfect_.
+#[derive(Copy, Clone, PartialEq)]
+#[repr(C)]
+pub struct Interval {
+    lower: Value,
+    upper: Value,
+}
+
+impl std::fmt::Debug for Interval {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
+        f.debug_tuple("")
+            .field(&self.lower)
+            .field(&self.upper)
+            .finish()
+    }
+}
+
+impl Interval {
+    pub const HALF: Self = Self::const_splat(0.5);
+    pub const MONE: Self = Self::const_splat(-1.0);
+    pub const NAN: Self = Self::const_splat(core::f32::NAN);
+    pub const ONE: Self = Self::const_splat(1.0);
+    pub const PI: Self = Self::const_splat(core::f32::consts::PI);
+    pub const ZERO: Self = Self::const_splat(0.0);
+
+    /// Builds a new interval
+    ///
+    /// There are two kinds of valid interval:
+    /// - `[lower, upper]` where `lower <= upper`
+    /// - `[NaN, NaN]`
+    ///
+    /// # Panics
+    /// Panics if the resulting interval would be invalid
+    #[inline]
+    pub fn new(lower: Value, upper: Value) -> Self {
+        assert!(
+            (upper.simd_ge(lower) | (lower.is_nan() & upper.is_nan())).all(),
+            "invalid interval [{lower:?}, {upper:?}]"
+        );
+        Self { lower, upper }
+    }
+
+    pub fn splat(value: f32) -> Interval {
+        Interval::new(Value::splat(value), Value::splat(value))
+    }
+
+    pub fn splat2(lower: f32, upper: f32) -> Interval {
+        Interval::new(Value::splat(lower), Value::splat(upper))
+    }
+
+    const fn const_splat(value: f32) -> Interval {
+        Interval {
+            lower: Value::splat(value),
+            upper: Value::splat(value),
+        }
+    }
+
+    /// Returns the lower bound of the interval
+    #[inline]
+    pub fn lower(&self) -> Value {
+        self.lower
+    }
+
+    /// Returns the upper bound of the interval
+    #[inline]
+    pub fn upper(&self) -> Value {
+        self.upper
+    }
+
+    /// Checks whether the given value is (strictly) contained in the interval
+    #[inline]
+    pub fn contains(&self, v: Value) -> Mask {
+        v.simd_ge(self.lower) & v.simd_le(self.upper)
+    }
+
+    /// Returns `true` if either bound of the interval is `NaN`
+    #[inline]
+    pub fn has_nan(&self) -> Mask {
+        self.lower.is_nan() | self.upper.is_nan()
+    }
+
+    /// Calculates the absolute value of the interval
+    #[inline]
+    pub fn abs(self) -> Self {
+        let llt = self.lower.simd_lt(VALUE_0);
+        let ugt = self.upper.simd_gt(VALUE_0);
+        let lower = llt.select(ugt.select(VALUE_0, -self.upper), self.lower);
+        let upper = llt.select(
+            ugt.select(self.upper.simd_max(-self.lower), -self.lower),
+            self.upper,
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Squares the interval
+    ///
+    /// Note that this has tighter bounds than multiplication, because we know
+    /// that both sides of the multiplication are the same value.
+    #[inline]
+    pub fn square(self) -> Self {
+        let ult = self.upper.simd_lt(VALUE_0);
+        let lgt = self.lower.simd_gt(VALUE_0);
+        let has_nan = self.has_nan();
+        let lower = ult.select(
+            self.upper * self.upper,
+            lgt.select(self.lower * self.lower, has_nan.select(VALUE_NAN, VALUE_0)),
+        );
+        let upper = ult.select(
+            self.lower * self.lower,
+            lgt.select(
+                self.upper * self.upper,
+                has_nan.select(VALUE_NAN, {
+                    let k = self.lower.abs().simd_max(self.upper.abs());
+                    k * k
+                }),
+            ),
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the sine of the interval
+    #[inline]
+    pub fn sin(self) -> Self {
+        let same_cycle = ((self.lower / VALUE_0) + VALUE_05)
+            .floor()
+            .simd_eq(((self.upper / VALUE_PI) + VALUE_05).floor());
+        let up = (((self.upper / VALUE_PI) + VALUE_05).floor()) % (VALUE_2);
+        let whole_cycle = (((self.upper / VALUE_PI) + VALUE_05).floor()
+            - ((self.lower / VALUE_PI) + VALUE_05).floor())
+        .simd_gt(VALUE_1);
+        let temp0 = self.lower.sin();
+        let temp1 = self.upper.sin();
+        let lower = self.has_nan().select(
+            VALUE_NAN,
+            ((!whole_cycle & (up.simd_eq(VALUE_1))) | same_cycle)
+                .select(temp0.simd_min(temp1), VALUE_M1),
+        );
+        let upper = self.has_nan().select(
+            VALUE_NAN,
+            ((!whole_cycle & (up.simd_eq(VALUE_0))) | same_cycle)
+                .select(temp0.simd_max(temp1), VALUE_1),
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the cosine of the interval
+    #[inline]
+    pub fn cos(self) -> Self {
+        let same_cycle = (self.lower / VALUE_0)
+            .floor()
+            .simd_eq((self.upper / VALUE_PI).floor());
+        let up = ((self.upper / VALUE_PI).floor()) % (VALUE_2);
+        let whole_cycle =
+            ((self.upper / VALUE_PI).floor() - (self.lower / VALUE_PI).floor()).simd_gt(VALUE_1);
+        let temp0 = self.lower.cos();
+        let temp1 = self.upper.cos();
+        let lower = self.has_nan().select(
+            VALUE_NAN,
+            ((!whole_cycle & (up.simd_eq(VALUE_0))) | same_cycle)
+                .select(temp0.simd_min(temp1), VALUE_M1),
+        );
+        let upper = self.has_nan().select(
+            VALUE_NAN,
+            ((!whole_cycle & (up.simd_eq(VALUE_1))) | same_cycle)
+                .select(temp0.simd_max(temp1), VALUE_1),
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the tangent of the interval
+    ///
+    /// Returns the `NAN` interval if the result contains a undefined point
+    #[inline]
+    pub fn tan(self) -> Self {
+        let size = self.upper - self.lower;
+        let lower_tmp = Value::from_array(self.lower.to_array().map(|f| f.tan()));
+        let upper_tmp = Value::from_array(self.upper.to_array().map(|f| f.tan()));
+        let lower =
+            (size.simd_lt(VALUE_PI) & upper_tmp.simd_ge(lower_tmp)).select(lower_tmp, VALUE_NAN);
+        let upper =
+            (size.simd_lt(VALUE_PI) & upper_tmp.simd_ge(lower_tmp)).select(upper_tmp, VALUE_NAN);
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the arcsine of the interval
+    ///
+    /// Returns the `NAN` interval if the input is invalid
+    #[inline]
+    pub fn asin(self) -> Self {
+        let lower = (self.lower.simd_lt(VALUE_M1) | self.upper.simd_gt(VALUE_1)).select(
+            VALUE_NAN,
+            Value::from_array(self.lower.to_array().map(|f| f.asin())),
+        );
+        let upper = (self.lower.simd_lt(VALUE_M1) | self.upper.simd_gt(VALUE_1)).select(
+            VALUE_NAN,
+            Value::from_array(self.upper.to_array().map(|f| f.asin())),
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the arccosine of the interval
+    ///
+    /// Returns the `NAN` interval if the input is invalid
+    #[inline]
+    pub fn acos(self) -> Self {
+        let lower = (self.lower.simd_lt(VALUE_M1) | self.upper.simd_gt(VALUE_1)).select(
+            VALUE_NAN,
+            Value::from_array(self.upper.to_array().map(|f| f.asin())),
+        );
+        let upper = (self.lower.simd_lt(VALUE_M1) | self.upper.simd_gt(VALUE_1)).select(
+            VALUE_NAN,
+            Value::from_array(self.lower.to_array().map(|f| f.asin())),
+        );
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the arctangent of the interval
+    #[inline]
+    pub fn atan(self) -> Self {
+        let lower = Value::from_array(self.lower.to_array().map(|f| f.asin()));
+        let upper = Value::from_array(self.upper.to_array().map(|f| f.asin()));
+        Interval::new(lower, upper)
+    }
+
+    /// Computes the exponent function applied to the interval
+    #[inline]
+    pub fn exp(self) -> Self {
+        Interval::new(self.lower.exp(), self.upper.exp())
+    }
+
+    /// Computes the natural log of the input interval
+    ///
+    /// Returns the `NAN` interval if the input contains zero
+    #[inline]
+    pub fn ln(self) -> Self {
+        if self.lower <= 0.0 {
+            f32::NAN.into()
+        } else {
+            Interval::new(self.lower.ln(), self.upper.ln())
+        }
+    }
+
+    /// Calculates the square root of the interval
+    ///
+    /// If the interval contains values below 0, returns a `NAN` interval.
+    #[inline]
+    pub fn sqrt(self) -> Self {
+        if self.lower < 0.0 {
+            f32::NAN.into()
+        } else {
+            Interval::new(self.lower.sqrt(), self.upper.sqrt())
+        }
+    }
+
+    /// Calculates the reciprocal of the interval
+    ///
+    /// If the interval includes 0, returns the `NAN` interval
+    #[inline]
+    pub fn recip(self) -> Self {
+        if self.lower > 0.0 || self.upper < 0.0 {
+            Interval::new(1.0 / self.upper, 1.0 / self.lower)
+        } else {
+            f32::NAN.into()
+        }
+    }
+
+    /// Calculates the minimum of two intervals
+    ///
+    /// Returns both the result and a [`Choice`] indicating whether one side is
+    /// always less than the other.
+    ///
+    /// If either side is `NAN`, returns the `NAN` interval and `Choice::Both`.
+    #[inline]
+    pub fn min_choice(self, rhs: Self) -> (Self, Choice) {
+        if self.has_nan() || rhs.has_nan() {
+            return (f32::NAN.into(), Choice::Both);
+        }
+        let choice = if self.upper < rhs.lower {
+            Choice::Left
+        } else if rhs.upper < self.lower {
+            Choice::Right
+        } else {
+            Choice::Both
+        };
+        (
+            Interval::new(self.lower.min(rhs.lower), self.upper.min(rhs.upper)),
+            choice,
+        )
+    }
+
+    /// Calculates the maximum of two intervals
+    ///
+    /// Returns both the result and a [`Choice`] indicating whether one side is
+    /// always greater than the other.
+    ///
+    /// If either side is `NAN`, returns the `NAN` interval and `Choice::Both`.
+    #[inline]
+    pub fn max_choice(self, rhs: Self) -> (Self, Choice) {
+        if self.has_nan() || rhs.has_nan() {
+            return (f32::NAN.into(), Choice::Both);
+        }
+        let choice = if self.lower > rhs.upper {
+            Choice::Left
+        } else if rhs.lower > self.upper {
+            Choice::Right
+        } else {
+            Choice::Both
+        };
+        (
+            Interval::new(self.lower.max(rhs.lower), self.upper.max(rhs.upper)),
+            choice,
+        )
+    }
+
+    /// Calculates the short-circuiting `AND` of two intervals
+    ///
+    /// Returns both the result and a [`Choice`] indicating whether one side is
+    /// always selected.  An unambiguous 0 in `self` selects itself; an
+    /// unambiguous 1 selects the opposite branch.
+    #[inline]
+    pub fn and_choice(self, rhs: Self) -> (Self, Choice) {
+        if self.has_nan() || rhs.has_nan() {
+            (f32::NAN.into(), Choice::Both)
+        } else if self.lower == 0.0 && self.upper == 0.0 {
+            (0.0.into(), Choice::Left)
+        } else if !self.contains(0.0) {
+            (rhs, Choice::Right)
+        } else {
+            // The output will either be the RHS or zero, so extend the interval
+            // to include zero in it.
+            (
+                Interval::new(rhs.lower.min(0.0), rhs.upper.max(0.0)),
+                Choice::Both,
+            )
+        }
+    }
+
+    /// Calculates the short-circuiting `OR` of two intervals
+    ///
+    /// Returns both the result and a [`Choice`] indicating whether one side is
+    /// always selected.  An unambiguous 0 in `self` selects the opposite
+    /// branch; an unambiguous 1 selects itself.
+    #[inline]
+    pub fn or_choice(self, rhs: Self) -> (Self, Choice) {
+        if self.has_nan() || rhs.has_nan() {
+            (f32::NAN.into(), Choice::Both)
+        } else if !self.contains(0.0) {
+            (self, Choice::Left)
+        } else if self.lower == 0.0 && self.upper == 0.0 {
+            (rhs, Choice::Right)
+        } else {
+            // The output could be anywhere in either interval
+            (
+                Interval::new(self.lower.min(rhs.lower), self.upper.max(rhs.upper)),
+                Choice::Both,
+            )
+        }
+    }
+
+    /// Returns the midpoint of the interval
+    #[inline]
+    pub fn midpoint(self) -> f32 {
+        (self.lower + self.upper) / 2.0
+    }
+
+    /// Splits the interval at the midpoint
+    ///
+    /// ```
+    /// # use fidget::types::Interval;
+    /// let a = Interval::new(0.0, 1.0);
+    /// let (lo, hi) = a.split();
+    /// assert_eq!(lo, Interval::new(0.0, 0.5));
+    /// assert_eq!(hi, Interval::new(0.5, 1.0));
+    /// ```
+    #[inline]
+    pub fn split(self) -> (Self, Self) {
+        let mid = self.midpoint();
+        (
+            Interval::new(self.lower, mid),
+            Interval::new(mid, self.upper),
+        )
+    }
+
+    /// Linear interpolation from `lower` to `upper`
+    ///
+    /// ```
+    /// # use fidget::types::Interval;
+    /// let a = Interval::new(0.0, 2.0);
+    /// assert_eq!(a.lerp(0.5), 1.0);
+    /// assert_eq!(a.lerp(0.75), 1.5);
+    /// assert_eq!(a.lerp(2.0), 4.0);
+    /// ```
+    #[inline]
+    pub fn lerp(self, frac: f32) -> f32 {
+        self.lower * (1.0 - frac) + self.upper * frac
+    }
+
+    /// Calculates the width of the interval
+    ///
+    /// ```
+    /// # use fidget::types::Interval;
+    /// let a = Interval::new(2.0, 3.0);
+    /// assert_eq!(a.width(), 1.0);
+    /// let b = Interval::new(2.0, 5.0);
+    /// assert_eq!(b.width(), 3.0);
+    /// ```
+    #[inline]
+    pub fn width(self) -> f32 {
+        self.upper - self.lower
+    }
+
+    /// Checks that the two values are roughly equal, panicking otherwise
+    #[cfg(test)]
+    pub(crate) fn compare_eq(&self, other: Self) {
+        let d = (self.lower - other.lower)
+            .abs()
+            .max((self.upper - other.upper).abs());
+        if d >= 1e-6 {
+            panic!("lhs != rhs ({self:?} != {other:?})");
+        }
+    }
+
+    /// Least non-negative remainder
+    #[inline]
+    pub fn rem_euclid(&self, other: Interval) -> Self {
+        // TODO optimize this more?
+        if self.has_nan() || other.has_nan() || other.contains(0.0) {
+            f32::NAN.into()
+        } else if other.lower == other.upper && other.lower > 0.0 {
+            let a = self.lower / other.lower;
+            let b = self.upper / other.lower;
+            if a != a.floor() && a.floor() == b.floor() {
+                Interval::new(
+                    self.lower.rem_euclid(other.lower),
+                    self.upper.rem_euclid(other.lower),
+                )
+            } else {
+                Interval::new(0.0, other.abs().upper())
+            }
+        } else {
+            Interval::new(0.0, other.abs().upper())
+        }
+    }
+
+    /// Largest value that is less-than-or-equal to this value
+    #[inline]
+    pub fn floor(&self) -> Self {
+        Interval::new(self.lower.floor(), self.upper.floor())
+    }
+
+    /// Smallest value that is greater-than-or-equal to this value
+    #[inline]
+    pub fn ceil(&self) -> Self {
+        Interval::new(self.lower.ceil(), self.upper.ceil())
+    }
+
+    /// Rounded value
+    #[inline]
+    pub fn round(&self) -> Self {
+        Interval::new(self.lower.round(), self.upper.round())
+    }
+
+    /// Four-quadrant arctangent
+    #[inline]
+    pub fn atan2(self, x: Self) -> Self {
+        if self.has_nan() || x.has_nan() {
+            f32::NAN.into()
+        } else {
+            // TODO optimize this further
+            Interval::new(-std::f32::consts::PI, std::f32::consts::PI)
+        }
+    }
+}
+
+impl std::fmt::Display for Interval {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        write!(f, "({}, {})", self.lower, self.upper)
+    }
+}
+
+impl From<[f32; 2]> for Interval {
+    #[inline]
+    fn from(i: [f32; 2]) -> Interval {
+        Interval::new(i[0], i[1])
+    }
+}
+
+impl From<f32> for Interval {
+    #[inline]
+    fn from(f: f32) -> Self {
+        Interval::new(f, f)
+    }
+}
+
+impl std::ops::Add<Interval> for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn add(self, rhs: Self) -> Self {
+        Interval::new(self.lower + rhs.lower, self.upper + rhs.upper)
+    }
+}
+
+impl std::ops::Mul<Interval> for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, rhs: Self) -> Self {
+        if self.has_nan() || rhs.has_nan() {
+            return f32::NAN.into();
+        }
+        let mut out = [0.0; 4];
+        let mut k = 0;
+        for i in [self.lower, self.upper] {
+            for j in [rhs.lower, rhs.upper] {
+                out[k] = i * j;
+                k += 1;
+            }
+        }
+        let mut lower = out[0];
+        let mut upper = out[0];
+        for &v in &out[1..] {
+            lower = lower.min(v);
+            upper = upper.max(v);
+        }
+        Interval::new(lower, upper)
+    }
+}
+
+impl std::ops::Mul<f32> for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn mul(self, rhs: f32) -> Self {
+        if self.has_nan() || rhs.is_nan() {
+            f32::NAN.into()
+        } else if rhs < 0.0 {
+            Interval::new(self.upper * rhs, self.lower * rhs)
+        } else {
+            Interval::new(self.lower * rhs, self.upper * rhs)
+        }
+    }
+}
+
+impl std::ops::Div<Interval> for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn div(self, rhs: Self) -> Self {
+        if self.has_nan() {
+            return f32::NAN.into();
+        }
+        if rhs.lower > 0.0 || rhs.upper < 0.0 {
+            let mut out = [0.0; 4];
+            let mut k = 0;
+            for i in [self.lower, self.upper] {
+                for j in [rhs.lower, rhs.upper] {
+                    out[k] = i / j;
+                    k += 1;
+                }
+            }
+            let mut lower = out[0];
+            let mut upper = out[0];
+            for &v in &out[1..] {
+                lower = lower.min(v);
+                upper = upper.max(v);
+            }
+            Interval::new(lower, upper)
+        } else {
+            f32::NAN.into()
+        }
+    }
+}
+
+impl std::ops::Sub<Interval> for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn sub(self, rhs: Self) -> Self {
+        Interval::new(self.lower - rhs.upper, self.upper - rhs.lower)
+    }
+}
+
+impl std::ops::Neg for Interval {
+    type Output = Self;
+
+    #[inline]
+    fn neg(self) -> Self {
+        Interval::new(-self.upper, -self.lower)
+    }
+}
+
+#[cfg(test)]
+mod test {
+    use super::*;
+
+    #[test]
+    fn test_interval() {
+        let a = Interval::new(0.0, 1.0);
+        let b = Interval::new(0.5, 1.5);
+        let (v, c) = a.min_choice(b);
+        assert_eq!(v, [0.0, 1.0].into());
+        assert_eq!(c, Choice::Both);
+    }
+}

src/main.rs

@@ -2,6 +2,16 @@
 #![feature(array_chunks)]
 #![feature(portable_simd)]
 
+const DYNAMIC_STATE: bool = true;
+const MSAA_ENABLE: bool = false;
+const MSAA_SAMPLES: u32 = 4;
+const SAMPLE_RATE_SHADING: Option<f32> = None; //1.0;
+
+const DUMP_SPV_TO_FILE: bool = false;
+const PIPELINE_CACHING: bool = false;
+
+const BYTECODE_IMITATES_GLSL: bool = false;
+
 use std::{
     error::Error,
     fs::{File, remove_file, rename},
@@ -12,13 +22,6 @@ use std::{
     time::Instant,
 };
 
-const DYNAMIC_STATE: bool = true;
-const MSAA_ENABLE: bool = false;
-const MSAA_SAMPLES: u32 = 4;
-const SAMPLE_RATE_SHADING: Option<f32> = None; //1.0;
-
-const MSAA_SAMPLES_ACTUAL: u32 = if MSAA_ENABLE { MSAA_SAMPLES } else { 1 };
-
 use bytemuck::{Pod, Zeroable};
 use egui_winit_vulkano::{Gui, GuiConfig};
 use foldhash::{HashMap, HashMapExt, HashSet};
@@ -94,7 +97,10 @@ mod ssa;
 mod instruction_set;
 
 mod interpreter;
-use crate::interpreter::Interpreter;
+mod interval;
+use crate::interpreter::{IntervalInterpreter, PointInterpreter};
+
+const MSAA_SAMPLES_ACTUAL: u32 = if MSAA_ENABLE { MSAA_SAMPLES } else { 1 };
 
 fn main() -> Result<(), impl Error> {
     CombinedLogger::init(vec![
@@ -626,7 +632,12 @@ impl ApplicationHandler for App {
             .flatten()
             .collect::<Vec<_>>();
 
-        File::create("out.spv").unwrap().write(&module_as_bytes).unwrap();
+        if DUMP_SPV_TO_FILE {
+            File::create("out.spv")
+                .unwrap()
+                .write(&module_as_bytes)
+                .unwrap();
+        }
 
         let smodule = unsafe {
             ::vulkano::shader::ShaderModule::new(
@@ -918,7 +929,12 @@ impl App {
                 .flatten()
                 .collect::<Vec<_>>();
 
-            File::create("out.spv").unwrap().write(&module_as_bytes).unwrap();
+            if DUMP_SPV_TO_FILE {
+                File::create("out.spv")
+                    .unwrap()
+                    .write(&module_as_bytes)
+                    .unwrap();
+            }
 
             let smodule = unsafe {
                 ::vulkano::shader::ShaderModule::new(
@@ -1162,13 +1178,16 @@ impl App {
                 10000.0, 10000.0, 10000.0, -10000.0, -10000.0, -10000.0, 0.0, 0.0,
             ]);
 
-            let mut interpreter = Interpreter::new(csg);
+            let mut interpreter = PointInterpreter::new(csg);
             let interpreter_out = interpreter
                 .scene(INTERPRET_INPUT_X, INTERPRET_INPUT_Y, INTERPRET_INPUT_Z)
-                .mul_add(INTERPRET_MUL, INTERPRET_ADD).to_array();
+                .mul_add(INTERPRET_MUL, INTERPRET_ADD)
+                .to_array();
 
-            let highest_corner = Vec3::new(interpreter_out[0], interpreter_out[1], interpreter_out[2]);
+            let highest_corner =
-            let lowest_corner = Vec3::new(interpreter_out[3], interpreter_out[4], interpreter_out[5]);
+                Vec3::new(interpreter_out[0], interpreter_out[1], interpreter_out[2]);
+            let lowest_corner =
+                Vec3::new(interpreter_out[3], interpreter_out[4], interpreter_out[5]);
 
             push_constants.dimensions = (highest_corner - lowest_corner).to_array().into();
             push_constants.lowest_corner = lowest_corner.into();
@@ -1512,6 +1531,7 @@ where
 }
 
 fn dump_pipeline_cache(cache: Arc<PipelineCache>) {
+    if PIPELINE_CACHING {
         if let Ok(data) = cache.get_data() {
             info!("Dumping pipeline cache...");
             let tmp_path = relpath("pipeline_cache.bin.tmp");
@@ -1526,10 +1546,12 @@ fn dump_pipeline_cache(cache: Arc<PipelineCache>) {
             }
             info!("Pipeline cache dump exit");
         }
+    }
 }
 
 fn get_pipeline_cache(device: Arc<Device>) -> Arc<PipelineCache> {
     let initial_data = {
+        if PIPELINE_CACHING {
             if let Ok(mut file) = File::open(relpath("pipeline_cache.bin")) {
                 let mut data = Vec::new();
                 if file.read_to_end(&mut data).is_ok() {
@@ -1541,6 +1563,9 @@ fn get_pipeline_cache(device: Arc<Device>) -> Arc<PipelineCache> {
             } else {
                 Vec::new()
             }
+        } else {
+            Vec::new()
+        }
     };
 
     unsafe {