# This file is part of OLMS (Open Library of Mathematical Sequences).
# Copyright (c) 2017: Peter Luschny. License is MIT.

module SwingingFactorial
using Primes

doc"""
Returns the accumulated product of an array.
"""
function ∏(A)
    function prod(a::Int, b::Int)
        n = b - a
        if n < 24
            p = BigInt(1)
            for k in a:b
                p *= A[k]
            end
            return p
        end
        m = div(a + b, 2)
        prod(a, m) * prod(m + 1, b)
    end
    prod(1, length(A))
end

const SwingOddpart = [1,1,1,3,3,15,5,35,35, 315, 63, 693, 231,
   3003, 429, 6435, 6435, 109395,12155,230945,46189,969969,
   88179,2028117, 676039,16900975,1300075,35102025,5014575,
   145422675,9694845,300540195,300540195]

doc"""
Computes the odd part of the swinging factorial ``n≀``. Cf. A163590.
"""
function swing_oddpart(n::Int)
    n < 33 && return BigInt(SwingOddpart[n+1])

    sqrtn = isqrt(n)
    factors = primes(div(n,2) + 1, n)
    r = primes(sqrtn + 1, div(n, 3))
    s = filter(x -> isodd(div(n, x)),r)
    append!(factors, s)

    for prime in primes(3, sqrtn)
        p, q = 1, n
        while true
            q = div(q, prime)
            q == 0 && break
            isodd(q) && (p *= prime)
        end
        p > 1 && push!(factors, p)
    end
    return ∏(factors)
end

doc"""
Computes the swinging factorial (a.k.a. swing numbers n≀). Cf. A056040.
"""
function swing(n::Int)
    sh = count_ones(div(n,2))
    swing_oddpart(n) << sh
end

const FactorialOddPart = [1, 1, 1, 3, 3, 15, 45, 315, 315, 2835, 14175, 155925,
    467775, 6081075, 42567525, 638512875, 638512875, 10854718875, 97692469875,
    1856156927625, 9280784638125, 194896477400625, 2143861251406875,
    49308808782358125, 147926426347074375, 3698160658676859375]

doc"""
Return the largest odd divisor of ``n!``. Cf. A049606.
"""
function factorial_oddpart(n::Int)
    n < length(FactorialOddPart) && return BigInt(FactorialOddPart[n+1])
    swing_oddpart(n)*(factorial_oddpart(div(n,2))^2)
end

doc"""
Return the factorial ``n! = 1*2*...*n``, which is the order of the
symmetric group S_n or the number of permutations of n letters. Cf. A000142.
"""
function factorial(n::Int)
    n < 0 && ArgumentError("Argument must be ≥ 0")
    sh = n - count_ones(n)
    factorial_oddpart(n) << sh
end

for n in 0:90
    println(factorial(n))
end

end # module