Implement inverseTotient

[Bodigrim]

This is still an early draft of a future framework for inversion of multiplicative functions, part of #60, implementing an algorithm along the lines of https://arxiv.org/pdf/1401.6054.pdf

It is fascinating how powerful polymorphism is. Here is how we can compute a list of inverses of the totient function at 120. These are all numbers n such that totient n == 120.

> inverseTotient 120 :: [Integer]
[155,310,183,366,225,450,175,350,231,462,143,286,244,372,396,308,248]

What if there are too many inverses to keep in memory and it is enough only count them? Just switch to another type:

> inverseTotient 120 :: Counter Integer
Counter 17

Interested in minimal or maximal preimage? No worries.

> inverseTotient 120 :: Compose Option Min Integer
Compose (Option (Just (Min 143)))
> inverseTotient 120 :: Compose Option Max Integer
Compose (Option (Just (Max 462)))

Two largest preimages? Easy as a pie.

> inverseTotient 120 :: Max2 Integer
Max2 (Just 450) (Just 462)

[rockbmb]

@Bodigrim this looks great, I'm taking a look at the paper as well. I don't think we have σ_k implemented in the library, do we?

[Bodigrim]

Actually we have, https://github.com/cartazio/arithmoi/blob/master/Math/NumberTheory/ArithmeticFunctions/Standard.hs#L112

[rockbmb]

@Bodigrim found what I was looking for, nevermind.

[Bodigrim]

It is getting tougher than I thought. First, there must be a bug, taking into account failing tests for sigma function. Second, I need more experiments to achieve a decent performance.

[Bodigrim]

OK, it seems to be ready. At least I was able to reproduce all results from the paper, including the inverse sum-of-divisors of 10^1000. @rockbmb @b-mehta and anyone else, please review.

I expect that there are dark corners, where implementation becomes quite convoluted and obscure for an external observer. Please report such places and I'll add clarifying comments. (Since I am deeply in the context, it is difficult for me to identify them myself.)

[rockbmb]

@Bodigrim

I was able to reproduce all results from the paper inverse sum-of-divisors of 10^1000.

Excellent work. I'll take a look, although I cannot guarantee any interesting insights, this code seems quite profound.

[rockbmb]

@Bodigrim I can't reproduce the 10^1000 result, how long did execution take for you? I presume that you didn't do it in ghci, but not even an executable did it for me.

[Bodigrim]

It took 4139 sec and 576 Mb or RAM.

[b-mehta]

Nice work! I'll take a close look in the next few days.

[Bodigrim]

@rockbmb thanks for valuable suggestions.

[rockbmb]

Took another look, couldn't find anything significant.

LGTM, if you don't have anything else to add or fix I suggest you merge it now since making this PR larger will make it more difficult to review, and also should this be rolled back for some reason, larger PRs that touch many things are more difficult to handle.

[b-mehta]

For interest, we can also compute the product of preimages:

data Prod a = Prod { count :: Word
                   , value :: a
                   }
                   deriving Show

instance Monoid a => Semiring (Prod a) where
  zero = Prod 0 mempty
  one = Prod 1 mempty
  (Prod n a) `plus` (Prod m b) = Prod (n + m) (a <> b)
  (Prod n a) `times` (Prod m b) = Prod (n * m) (stimes m a <> stimes n b)

then

>>> inverseTotient (Prod 1 . Product) 120             
Prod {count = 17, value = Product {getProduct = 239173504760939546437745014990080000000000}}

matches the expected value, and counts the number of preimages as a bonus. (It is not hard to check this forms a commutative semiring, and that the map from P_fin(Z) is a weak homomorphism: note that disjointness is required here!).

Another semiring can be given by wrapping Set to ensure it contains no more than a fixed number of elements, to get the first n preimages:

newtype FirstThree a = Three (S.Set a)
  deriving Show

instance (Ord a, Monoid a) => Semiring (FirstThree a) where
  zero = Three zero
  one = Three one
  (Three a) `plus` (Three b) = Three (S.take 3 (a `plus` b))
  (Three a) `times` (Three b) = Three (S.take 3 (a `times` b))

then

>>> inverseSigma (Three . S.singleton . Product) 120
Three (fromList [Product {getProduct = 54},Product {getProduct = 56},Product {getProduct = 87}])

[Bodigrim]

@rockbmb @b-mehta thanks for reviewing!

[Bodigrim]

@b-mehta Prod and FirstThree are very cunning, cool stuff :)

[rockbmb]

@b-mehta using your Data.Ord.Down tip with the FirstThree datatype means we can also find the largest N preimages, by doing newtype LastThree a = Three (S.Set (Down a)). This is really quite cunning.

[b-mehta]

Absolutely! I was tempted to use the trick from the Chudnovsky tests to generalise to First n also...