# Performance advice needed

**URL:** <https://discourse.julialang.org/t/performance-advice-needed/33467>\
**Category:** Performance\
**Tags:** fftw\
**Created:** [January 16, 2020, 11:58pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467 "2020-01-16T23:58:53Z")\
**Posts on this page:** 18\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 16, 2020, 11:58pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/1 "2020-01-16T23:58:53Z")

</div>

Hi,

I am working on a function with calculation that have structure similar to FFT. I’m a bit disappointed, cause I hoped it will be 10x slower than FFTW but benchmarks show its 100x slower. i tried to put @simd and @inbounds but it have minimal effect. Could you tell me if there is something obviously wrong? I don’t ask about algorithm analysis, but just the code quality. I am not interested in adding multithreading, cause this package is just a part of a bigger thing, and parallelisation will be potentially done on a higher level.

The code is in `fnt!` function.

[https://github.com/jakubwro/NumberTheoreticTransforms.jl/blob/ed0e99694da5a64a47ff301a9e9c8ea0ae2737d7/src/fnt.jl#L56](https://github.com/jakubwro/NumberTheoreticTransforms.jl/blob/ed0e99694da5a64a47ff301a9e9c8ea0ae2737d7/src/fnt.jl#L56)

Example benchmark code:

```julia
using NumberTheoreticTransforms, FFTW

const x = mod.(rand(Int, 4096), 65537);
@btime fnt($x, $169, $65537); # 6.306 ms (4 allocations: 32.42 KiB)
@btime fft($x); #61.354 μs (53 allocations: 131.02 KiB)

const x2 = mod.(rand(Int, 8192), 65537);
@btime fnt($x2, $225, $65537); #15.061 ms (4 allocations: 64.45 KiB)
@btime fft($x2); #130.380 μs (53 allocations: 259.02 KiB)

```

---

<div class="post-metadata">

**Author:** ![rdeits](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rdeits/32/286_2.png) [@rdeits](https://discourse.julialang.org/u/rdeits)\
**Post date:** [January 17, 2020, 12:13am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/2 "2020-01-17T00:13:17Z")

</div>

```julia
for M in 2 .^ [0:logN-1;]

```

This line is doubly inefficient: it first constructs `0:logN`, which is cheap and fine, but then it collects that into a new vector, then it computes an entire new vector of `2 ^ ...`. That’s two extra vector allocations that you don’t need.

Instead, you could just iterate with something like:

```julia
for i in 0:logN 
  x = 2^i
  ...

```

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 12:19am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/3 "2020-01-17T00:19:40Z")

</div>

No performance change, this vector is really short. But thanks, I didn’t like this line anyway 🙂

---

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [January 17, 2020, 12:40am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/4 "2020-01-17T00:40:15Z")

</div>

Did you try loading and storing `x[i]` and `x[j]` only once per iteration (by assigning them to local variables)?

---

<div class="post-metadata">

**Author:** ![tkf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkf/32/17635_2.png) [@tkf](https://discourse.julialang.org/u/tkf)\
**Post date:** [January 17, 2020, 12:42am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/5 "2020-01-17T00:42:15Z")

</div>

Also, can’t `powermod` be pulled out of the inner loop?

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 17, 2020, 6:21am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/6 "2020-01-17T06:21:06Z")

</div>

@tkf is right, `powermod` should be avoided. `Mod` itself is rather costly operation, and `powermod` is even more expensive.  
You can remove `powermod` from inner loop and iteratively calculate it, which gives the necessary speed up.

```julia
function fnt2!(x::Array{T, 1}, g::T, q::T) where {T<:Integer}
    N = length(x)
    @assert ispow2(N)
    @assert isfermat(q)

    radix2sort!(x)

    logN = log2(N) |> ceil |> T

    for M in 2 .^ [0:logN-1;] # TODO: not very readable
        interval = 2M
        p = div(N, interval)
        gp = powermod(g, p, q)
        W = 1
        for m in 1:M
            for i in m:interval:N
               j = i + M
               Wxj = W * x[j]
               x[i], x[j] = x[i] + Wxj, x[i] - Wxj
               x[i] = mod(x[i], q)
               x[j] = mod(x[j], q)
            end
            W = mod(W*gp, q)
        end
    end

    return x
end

function fnt2(x::Array{T}, g::T, q::T) where {T<:Integer}
    return fnt2!(copy(x), g, q)
end

# Sanity check
const x = mod.(rand(Int, 4096), 65537);
all(fnt2(x, 169, 65537) .== fnt(x, 169, 65537)) # true

@btime fnt($x, $225, $65537) # 6.923 ms (4 allocations: 32.42 KiB)
@btime fnt2($x, $225, $65537) # 431.841 μs (4 allocations: 32.42 KiB)

```

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 17, 2020, 6:30am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/7 "2020-01-17T06:30:37Z")

</div>

You can get additional speed up by turning 2 `mod` operations in the inner loop to one

```julia
function fnt3!(x::Array{T, 1}, g::T, q::T) where {T<:Integer}
    N = length(x)
    @assert ispow2(N)
    @assert isfermat(q)

    radix2sort!(x)

    logN = log2(N) |> ceil |> T

    for M in 2 .^ [0:logN-1;] # TODO: not very readable
        interval = 2M
        p = div(N, interval)
        gp = powermod(g, p, q)
        W = 1
        for m in 1:M
            for i in m:interval:N
               j = i + M
               Wxj = mod(W * x[j], q)
               x[i], x[j] = x[i] + Wxj, x[i] - Wxj + q
               x[i] = x[i] >= q ? x[i] - q : x[i]
               x[j] = x[j] >= q ? x[j] - q : x[j]
            end
            W = mod(W*gp, q)
        end
    end

    return x
end

function fnt3(x::Array{T}, g::T, q::T) where {T<:Integer}
    return fnt3!(copy(x), g, q)
end

# Sanity check
const x = mod.(rand(Int, 4096), 65537)
all(fnt3(x, 169, 65537) .== fnt(x, 169, 65537)) # true

@btime fnt3($x, $169, $65537) # 245.112 μs (4 allocations: 32.42 KiB)

```

EDIT: Changed slightly definition of `x[j]`, it doesn’t affect performance, but initial version introduced bug for unsigned integers, because

```julia
mod(UInt16(7) - UInt16(9), UInt16(11)) # 0x0007
mod(UInt16(7) - UInt16(9) + UInt16(11), UInt16(11)) # 0x0009

# and also
UInt16(7) - UInt16(9) < 0 # false

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 17, 2020, 7:41am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/8 "2020-01-17T07:41:29Z")

</div>

The next performance improvement would be to calculate q^-1 (in the group sense), so your mood operations can be done as `x[j]-x[j]q^-1`.

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 10:32am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/9 "2020-01-17T10:32:57Z")

</div>

@Skoffer, it’s brilliant, all unit tests passed.

```julia
julia> @btime fnt($x2, $225, $65537); #15.061 ms (4 allocations: 64.45 KiB)
  644.937 μs (2 allocations: 64.08 KiB)

```

@tkf, after storing x[i] and x[j] in variables

```julia
julia> @btime fnt($x2, $225, $65537); #15.061 ms (4 allocations: 64.45 KiB)
  593.164 μs (2 allocations: 64.08 KiB)

```

I also added @inbounds to outer loop:

```julia
julia> @btime fnt($x2, $225, $65537); #15.061 ms (4 allocations: 64.45 KiB)
  574.171 μs (2 allocations: 64.08 KiB)

```

Adding @simd has no effect, probably cause my laptop is from 2013. I will check tomorrow on a desktop.

Summing up: 26x speedup, and about 5x slower than FFTW. Amazing, I had no idea that mod operations are so costly.

I assume I can commit your suggestions to my MIT licensed lib, fine?

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 10:35am UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/10 "2020-01-17T10:35:33Z")

</div>

@Oscar_Smith, I am not sure I understand. q is the modulus, so q == 0 mod q, so it has no inverse.

---

<div class="post-metadata">

**Author:** ![Skoffer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skoffer/32/378_2.png) [@Skoffer](https://discourse.julialang.org/u/Skoffer)\
**Post date:** [January 17, 2020, 3:03pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/11 "2020-01-17T15:03:06Z")

</div>

Actually you can get additional boost by exploiting the fact, that q is fermat number. This is little more involved, of course. Idea is the following:

1. We need to find `a mod q` where `a < q^2` and `q = 2^2^p + 1` where p some integer.
2. `a` can be presented in the form `a = n*q + a'`, here `a' = a mod q` and n some integer.
3. `a` can be further presented in the form `a = n*q + a' = n*(q - 1) + n + a'`, so `n + a' = a mod (q - 1) = a & 2^2^p - 1`, here we used the fact that `mod 2^l` in binary representation is just logical and with `2^j - 1`.
4. `n` can be found as `a div q - 1 = a div 2^2^p = a >>> 2^p` where we have used the fact that in binary representation division by power of 2 equals to corresponding shift.

Together it gives us this nice function

```julia
function fermat_mod(a::T, q::T) where T <: Integer
    x = a & (q - T(2)) - a >>> trailing_zeros(q - T(1)) + q
    x = x >= q ? x - q : x
end

```

in this function it is assumed that `q` is fermat and `a < q^2`.

Adding it into `fnt!` function yields us

```julia
function fnt2!(x::Array{T, 1}, g::T, q::T) where {T<:Integer}
    N = length(x)
    @assert ispow2(N)
    @assert isfermat(q)

    radix2sort!(x)

    logN = log2(N) |> ceil |> T

    for M in 2 .^ [0:logN-1;] # TODO: not very readable
        interval = 2M
        p = div(N, interval)
        gp = powermod(g, p, q)
        W = 1
        for m in 1:M
            for i in m:interval:N
               j = i + M
               Wxj = W * x[j]
               Wxj = Wxj & (q - T(2)) - Wxj >>> trailing_zeros(q - T(1)) + q
               Wxj = Wxj >= q ? Wxj + q : Wxj

               x[i], x[j] = x[i] + Wxj, x[i] - Wxj + q
               x[i] = x[i] >= q ? x[i] - q : x[i]
               x[j] = x[j] >= q ? x[j] - q : x[j]
            end
            W = mod(W*gp, q)
        end
    end

    return x
end

```

And benchmark

```julia
@btime fnt2($x, $169, $65537) # 176.820 μs (4 allocations: 32.42 KiB) 

```

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 3:45pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/12 "2020-01-17T15:45:41Z")

</div>

Excellent! You will beat FFTW till tomorrow 🙂

@btime fnt($x2, $225, $65537);  
600.204 μs (4 allocations: 64.45 KiB)

@btime fnt3($x2, $225, $65537);  
451.002 μs (4 allocations: 64.45 KiB)

@btime fft($x2); #130.380 μs (53 allocations: 259.02 KiB)  
130.582 μs (53 allocations: 259.02 KiB)

451/130 # fnt / fft  
3.4692307692307693

Your reasoning is correct, but there were a minor typo that caused errors for inputs close to q.

`Wxj = Wxj >= q ? Wxj + q : Wxj` → `Wxj = Wxj >= q ? Wxj - q : Wxj`

This is performance I could not even imagine when writing original post.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 17, 2020, 4:31pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/13 "2020-01-17T16:31:22Z")

</div>

Sorry for the confusion, I was typing on my phone while tired. I’m having a slightly hard time explaining what I mean, but [https://gmplib.org/~tege/divcnst-pldi94.pdf](https://gmplib.org/~tege/divcnst-pldi94.pdf) does a good job of it. The TLDR is the group I was talking about was the `Int64`s under multiplication. Since modulo can be computed cheaply once you have division, this should speed up the repeated mod q operations.

---

<div class="post-metadata">

**Author:** ![kristoffer.carlsson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kristoffer.carlsson/32/22_2.png) [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)\
**Post date:** [January 17, 2020, 4:43pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/14 "2020-01-17T16:43:24Z")

</div>

> [@Oscar\_Smith](#):
>
> Sorry for the confusion, I was typing on my phone while tired. I’m having a slightly hard time explaining what I mean, but [https://gmplib.org/~tege/divcnst-pldi94.pdf](https://gmplib.org/~tege/divcnst-pldi94.pdf) does a good job of it. The TLDR is the group I was talking about was the `Int64` s under multiplication. Since modulo can be computed cheaply once you have division, this should speed up the repeated mod q operations.

Note that Julia already has something like this in [https://github.com/JuliaLang/julia/blob/master/base/multinverses.jl](https://github.com/JuliaLang/julia/blob/master/base/multinverses.jl).

```julia
julia> smi = Base.multiplicativeinverse(3432)
Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64}(3432, 5503923639708211205, 0, 0x0a)

julia> div(21313221, smi)
6210

julia> div(21313221, 3432)
6210

```

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 17, 2020, 4:44pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/15 "2020-01-17T16:44:23Z")

</div>

Good to know, that makes this much easier to actually use.

---

<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:** [January 17, 2020, 5:18pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/16 "2020-01-17T17:18:49Z")

</div>

> [@Jakub\_Wronowski](#):
>
> `@btime fft($x2); #130.380 μs (53 allocations: 259.02 KiB) 130.582 μs (53 allocations: 259.02 KiB)`

Note that by calling `fft(x2)` you are not really seeing FFTW’s full speed. Try:

```julia
p = plan_fft(x2, flags=FFTW.PATIENT)
@btime $p * $x2

```

or with pre-allocated output and pre-allocated conversion of `x2` to the type that FFTW uses:

```julia
x2c = ComplexF64.(x2); y2c = copy(x2c);
@btime mul!($y2c, $p, $x2c);

```

or taking advantage of the fact that the input is real:

```julia
pr = plan_rfft(x2, flags=FFTW.PATIENT);
@btime $pr * $x2;

```

plus pre-allocated input/output:

```julia
x2f = Float64.(x2); y2cr = Array{ComplexF64}(undef, length(x2f)÷2+1);
@btime mul!($y2cr, $pr, $x2f);

```

With all of the tricks, for the final result I get about a factor of 5–8 faster than `fft(x2)`.

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 5:39pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/17 "2020-01-17T17:39:22Z")

</div>

Thanks for explaining, I gave FFT just as an example of algorithm that should have similar complexity to see how slow we are with FNT. I benchmarked this whole plan creation and preallocation and it’s 211 μs, so I understand it will be faster if we have a lot of vectors to transform, but surely you are right here.  
I added again suggestions from @tkf, so now fnt() runs 382 μs (451 μs was before). I think this is amazing result for pure Julia implementation.  
BTW, are you aware of any pure Julia FFT implementations? Maybe I should compare timings also with them.

---

<div class="post-metadata">

**Author:** ![Jakub\_Wronowski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jakub_wronowski/32/204030_2.png) [@Jakub\_Wronowski](https://discourse.julialang.org/u/Jakub_Wronowski)\
**Post date:** [January 17, 2020, 5:55pm UTC](https://discourse.julialang.org/t/performance-advice-needed/33467/18 "2020-01-17T17:55:02Z")

</div>

Thank you for sharing, it looks very interesting, but @Skoffer eliminated all mod operations, except ` W = mod(W * gp, q)`, which is executed just few times. I will try to use it anyway on the original source I posted.

UPDATE: as I promissed I did testing with those multiplicative inverses. I did also some further improvements which gave:

- 349.411 μs for @Skoffer’s solution
- 385.504 μs for @Oscar_Smith’s solution
- 596.341 μs for using `Base.mod`

Also this multiplicative inverse works just for machine word sized integers. This FNT approach is really only useful with BigInts to get a very high precision in number representation so this is a deal breaker for me.