# Help with speeding up this code

**URL:** <https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407>\
**Category:** Performance\
**Created:** [July 9, 2023, 5:17pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407 "2023-07-09T17:17:16Z")\
**Posts on this page:** 12\
**Page:** 2

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [July 9, 2023, 10:04pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/21 "2023-07-09T22:04:18Z")

</div>

If not all performance needs to be squeezed out, just memoizing the straight-forward implementation might already be enough:

```julia
using Memoize

@memoize Dict function facm(n::T) where {T <: Integer}
    if iszero(n)
        one(T)
    else
        n * facm(n - 1)
    end
end

@memoize Dict function loopm(n::T, k::T) where {T <: Integer}
    if n < zero(T) || k <= zero(T)
        error("Invalid values")
    elseif n > zero(T) && n < k
        error("Invalid values")
    end

    if k == one(T) || n == zero(T)
        return one(T)
    end
    cnt = zero(T)
    for j = 1:(n÷k)
        m = min(k-1, n-j*k)
        if m == zero(T)
             m = one(T)
        end
        recur = sum(loopm(n-j*k, t) for t = 1:m)
        cnt += recur*(facm(n) ÷ facm(n-j*k)) ÷ (facm(j)*(k^j))
    end
    return cnt
end

function main(n)
    [loopm(n,i) / facm(n) for i = 1:n]
end

```

---

<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:** [July 9, 2023, 11:17pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/22 "2023-07-09T23:17:35Z")

</div>

> [@Mark\_Szlazak](#):
>
> ```julia
> if k == 1 || n == 0
> return 1.0
> end
> 
> ```

My generic advice here is that if you want recursion to be fast (in any language), you should ideally enlarge the base case (“coarsen” the recursion) — see [this comment](https://discourse.julialang.org/t/recursive-call-vs-while-loop/7723/18).

(That being said, maybe with memo-ization here you won’t call the trivial base cases enough to matter?)

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [July 9, 2023, 11:56pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/23 "2023-07-09T23:56:32Z")

</div>

Did a little refactoring, as gdalle said, to make the code totally unrecognizable. Specifically, used the semantic meaning of the `loops` values to find a more simply calculated expression.  
It goes like this (using memoization also):

```julia
# returns number of perms of 1..n with largest cycle of size k
function loop!(M,CM,n,k)
    n < k && return 0.0 # base and trivial cases
    n == 0 && k == 0 && return 1.0
    k <= 0 && return 0.0
    n <= 0 && return 0.0
    M[n,k] >= 0.0 && return M[n,k] # memoization
    if k == 1
        M[n, k] = 1.0
        return 1.0
    end
    factor = 1.0 # non base case calculation
    res = 0.0
    for i in 1:k-1
        res += factor * loop!(M,CM,n-i,k)
        factor *= n-i
    end
    if n>k
        j = k
        while j>1 && CM[n-k,j] < 0.0
            j -= 1
        end
        for i in j:k # second memoization for M cumsum
            CM[n-k, i] = (i>1 ? CM[n-k, i-1] : 0.0) + loop!(M,CM,n-k, i)
        end    
        res += factor * CM[n-k,k]
    else
        res += factor
    end
    M[n,k] = res
    return res
end

```

To use this function, the memoization matrices need to be initialized and passed as parameters, as follows:

```julia
julia> const M = fill(-1.0,100,100);

julia> const CM = fill(-1.0,100,100);

julia> @time [loop!(M, CM, 100, i) for i in 1:100];
  0.027222 seconds (23.15 k allocations: 1.561 MiB, 98.51% compilation time)

```

The timing is not very informative, ~~as I’m not including any other time to compare with~~. But it seems at least \>1000x speedup.

UPDATE: To benchmark, adding a `main` as in previous posts:

```julia
function main(n)
    M = fill(-1.0, n,n)
    CM = fill(-1.0, n,n)
    y = Vector{Float64}(undef, n)
    fac = cumprod(float.(1:n))
    for i = 1:n
        y[i] = loop!(M, CM, n, i) / fac[n]
    end
    return y
end

```

with this `main` and comparing to one of the original `oldmain`s:

```julia
julia> @time oldmain(100);
  5.305958 seconds (4 allocations: 3.500 KiB)

julia> @time main(100);
  0.001171 seconds (7 allocations: 158.969 KiB)

```

( 5000x speedup )

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [July 10, 2023, 7:34am UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/24 "2023-07-10T07:34:22Z")

</div>

Since you kept the floats this might come at the cost of numerical errors that are hard to quantify. I’d love to give IntervalArithmetic a go here, to figure out how much precision we lose when giving up on BigInt

---

<div class="post-metadata">

**Author:** ![Mark\_Szlazak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mark_szlazak/32/51393_2.png) [@Mark\_Szlazak](https://discourse.julialang.org/u/Mark_Szlazak)\
**Post date:** [July 10, 2023, 5:11pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/25 "2023-07-10T17:11:09Z")

</div>

I tried Dan’s code at n=100 and it is fastest but that is not the case when I tried it with n=200. The code of gdalle with better memoization near the top of this thread runs a lot faster at n=200. Here is what I measured:

Dan(200):  
1004.275557 seconds (61.17 k allocations: 4.950 MiB, 0.01% compilation time)

gdale(200):  
2.237682 seconds (23.36 M allocations: 488.616 MiB, 25.59% gc time, 5.58% compilation time)

Also at n=100, the bar graph with Dan’s code is what it should be. However at n=200 for Dan’s code, I only get the first couple of bars in bar graph and nothing more. Many similar warnings to the following accompany the bar graph:

┌ Warning: Indices Base.OneTo(200) of attribute `seriescolor` does not match data indices 1:4.  
└ @ Plots C:\Users\Mark.julia\packages\Plots\QWCV6\src\utils.jl:141  
┌ Info: Data contains NaNs or missing values, and indices of `seriescolor` vector do not match data indices.  
│ If you intend elements of `seriescolor` to apply to individual NaN-separated segments in the data,  
│ pass each segment in a separate vector instead, and use a row vector for `seriescolor`. Legend entries  
│ may be suppressed by passing an empty label.  
│ For example,  
└ plot([1:2,1:3], [[4,5],[3,4,5]], label=[“y” “”], seriescolor=[1 2])

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [July 10, 2023, 8:05pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/26 "2023-07-10T20:05:24Z")

</div>

> [@Mark\_Szlazak](#):
>
> The code of gdalle with better memoization near the top of this thread runs a lot faster at n=200.

That’s weird. Maybe when we memoize too much we end up spending a lot of time in the memory access

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [July 10, 2023, 8:44pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/27 "2023-07-10T20:44:37Z")

</div>

There seems to be some overflow (or other trouble) to produce NaNs. Also, my version uses a lot more memory (it says 5GB) which could potentially exhaust the memory on the machine which causes serious slowdown.

One workaround would be to use 128 float numbers. This can be found in DoubleFloats package ([GitHub - JuliaMath/DoubleFloats.jl: math with more good bits](https://github.com/JuliaMath/DoubleFloats.jl) ).

It’s worth finding the first n for which Dan(n) starts misbehaving (for easier debugging).

ADDENDUM:

The maximum number representable in Float64 is **1.7E+308**.  
While 200! is **7.9E+374**. So that’s a problem.

The DoubleFloat solution should work for that. As for the speed, it may be related, as NaN would mess up the memoization conditions.

---

<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:** [July 10, 2023, 8:53pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/28 "2023-07-10T20:53:37Z")

</div>

doublefloats won’t help here. they don’t have a wider range. you might want float128 instead.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [July 10, 2023, 8:58pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/29 "2023-07-10T20:58:08Z")

</div>

> [@Oscar\_Smith](#):
>
> you might want float128 instead

Thanks. I was just quickly googling longer floats. Without a REPL open.  
Float128 is provided by the Quadmath package.  
Example:

```julia
julia> using Quadmath

julia> prod([Float128(i) for i in 1:200])
7.88657867364790503552363213932184275e+374

```

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [July 10, 2023, 9:11pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/30 "2023-07-10T21:11:46Z")

</div>

The fixed up code works okay for `n=200`:

```julia
using Quadmath

# returns number of perms of 1..n with largest cycle of size k
function loop!(M::Matrix{T},CM::Matrix{T},n::Integer,k::Integer) where {T}
    n < k && return zero(T) # base and trivial cases
    n == 0 && k == 0 && return one(T)
    k <= 0 && return zero(T)
    n <= 0 && return zero(T)
    M[n,k] >= 0 && return M[n,k] # memoization
    if k == 1
        M[n, k] = one(T)
        return one(T)
    end
    factor = one(T) # non base case calculation
    res = zero(T)
    for i in 1:k-1
        res += factor * loop!(M,CM,n-i,k)
        factor *= n-i
    end
    if n>k
        j = k
        while j>1 && CM[n-k,j] < 0
            j -= 1
        end
        for i in j:k # second memoization for M cumsum
            CM[n-k, i] = (i>1 ? CM[n-k, i-1] : zero(T)) + loop!(M,CM,n-k, i)
        end    
        res += factor * CM[n-k,k]
    else
        res += factor
    end
    M[n,k] = res
    return res
end

function main(::Type{T}, n) where {T<:Real}
    M = fill(-one(T), n,n)
    CM = fill(-one(T), n,n)
    y = Vector{T}(undef, n)
    fac = cumprod(T.(1:n))
    for i = 1:n
        y[i] = loop!(M, CM, n, i) / fac[n]
    end
    return y
end

main(n) = main(Float64, n)

```

With these definitions:

```julia
julia> main(10)
10-element Vector{Float64}:
 2.755731922398589e-7
 0.00261656746031746
...
julia> main(Float128, 200)
200-element Vector{Float128}:
 1.26797695348096242175301637107480335e-375
 4.65660197354737973222078346996522418e-183
...

```

---

<div class="post-metadata">

**Author:** ![Mark\_Szlazak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mark_szlazak/32/51393_2.png) [@Mark\_Szlazak](https://discourse.julialang.org/u/Mark_Szlazak)\
**Post date:** [July 10, 2023, 10:45pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/31 "2023-07-10T22:45:29Z")

</div>

Dan(200):  
0.200233 seconds (76.84 k allocations: 6.629 MiB, 43.89% compilation time)

---

<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:** [July 26, 2023, 8:37pm UTC](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407/33 "2023-07-26T20:37:11Z")

</div>

Split this topic on July 26, 2023, 8:37pm UTC

3 posts were split to a new topic: [Question about sum(x → …, 1:n) notation](https://discourse.julialang.org/t/question-about-sum-x-1-n-notation/102132)

[Previous page](https://discourse.julialang.org/t/help-with-speeding-up-this-code/101407.md?page=1)