{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Math.NumberTheory.Recurrences.Bilinear
(
binomial
, binomialRotated
, binomialLine
, binomialDiagonal
, binomialFactors
, stirling1
, stirling2
, lah
, eulerian1
, eulerian2
, bernoulli
, euler
, eulerPolyAt1
, faulhaberPoly
) where
import Data.Euclidean (GcdDomain(..))
import Data.List (scanl', zipWith4)
import Data.Maybe
import Data.Ratio
import Data.Semiring (Semiring(..))
import Numeric.Natural
import Math.NumberTheory.Recurrences.Linear (factorial)
import Math.NumberTheory.Primes
binomial :: Semiring a => [[a]]
binomial :: forall a. Semiring a => [[a]]
binomial = forall a. (a -> a) -> a -> [a]
iterate (\[a]
l -> forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Semiring a => a -> a -> a
plus ([a]
l forall a. [a] -> [a] -> [a]
++ [forall a. Semiring a => a
zero]) (forall a. Semiring a => a
zero forall a. a -> [a] -> [a]
: [a]
l)) [forall a. Semiring a => a
one]
{-# SPECIALIZE binomial :: [[Int]] #-}
{-# SPECIALIZE binomial :: [[Word]] #-}
{-# SPECIALIZE binomial :: [[Integer]] #-}
{-# SPECIALIZE binomial :: [[Natural]] #-}
binomialRotated :: Semiring a => [[a]]
binomialRotated :: forall a. Semiring a => [[a]]
binomialRotated = forall a. (a -> a) -> a -> [a]
iterate (forall a. [a] -> [a]
tail forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' forall a. Semiring a => a -> a -> a
plus forall a. Semiring a => a
zero) (forall a. a -> [a]
repeat forall a. Semiring a => a
one)
{-# SPECIALIZE binomialRotated :: [[Int]] #-}
{-# SPECIALIZE binomialRotated :: [[Word]] #-}
{-# SPECIALIZE binomialRotated :: [[Integer]] #-}
{-# SPECIALIZE binomialRotated :: [[Natural]] #-}
binomialLine :: (Enum a, GcdDomain a) => a -> [a]
binomialLine :: forall a. (Enum a, GcdDomain a) => a -> [a]
binomialLine a
n = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl'
(\a
x (a
k, a
nk1) -> forall a. HasCallStack => Maybe a -> a
fromJust forall a b. (a -> b) -> a -> b
$ (a
x forall a. Semiring a => a -> a -> a
`times` a
nk1) forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
k)
forall a. Semiring a => a
one
(forall a b. [a] -> [b] -> [(a, b)]
zip [forall a. Semiring a => a
one..a
n] [a
n, forall a. Enum a => a -> a
pred a
n..forall a. Semiring a => a
one])
{-# SPECIALIZE binomialLine :: Int -> [Int] #-}
{-# SPECIALIZE binomialLine :: Word -> [Word] #-}
{-# SPECIALIZE binomialLine :: Integer -> [Integer] #-}
{-# SPECIALIZE binomialLine :: Natural -> [Natural] #-}
binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a]
binomialDiagonal :: forall a. (Enum a, GcdDomain a) => a -> [a]
binomialDiagonal a
n = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl'
(\a
x a
k -> forall a. HasCallStack => Maybe a -> a
fromJust (a
x forall a. Semiring a => a -> a -> a
`times` (a
n forall a. Semiring a => a -> a -> a
`plus` a
k) forall a. GcdDomain a => a -> a -> Maybe a
`divide` a
k))
forall a. Semiring a => a
one
[forall a. Semiring a => a
one..]
{-# SPECIALIZE binomialDiagonal :: Int -> [Int] #-}
{-# SPECIALIZE binomialDiagonal :: Word -> [Word] #-}
{-# SPECIALIZE binomialDiagonal :: Integer -> [Integer] #-}
{-# SPECIALIZE binomialDiagonal :: Natural -> [Natural] #-}
binomialFactors :: Word -> Word -> [(Prime Word, Word)]
binomialFactors :: Word -> Word -> [(Prime Word, Word)]
binomialFactors Word
n Word
k
| Word
n forall a. Ord a => a -> a -> Bool
< Word
2
= []
| Bool
otherwise
= forall a. (a -> Bool) -> [a] -> [a]
filter ((forall a. Eq a => a -> a -> Bool
/= Word
0) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> b
snd)
forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (\Prime Word
p -> (Prime Word
p, Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) Word
n forall a. Num a => a -> a -> a
- Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) (Word
n forall a. Num a => a -> a -> a
- Word
k) forall a. Num a => a -> a -> a
- Word -> Word -> Word
mult (forall a. Prime a -> a
unPrime Prime Word
p) Word
k))
[forall a. Bounded a => a
minBound .. forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
precPrime Word
n]
where
mult :: Word -> Word -> Word
mult :: Word -> Word -> Word
mult Word
p Word
m = Word -> Word -> Word
go Word
mp Word
mp
where
mp :: Word
mp = Word
m forall a. Integral a => a -> a -> a
`quot` Word
p
go :: Word -> Word -> Word
go !Word
acc !Word
x
| Word
x forall a. Ord a => a -> a -> Bool
>= Word
p = let xp :: Word
xp = Word
x forall a. Integral a => a -> a -> a
`quot` Word
p in Word -> Word -> Word
go (Word
acc forall a. Num a => a -> a -> a
+ Word
xp) Word
xp
| Bool
otherwise = Word
acc
stirling1 :: (Num a, Enum a) => [[a]]
stirling1 :: forall a. (Num a, Enum a) => [[a]]
stirling1 = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall {a}. (Num a, Enum a) => [a] -> a -> [a]
f [a
1] [a
0..]
where
f :: [a] -> a -> [a]
f [a]
xs a
n = a
0 forall a. a -> [a] -> [a]
: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail (\a
_ a
x a
y -> a
x forall a. Num a => a -> a -> a
+ a
n forall a. Num a => a -> a -> a
* a
y) a
1 [a]
xs a
0
{-# SPECIALIZE stirling1 :: [[Int]] #-}
{-# SPECIALIZE stirling1 :: [[Word]] #-}
{-# SPECIALIZE stirling1 :: [[Integer]] #-}
{-# SPECIALIZE stirling1 :: [[Natural]] #-}
stirling2 :: (Num a, Enum a) => [[a]]
stirling2 :: forall a. (Num a, Enum a) => [[a]]
stirling2 = forall a. (a -> a) -> a -> [a]
iterate forall {a}. (Num a, Enum a) => [a] -> [a]
f [a
1]
where
f :: [a] -> [a]
f [a]
xs = a
0 forall a. a -> [a] -> [a]
: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail (\a
k a
x a
y -> a
x forall a. Num a => a -> a -> a
+ a
k forall a. Num a => a -> a -> a
* a
y) a
1 [a]
xs a
0
{-# SPECIALIZE stirling2 :: [[Int]] #-}
{-# SPECIALIZE stirling2 :: [[Word]] #-}
{-# SPECIALIZE stirling2 :: [[Integer]] #-}
{-# SPECIALIZE stirling2 :: [[Natural]] #-}
lah :: Integral a => [[a]]
lah :: forall a. Integral a => [[a]]
lah = forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall {b}. Integral b => b -> b -> [b]
f (forall a. [a] -> [a]
tail forall a. (Num a, Enum a) => [a]
factorial) [a
1..]
where
f :: b -> b -> [b]
f b
nf b
n = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl (\b
x b
k -> b
x forall a. Num a => a -> a -> a
* (b
n forall a. Num a => a -> a -> a
- b
k) forall a. Integral a => a -> a -> a
`div` (b
k forall a. Num a => a -> a -> a
* (b
k forall a. Num a => a -> a -> a
+ b
1))) b
nf [b
1..b
nforall a. Num a => a -> a -> a
-b
1]
{-# SPECIALIZE lah :: [[Int]] #-}
{-# SPECIALIZE lah :: [[Word]] #-}
{-# SPECIALIZE lah :: [[Integer]] #-}
{-# SPECIALIZE lah :: [[Natural]] #-}
eulerian1 :: (Num a, Enum a) => [[a]]
eulerian1 :: forall a. (Num a, Enum a) => [[a]]
eulerian1 = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall {a}. (Num a, Enum a) => [a] -> a -> [a]
f [] [a
1..]
where
f :: [a] -> a -> [a]
f [a]
xs a
n = a
1 forall a. a -> [a] -> [a]
: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail (\a
k a
x a
y -> (a
n forall a. Num a => a -> a -> a
- a
k) forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ (a
k forall a. Num a => a -> a -> a
+ a
1) forall a. Num a => a -> a -> a
* a
y) a
1 [a]
xs a
0
{-# SPECIALIZE eulerian1 :: [[Int]] #-}
{-# SPECIALIZE eulerian1 :: [[Word]] #-}
{-# SPECIALIZE eulerian1 :: [[Integer]] #-}
{-# SPECIALIZE eulerian1 :: [[Natural]] #-}
eulerian2 :: (Num a, Enum a) => [[a]]
eulerian2 :: forall a. (Num a, Enum a) => [[a]]
eulerian2 = forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall {a}. (Num a, Enum a) => [a] -> a -> [a]
f [] [a
1..]
where
f :: [a] -> a -> [a]
f [a]
xs a
n = a
1 forall a. a -> [a] -> [a]
: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail (\a
k a
x a
y -> (a
2 forall a. Num a => a -> a -> a
* a
n forall a. Num a => a -> a -> a
- a
k forall a. Num a => a -> a -> a
- a
1) forall a. Num a => a -> a -> a
* a
x forall a. Num a => a -> a -> a
+ (a
k forall a. Num a => a -> a -> a
+ a
1) forall a. Num a => a -> a -> a
* a
y) a
1 [a]
xs a
0
{-# SPECIALIZE eulerian2 :: [[Int]] #-}
{-# SPECIALIZE eulerian2 :: [[Word]] #-}
{-# SPECIALIZE eulerian2 :: [[Integer]] #-}
{-# SPECIALIZE eulerian2 :: [[Natural]] #-}
bernoulli :: Integral a => [Ratio a]
bernoulli :: forall a. Integral a => [Ratio a]
bernoulli = forall a.
Integral a =>
([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
helperForBEEP forall a. a -> a
id (forall a b. (a -> b) -> [a] -> [b]
map forall a. Fractional a => a -> a
recip [Ratio a
1..])
{-# SPECIALIZE bernoulli :: [Ratio Int] #-}
{-# SPECIALIZE bernoulli :: [Rational] #-}
faulhaberPoly :: (GcdDomain a, Integral a) => Int -> [Ratio a]
faulhaberPoly :: forall a. (GcdDomain a, Integral a) => Int -> [Ratio a]
faulhaberPoly Int
p
= forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) ((Ratio a
0forall a. a -> [a] -> [a]
:)
forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse
forall a b. (a -> b) -> a -> b
$ forall a. Int -> [a] -> [a]
take (Int
p forall a. Num a => a -> a -> a
+ Int
1) forall a. Integral a => [Ratio a]
bernoulli)
forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a. Integral a => a -> a -> Ratio a
% (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
pforall a. Num a => a -> a -> a
+a
1))
forall a b. (a -> b) -> a -> b
$ forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith forall a. Num a => a -> a -> a
(*) (forall a. (a -> a) -> a -> [a]
iterate forall a. Num a => a -> a
negate (if forall a. Integral a => a -> Bool
odd Int
p then a
1 else -a
1))
forall a b. (a -> b) -> a -> b
$ forall a. Semiring a => [[a]]
binomial forall a. [a] -> Int -> a
!! (Int
pforall a. Num a => a -> a -> a
+Int
1)
euler' :: forall a . Integral a => [Ratio a]
euler' :: forall a. Integral a => [Ratio a]
euler' = forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a.
Integral a =>
([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
helperForBEEP forall a. [a] -> [a]
tail [Ratio a]
as
where
as :: [Ratio a]
as :: [Ratio a]
as = forall a b c d. (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
zipWith3
(\a
sgn a
frac a
ones -> (a
sgn forall a. Num a => a -> a -> a
* a
ones) forall a. Integral a => a -> a -> Ratio a
% a
frac)
(forall a. [a] -> [a]
cycle [a
1, a
1, a
1, a
1, -a
1, -a
1, -a
1, -a
1])
(forall a. [a] -> [a]
dups (forall a. (a -> a) -> a -> [a]
iterate (a
2 forall a. Num a => a -> a -> a
*) a
1))
(forall a. [a] -> [a]
cycle [a
1, a
1, a
1, a
0])
dups :: forall x . [x] -> [x]
dups :: forall a. [a] -> [a]
dups = forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\x
n [x]
list -> x
n forall a. a -> [a] -> [a]
: x
n forall a. a -> [a] -> [a]
: [x]
list) []
{-# SPECIALIZE euler' :: [Ratio Int] #-}
{-# SPECIALIZE euler' :: [Rational] #-}
euler :: forall a . Integral a => [a]
euler :: forall a. Integral a => [a]
euler = forall a b. (a -> b) -> [a] -> [b]
map forall a. Ratio a -> a
numerator forall a. Integral a => [Ratio a]
euler'
eulerPolyAt1 :: forall a . Integral a => [Ratio a]
eulerPolyAt1 :: forall a. Integral a => [Ratio a]
eulerPolyAt1 = forall a. [a] -> [a]
tail forall a b. (a -> b) -> a -> b
$ forall a.
Integral a =>
([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
helperForBEEP forall a. [a] -> [a]
tail (forall a b. (a -> b) -> [a] -> [b]
map forall a. Fractional a => a -> a
recip (forall a. (a -> a) -> a -> [a]
iterate (Ratio a
2 forall a. Num a => a -> a -> a
*) Ratio a
1))
{-# SPECIALIZE eulerPolyAt1 :: [Ratio Int] #-}
{-# SPECIALIZE eulerPolyAt1 :: [Rational] #-}
zipIndexedListWithTail :: Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail :: forall b a. Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b]
zipIndexedListWithTail b -> a -> a -> b
f b
n [a]
as a
a = case [a]
as of
[] -> []
(a
x : [a]
xs) -> b -> a -> [a] -> [b]
go b
n a
x [a]
xs
where
go :: b -> a -> [a] -> [b]
go b
m a
y [a]
ys = case [a]
ys of
[] -> let v :: b
v = b -> a -> a -> b
f b
m a
y a
a in [b
v]
(a
z : [a]
zs) -> let v :: b
v = b -> a -> a -> b
f b
m a
y a
z in (b
v forall a. a -> [a] -> [a]
: b -> a -> [a] -> [b]
go (forall a. Enum a => a -> a
succ b
m) a
z [a]
zs)
{-# INLINE zipIndexedListWithTail #-}
helperForBEEP :: Integral a => ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
helperForBEEP :: forall a.
Integral a =>
([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a]
helperForBEEP [Ratio a] -> [Ratio a]
g [Ratio a]
xs = forall a b. (a -> b) -> [a] -> [b]
map ([Ratio a] -> Ratio a
f forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Ratio a] -> [Ratio a]
g) forall a. (Num a, Enum a) => [[a]]
stirling2
where
f :: [Ratio a] -> Ratio a
f = forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b c d e.
(a -> b -> c -> d -> e) -> [a] -> [b] -> [c] -> [d] -> [e]
zipWith4 (\Ratio a
sgn Ratio a
fact Ratio a
x Ratio a
stir -> Ratio a
sgn forall a. Num a => a -> a -> a
* Ratio a
fact forall a. Num a => a -> a -> a
* Ratio a
x forall a. Num a => a -> a -> a
* Ratio a
stir) (forall a. [a] -> [a]
cycle [Ratio a
1, -Ratio a
1]) forall a. (Num a, Enum a) => [a]
factorial [Ratio a]
xs
{-# SPECIALIZE helperForBEEP :: ([Ratio Int] -> [Ratio Int]) -> [Ratio Int] -> [Ratio Int] #-}
{-# SPECIALIZE helperForBEEP :: ([Rational] -> [Rational]) -> [Rational] -> [Rational] #-}