### Author:  [@linwaytin](https://discourse.julialang.org/u/linwaytin)
#### Post date: [May 7, 2026, 9:08am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/1 "2026-05-07T09:08:43Z")
This is an old post but it seems it hasn’t been discussed yet.
> **[Making Julia as Fast as C++ · FLOW Lab](https://flow.byu.edu/posts/julia-c++)**
An interesting observation the author found is that non-integer powers in Julia are fairly slow at that time, compared to C++. I am not sure if this is really the case. He mentioned removing those operations speeds up the Julia code by more than 10 times. This seems too much, but maybe at that time there were some problems in the non-integer powers.
Note that this post was in 2019 so kind of old, but I think Julia 1.0 had been released at that time.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 7, 2026, 9:37am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/2 "2026-05-07T09:37:04Z")
> [@linwaytin](#):
>
> it seems it hasn’t been discussed yet
A lot of things changed then (really, too many to enumerate), but the basic strategies remain relevant.
That said, it is unclear _what_ you want to discuss about this blog post. Asking a specific question would focus the discussion.
> [@linwaytin](#):
>
> maybe at that time there were some problems in the non-integer powers
Cf
>
>
> LLVM will lower an inverse sqrt to \`rsqrt\` instruction in fast mode, e.g.: https…://godbolt.org/z/wMLxq1
>
> However, \[\`sqrt\_fast\`\](https://github.com/JuliaLang/julia/blob/425c53d33cfe38e57051f784b4260ca0b799ad9a/base/fastmath.jl#L280) is missing the \`fast\` flag, preventing LLVM's lowering.
>
> On Skylake+ this less important, however for older architectures using \`rsqrt\` can be ~2-3x as fast, e.g. on Ivy Bridge, 6 vs. 14 cycles.
which was closed 5 years ago.
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [May 7, 2026, 10:02am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/3 "2026-05-07T10:02:18Z")
Really interesting, indeed its a bit old and I think they wanted to keep things standard julia but with a little of packages you can beat the pythonic version by 97x without changing much :
```julia-auto
using ConcreteStructs,LinearAlgebra, StaticArrays, BenchmarkTools
function generate_particles(PType, n, lambda; l=1, Gamma=SA[0.0,0.0,0.0])
sigma = l/(n-1)*lambda
particles = fill(zero(PType), n^3)
xs = range(0, stop=l, length=n)
# Adds particles in a regular lattice
ind = 1
for k in 1:n
for j in 1:n
for i in 1:n
X = SA[xs[i], xs[j], xs[k]]
particles[ind] = PType(X, Gamma, sigma)
ind += 1
end
end
end
return particles
end
"""
This is a particle struct made up of ambiguous (unspecified) types
"""
@concrete struct ParticleAmbiguous
# User inputs
X # Position
Gamma # Vectorial circulation
sigma # Smoothing radius
# Properties
U # Velocity at particle
J # Jacobian at particle J[i,j]=dUi/dxj
end
ParticleAmbiguous(X, Gamma, sigma; U=zeros(3), J=zeros(3,3)) = ParticleAmbiguous(X, Gamma, sigma, U, J )
# Empty initializer
Base.zero(::Type{<:ParticleAmbiguous}) = ParticleAmbiguous(SA[0.0,0.0,0.0], SA[0.0,0.0,0.0], 0.0)
# Winckelmans regularizing kernel
g_wnk(r) = r^3 * (r^2 + 2.5) / (r^2 + 1)^2.5
dgdr_wnk(r) = 7.5 * r^2 / (r^2 + 1)^3.5
const const4 = 1/(4*pi)
"""
Pythonic programming approach
"""
@views @fastmath function P2P_pythonic(particles, g, dgdr)
for Pi in particles
@inbounds for Pj in particles
dX = Pi.X - Pj.X
r = norm(dX)
if r != 0
gsgm = g(r / Pj.sigma)
dgsgmdr = dgdr(r / Pj.sigma)
crss = cross(-const4 * (dX/r^3), Pj.Gamma)
Pi.U .+= gsgm * crss
@inbounds @simd for j in 1:3
Pi.J[:, j] += ( dX[j] / (Pj.sigma*r) * (dgsgmdr*crss) -
gsgm * (3*dX[j]/r^2) * crss -
gsgm * (const4/r^3) * cross(SVector{3}(i==j for i in 1:3), Pj.Gamma) )
end
end
end
end
end
particles_amb = generate_particles(ParticleAmbiguous, 6, 1.0)
@btime P2P_pythonic($particles_amb, $g_wnk, $dgdr_wnk) # 2ms
```
I don’t know how it compares with the c++ if you want to try go ahead
---
### Author:  [@linwaytin](https://discourse.julialang.org/u/linwaytin)
#### Post date: [May 7, 2026, 10:37am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/4 "2026-05-07T10:37:42Z")
Thanks. It is useful to know that now Julia can easily generate code which has similar performance to C++.
And I didn’t know the package ConcreteStructs. Thanks a lot.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 7, 2026, 10:42am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/5 "2026-05-07T10:42:24Z")
> [@linwaytin](#):
>
> the package ConcreteStructs
It is a nice convenience package, but just makes it easier to write idiomatic Julia code, not the source of any speed improvement _per se_. The key is StaticArrays, which does the unrolling etc.
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [May 7, 2026, 10:51am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/6 "2026-05-07T10:51:13Z")
actually without it all the component would be understood as Any, but I think you meant you could just parametrize everywhere which is true
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [May 7, 2026, 10:54am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/7 "2026-05-07T10:54:36Z")
> [@yolhan\_mannes](#):
>
> I think you meant you could just parametrize everywhere which is true
Yes. It is of course tedious.
---
### Author:  [@DanielSussman](https://discourse.julialang.org/u/DanielSussman)
#### Post date: [May 7, 2026, 3:05pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/8 "2026-05-07T15:05:50Z")
Building on Tamas\_Papp point: The basic optimization strategies remain relevant regardless of the language. Tools like ConcreteStructs and StaticArrays are fantastic because they make it much easier to write idiomatic code that is fast, by giving the compiler the type stability and memory layout it needs to do its job.
I think the original post also highlights how hard it is to do these cross-language comparisons well: it’s easy to measure the difference between specific implementations of algorithms, but that is somewhat separate from the question of how close to language / hardware limits those implementations actually are.
For instance, even from a quick glance at the C++ code used in that post, there are obvious low-hanging optimizations to be made. With only a small amount of work one can get benchmarking results like:
```julia-auto
Activating project at `~/test/post-juliavscpp/code`
Benchmarking Original Julia...
Benchmarking Better Julia...
Benchmarking Original C++...
Benchmarking Better C++...
=================================================================
Version | Time (ms) | Speedup
-----------------------------------------------------------------
Julia (final post version) | 2.961 | 1.00x
Julia (yolhan_mannes version) | 2.354 | 1.26x
C++ (post version) | 1.210 | 2.45x
C++ (better version) | 0.400 | 7.40x
```
(All the usual caveats: precise profiling is hard, this is on my hardware, etc etc.)
I’m not trying to declare a “winner,” – I’ve been tremendously enjoying my experience with Julia, but I’ve been using it for almost an order of magnitude fewer years than I’ve been writing code in C/C++. I suspect that an expert could look at the Julia implementations and generate a yet-better version. I guess I just want to reinforce that the benefits of these posts probably lie in the broad principles for optimization within and across languages, but a lot of the specific details of the “X is slower/faster/as-fast-as Y!” discourse is often more about developer optimization effort.
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [May 7, 2026, 3:40pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/9 "2026-05-07T15:40:13Z")
I agree, its always skill issues, what I like about julia is that you can go from the straigntforward but insanly slow code to a very well optimised code in the same language, it also mean julia learning curve can feel overwhelming but still, i don’t know a lot of language that can do that within the same language.
Also if not happy just copy paste the llvm from clang and here is your perf (like it was easy) 😂
---
### Author:  [@linwaytin](https://discourse.julialang.org/u/linwaytin)
#### Post date: [May 7, 2026, 3:51pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/10 "2026-05-07T15:51:37Z")
The benchmark you show is interesting. Now your better C++ code is not only faster than the Julia version by at least 5 times, it also faster than the old C++ version by 3 times. What did you do to make this?
---
### Author:  [@DanielSussman](https://discourse.julialang.org/u/DanielSussman)
#### Post date: [May 7, 2026, 8:26pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/11 "2026-05-07T20:26:26Z")
Just a grab bag of what I think are pretty standard techniques (and work across languages): hoist invariants out of loops (and other applications of the “don’t compute the same thing twice” principle), use a SoA layout of data instead of AoS, etc.
For this particular code, it’s also worth noticing that the combination of `std::pow(...,2.5)` and compiler flag `-o3` is a pessimization (even if you’re also using the -ffast-math flag, depending on the compiler)
---
### Author:  [@linwaytin](https://discourse.julialang.org/u/linwaytin)
#### Post date: [May 7, 2026, 8:47pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/12 "2026-05-07T20:47:20Z")
Awesome, so maybe the same optimization can be applied to the Julia code.
What do you mean by “the combination of `std::pow(...,2.5)` and compiler flag `-o3` is a pessimization”?
---
### Author:  [@yolhan\_mannes](https://discourse.julialang.org/u/yolhan_mannes)
#### Post date: [May 7, 2026, 8:51pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/13 "2026-05-07T20:51:25Z")
yes I thought I could go simple
```julia-auto
AA = StructArray(particles_amb)
@btime P2P_pythonic($AA, $g_wnk, $dgdr_wnk)
```
but its not that easy, you will actually need to change the layout a lot because this will do vector of vectors and we don’t win anything.
A way to cheat but it becomes unfair is
```julia-auto
@views @fastmath function P2P_pythonic(particles, g, dgdr)
Threads.@threads for Pi in particles
@inbounds for Pj in particles
dX = Pi.X - Pj.X
r = norm(dX)
if r != 0
gsgm = g(r / Pj.sigma)
dgsgmdr = dgdr(r / Pj.sigma)
crss = cross(-const4 * (dX/r^3), Pj.Gamma)
Pi.U .+= gsgm * crss
@inbounds @simd for j in 1:3
Pi.J[:, j] += ( dX[j] / (Pj.sigma*r) * (dgsgmdr*crss) -
gsgm * (3*dX[j]/r^2) * crss -
gsgm * (const4/r^3) * cross(SVector{3}(i==j for i in 1:3), Pj.Gamma) )
end
end
end
end
end
```
and I get the 0.4ms of the c++ but I guess the c++ is single threaded
---
### Author:  [@aryavorskiy](https://discourse.julialang.org/u/aryavorskiy)
#### Post date: [May 7, 2026, 10:02pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/14 "2026-05-07T22:02:13Z")
There are some small tweaks you can do to make the thing even faster
- Use mutable `MVector`s instead of usual ones:
```julia-auto
ParticleAmbiguous(X, Gamma, sigma; U=@MArray(zeros(3)), J=@MArray(zeros(3,3))) = ParticleAmbiguous(X, Gamma, sigma, U, J)
```
- Remove `@views` - you do not need them anymore, since `MVector`s are stack-allocated
- You can also remove `@inbounds` and `@simd` on the inner loop for the sake of cleanness - it is too short to be properly vectorized, and the `@inbounds` is redundant because the `@inbounds` on the outer loop already does the job
On my machine this reduced the runtime from `2.13 ms` to `1.96 ms` (Apple M2 pro).
As the author of the post stated, there is not much else we can do, since the power function is the heaviest workload left. However we can do this:
```julia-auto
_pow2d5(x) = x^2 * sqrt(x) # x^2.5
_pow3d5(x) = x^3 * sqrt(x) # x^3.5
g_wnk_fast(r) = r^3 * (r^2 + 2.5) / _pow2d5(r^2 + 1)
dgdr_wnk_fast(r) = 7.5 * r^2 / _pow3d5(r^2 + 1)
```
This yields `0.67 ms` on my machine - close enough. Seriously though, generally it is not a very good idea, you may encounter precision loss
---
### Author:  [@DanielSussman](https://discourse.julialang.org/u/DanielSussman)
#### Post date: [May 7, 2026, 10:04pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/15 "2026-05-07T22:04:25Z")
The `std::pow` function, in general, is required to do extra work (handle subnormals, deal with zero and negatives, set `errno`), and here we can get away with using the optimized `std::sqrt` function.
Using only `-O3` compilers are not allowed to replace `pow(x, 2.5)` with `x*x*sqrt(x)`… Modern gcc and clang are pretty good about making this optimization if you also use the `-ffast-math` flag, but the way this interacts with vectorized versions of functions is tricky.
(Also, yes – I wrote a single-threaded cpp version, to keep things easy)
---
### Author:  [@DanielSussman](https://discourse.julialang.org/u/DanielSussman)
#### Post date: [May 7, 2026, 10:30pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/16 "2026-05-07T22:30:09Z")
Is the concern about precision correct? I’m always happy to be corrected, but I was under the impression that sqrt is guaranteed to be rounded correctly to 0.5ULP, whereas exponentiation in general isn’t… In this case I would expect the replacement of the pow with sqrt to be a change that is both safe and fast.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 7, 2026, 10:33pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/17 "2026-05-07T22:33:27Z")
They will be about roughly equally accurate.
---
### Author:  [@aryavorskiy](https://discourse.julialang.org/u/aryavorskiy)
#### Post date: [May 7, 2026, 10:46pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/18 "2026-05-07T22:46:44Z")
To be honest, the code just slightly reeked to me, I was not completely sure. However, we can run a test using `BigFloat`s:
```julia-auto
using Statistics
A = rand(10000) .+ 4
Ab = big.(A)
res = @. sqrt(Ab ^ 7) # let's consider this ideally precise
std(@. res - A^3.5) # 9.5e-15
std(@. res - A^3 * sqrt(A)) # 2.0e-14
```
2X precision difference 😀
---
### Author:  [@mahmah](https://discourse.julialang.org/u/mahmah)
#### Post date: [May 8, 2026, 4:20am UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/19 "2026-05-08T04:20:22Z")
So which of these fixes and suggestions are still valid these days?
---
### Author:  [@DanielSussman](https://discourse.julialang.org/u/DanielSussman)
#### Post date: [May 8, 2026, 1:03pm UTC](https://discourse.julialang.org/t/making-julia-as-fast-as-c/137016/20 "2026-05-08T13:03:49Z")
We’re outside of my area of expertise, so (a) happy to have an excuse to play around and learn a bit and (b) would love to hear from someone on this forum that knows more! My understanding was that `sqrt` itself is required to be accurate by IEEE to 0.5 ULP, and that as a transcendental function `pow` is recommended but not mandated to be accurate to the same level (and, hence, whether it is or not is is implementation-defined).
I like the idea of using BigFloat as an oracle to easily measure how accurate different versions are. At a minimum, though, I think we should be measuring in ULPs relative to the answer. Something like:
```julia-auto
using Statistics
function ulp_test(f, domain; n_samples=100000)
a, b = extrema(domain)
x = a .+ rand(Float64, n_samples) .* (b - a)
vals = f.(x)
big_float_x = BigFloat.(x)
exact_vals = f.(big_float_x) # "exact"
ulp_scale = eps.(Float64.(exact_vals)) # measure ULP relative to answer
ulp_errors = (vals .- exact_vals) ./ ulp_scale
return (
std_ulp = Float64(std(ulp_errors)),
max_abs_ulp = Float64(maximum(abs, ulp_errors))
)
end
```
Running a few tests:
```julia-auto
julia> ulp_test(x->sqrt(x),(0.0,1.0))
(std_ulp = 0.28975401566751907, max_abs_ulp = 0.49999874016920975)
julia> ulp_test(x->x^pi,(0.0,1.0))
(std_ulp = 0.875775912554217, max_abs_ulp = 9.312060290475852)
julia> ulp_test(x->x^3.1415926,(0.0,1.0))
(std_ulp = 0.2892737541146482, max_abs_ulp = 0.5210410859100494)
```
So far, as expected: the `sqrt` function correctly rounds, and exponentiation doesn’t have to according to this (probably not-quite-correct) function. On the other hand, switching to exponentiation by a quantity which is more clearly a float looks basically fine.
Finally, we can look at the question at hand:
```julia-auto
julia> ulp_test(x->x*x*x*sqrt(x),(0.0,100.0))
(std_ulp = 0.6003549004331392, max_abs_ulp = 2.5012293636064946)
julia> ulp_test(x->x^3.5,(0.0,100.0))
(std_ulp = 0.2887374737640995, max_abs_ulp = 0.5205982181157858)
```
Quite clear: in the first (faster) version, the ULP errors accumulate over all of the internal operations used to build up an expressions, as one would expect. In the second version, Julia is clearly using a slower version of exponentiation that maintains accuracy over the domain.
Whether this difference in accuracy is relevant is, of course, domain-specific.