def sumpow(n, p=1):
    """sumpow(stop_number[, power]) -> int

    Summation from 1**p to n**p, calculate by Bernulli/Faulgaber's formula.

    >>> sumpow(10)
    55
    >>> sumpow(10, 2)
    385

    Step of each seq is 1. If you need the summation of odd or even,
    Do some mathematical proove that gives, for example,

        summation(n if n is even) =      2 + 4 + 6 + ... + n
                                  = 2 * (1 + 2 + 3 + ... + n/2)
                                  = 2 * summation(n/2)

    The start number is also 1. You may now figure out that,

        summation([5..10]) = summation(10) - summation(5-1)

    Avalable power now implement only for p=[0..4]. TODO

    Complexity is O(p log p), (considered power done in O(log p) time).
    While sum(i**p for i in range(1, n+1)) required O(n log p) time.

    see also: <http://en.wikipedia.org/wiki/Bernoulli_number>"""
    if p == 0:
        return n
    elif p == 1:
        return n * (n + 1) // 2
    elif p == 2:
        return n * (n + 1) * (2 * n + 1) // 6
    elif p == 3:
        return (n * (n + 1) // 2) ** 2
    elif p == 4:
        return (6 * n ** 5 + 15 * n ** 4 + 10 * n ** 3 - n) // 30
    else:
        raise NotImplementedError
