# Should Julia have rsqrt? And e.g. support RSQRTSS

**URL:** <https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969>\
**Category:** Internals & Design\
**Tags:** numerics\
**Created:** [February 16, 2025, 6:56pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969 "2025-02-16T18:56:42Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [February 16, 2025, 6:56pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/1 "2025-02-16T18:56:42Z")

</div>

Neither Python (Numpy) nor C++ have rsqrt, but pytorch and tensorflow have it (I didn’t confirm in Flux.jl or Lux.jl or look carefully, though CUDA.jl/Nvidia and JavaScript stdlib have for single and double, and AMD/HIP has too), and CPUs have instructions for. Zig has this as no planned (likely since they optimize better already, see there for Gotbolt): [rsqrt as builtin · Issue #19302 · ziglang/zig · GitHub](https://github.com/ziglang/zig/issues/19302) and Rust has an open issue: [Add rsqrt method to Float trait · Issue #1 · rust-num/num-traits · GitHub](https://github.com/rust-num/num-traits/issues/1)

1/sqrt(x) is common enough, that you want it to be fast, by default, and even to use fastest approximate result. It’s probably the only CPU (FPU) instruction I can think of not well-supported in Julia.

At first looking at if it’s supported, it seemingly wasn’t in Julia, i.e. code\_native was verbose, though is seems supported, maybe in a good ~~best~~ possible way, for full accuracy (and bit-identical results across Intel and AMD), but not for speed.

Some background:

```julia-auto
julia> sqrt(2)
1.4142135623730951

```

seems right, the only correct result (except for infinite-length… irrational), but I can argue for 1.4 accurate enough, at least no longer than 1.4142135f0 casted back to Float64 if wanted. It depends on the accuracy of the number 2, if one significant digit, it’s [1.5, 2.5]. sqrt(2.0) is still only 1.4, or even 1.396, sqrt(2.00) is 1.41.

So I actually think `sqrt` of Float64, should return Float32, but that ship has sailed (at least for now), but for `rsqrt` the implied `sqrt` there can return Float32, while the final seemingly needs to return Float64, because of values in [0, 1], as done in benchmark code below, actually only values in [0, 1f-46] otherwise overflow).

```julia-auto
julia> Float32(Float64(1/sqrt(1f-46)))
Inf32

```

> <https://stackoverflow.com/questions/15175654/sqrt-vs-rsqrt-vs-sse-mm-rsqrt-ps-benchmark>

> **My conclusion is not it is not worth bothering with SSE2 unless we make calculations on no less than 4 variables. (Maybe this applies to only rsqrt here but it is an expensive calculation (it also includes multiple multiplications) so it probably applies to other calculations too)**
> 
> **Also `sqrt(x)` is faster than `x*rsqrt(x)` with two iterations, and x\*rsqrt(x) with one iteration is too inaccurate for distance calculation.**

> <https://stackoverflow.com/questions/58614226/is-there-a-c-function-that-returns-exactly-the-value-of-the-built-in-cpu-opera>

> There is no function in the standard library that does this, but your compiler might optimize the expression 1 / sqrt(value) such that it does emit the RSQRTSS instruction.

[It doesn’t for Julia, even with fastmath.]

> **[rr Trace Portability: Diverging Behavior of RSQRTSS in AMD vs Intel](https://robert.ocallahan.org/2021/09/rr-trace-portability-diverging-behavior.html)**

> Unfortunately, today I discovered a new difference between AMD and Intel: the [RSQRTSS](https://www.felixcloutier.com/x86/rsqrtss) instruction. Perhaps this is unsurprising, since it is described as: “computes an _approximate_ reciprocal of the square root of the low single-precision floating-point value in the source operand” (emphasis mine). A simple test program:

```julia-auto
..
asm ("rsqrtss %1,%0" : "=x"(out) : "x"(in));
..

On Intel Skylake I get

out = 3d7ff000, float = 0.062485

On AMD Rome I get

out = 3d7ff800, float = 0.062492

[Neither values are correct, or within [prevfloat, nextfloat] of the correctly computed 0.0625, unless 256 has only 3 or 4 significant decimal digits, i.e. going down to a range covering 256.95f0 is enough.]

Intel's result just stays within the documented 1.5 x 2-12 relative error bound. (Seems unfortunate given that the exact reciprocal square root of 256 is so easily computed to 0.0625, but whatever...) 

```

https://embed.reddit.com/r/programming/comments/12tys7d/revisiting_the_fast_inverse_square_root_is_it/?embed=true&ref_source=embed&ref=share

@foobar_lv, @ChrisRackauckas

> [@Help speeding up a test progam from a book](https://discourse.julialang.org/t/help-speeding-up-a-test-progam-from-a-book/6383/16):
>
> What you really want to do is emit the juicy RSQRTSS and RCPSS instructions (on x86) and their vectorized counterparts. In terms of speedup, single precision is chump change compared to these beautiful bastards (compute approximate reciprocal square-root and approx reciprocal of floats for almost the same price as multipliation, as of [http://www.agner.org/optimize/instruction\_tables.pdf](http://www.agner.org/optimize/instruction_tables.pdf)). Which brings me to the question: Is there a sensible way of accessing these instructions in Julia? Like @r…

> **[GitHub - singularitti/FastInverseSqrt.jl: A pure Julia implementation of the fast inverse...](https://github.com/singularitti/FastInverseSqrt.jl)**
>
> A pure Julia implementation of the fast inverse square root algorithm

> **[Language Marketing for Julia compared to Fortran](https://fortran-lang.discourse.group/t/language-marketing-for-julia-compared-to-fortran/3309/65)**
>
> @mecej4 yes, exactly. That’s why I think this “numerical benchmark” is not a good benchmark. So the only other numerical benchmark in there is the n-body benchmark, which I think actually is not bad. It’s a bit specialized (mostly how quickly you...

```julia-auto
    # When called with @fastmath, 1 / sqrt(x::Float32) will use
    # SSE single-precision rsqrt approximation followed by a single
    # iteration of the Newton-Raphson method. This, followed by an
    # additional double-precision Newton-Raphson iteration gives
    # sufficient precision for this problem and is significantly
    # faster than double-precision division and sqrt.
    rd = @ntuple 10 k-> @fastmath Float64(1 / sqrt(Float32(dsq[k])))

```

```julia-auto
julia> @code_lowered rsqrt(10.0f0)
CodeInfo(
1 ─ %1 = Base.FastMath
│ %2 = Base.getproperty(%1, :div_fast)
│ %3 = Base.FastMath
│ %4 = Base.getproperty(%3, :sqrt_fast)
│ %5 = Main.Float32(Main.x)
│ %6 = (%4)(%5)
│ %7 = (%2)(1, %6)
│ %8 = Main.Float64(%7)
└── return %8
)

```

One use of `rsqrt` if in Diff Attension for DIFF Transformer (seem my Julia code translated from the paper there):

> [@Community Interest Check: LLMs from Scratch in Pure Julia](https://discourse.julialang.org/t/community-interest-check-llms-from-scratch-in-pure-julia/121796/8):
>
> Consider implementing DIFF Transformer from October paper (it seems better across the board), I started translating to Julia: function DiffAttn(X, W\_q, W\_k, W\_v, λ) Q1, Q2 = split(X \* W\_q) K1, K2 = split(X \* W\_k) V = X \* W\_v # Qi, Ki: [b, n, d]; V: [b, n, 2d] s = 1 / sqrt(d) # torch.rsqrt A1 = Q1 \* K1.transpose(−1, −2) \* s A2 = Q2 \* K2.transpose(−1, −2) \* s return (softmax(A1) − λ \* softmax(A2)) \* V end I leave in the Python pseudocode as is: def MultiHead(X, W\_q, W\_k, W\_v, W…

[Offtopic, but also curious, does Julia ever emit the “[FPATAN — Partial Arctangent](https://www.felixcloutier.com/x86/fpatan)” instruction? Seemingly not. Couldn’t confirm maybe behind the `call rax` I only see in the assembly. Is there a good way to know and/or get full inlined assembly?  
[x86 FPATAN instruction · Issue #19330 · ziglang/zig · GitHub](https://github.com/ziglang/zig/issues/19330)]

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [February 16, 2025, 10:07pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/2 "2025-02-16T22:07:07Z")

</div>

Already tracked among the issues on Github:

- [IEEE-754 recommended functions · Issue #6148 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/6148)
- [We need to provide these IEEE754-2019 functions · Issue #32877 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/32877)

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [February 17, 2025, 1:29am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/3 "2025-02-17T01:29:35Z")

</div>

> [@Palli](#):
>
> but pytorch and tensorflow have it

That’s likely coming from Stable HLO: [StableHLO Specification &nbsp;|&nbsp; OpenXLA Project](https://openxla.org/stablehlo/spec#rsqrt). But it doesn’t look like LLVM has a dedicated intrinsic: [LLVM Language Reference Manual — LLVM 21.0.0git documentation](https://llvm.org/docs/LangRef.html)

---

<div class="post-metadata">

**Author:** ![wsmoses](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wsmoses/32/26497_2.png) [@wsmoses](https://discourse.julialang.org/u/wsmoses)\
**Post date:** [February 17, 2025, 1:31am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/4 "2025-02-17T01:31:08Z")

</div>

Fun fact Lux.jl can use that automatically via Reactant.jl (which will emit to use rsqrt where appropriate)

cc @avik-pal

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [February 17, 2025, 12:29pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/5 "2025-02-17T12:29:20Z")

</div>

It’s good to see on the TODO list there for IEEE-754 recommended functions.

But then a more controversial suggestion. What should be the return type? I suggest Float32 for Float64 input, and any input really, e.g. for Float16 too, except BigFloat, and I suppose bfloat should return the input type then bfloat.

This is the fastest option, can access the fastest (single-precision) assembly instruction.

Since it’s a function, we could have rsqrt(x; accurate=false) to call for same return type, with the default overridden. I don’t know if there’s president for return type to depend on the keyword argument not input type. I don’t think IEEE specifies that Float64 must return Float64.

There is precedent to NOT just rely on input type, e.g. 1/3 is Float64 not an Int or a Rational.

> Should all of these go into `Base`?

It’s technically breaking if exported. But I rather want to have this in Base, for discoverability and ease of use to access a single assembly instruction… and doing `Base.rsqrt` is awkward, so a POC would be nice and PkgEval ASAP.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [February 17, 2025, 1:28pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/6 "2025-02-17T13:28:16Z")

</div>

> [@Palli](#):
>
> I don’t think IEEE specifies that Float64 must return Float64.

It depends on the operation. Of the many terms Clause 5 of IEEE 754 uses, one is _formatOf_, which “specifies the floating-point destination _format_, which might be different from the floating-point operands’ formats”. Clause 5 specifies _formatOf_- **squareRoot** (_source1_), and while subclause 9.2 does not use these terms for operations like **rSqrt** , we can infer the same thing from its definition 1/√(x).

However in practice, people do expect matching input and output numeric types when feasible. This is true of Julia’s `sqrt`:

```julia
julia> typeof.(sqrt.(one.([Float16 Float32 Float64 BigFloat])))
1×4 Matrix{DataType}:
 Float16 Float32 Float64 BigFloat

```

so `rsqrt` must match to conform to its IEEE definition.

> [@Palli](#):
>
> This is the fastest option, can access the fastest (single-precision) assembly instruction.

If this is your intent, we can simply support `rsqrt(::Float32)` without violating the expectation of matching types, and it’s simple and typical to manually convert scalars or collections to `Float32` for performance. Consider the possibility that people need `rsqrt` to support another format like double precision…

> [@Palli](#):
>
> Float32 for Float64 input, and any input really, e.g. for Float16…

then we’d need a breaking change or a deprecation in favor of a new function.

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [February 17, 2025, 3:15pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/7 "2025-02-17T15:15:45Z")

</div>

> [@Palli](#):
>
> I don’t know if there’s president for return type to depend on the keyword argument not input type.

That’d make the function type-unstable, since keyword arguments (let alone their _values_ rather than their types) don’t participate in dispatch.

> [@Palli](#):
>
> There is precedent to NOT just rely on input type, e.g. 1/3 is Float64 not an Int or a Rational.

That comparison doesn’t make much sense to me. Contrary to other arithmetic operations (apart from un-/signedness and overflow details), (non-integer) division isn’t closed over integer numbers, unless you want to throw an error for `1/3` the only practical option is to always return a floating point number for integer input. But here you’re suggesting not respecting the _ **precision** _ of the input, which is a big gotcha.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 17, 2025, 11:32pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/8 "2025-02-17T23:32:46Z")

</div>

Furthermore, `/(::Int, ::Int)` is type stable, so the output type does _indeed_ rely just on the input types.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [February 18, 2025, 1:37pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/9 "2025-02-18T13:37:37Z")

</div>

Sometimes we want most accurate, or appropriate, same type back, and sometimes we just want access to the fastest method.

Most of the time not all of the 53 bits in Float64 are significant, and we just use Float64 out of habit. If only slightly less than half of them are then it makes sense to do calculation/output in Float32. Or top of that, the RSQRTSS only gives “about 11 significant bits of precision”, and it depends on Intel or AMD.

I find it intriguing `rsqrt(n)` (with its assembly instruction) is actually the faster way to compute `sqrt(n)` by then multiplying by n, and the combination is faster, another reasons to have rsqrt, and as fast as possible.

> **[How do computers calculate square roots?](https://www.quora.com/How-do-programming-languages-and-calculators-perform-square-root-calculations/answer/Andrew-Bromage?no_redirect=1)**
>
> Answer (1 of 21): The Babylonian Method.
> 
> Take the number we are trying to find the square root of, a. The computer plugs this number into the formula y=0.5(x+a/x). Where x is any positive number. Then y becomes the new value for x.
> 
> (It is important...

> **Newton-Raphson, take two**
> 
> Now comes the first key insight: Instead of computing the square root, compute the reciprocal of the square root. […] It turns out that this is a far easier number to compute, and if you need the square root, multiply this number by n and you are done.

> **[In computer programming, why do you avoid the square root?](https://www.quora.com/In-computer-programming-why-do-you-avoid-the-square-root/answer/Andrew-Bromage)**
>
> Andrew Bromage's answer: If you need to calculate a square root, you should not avoid the square root operation. If you are unsure whether or not you need to calculate a square root, and want to know what the alternatives are, then you need to...

> I’m going to use [Agner Fog’s instruction tables](http://www.agner.org/optimize/) as my source […] So here are the relative costs of some basic floating-point operations. Note that I am only going to deal with non-vector operations for simplicity.
> 
> - […] Double precision floating point add (ADDSD/SUBSD): latency 3, throughput 1
> - Single precision floating point multiply (MULSS): latency 5, throughput 0.5
> - Double precision floating point multiply (MULSD): latency 5, throughput 0.5
> - Single precision floating point divide (DIVSS): latency 10–13, throughput 7
> - Double precision floating point divide (DIVSD): latency 10–20, throughput 8–14
> - […] Single precision floating point square root (SQRTSS): latency 11, throughput 7
> - Double precision floating point square root (SQRTSD): latency 16, throughput 8–14
> - Single precision approximate reciprocal square root (about 11 significant bits of precision, RSQRTSS): latency 3, throughput 1

So `sqrt` can be calculated in 3+5=8 cycles vs 11 (Float32) or 16 (Float64) (yes, throughput is different, and accuracy is less).

It felt like we needed only half as many digits in the result, since **i** sqrt gives always half as many binary digits, but it actually seems I was wrong on only half as many significant digits:

> **[StandardsofAccuracy.pdf](https://www.soest.hawaii.edu/martel/Courses/GG303/StandardsofAccuracy.pdf)**
>
> 829.23 KB

> The number of significant figures in a square root is equal to the number of significant figures in its square

[https://proofwiki.org/wiki/Number\_of\_Significant\_Figures\_in\_Result\_of\_Square\_Root](https://proofwiki.org/wiki/Number_of_Significant_Figures_in_Result_of_Square_Root)

> <https://math.stackexchange.com/questions/1120419/how-to-determine-significant-figures-involving-radicals-and-exponents>

[https://www.reddit.com/r/learnmath/comments/161qt56/how\_do\_sig\_figs\_apply\_to\_exponents\_and\_square/](https://www.reddit.com/r/learnmath/comments/161qt56/how_do_sig_figs_apply_to_exponents_and_square/)

[https://brainly.in/question/26701421](https://brainly.in/question/26701421)[math - How many Binary Digits needed for X Decimal Digits of Square Root of 2 - Stack Overflow](https://stackoverflow.com/questions/61730548/how-many-binary-digits-needed-for-x-decimal-digits-of-square-root-of-2)[3.18: Significant Figures in Multiplication and Division - Chemistry LibreTexts](https://chem.libretexts.org/Bookshelves/Introductory_Chemistry/Introductory_Chemistry_(CK-12)/03%3A_Measurements/3.18%3A_Significant_Figures_in_Multiplication_and_Division)

> [@Benny](#):
>
> > [@Palli](#):
> >
> > Float32 for Float64 input, and any input really, e.g. for Float16…
> 
> then we’d need a breaking change or a deprecation in favor of a new function.

It wouldn’t be a breaking change for a new (rsqrt) function (I wasn’t proposing for ow, to change others, e.g. sqrt), to not promise more than e.g. 11 significant bits by default. It would be unusual, and maybe it needs to be a non-default. I think we’re back to suggesting that, same output type, fully rounded, and non-default for lesser accuracy by keyword argument (still returned in same type likely, I’m though unconvinced the function needs to be type-stable).

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [February 18, 2025, 2:45pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/10 "2025-02-18T14:45:41Z")

</div>

> [@Palli](#):
>
> It wouldn’t be a breaking change for a new (rsqrt) function (I wasn’t proposing for ow, to change others, e.g. sqrt), to not promise more than e.g. 11 significant bits by default.

That’s not a breaking change because it’s not a change, just API. However, say some platform comes out with a `rsqrt` instruction for `Float64`. People would naturally expect `rsqrt(::Float64)` to use that instruction and return `Float64`. _That_ would be a breaking change, and it could have been prevented by not using `RSQRTSS` for other input types other than what its name suggests. Alternatively we could make a `rsqrtss` or `rsqrtFloat32` function just for the clarity, though the implicit type conversion still seems unnecessary.

> [@Palli](#):
>
> I’m though unconvinced the function needs to be type-stable

We’d definitely want that; it’s contradictory to introduce type instability overheads if we’re trying to use a hardware instruction for better performance. Besides, you don’t need a keyword argument for a default value, a trailing positional argument can also serve that purpose AND allows multiple methods.

> [@giordano](#):
>
> That’d make the function type-unstable, since keyword arguments (let alone their _values_ rather than their types) don’t participate in dispatch.

Possibly, but the reasoning is incorrect. Keyword arguments not participating in dispatch means _we_ can’t specialize the wider function, that is write multiple methods that vary only the keyword arguments. The _compiler_ still specializes a method over its keyword argument’s types. The issue is how it happens. Currently, keyworded calls are implemented by the keyword arguments going into a `NamedTuple` before being passed to an underlying method with only positional arguments.

- Specializing over `accurate=false` only gives us one `Bool` input type for two possible output types: `Float32` or the input value’s type. That’s type-unstable, whether `accurate` is positional or a keyword makes no difference.
- If we provide a conversion type instead, the `NamedTuple` and thus the underlying method only sees `DataType`. `::Type{T}` annotations of the keyword argument do not reach the underlying method and can’t introduce type stability as it could for positional arguments.
- If we provide a value of the target type instead, say its `zero`, then its type makes it into the `NamedTuple` and the underlying method can be compiled for it.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 19, 2025, 2:06pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/11 "2025-02-19T14:06:37Z")

</div>

> [@Palli](#):
>
> Most of the time not all of the 53 bits in Float64 are significant, and we just use Float64 out of habit.

The lower the precision, the harder one has to think about numerical error. Most people use 64-bit float because a lot of the time (but of course not always) it gives reasonable precision without explicitly thinking about the numerical analysis of an algorithm too hard. It is much easier to run into issues with 32-bit float.

Regarding `RSQRTSS`: if it is ever included in Base, it would be best to highlight that it is an _approximation_, starting with the name. `approximate_rsqrt` or something like that would be best. Initially, a `::Float32` method would be sufficient, which would required that users convert explicitly. It does not really make sense to have a generic method for something so specific to a single type.

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [February 19, 2025, 9:24pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/12 "2025-02-19T21:24:56Z")

</div>

> [@Palli](#):
>
> Most of the time not all of the 53 bits in Float64 are significant, and we just use Float64 out of habit. If only slightly less than half of them are then it makes sense to do calculation/output in Float32. Or top of that, the RSQRTSS only gives “about 11 significant bits of precision”, and it depends on Intel or AMD.

I don’t understand who “we” are in this story. `rsqrt` can be used for many things, like normalizing vectors. When your computation goes on for days and weeks, the use of `Float64` is not a mere habit, it is a necessity to keep the errors at a manageable level.

There are applications where it’s sufficient to have an approximate square root. Some RISC architectures let the user (i.e. the compiler) do the iterations of some Newton-like method (or Babylonian method, or Henon’s method, whatever, it’s all the same), and two or three iterations might be sufficient in some cases. In other applications, the highest possible precision is needed. I’ve seen older fortran programs in biotech doing `REAL*16` computations (emulated) for this purpose.

In case the rsqrt instruction has 32-bit precision, an `rsqrt(::Float64)` should anyway return a `Float64`, but properly documented that it might have lower precision. This is important if/when hardware in the future supports 64-bits in its ISA.

---

<div class="post-metadata">

**Author:** ![barnett567](https://avatars.discourse-cdn.com/v4/letter/b/e0b2c6/32.png) [@barnett567](https://discourse.julialang.org/u/barnett567)\
**Post date:** [February 20, 2025, 3:34am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/13 "2025-02-20T03:34:52Z")

</div>

I’d just like to clarify that, despite some wrong claims in this thread, `sqrt` (and `rsqrt`) are _stable_ operations. Formally: the _relative condition numbers_ \kappa of the evaluation problem is {\cal O}(1). Thus it is meaningful to return `Float64` precision for `Float64` input, and indeed using a couple of Newton iterations as usual for the sqrt should get you to \<1 ULP from the truth. Of course, there is also a place for fast more approximate algorithms. A great place to learn about this and the concept of _backward stability_ is Lectures 12-15 of Trefethen & Bau’s book on Numerical Linear Algebra.

---

<div class="post-metadata">

**Author:** ![Palli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palli/32/3380_2.png) [@Palli](https://discourse.julialang.org/u/Palli)\
**Post date:** [May 24, 2026, 1:54am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/14 "2026-05-24T01:54:40Z")

</div>

> [@sgaure](#):
>
> In case the rsqrt instruction has 32-bit precision, an `rsqrt(::Float64)` should anyway return a `Float64`, but properly documented that it might have lower precision.

Exactly. I’m not so opposed to the quick conversion to the larger type.

> [@barnett567](#):
>
> it is meaningful to return `Float64` precision for `Float64` input, and indeed using a couple of Newton iterations as usual for the sqrt should get you to \<1 ULP from the truth.

I think we can’t do that, i.e. we want access to the (approximate) float/Float64 [RSQRTPS — Compute Reciprocals of Square Roots of Packed Single Precision Floating-PointValues](https://www.felixcloutier.com/x86/rsqrtps)

If we use the instruction for Float64, we can get the double precision precision, but not easily the range (without branchy code, why I think it should be opt-in, for Float64). It seems to me there’s no corresponding Float64 instruction since this is an approximation instruction anyway. It would then be RSQRPD but I don’t find with -D ending.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 24, 2026, 2:10am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/15 "2026-05-24T02:10:03Z")

</div>

> [@Palli](#):
>
> But then a more controversial suggestion. What should be the return type? I suggest Float32 for Float64 input, and any input really, e.g. for Float16 too, except BigFloat, and I suppose bfloat should return the input type then bfloat.

This is not reasonable to me (even for a crude `rsqrt` approximation that only has a few accurate digits) because the exponent ranges are different, so many inputs would overflow to `Inf` or underflow to `0.0`. For example, `rsqrt(1e-80)` would overflow to `Inf32` if the answer were returned in `Float32` precision.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [May 26, 2026, 10:39am UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/16 "2026-05-26T10:39:13Z")

</div>

> [@stevengj](#):
>
> For example, `rsqrt(1e-80)` would overflow to `Inf32` if the answer were returned in `Float32` precision.

But returning `Inf64` is totally fine and OK, isn’t it? (convert input to Float32, RSQRTPS, convert to Float64, return, check that LLVM can optimize out the round-trip if the return-value is immediately converted to Float32)

The docstring should anyways have the form of

> `rsqrt(x)`  
> Compute a rough approximation of inv(sqrt(x)).
> 
> The context is that many architectures have very fast hardware support for this function, but make only limited promises about the actual result and reproducibility on different hardware of the same architecture. This function exists in order to save you the hassle of writing inline-assembly, and give you portability in the sense of reasonable fallbacks existing if there is no hardware support. But in order to know what is actually happening, please consult the `code_native`, your processor manual and appropriate blog posts that reverse engineered the precise behavior on your specific machine.
> 
> Use this instruction only if bitwise reproducibility between different CPUs don’t matter to you, and if your algorithm can deal with the limited precision and dynamic range.

I had to look this up, but you at least get the same result on P and E cores on the same CPU (I think that is a mere observation about existing silicon, and not a promise of a spec).

---

<div class="post-metadata">

**Author:** ![barnett567](https://avatars.discourse-cdn.com/v4/letter/b/e0b2c6/32.png) [@barnett567](https://discourse.julialang.org/u/barnett567)\
**Post date:** [May 26, 2026, 7:00pm UTC](https://discourse.julialang.org/t/should-julia-have-rsqrt-and-e-g-support-rsqrtss/125969/17 "2026-05-26T19:00:05Z")

</div>

Ok, sorry, my comment was about the achievable error in the true sqrt and 1/sqrt operations, not this approximate operation. I now see what the goal of the question was: wrapping the Float32 SIMD instrinsics RSQRTSS, RSQRTPS, etc, and nothing more. As approximate operations these only give 2^{-11} relative accuracy, and map Float32-\>Float32. They are for experts only, since they need to be followed up with a couple of Newton iters to get close to machine accuracy (2^{-23} in Float32, or 2^{-53} if cast to Float64 and perform the iters in Float64). I leave it for others to say if exposure of this (potentially dangerous) SIMD operation belongs in a SIMD library vs Base.