Add a 'partition' recurrence #115

rockbmb · 2018-07-28T00:59:27Z

Closes #107.

rockbmb · 2018-07-28T02:38:22Z

Benchmarks for Math.NumberTheory.partition as of commit 5ff9a29:

benchmarked Recurrencies/Partition function/partition/!!100
time                 237.4 ns   (203.8 ns .. 261.8 ns)
                     0.923 R²   (0.897 R² .. 0.953 R²)
mean                 268.6 ns   (256.3 ns .. 282.2 ns)
std dev              45.85 ns   (39.81 ns .. 54.31 ns)
variance introduced by outliers: 83% (severely inflated)

benchmarked Recurrencies/Partition function/partition/!!1000
time                 3.403 μs   (3.387 μs .. 3.419 μs)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 3.446 μs   (3.437 μs .. 3.457 μs)
std dev              33.26 ns   (26.85 ns .. 44.39 ns)

benchmarked Recurrencies/Partition function/partition/!!10000
time                 44.19 μs   (43.89 μs .. 44.51 μs)
                     1.000 R²   (1.000 R² .. 1.000 R²)
mean                 44.07 μs   (43.98 μs .. 44.17 μs)
std dev              325.4 ns   (271.7 ns .. 417.7 ns)

Bodigrim · 2018-07-28T11:23:01Z

Math/NumberTheory/Recurrencies.hs

+--
+-- Note: @tail@ is applied to @pents@ because otherwise the calculation of
+-- @p(n)@ would involve a duplicated @p(n-1)@ term (see the above example).
+partition :: forall a . Integral a => [a]


It is valuable do not impose Integral a constraint here. Can we do with Num a only?

The reason is that since partition grows very fast, we are often interested not in the value of partition itself, but rather, for instance, in partition(x) mod 100. It is very handy to compute it as partition !! x :: Mod 100. But Mod is not Integral, it has only Num interface.

Bodigrim · 2018-07-28T11:24:53Z

Math/NumberTheory/Recurrencies.hs

+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE RankNTypes #-}
+
+module Math.NumberTheory.Recurrencies


Let us move this stuff to a new Math.NumberTheory.Recurrencies.Pentagonal and use Math.NumberTheory.Recurrencies to re-export Linear, Bilinear and Pentagonal.

Bodigrim · 2018-07-28T12:25:55Z

Math/NumberTheory/Recurrencies.hs

+-- >>> take 10 pents
+-- [0, 1, 2, 5, 7, 12 ,15, 22, 26, 35]
+pents :: Integral a => [a]
+pents = map pent pentIndexes


If we need only consecutive applications of pent, then we can speed up computations using identities

map pent [0..] = scanl (\acc n -> acc + 3 * n - 2) 0 [1..] map pent [0,-1..] = scanl (\acc n -> acc + 3 * n - 1) 0 [1..]

This implementation avoids both multiplication and division and requires Num a only.

The first identity holds, but the second one does not seem to:

[0,2,7,15,26,40,57,77,100,126]

where the first few generalized pentagonal numbers are [0, 1, 2, 5, 7, 12, 15, 22, 26, 35, 40, 51, 57, 70, 77, 92, 100, 117, 126]. I would bet only those with positive indexes are being calculated here.

> take 10 $ scanl (\acc n -> acc + 3 * n - 1) (-1) [-1, -2..] [-1,-5,-12,-22,-35,-51,-70,-92,-117,-145]

This means that

map pents [0..] = interleave (scanl (\acc n -> acc + 3 * n - 1) 0 [1..]) (map abs $ scanl (\acc n -> acc + 3 * n - 1) (-1) [-1, -2..])

where interleave :: [a] -> [a] -> [a] weaves two lists together, element by element, starting with the left argument list.

The first identity holds, but the second one does not seem to

The second identity looks good to me:

> take 10 (map pent [0,-1..]) [0,2,7,15,26,40,57,77,100,126] > take 10 (scanl (\acc n -> acc + 3 * n - 1) 0 [1..]) [0,2,7,15,26,40,57,77,100,126]

@Bodigrim My apologies, I misread the left hand side of the second identity...

@Bodigrim I just noticed, won't a Enum a be necessary for pents in any case?

You are right about Enum a, my fault.

No problem, although I made another mistake in thinking it a problem. (Enum a, Num a, Ord a) is necessary for Integral a, but not sufficient. I distractedly assumed it was, and that the former implied the latter.

In fact, Mod m has instances for (Enum a, Num a, Ord a), but none for Integral a, which means having those 3 constraints on partition is just fine. Will push a fix.

@Bodigrim

'pents' has been refactored per @Bodigrim's remarks.

Bodigrim · 2018-07-28T21:17:59Z

Math/NumberTheory/Recurrencies.hs

+-- [1, 2, -3, -4, 5, 6]
+pentagonalSigns :: Num a => [a] -> [a]
+pentagonalSigns = helper True
+ where


Can we express this manual recursion via some list combinators? Looks like zipWith and cycle may be applicable.

Bodigrim · 2018-07-28T21:18:19Z

Math/NumberTheory/Recurrencies.hs

+ interleave _ _ = []
+
+-- | When calculating the @n@-th partition number @p(n)@ using the sum
+-- @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(11) + ...@, the signs of each


Typo, should be p(n-11).

Bodigrim · 2018-07-28T21:18:30Z

Math/NumberTheory/Recurrencies.hs

+ helper _ _ = []
+
+-- @partition !! n@ calculates the @n@-th partition number:
+-- @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(11) + ...@, where @p(0) = 1@


Typo, should be p(n-11).

Bodigrim · 2018-07-28T21:21:12Z

Math/NumberTheory/Recurrencies.hs

+--
+-- Note: @tail@ is applied to @pents@ because otherwise the calculation of
+-- @p(n)@ would involve a duplicated @p(n-1)@ term (see the above example).
+partition :: forall a . (Enum a, Num a, Ord a) => [a]


Please write a test checking that, for example, partition :: [Mod 10] coincides with partition :: [Integer] modulo 10.

'partition' and 'pents' are now found in 'Math.NumberTheory.Recurrencies.Pentagonal', which is reexported by the reintroduced 'Math.NumberTheory.Recurrencies' module.

Also fixes the name of a test in 'Math.NumberTheory.Recurrencies.PentagonalTests'.

rockbmb · 2018-07-29T02:04:32Z

@Bodigrim can you try doing take 10 (Math.NumberTheory.Recurrencies.partition) :: [Math.NumberTheory.Moduli.Mod 10] in a cabal new-repl session?

It calculates the head of the list, and that's as far as it goes. Consumed 12GB of RAM before I killed the session.

rockbmb · 2018-07-29T05:47:04Z

test-suite/Math/NumberTheory/Recurrencies/PentagonalTests.hs

+ p = case someNatVal n of
+ SomeNat (_ :: Proxy t) -> t
+ driver :: KnownNat p => Bool
+ driver = take m' (partition :: KnownNat n => [Mod p]) ==


@Bodigrim I haven't made this work yet, but the idea is present in this and the next lines. I'm not sure how I'll reconcile the [SomeMod] produced by helper with the [Mod n] that partition return.

We can either apply fromInteger :: Integer -> Mod p to partition :: [Integer] and compare with partition :: [Mod p], or on contrary extract integers from partition :: [Mod p] by getVal and compare with map (flip mod p) partition :: [Integer].

Or wrap everything up to SomeMod, employing its Eq instance, applying SomeMod constructor to Mod p and flip modulo p to Integer.

or on contrary extract integers from partition :: [Mod p] by getVal and compare with map (flip mod p) partition :: [Integer]

I'll need to somehow reflect the Natural passed to the test as an argument into a Nat. someNatVal does :: Natural -> SomeNat, but I'm not sure how to do SomeNat -> Nat, but at the type level.

Actually you cannot pass any term of type Nat to partitionProperty2, because Nat is an uninhabited type (but "inhabited" kind). So we should go other way around: pass Natural to partitionProperty and reflect it on type level as Nat, e. g., following pattern from modulo.

I did it just before you posted this, now I'm fixing the fact that the test will fail because smallcheck: spec: internal error: Unable to commit 1048576 bytes of memory. This is related to below's comments about using Map Int a in partition. I'll commit the corrected test for now.

Bodigrim · 2018-07-29T09:14:54Z

Math/NumberTheory/Recurrencies/Pentagonal.hs

+ let n' = (sum .
+ pentagonalSigns .
+ map (\m -> dict M.! (n - m)) .
+ takeWhile (\m -> n - m >= 0) .


@Bodigrim can you try doing take 10 (Math.NumberTheory.Recurrencies.partition) :: [Math.NumberTheory.Moduli.Mod 10] in a cabal new-repl session?

It calculates the head of the list, and that's as far as it goes. Consumed 12GB of RAM before I killed the session.

I can reproduce the issue. AFAIU it happens because the condition n - m >= 0 :: [Mod m] is always true in modular arithmetic. Can we actually employ Map Int a or even IntMap a for caching here?

Can we actually employ Map Int a or even IntMap a for caching here?

Are you asking if it is possible, or if we should it at all?

My impression is that it is possible and desirable. maxBound :: Int is big enough to store indices, so no worries about possible overflow here.

Even after using Map Int a, doing partition :: [Mod 10] won't calculate more than 7 terms. I think it's because of the following line, tail) pents. By asking for partition with type [Mod 10], what happens is that pents also has that type, and this spoils the remainder of the calculation. The following takeWhile won't stop where it's supposed to, and so on.

I'll try giving pents type [Int] and see what happens.

rockbmb · 2018-07-30T02:48:24Z

@Bodigrim Travis CI's Haddock builds tell me:

50% (  2 /  4) in 'Math.NumberTheory.Recurrencies.Pentagonal'
  Missing documentation for:
    Module header

The module header part is the one that has the Copyright:/Maintainer fields, right? Is there anything in particular I need to put there?

rockbmb · 2018-07-30T04:13:22Z

Benchmark with Map:

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!1000
time                 8.449 μs   (7.654 μs .. 9.516 μs)
                     0.919 R²   (0.883 R² .. 0.952 R²)
mean                 8.802 μs   (8.394 μs .. 9.332 μs)
std dev              1.450 μs   (1.187 μs .. 1.838 μs)
variance introduced by outliers: 81% (severely inflated)

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!10000
time                 74.68 μs   (66.34 μs .. 84.05 μs)
                     0.888 R²   (0.822 R² .. 0.943 R²)
mean                 97.04 μs   (91.75 μs .. 106.8 μs)
std dev              21.70 μs   (15.41 μs .. 34.40 μs)
variance introduced by outliers: 89% (severely inflated)

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!100000
time                 3.629 ms   (3.455 ms .. 3.833 ms)
                     0.981 R²   (0.961 R² .. 0.994 R²)
mean                 3.370 ms   (3.294 ms .. 3.461 ms)
std dev              241.9 μs   (184.1 μs .. 371.2 μs)
variance introduced by outliers: 37% (moderately inflated)

Benchmark criterion: FINISH

Benchmark with IntMap:

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!1000
time                 4.951 μs   (4.893 μs .. 5.000 μs)
                     0.999 R²   (0.999 R² .. 1.000 R²)
mean                 4.969 μs   (4.941 μs .. 5.009 μs)
std dev              114.5 ns   (89.21 ns .. 157.1 ns)

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!10000
time                 59.54 μs   (58.64 μs .. 60.56 μs)
                     0.998 R²   (0.996 R² .. 1.000 R²)
mean                 60.45 μs   (60.03 μs .. 61.42 μs)
std dev              1.987 μs   (1.065 μs .. 3.635 μs)
variance introduced by outliers: 16% (moderately inflated)

benchmarked Recurrencies/Pentagonal/Partition function/partition/!!100000
time                 3.019 ms   (2.993 ms .. 3.041 ms)
                     1.000 R²   (0.999 R² .. 1.000 R²)
mean                 3.023 ms   (3.009 ms .. 3.035 ms)
std dev              38.90 μs   (29.73 μs .. 56.17 μs)

Benchmark criterion: FINISH

rockbmb · 2018-07-30T04:30:19Z

@Bodigrim Inside partition, do you have any preference for

go :: forall a . (Enum a, Num a, Ord a) => IM.IntMap a -> Int -> [a]
    go dict !n =
        let intN = fromEnum n
            n' = (sum .
                  pentagonalSigns .
                  map (\m -> dict IM.! (intN - m)) .
                  takeWhile (\m -> intN >= m) .
                  tail) (pents :: [Int])
            dict' = IM.insert intN n' dict
        in n' : go dict' (n+1)

or

go :: forall a . (Enum a, Num a, Ord a) => IM.IntMap a -> [Int] -> [a]
    go dict (n : ns) =
        let intN = fromEnum n
            n' = (sum .
                  pentagonalSigns .
                  map (\m -> dict IM.! (intN - m)) .
                  takeWhile (\m -> intN >= m) .
                  tail) (pents :: [Int])
            dict' = IM.insert intN n' dict
        in n' : go dict' ns
    go _ _ = []

?

Bodigrim · 2018-07-30T21:25:10Z

Math/NumberTheory/Recurrencies/Pentagonal.hs

+--
+-- Note: @tail@ is applied to @pents@ because otherwise the calculation of
+-- @p(n)@ would involve a duplicated @p(n-1)@ term (see the above example).
+partition :: forall a . (Enum a, Num a, Ord a) => [a]


Enum and Ord constraints seem to be redundant.

Bodigrim · 2018-07-30T21:25:26Z

Math/NumberTheory/Recurrencies/Pentagonal.hs

+partition :: forall a . (Enum a, Num a, Ord a) => [a]
+partition = 1 : go (IM.singleton 0 1) [1..]
+ where
+ go :: forall a . (Enum a, Num a, Ord a) => IM.IntMap a -> [Int] -> [a]


Enum and Ord constraints seem to be redundant.

Bodigrim · 2018-07-30T21:26:18Z

Math/NumberTheory/Recurrencies/Pentagonal.hs

+ where
+ go :: forall a . (Enum a, Num a, Ord a) => IM.IntMap a -> [Int] -> [a]
+ go dict (n : ns) =
+ let intN = fromEnum n


AFAIU both n and intN have type Int, so this cast is redundant.

Bodigrim · 2018-07-30T21:30:27Z

test-suite/Math/NumberTheory/Recurrencies/PentagonalTests.hs

+ pentagonalSigns .
+ map (\m -> partition' (n - m)) .
+ takeWhile (\m -> n - m >= 0) .
+ tail $ pents)


If your intention was to compare partition against the reference implementation, it is worth to copy-paste reference definitions of pents and pentagonalSigns to tests as well. Otherwise, if one breaks, for instance, pents, this test would not be able to catch it.

It is also worth to check that first, like, 20 values of partition match the baseline. Brutal, but practical.

@Bodigrim it is probably not necessary to copy-paste pents as pentagonalNumbersProperty1 in this module checks ~~that~~ if the sequence is correct.

Bodigrim · 2018-07-30T21:31:58Z

@Bodigrim Inside partition, do you have any preference for...

I'd prefer the first snippet.

Redundant 'Enum, Ord' constraints have been removed from 'partition'.

A note about the maximum effective index of `partition`'s index has also been added.

rockbmb · 2018-07-31T02:18:25Z

@Bodigrim The PR should be complete now.

Bodigrim · 2018-07-31T17:33:31Z

Merged. Thanks for your efforts and contribution!

rockbmb added 2 commits July 27, 2018 05:47

[Bodigrim#107] Add 'partition' function

c558fae

[Bodigrim#107] Refactor, add proper comments

5ff9a29

[Bodigrim#107] Use a 'Map' to memoize previous values of 'partition'

ed802a8

Bodigrim reviewed Jul 28, 2018

View reviewed changes

rockbmb added 2 commits July 28, 2018 19:54

[Bodigrim#107] Refactor 'pents'

3922eef

'pents' has been refactored per @Bodigrim's remarks.

[Bodigrim#107] Remove 'Integral' constraints

279ac8d

Bodigrim reviewed Jul 28, 2018

View reviewed changes

rockbmb added 3 commits July 28, 2018 22:11

[Bodigrim#107] Move 'partition' and its helpers into a separate module

f150c01

'partition' and 'pents' are now found in 'Math.NumberTheory.Recurrencies.Pentagonal', which is reexported by the reintroduced 'Math.NumberTheory.Recurrencies' module.

[Bodigrim#107] Refactor 'pentagonalSigns'

f7100f4

[Bodigrim#107] Fix comments involving partition's definition

e2f61c6

Also fixes the name of a test in 'Math.NumberTheory.Recurrencies.PentagonalTests'.

rockbmb force-pushed the partition-function branch from e81a619 to e2f61c6 Compare July 29, 2018 00:48

rockbmb commented Jul 29, 2018

View reviewed changes

[Bodigrim#107][WIP] Add residues modulo test to 'partition

44d690a

rockbmb force-pushed the partition-function branch from 09858ad to 44d690a Compare July 29, 2018 05:51

Bodigrim reviewed Jul 29, 2018

View reviewed changes

rockbmb force-pushed the partition-function branch 2 times, most recently from 3bfa707 to 44d690a Compare July 29, 2018 17:13

rockbmb added 2 commits July 29, 2018 17:14

Change newly created file's line endings to LF

ed8979a

[Bodigrim#107] Fix partition's residues modulo test

5fa37ea

rockbmb force-pushed the partition-function branch 2 times, most recently from 7047f84 to 8d619e5 Compare July 30, 2018 02:46

rockbmb force-pushed the partition-function branch from 8d619e5 to 2d89c42 Compare July 30, 2018 03:38

rockbmb force-pushed the partition-function branch from e4fbc32 to db8555d Compare July 30, 2018 04:35

rockbmb force-pushed the partition-function branch from db8555d to c7dadbd Compare July 30, 2018 06:15

Bodigrim reviewed Jul 30, 2018

View reviewed changes

rockbmb added 2 commits July 30, 2018 22:24

[Bodigrim#107] Correct partition :: [Mod n]'s behavior

4ffd44a

[Bodigrim#107] Use 'IntMap' instead of 'Map' in 'partition'

cd69906

rockbmb force-pushed the partition-function branch from bbbd45c to 38d8ed5 Compare July 30, 2018 22:24

[Bodigrim#107] Refactor 'partition' and its tests

b8f767b

Redundant 'Enum, Ord' constraints have been removed from 'partition'.

rockbmb force-pushed the partition-function branch 2 times, most recently from 2a041f4 to 988acfb Compare July 31, 2018 01:22

[Bodigrim#107] Add module header documentation

4cfdfef

A note about the maximum effective index of `partition`'s index has also been added.

rockbmb force-pushed the partition-function branch from 988acfb to 4cfdfef Compare July 31, 2018 02:12

Bodigrim merged commit 06b399d into Bodigrim:master Jul 31, 2018

rockbmb deleted the partition-function branch July 31, 2018 17:44

Add a 'partition' recurrence #115

Add a 'partition' recurrence #115

Conversation

rockbmb commented Jul 28, 2018 • edited Loading

rockbmb commented Jul 28, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Bodigrim Jul 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rockbmb Jul 28, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rockbmb commented Jul 29, 2018

rockbmb Jul 29, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rockbmb Jul 29, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rockbmb commented Jul 30, 2018

rockbmb commented Jul 30, 2018

rockbmb commented Jul 30, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Bodigrim Jul 30, 2018 • edited Loading

Choose a reason for hiding this comment

rockbmb Jul 30, 2018 • edited Loading

Choose a reason for hiding this comment

Bodigrim commented Jul 30, 2018

rockbmb commented Jul 31, 2018

Bodigrim commented Jul 31, 2018

rockbmb commented Jul 28, 2018 •

edited

Loading

Bodigrim Jul 28, 2018 •

edited

Loading

rockbmb Jul 28, 2018 •

edited

Loading

rockbmb Jul 29, 2018 •

edited

Loading

rockbmb Jul 29, 2018 •

edited

Loading

Bodigrim Jul 30, 2018 •

edited

Loading

rockbmb Jul 30, 2018 •

edited

Loading