module SeriesIterators

export iterm, icumsum
export bernoulli

using IterTools
using Base.Iterators

import Base: iterate, IteratorSize, length, size

"""
    iterm(g, ind)

Assume a sequence of terms a[k] can be generated by

    a[i0], b = g(i0)
    a[i], b =  g(i, b) for subsequent indices

The `iterm(g, ind)` iterates the sequence, when `ind` is the
iterator of the indices; first index is `i0`.
"""
iterm(gen, it) = TermIterator(gen, it)
struct TermIterator{F,I}
    generate::F
    indices::I
end
function iterate(t::TermIterator)
    vs = iterate(t.indices)
    k, s = @IterTools.ifsomething vs
    v = t.generate(k)
    v[1], (s, v)
end
function iterate(t::TermIterator, s)
    sx, p = s
    vs = iterate(t.indices, sx)
    k, s = @IterTools.ifsomething vs
    v = t.generate(k, p)
    v[1], (s, v)
end
IteratorSize(::Type{<:TermIterator{<:Any,I}}) where I = IteratorSize(I)
length(t::TermIterator) = length(t.indices)
size(t::TermIterator) = size(t.indices)

function gen_bernoulli(::Type{T}) where T<:Integer
    B = Vector{Rational{T}}(undef, 2)
    B[1] = 1
    B[2] = -1 // 2
    Z = zero(Rational{T})
    function _bernoulli(n::Integer)
        b = B
        n <= 1 && return b[n+1]
        isodd(n) && return Z
        l = length(b)
        if l * 2 - 2 <= n
            ext = min(1_000_000 ÷ l^2, 40)
            m = max(n ÷ 2 + 2, l + ext)
            n0 = (length(b) - 1 ) * 2
            sizehint!(b, ((m + 80) >> 6) << 6)
            resize!(b, m)
            _bernoulli!(b, n0, (m - 1) * 2 - 1, n)
        end
        b[n÷2 + 2]
    end
    _bernoulli
end

function _bernoulli!(B::Vector{Rational{T}}, n0::Integer, n::Integer, nn::Integer) where T
    lastx = 0
    try
        for m = n0:2:n
            sum = -B[1] // T(m + 1)
            sum -= B[2]
            mok = T((m * (m - 1)) ÷ 2)
            for k = 2:2:m-1
                sum -= B[k÷2 + 2] // T(m - k + 1) * mok
                mok *= T(m - k) // T(k + 1)
                mok *= T(m - k - 1) // T(k + 2)
            end
            lastx = m ÷ 2 + 2
            B[lastx] = sum
        end
    catch
        if lastx > 0
            resize!(B, lastx)
        end
        if lastx < nn ÷ 2 + 2
            rethrow()
        end
    end
    nothing
end
"""
    bernoulli(n)

Return Bernoulli number `B_n`. `B_1 = -1/2`.
"""
const bernoulli = gen_bernoulli(BigInt)

struct CumSum{T}
    itr::T
end

"""
    icumsum(itr)

Return iterator, which lazily delivers the same results as `cumsum(itr)`.
"""
icumsum(itr) = CumSum(itr)
function iterate(c::CumSum)
    v, s = @IterTools.ifsomething iterate(c.itr)
    v, (s, v)
end
function iterate(c::CumSum, (s, p))
    v, s = @IterTools.ifsomething iterate(c.itr, s)
    v = p + v
    v, (s, v)
end
IteratorSize(::Type{<:CumSum{T}}) where T = IteratorSize(T)
length(c::CumSum) = length(c.itr)
size(c::CumSum) = size(c.itr)

end
