################################################################################
# Quantum analogs
#
# Copyright (C) 2020 Ulrich Thiel, ulthiel.com/math
################################################################################

export quantum_integer, quantum_factorial, quantum_binomial

@doc raw"""
	quantum_integer(n::Int, q::RingElem)

Let ``n âˆˆ â„¤`` and let ``â„š(ğ��ª)`` be the fraction field of the polynomial ring ``â„¤[ğ��ª]``. The **quantum integer** ``[n]_ğ��ª âˆˆ â„š(ğ��ª)`` of ``n`` is defined as
```math
[n]_ğ��ª â‰” \frac{ğ��ª^n-1}{ğ��ª-1} \;.
```
We have
```math
[n]_ğ��ª = \sum_{i=0}^{n-1} ğ��ª^i âˆˆ â„¤[ğ��ª] \quad \text{if } n â‰¥ 0
```
and
```math
[n]_ğ��ª = -ğ��ª^{n} [-n]_ğ��ª \quad \text{for any } n âˆˆ â„¤ \;,
```
hence
```math
[n]_ğ��ª = - \sum_{i=0}^{-n-1} ğ��ª^{n+i} âˆˆ â„¤[ğ��ª^{-1}] \quad \text{ if } n < 0 \;.
```
This shows in particular that actually
```math
[n]_ğ��ª âˆˆ â„¤[ğ��ª,ğ��ª^{-1}] âŠ‚ â„š(ğ��ª) \quad \text{ for any } n âˆˆ â„¤ \;.
```
Now, for an element ``q`` of a ring ``R`` we define ``[n]_q âˆˆ R`` as the specialization of ``[n]_ğ��ª`` in ``q`` using the two equations aboveâ€”assuming that ``q`` is invertible in ``R`` if ``n<0``. Note that for ``q=1`` we obtain
```math
[n]_1 = n \quad \text{for for any } n âˆˆ â„¤ \;,
```
so the quantum integers are "deformations" of the usual integers.

If the argument ``q`` is not specified, then ``[n]_ğ��ª`` is returned as an element of ``â„¤[ğ��ª,ğ��ª^{-1}]``. If ``q`` is a Julia integer, then it is taken as an element of ```ZZ``` if ``n â‰¥ 0`` or ``q = Â± 1``, otherwise it is taken as an element of ```QQ```.

# Examples
```julia-repl
julia> quantum_integer(3)
q^2 + q + 1
julia> quantum_integer(-3)
-q^-1 - q^-2 - q^-3
julia> quantum_integer(-3,2)
-7//8
```

# References
1. Conrad, K. (2000). A q-analogue of Mahler expansions. I. *Adv. Math., 153*(2), 185--230. [https://doi.org/10.1006/aima.1999.1890](https://doi.org/10.1006/aima.1999.1890)
2. Kac, V. & Cheung, P. (2002). *Quantum calculus*. Springer-Verlag, New York.
3. Wikipedia, [Q-analog](https://en.wikipedia.org/wiki/Q-analog).
"""
function quantum_integer(n::Int, q::RingElem)

	R = parent(q)
	if isone(q)
		return R(n)
	else
		z = zero(R)
		if n >= 0
			for i=0:n-1
				z += q^i
			end
		else
			for i=0:-n-1
				z -= q^(n+i)
			end
		end
		return z
	end
end

function quantum_integer(n::Int, q::Int)
	if n >= 0 || q == 1 || q == -1
		return quantum_integer(n,ZZ(q))
	else
		return quantum_integer(n,QQ(q))
	end
end

function quantum_integer(n::Int)
	R,q = LaurentPolynomialRing(ZZ, "q")
	return quantum_integer(n,q)
end


@doc raw"""
	quantum_factorial(n::Int, q::RingElem)

For a non-negative integer ``n`` and an element ``q`` of a ring ``R`` the **quantum factorial** ``[n]_q! âˆˆ R`` is defined as
```math
[n]_q! â‰” [1]_q â‹… â€¦ â‹… [n]_q âˆˆ R \;.
```
Note that for ``q=1`` we obtain
```math
[n]_1! = n! \quad \text{ for all } n âˆˆ â„¤ \;,
```
hence the quantum factorial is a "deformation" of the usual factorial.

# Examples
```julia-repl
julia> quantum_factorial(3)
q^3 + 2*q^2 + 2*q + 1
```
"""
function quantum_factorial(n::Int, q::RingElem)

	n >= 0 || throw(ArgumentError("n â‰¥ 0 required"))

	R = parent(q)
	if isone(q)
		return R(factorial(n))
	else
		z = one(R)
		for i=1:n
			z *= quantum_integer(i,q)
		end
		return z
	end
end

function quantum_factorial(n::Int, q::Int)
	return quantum_factorial(n,ZZ(q))
end

function quantum_factorial(n::Int)
	R,q = LaurentPolynomialRing(ZZ, "q")
	return quantum_factorial(n,q)
end

@doc raw"""
	quantum_binomial(n::Int, q::RingElem)

Let ``k`` be a non-negative integer and let ``n âˆˆ â„¤``. The **quantum binomial** ``\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª \in â„š(ğ��ª)`` is defined as
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª â‰” \frac{[n]_ğ��ª!}{[k]_ğ��ª! [n-k]_ğ��ª!} = \frac{[n]_ğ��ª [n-1]_ğ��ªâ‹… â€¦ â‹… [n-k+1]_ğ��ª}{[k]_ğ��ª!}
```
Note that the first expression is only defined for ``n â‰¥ k`` since the quantum factorials are only defined for non-negative integersâ€”but the second  expression is well-defined for all ``n âˆˆ â„¤`` and is used for the definition. In Corad (2000) it is shown that
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª = \sum_{i=0}^{n-k} q^i \begin{bmatrix} i+k-1 \\ k-1 \end{bmatrix}_ğ��ª \quad \text{if } n â‰¥ k > 0 \;.
```
Since
```math
\begin{bmatrix} n \\ 0 \end{bmatrix}_ğ��ª = 1 \quad \text{for all } n âˆˆ â„¤
```
and
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª = 0 \quad \text{if } 0 â‰¤ n < k \;,
```
it follows inductively that
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª âˆˆ â„¤[ğ��ª] \quad \text{if } n â‰¥ 0 \;.
```
For all ``n âˆˆ â„¤`` we have the relation
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª = (-1)^k ğ��ª^{-k(k-1)/2+kn} \begin{bmatrix} k-n-1 \\ k \end{bmatrix}_ğ��ª \;,
```
which shows that
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª âˆˆ â„¤[ğ��ª^{-1}] \quad \text{if } n < 0 \;.
```
In particular,
```math
\begin{bmatrix} n \\ k \end{bmatrix}_ğ��ª âˆˆ â„¤[ğ��ª,ğ��ª^{-1}] \quad \text{for all } n âˆˆ â„¤ \;.
```
Now, for an element ``q`` of a ring ``R`` we define ``\begin{bmatrix} n \\ k \end{bmatrix}_q`` as the specialization of ``\begin{bmatrix} n \\ k \end{bmatrix}_{\mathbf{q}}`` in ``q``, where ``q`` is assumed to be invertible in ``R`` if ``n < 0``.

Note that for ``q=1`` we obtain
```math
\begin{bmatrix} n \\ k \end{bmatrix}_1 = {n \choose k} \;,
```
hence the quantum binomial coefficient is a "deformation" of the usual binomial coefficient.

# Examples
```julia-repl
julia> quantum_binomial(4,2)
q^4 + q^3 + 2*q^2 + q + 1
julia> quantum_binomial(19,5,-1)
36
julia> K,i = CyclotomicField(4);
julia> quantum_binomial(17,10,i)
0
```

# References
1. Conrad, K. (2000). A q-analogue of Mahler expansions. I. *Adv. Math., 153*(2), 185--230. [https://doi.org/10.1006/aima.1999.1890](https://doi.org/10.1006/aima.1999.1890)
2. Wikipedia, [Gaussian binomial coefficient](https://en.wikipedia.org/wiki/Gaussian_binomial_coefficient)
"""
function quantum_binomial(n::Int, k::Int, q::RingElem)

	k >= 0 || throw(ArgumentError("k â‰¥ 0 required"))

	R = parent(q)
	if isone(q)
		return R(binomial(n,k))
	elseif k == 0
		return one(R)
	elseif k == 1
		return quantum_integer(n,q)
	elseif n >= 0
		if n < k
			return zero(R)
		else
			z = zero(R)
			for i=0:n-k
				z += q^i * quantum_binomial(i+k-1,k-1,q)
			end
			return z
		end
	elseif n<0
		return (-1)^k * q^(div(-k*(k-1),2) + k*n) * quantum_binomial(k-n-1,k,q)
	end

end

function quantum_binomial(n::Int, k::Int, q::Int)
	if n > 0
		return quantum_binomial(n,k,ZZ(q))
	else
		return quantum_binomial(n,k,QQ(q))
	end
end

function quantum_binomial(n::Int, k::Int)
	R,q = LaurentPolynomialRing(ZZ, "q")
	return quantum_binomial(n,k,q)
end
