# 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:** 1\
**Showing post:** 11

<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) 

```

---

_[View the full topic](https://discourse.julialang.org/t/performance-advice-needed/33467)._