arithmoi-0.12.1.0: Efficient basic number-theoretic functions.

Copyright | (c) 2016 Andrew Lelechenko
License | MIT
Maintainer | Andrew Lelechenko <andrew.lelechenko@gmail.com>
Safe Haskell | Safe-Inferred
Language | Haskell2010

Math.NumberTheory.Recurrences.Bilinear

Contents

- Pascal triangle
- Other recurrences

Description

Bilinear recurrent sequences and Bernoulli numbers, roughly covering Ch. 5-6 of Concrete Mathematics by R. L. Graham, D. E. Knuth and O. Patashnik.

Note on memory leaks and memoization. Top-level definitions in this module are polymorphic, so the results of computations are not retained in memory. Make them monomorphic to take advantages of memoization. Compare

    >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,385,512 bytes)
    1
    >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,616 bytes)
    1

against

    >>> let binomial' = binomial :: [[Integer]]
    >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,696 bytes)
    1
    >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 391,152 bytes)
    1

Synopsis

- binomial :: Semiring a => [[a]]
- binomialRotated :: Semiring a => [[a]]
- binomialLine :: (Enum a, GcdDomain a) => a -> [a]
- binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a]
- binomialFactors :: Word -> Word -> [(Prime Word, Word)]
- stirling1 :: (Num a, Enum a) => [[a]]
- stirling2 :: (Num a, Enum a) => [[a]]
- lah :: Integral a => [[a]]
- eulerian1 :: (Num a, Enum a) => [[a]]
- eulerian2 :: (Num a, Enum a) => [[a]]
- bernoulli :: Integral a => [Ratio a]
- euler :: forall a. Integral a => [a]
- eulerPolyAt1 :: forall a. Integral a => [Ratio a]
- faulhaberPoly :: (GcdDomain a, Integral a) => Int -> [Ratio a]

# Pascal triangle

binomial :: Semiring a => [[a]] Source #

Infinite zero-based table of binomial coefficients (also known as Pascal triangle).

    binomial !! n !! k == n! / k! / (n - k)!

Note that `binomial` !! n !! k is asymptotically slower than `binomialLine` n !! k, but imposes only `Semiring` constraint.

    >>> take 6 binomial :: [[Int]]
    [[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1],[1,5,10,10,5,1]]

binomialRotated :: Semiring a => [[a]] Source #

Pascal triangle, rotated by 45 degrees.

    binomialRotated !! n !! k == (n + k)! / n! / k! == binomial !! (n + k) !! k

Note that `binomialRotated` !! n !! k is asymptotically slower than `binomialDiagonal` n !! k, but imposes only `Semiring` constraint.

    >>> take 6 (map (take 6) binomialRotated) :: [[Int]]
    [[1,1,1,1,1,1],[1,2,3,4,5,6],[1,3,6,10,15,21],[1,4,10,20,35,56],[1,5,15,35,70,126],[1,6,21,56,126,252]]

binomialLine :: (Enum a, GcdDomain a) => a -> [a] Source #

The n-th (zero-based) line of `binomial` (and the n-th diagonal of `binomialRotated`).

    >>> binomialLine 5
    [1,5,10,10,5,1]

binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a] Source #

The n-th (zero-based) diagonal of `binomial` (and the n-th line of `binomialRotated`).

    >>> take 6 (binomialDiagonal 5)
    [1,6,21,56,126,252]

binomialFactors :: Word -> Word -> [(Prime Word, Word)] Source #

Prime factors of a binomial coefficient.

    binomialFactors n k == factorise (binomial !! n !! k)

    >>> binomialFactors 10 4
    [(Prime 2,1),(Prime 3,1),(Prime 5,1),(Prime 7,1)]

# Other recurrences

stirling1 :: (Num a, Enum a) => [[a]] Source #

Infinite zero-based table of Stirling numbers of the first kind.

    >>> take 5 (map (take 5) stirling1)
    [[1],[0,1],[0,1,1],[0,2,3,1],[0,6,11,6,1]]

Complexity: `stirling1 !! n !! k` is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks `stirling1 !! i !! j` for `0 <= i <= n` and `max(0, k - n + i) <= j <= k`.

One could also consider `unsignedStirling1st` from combinat package to compute stand-alone values.

stirling2 :: (Num a, Enum a) => [[a]] Source #

Infinite zero-based table of Stirling numbers of the second kind.

    >>> take 5 (map (take 5) stirling2)
    [[1],[0,1],[0,1,1],[0,1,3,1],[0,1,7,6,1]]

Complexity: `stirling2 !! n !! k` is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks `stirling2 !! i !! j` for `0 <= i <= n` and `max(0, k - n + i) <= j <= k`.

One could also consider `stirling2nd` from combinat package to compute stand-alone values.

lah :: Integral a => [[a]] Source #

Infinite one-based table of Lah numbers. `lah !! n !! k` equals to lah(n + 1, k + 1).

    >>> take 5 (map (take 5) lah)
    [[1],[2,1],[6,6,1],[24,36,12,1],[120,240,120,20,1]]

Complexity: `lah !! n !! k` is O(n ln n) bits long, its computation takes O(k n ln n) time and forces thunks `lah !! n !! i` for `0 <= i <= k`.

eulerian1 :: (Num a, Enum a) => [[a]] Source #

Infinite zero-based table of Eulerian numbers of the first kind.

    >>> take 5 (map (take 5) eulerian1)
    [[],[1],[1,1],[1,4,1],[1,11,11,1]]

Complexity: `eulerian1 !! n !! k` is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks `eulerian1 !! i !! j` for `0 <= i <= n` and `max(0, k - n + i) <= j <= k`.

eulerian2 :: (Num a, Enum a) => [[a]] Source #

Infinite zero-based table of Eulerian numbers of the second kind.

    >>> take 5 (map (take 5) eulerian2)
    [[],[1],[1,2],[1,8,6],[1,22,58,24]]

Complexity: `eulerian2 !! n !! k` is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks `eulerian2 !! i !! j` for `0 <= i <= n` and `max(0, k - n + i) <= j <= k`.

bernoulli :: Integral a => [Ratio a] Source #

Infinite zero-based sequence of Bernoulli numbers, computed via connection with `stirling2`.

    >>> take 5 bernoulli
    [1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30]

Complexity: `bernoulli !! n` is O(n ln n) bits long, its computation takes O(n^3 ln n) time and forces thunks `stirling2 !! i !! j` for `0 <= i <= n` and `0 <= j <= i`.

One could also consider `bernoulli` from combinat package to compute stand-alone values.

euler :: forall a. Integral a => [a] Source #

The same sequence as `euler'`, but with type `[a]` instead of `[Ratio a]` as the denominators in `euler'` are always `1`.

    >>> take 10 euler :: [Integer]
    [1,0,-1,0,5,0,-61,0,1385,0]

eulerPolyAt1 :: forall a. Integral a => [Ratio a] Source #

Infinite zero-based list of the `n`-th order Euler polynomials evaluated at `1`. The algorithm used was derived from Algorithms for Bernoulli numbers and Euler numbers by Kwang-Wu Chen, third formula of the Corollary in page 7. Element-by-element division of sequences A1986631 and A006519 in OEIS.

    >>> take 10 eulerPolyAt1 :: [Rational]
    [1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2]

faulhaberPoly :: (GcdDomain a, Integral a) => Int -> [Ratio a] Source #

Faulhaber's formula.

    >>> sum (map (^ 10) [0..100])
    959924142434241924250
    >>> sum $ zipWith (*) (faulhaberPoly 10) (iterate (* 100) 1)
    959924142434241924250 % 1