fastsin

Let's build the fastest accurate sinf / cosf in the world.

Not the fastest — CUDA already has __sinf, one hardware instruction, and nothing will beat it. Not the most accurate either — that is what sinf is for. The interesting target is the corner nobody occupies: sinf's accuracy at a speed close to the intrinsic's, across the whole float range.

Today, on six GPUs across four architectures, fastsin::sin is 6–8× faster than sinf at |x| = 10⁶ with identical error. That is a starting point, not a finish line. Contributions that push either axis are what this repo is for.

#include "fastsin.cuh"

__global__ void k(const float* x, float* y, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) y[i] = fastsin::sin(x[i]);   // drop-in for sinf
}

Header-only, MIT, C++17, sm_75 and up. Copy include/ and include the header.

# -arch must match your GPU: it selects the tuning constants, not just the SASS
ARCH=sm_$(nvidia-smi --query-gpu=compute_cap --format=csv,noheader | head -1 | tr -d '.')
nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude your_code.cu

# check what your build actually did
nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude -o config bench/config.cu && ./config

	valid |x|	abs error
sin cos sincos	all	~1.1e−07	start here — dispatches on |x|
sin_warp cos_warp sincos_warp	all	~1.1e−07	warp-uniform; needs a full converged warp
sin_small cos_small sincos_small	< 5e6	~9e−08	fastest accurate path. Silently wrong past 5e6
sin_big cos_big sincos_big	all	~1.1e−07	branch-free, no dispatch
sin_rough cos_rough sincos_rough	< 1e3	~4e−07	__sinf-class accuracy, 1000× wider range

-use_fast_math silently rewrites every sinf in your code into __sinf — error goes from 8.7e−08 to 0.14 at |x| = 10⁶ with no diagnostic. fastsin is unaffected, so it also works as a way to keep accurate trig in a build that needs fast-math elsewhere. bench/config.cu detects the flag.

---

References

The algorithms are not new. These are the sources they come from.

Argument reduction — used directly

• Payne, M. H. and Hanek, R. N. Radian Reduction for Trigonometric Functions. ACM SIGNUM Newsletter 18(1), 1983, 19–24. — the infinite-precision reduction sin_big implements.
• Cody, W. J. and Waite, W. Software Manual for the Elementary Functions. Prentice-Hall, 1980. — the splitting of π/2 into few-significant-bit pieces that sin_small uses.
• Ng, K. C. Argument Reduction for Huge Arguments: Good to the Last Bit. SunPro technical report, 1992. — the practical treatment of how many bits of 2/π you actually need, which bench/exhaustive.cu re-derives by brute force for binary32.

Background and error analysis

• Muller, J.-M. Elementary Functions: Algorithms and Implementation. Birkhäuser, 3rd ed., 2016. — the standard reference for everything here.
• Brisebarre, N., Defour, D., Kornerup, P., Muller, J.-M. and Revol, N. A New Range-Reduction Algorithm. IEEE Transactions on Computers 54(3), 2005, 331–339.
• Boldo, S., Daumas, M. and Li, R.-C. Formally Verified Argument Reduction with a Fused Multiply-Add. IEEE Transactions on Computers 58(8), 2009, 1139–1145.
• Daumas, M., Mazenc, C., Merrheim, X. and Muller, J.-M. Modular Range Reduction: A New Algorithm for Fast and Accurate Computation of the Elementary Functions. Journal of Universal Computer Science 1(3), 1995.
• Goldberg, D. What Every Computer Scientist Should Know About Floating-Point Arithmetic. ACM Computing Surveys 23(1), 1991, 5–48.

Vendor documentation

• CUDA C++ Programming Guide, Appendix: Mathematical Functions — the documented 2 ULP bound for sinf that everything here is measured against.

---

Contributing

The ask is simple: show me a measurement, and I will merge it.

A contribution here is not "I think this is faster." It is: at this range of |x|, this function beats the current one by this much on speed or accuracy, and here is the run that proves it. Bring that and it goes in.

What a proposal looks like

1. Which entry point, or a new one.
2. Which range of |x| it wins in. Everything in this library is a trade that holds over some interval and fails outside it; saying where is not a caveat, it is the result.
3. How much, on which axis — throughput, absolute error, or both.
4. The output of the harness, not a summary of it.

The harness

ARCH=sm_$(nvidia-smi --query-gpu=compute_cap --format=csv,noheader | head -1 | tr -d '.')

# throughput and error across twenty decades of |x| -- this is the evidence
nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude -o sweep bench/sweep.cu
./sweep > results/sweep/YourGPU.csv
python bench/plot.py                    # regenerates docs/curves.svg

# correctness gates -- both must pass, no exceptions
nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude -o exhaustive bench/exhaustive.cu
./exhaustive                            # all 2^31 positive floats vs double

nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude -o test_all bench/test_all.cu
./test_all                              # 96M samples + inf, NaN, subnormals

# optional: the same inputs judged in ULP instead, which is the metric this
# library deliberately does not hold. Diagnostic, not a gate.
nvcc -O3 -arch=$ARCH -std=c++17 -Iinclude -o ulp bench/ulp.cu && ./ulp

The bar

• exhaustive reports PASS and test_all reports SPECIAL_SUMMARY mismatches=0 verdict=PASS.
• No existing entry point gets slower or less accurate. If a change trades one for the other, it becomes a new entry point rather than replacing one.
• A tuning constant that differs per architecture goes in include/fastsin_tuning.cuh with the measurement in a comment. That file has no unsourced numbers in it and should stay that way.
• Negative results are welcome as issues. Several things that obviously should have worked did not, and knowing which is worth as much as the code.

Where the room is

Concretely, what I would most like to see:

• Blackwell. sm_100 / sm_120 compile and are correct, but fall through to conservative defaults. Someone with the hardware running bench/frontier.cu and adding a measured branch would close the last gap in the architecture table.
• sin_big past 154 Gop/s. It is still 1.7× behind sin_small on an RTX 4060. bench/ablate.cu prices each remaining piece.
• A cheaper dispatcher. nvcc if-converts the |x| branch and runs both paths; __noinline__ measured worse. There may be a third way.
• Other precisions. __half, __nv_bfloat16, and double all have the same shape of problem and none of them are here.
• Other vendors. HIP and SYCL ports.

---

Where it stands

RTX 4060. Top: throughput. Bottom: absolute error. The amber field is |x| ≥ 105615, where CUDA's sinf leaves its fast path — a constant measured to one ULP and identical on every GPU tested. The grey field is worse than 2 ULP, the bound CUDA documents for sinf.

At |x| = 10⁶, same error in every row (8.45e−08):

GPU	arch	sinf	fastsin::sin
Tesla T4	sm_75	22.9	138.1	6.02×
A100-40GB	sm_80	65.2	441.7	6.77×
A10	sm_86	46.9	335.6	7.16×
L4	sm_89	43.9	295.7	6.73×
RTX 4060	sm_89	27.4	219.1	8.01×
H100	sm_90	112.6	882.5	7.84×

Gop/s, one evaluation per op. This is not a speed-for-accuracy trade — the error columns are identical to three significant figures.

Verified exhaustively, not sampled. bench/exhaustive.cu evaluates every one of the 2³¹ positive floats against a double-precision reference; negatives mirror exactly, since sin is odd and cos is even.

SCAN sin_big   max_abs_err=1.1042e-07   PASS
SCAN cos_big   max_abs_err=1.0948e-07   PASS

SPECIAL_SUMMARY mismatches=0 verdict=PASS on all six GPUs above.

Not done yet

• Blackwell is untuned.

• Accuracy is absolute, not ULP — and the gap is not small. Measured over 96M inputs by bench/ulp.cu:

  	max abs error	max ULP
  sinf	8.81e−08	1.49
  sin_big	1.09e−07	3,079,907

  Near a zero of sin the correctly rounded result is tiny, so its ULP is tiny, and a reduction that is comfortably within 1e−07 absolute is millions of ULP out in relative terms. sinf holds its documented bound there and this does not. If you need relative accuracy near multiples of π, use sinf.

• sin_small returns a wrong number past |x| = 5e6 with no warning. Use fastsin::sin unless you know your range.

• Only NVIDIA, only float.

---

Support

If this saved you time, buy me a coffee.

Every hour of GPU rental for the cross-architecture sweeps comes out of pocket — the numbers in the table above are six cloud GPUs, re-run from scratch on every change to the reduction. Coffee keeps the sweeps running and the next architecture measured instead of guessed.

License

MIT. See LICENSE.