3
\$\begingroup\$

Over the past month or so, I’ve been trying to create an extremely fast, platform agnostic, auto-vectorizing sin/cos function for fun. I initially started with sleef-rs’s fast sin function, and decoupled it from the library, made it more extensible, and optimized it further.

During my optimization process, I’ve been closely watching the generated assembly and llvm-mca on multiple platforms to make sure progressions on some platforms didn’t result in regressions on others. Right now, this is faster than SLEEF and Intel SVML, however, those are typically used for more accurate results.

Here’s the code currently:

#![feature(core_intrinsics)]
#![no_std]

use core::intrinsics::*;

/// Inputs valid between [-2^23, 2^23].
/// Precision can set between 0 and 3, with 0 being the fastest and least
/// precise, and 3 being the slowest and most precise.
/// If COS is set to true, the period is offset by PI/2.
///
/// As the inputs get further from 0, the accuracy gets continuously worse
/// due to nature of the fast range reduction.
///
/// This function should auto vectorize under LLVM with opt-level=3.
///
/// The coefficient constants were derived from the constants defined here:
/// https://publik-void.github.io/sin-cos-approximations/#_cos_abs_error_minimized_degree_6
#[inline]
pub unsafe fn sin_fast_approx<const PRECISION: usize, const COS: bool>(x: f32) -> f32 {
    let pi_multiples = fadd_fast(
        fmul_fast(x, core::f32::consts::FRAC_1_PI),
        if COS { 0.0_f32 } else { -0.5_f32 },
    );
    let rounded_multiples = nearbyintf32(pi_multiples);
    let pi_fraction = pi_multiples - rounded_multiples;
    let fraction_squared = pi_fraction * pi_fraction;

    let coeffs = {
        const COEFF_TABLE: [f32; 14] = [
            -4.0_f32,
            0.9719952_f32,
            3.5838444_f32,
            -4.8911867_f32,
            0.99940324_f32,
            -1.2221271_f32,
            4.0412836_f32,
            -4.933938_f32,
            0.9999933_f32,
            0.2196968_f32,
            -1.3318802_f32,
            4.058412_f32,
            -4.934793_f32,
            0.99999994_f32,
        ];

        let shifted_degree = PRECISION + 1;
        let slice_start = (((shifted_degree * shifted_degree) + shifted_degree) / 2) - 1;
        let slice_end = slice_start + PRECISION + 2;
        &COEFF_TABLE[slice_start..slice_end]
    };

    let mut output = coeffs[0];
    for i in 1..coeffs.len() {
        output = fadd_fast(fmul_fast(fraction_squared, output), coeffs[i]);
    }

    let parity_sign = (rounded_multiples.to_int_unchecked::<i32>() as u32) << 31_u32;
    f32::from_bits(output.to_bits() ^ parity_sign)
}

Code and built-in benchmarks on Compiler Explorer (Keep in mind that the hardware used to run the benchmarks is not consistent between tests. I would advise looking at llvm-mca to check the current architecture to interpret results.)

The main thing I’ve been watching is the separate truncation and rounding, attempting to combine them if possible. I’ve tried adding and subtracting 2^23, shifting the parity bit in between, but that was only faster on intel processors. I tried bit-shifting things around to combine parts of the operations, but that was slower across the board.

What I was able to do was, using platform-specific instructions, merge the round and int conversion into a single vcvtps2pd, and converting it back with a vcvtpd2ps, but I can’t seem to get LLVM to generate it for any situations. If I were able to get this to generate, it would be a pretty large performance gain for Intel CPUs that have a vroundps instruction with a 6-8 cycle latency.

Are there any ways to make this faster, more accurate (without sacrificing performance), or cleaner?

\$\endgroup\$
2
  • \$\begingroup\$ Hey, welcome to CodeReview, your question is totally valid as it is, but I would suggest adding an automatic benchmark and test so that reviewers can test their ideas more efficiently. \$\endgroup\$
    – Caridorc
    Mar 23, 2023 at 2:23
  • \$\begingroup\$ @Caridorc Just added some benchmarks for performance and error that I forgot to put into the Compiler Explorer link initially. \$\endgroup\$
    – burgerdude
    Mar 23, 2023 at 3:17

1 Answer 1

1
\$\begingroup\$

I’ll defer to some of the greats and do this old-school. The Apollo 11 Guidance Computer approximation of sine (whose authors included Margaret Hamilton), calculated a three-term polynomial, .7853134·x - .3216147·x³ + 0.036551·x⁵ (close but not identical to a Taylor series approximation), This ran on a computer less powerful than your USB port, and was accurate enough to get a spaceship to the moon.

But that was well before my time. And at least she had a FPU! In the ’80s, it was commen to keep a table of one-eighth the unit circle, and rotate or reflect a lookup into that table to the actual sin or cos value requested using subrraction and/or reversing the sign.

For kicks, I went ahead and implemented a solution like that, which keeps a const table of the sin values in the first octant, maps all other values to the first octant using trigonometric identities, and rounds to the nearest discrete entry in the table.

I doubt this will actually be faster, since it has unpredictable branches, and a larger table with acceptable accuracy would incur a lot of cache misses. Also, Rust does not support floating-point math in const or static computations, so you would either need to generate a very large array expression in the source, or use lazy_static. Like I said, for kicks.

In theory, though, since the 8:1 mapping takes out the sine bit and two bits of precision, you could cover any f32 to the limit of its 24 bits of mantissa with “only” 4 Mi entries, using “only” 16 MiB of memory. Wonder what Margaret Hamilton would say about that.

And what programmers on computers without FPUs actually did back then was compute a fixed-precision table that used only integer math. If you could do that, you could actually turn this into aconst fn.

use std::f64::consts::TAU;

const SIN_TABLE_ENTRIES : usize = 32;
const SIN_OCTANT : [f64; SIN_TABLE_ENTRIES + 1] = [ 
    0.0,  0.024541228522912288, 0.049067674327418015, 0.07356456359966743,
    0.0980171403295606, 0.1224106751992162, 0.14673047445536175, 0.17096188876030122,
    0.19509032201612825, 0.2191012401568698, 0.24298017990326387, 0.26671275747489837,
    0.29028467725446233, 0.3136817403988915, 0.33688985339222005, 0.3598950365349881,
    0.3826834323650898, 0.40524131400498986, 0.4275550934302821, 0.44961132965460654,
    0.47139673682599764, 0.49289819222978404, 0.5141027441932217, 0.5349976198870972,
    0.5555702330196022, 0.5758081914178453, 0.5956993044924334, 0.6152315905806268,
    0.6343932841636455, 0.6531728429537768, 0.6715589548470183, 0.6895405447370668,
    0.0 // THis entry is a workaround for calculating a modulus of SIN_TABLE_ENTRIES instead of 0,
 ];


pub fn sin_approx(theta : f64) -> f64 {
    let approx_gradient = (theta*8.0*SIN_TABLE_ENTRIES as f64/TAU).round() as i64 % (SIN_TABLE_ENTRIES*8) as i64;
    let normalized = if approx_gradient < 0 {
            (approx_gradient + 8*SIN_TABLE_ENTRIES as i64) as usize 
        } else {
            approx_gradient as usize
        };

    if normalized < SIN_TABLE_ENTRIES {
        SIN_OCTANT[normalized]
    } else if normalized < 2*SIN_TABLE_ENTRIES {
        let complement = SIN_OCTANT[2*SIN_TABLE_ENTRIES - normalized];
        (1.0 - complement*complement).sqrt()
    } else if normalized < 3*SIN_TABLE_ENTRIES {
        let complement = SIN_OCTANT[normalized-2*SIN_TABLE_ENTRIES];
        (1.0 - complement*complement).sqrt()
    } else if normalized < 4*SIN_TABLE_ENTRIES {
        SIN_OCTANT[4*SIN_TABLE_ENTRIES - normalized]
    } else if normalized < 5*SIN_TABLE_ENTRIES {
        -SIN_OCTANT[normalized - 4*SIN_TABLE_ENTRIES]
    } else if normalized < 6*SIN_TABLE_ENTRIES {
        let complement = SIN_OCTANT[6*SIN_TABLE_ENTRIES-normalized];
        -(1.0 - complement*complement).sqrt()
    } else if normalized < 7*SIN_TABLE_ENTRIES {
        let complement = SIN_OCTANT[normalized-6*SIN_TABLE_ENTRIES];
        -(1.0 - complement*complement).sqrt()
    } else if normalized < 8*SIN_TABLE_ENTRIES {
        -SIN_OCTANT[8*SIN_TABLE_ENTRIES-normalized]
    } else {
        panic!("The normalized index was not normalized!")
    }
}

And a quick test:

pub fn main() {
    (-24 as i32..=24).map(move|n|{n as f64 * TAU/12.0})
                     .map(move|x|{ ( x, x.sin(), sin_approx(x) ) })
                     .for_each(move|(theta, x1, x2)| {
                          println!("sin {} ≈ {} ≈ {}", theta, x1, x2); } )
}

Update

This answer was more a different solution than an actual code review, so which I hope you don’t think is typical of this site. I’d wanted to investigate a few things about your answer, did not originally have time to finish, and what I found was, LGTM, ship it.

The one thing I might recommend is that it would be cleaner, in my opinion, if you wrote the calculation of the coefficent arrays as something like a match PRECISION block ending in _ => !unreachable(). That makes a bug where someone sets PRECISION to 5 fail fast, with a clear reason in the code why it failed. However, there is already a comment not far above explaining what the allowed values of PRECISION are. It would also be good to make the coefficients const, and const arrays of floating-point numbers are possible, just not const expressions that use floating-point operations. but as you know, Rust stops you from using the outer function’s const attribute.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ I'll accept this as the answer because I took your idea and decided to make the coefficient table use a match expression, which is much cleaner and includes an unreachable!() for invalid values. However, it would be nice to squeeze a bit more performance out of this function if possible. \$\endgroup\$
    – burgerdude
    Apr 2, 2023 at 23:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.