# 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