Hypot op

[alexneufeld]

I noticed that $\sqrt(x^2+y^2)$ is a common idiom in SDFs.

The wikipedia page on Pythagorean addition has this fun quote:

Donald Knuth has written that "Most of the square root operations in computer programs could probably be avoided if [Pythagorean addition] were more widely available, because people seem to want square roots primarily when they are computing distances."

In theory, adding a Hypot operator should make some instruction tapes shorter and therefore faster to compute. Example trees for the dreaded box sdf:

before:

after:

So I made a prototype branch with just enough code to test a JIT float-slice evaluator

    fn build_hypot(&mut self, out_reg: u8, lhs_reg: u8, rhs_reg: u8) {
        dynasm!(self.0.ops
            ; vmulps ymm0, Ry(reg(lhs_reg)), Ry(reg(lhs_reg))
            ; vmulps ymm1, Ry(reg(rhs_reg)), Ry(reg(rhs_reg))
            // unlike a real hypot implementation, this intermediate add can cause +inf results
            ; vaddps ymm2, ymm0, ymm1
            ; vsqrtps Ry(reg(out_reg)), ymm2
        );
    }

The results are not promising, even for a favorable benchmark of nothing but hypotenuse calculations:

 sqrt: eval'd 10000000 points in 3.0443s
 hypot: eval'd 10000000 points in 2.9726s

benchmark code

use  fidget::{context::Tree, jit::JitShape};
use fidget::shape::EzShape;
use rand::RngExt;
use std::time;


fn rec_hypot(depth: usize, maxdepth: usize) -> Tree {
    if depth == maxdepth {
        let mut rng = rand::rng();
        let ax1i: usize = rng.random_range(0..=2);
        let ax2i: usize = rng.random_range(0..=2);
        let ax1 = match ax1i {0 => Tree::x(), 1 => Tree::y(), 2 => Tree::z(), _ => unreachable!()};
        let ax2 = match ax2i {0 => Tree::x(), 1 => Tree::y(), 2 => Tree::z(), _ => unreachable!()};
        let ax1c: f32 = rng.random_range(-10.0..10.0);
        let ax2c: f32 = rng.random_range(-10.0..10.0);
        return (ax1 - ax1c).hypot(ax2 - ax2c);
    } else {
        let rhs = rec_hypot(depth + 1, maxdepth);
        let lhs = rec_hypot(depth + 1, maxdepth);
        return lhs.hypot(rhs);
    }
}

fn rec_sqrt(depth: usize, maxdepth: usize) -> Tree {
    if depth == maxdepth {
        let mut rng = rand::rng();
        let ax1i: usize = rng.random_range(0..=2);
        let ax2i: usize = rng.random_range(0..=2);
        let ax1 = match ax1i {0 => Tree::x(), 1 => Tree::y(), 2 => Tree::z(), _ => unreachable!()};
        let ax2 = match ax2i {0 => Tree::x(), 1 => Tree::y(), 2 => Tree::z(), _ => unreachable!()};
        let ax1c: f32 = rng.random_range(-10.0..10.0);
        let ax2c: f32 = rng.random_range(-10.0..10.0);
        return ((ax1 - ax1c).square() + (ax2 - ax2c).square()).sqrt();
    } else {
        let rhs = rec_sqrt(depth + 1, maxdepth);
        let lhs = rec_sqrt(depth + 1, maxdepth);
        return (lhs.square()+rhs.square()).sqrt();
    }
}

fn main() {
    let maxdepth=10;
    let h = rec_hypot(0, maxdepth);
    let old = rec_sqrt(0, maxdepth);
    let mut rng = rand::rng();
    let n = 10_000_000;
    let xvals: Vec<f32> = (0..n).map(|_| rng.random_range(-10.0..10.0)).collect();
    let yvals: Vec<f32> = (0..n).map(|_| rng.random_range(-10.0..10.0)).collect();
    let zvals: Vec<f32> = (0..n).map(|_| rng.random_range(-10.0..10.0)).collect();

    let shape = JitShape::from(old);
    let mut eval = JitShape::new_float_slice_eval();
    let tape = shape.ez_float_slice_tape();
    let stopwatch = time::Instant::now();
    let _out2 = eval.eval(&tape, &xvals, &yvals, &zvals).unwrap();
    println!("sqrt: eval'd {} points in {:.4}s", n, stopwatch.elapsed().as_secs_f32());

    let shape = JitShape::from(h);
    let mut eval = JitShape::new_float_slice_eval();
    let tape = shape.ez_float_slice_tape();
    let stopwatch = time::Instant::now();
    let _out1 = eval.eval(&tape, &xvals, &yvals, &zvals).unwrap();
    println!("hypot: eval'd {} points in {:.4}s", n, stopwatch.elapsed().as_secs_f32());

    // example output:
    // sqrt: eval'd 10000000 points in 3.0443s
    // hypot: eval'd 10000000 points in 2.9726s

}

This was a just-for-fun experiment, but maybe also a useful datapoint. I had some fun poking at fidgets internals at least :)

[mkeeter]

Thanks for the report!

I had done similar experiments, and came to a similar conclusion – it seems like it should help, but I didn't see any real change in practice. The JIT performance doesn't specifically surprise me, though: I think that the JIT would generate nigh-identical machine code for the manual hypot operator, since the build_square / build_add / build_sqrt builders emit the same vmulps / vaddps / vsqrtps instructions.

(my implementation was at the tape level, and I just pushed the op-fusion branch in case you're interested)

It may still be worthwhile for the GPU implementation, where tape space is at more of a premium...