Implement inverseTotient #142

Bodigrim · 2018-09-28T22:57:36Z

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 · 2018-09-29T17:23:21Z

@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 · 2018-09-29T18:02:44Z

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

rockbmb

How did you profile the function? I don't think you uploaded the Inverse.hs file that your executable used.

rockbmb · 2018-10-24T21:48:43Z

@Bodigrim found what I was looking for, nevermind.

Bodigrim · 2018-12-04T00:14:37Z

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.

…licit unit leading coefficient

Bodigrim · 2018-12-06T23:58:04Z

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 · 2018-12-07T01:06:33Z

@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 · 2018-12-07T18:40:51Z

@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 · 2018-12-07T20:04:09Z

It took 4139 sec and 576 Mb or RAM.

b-mehta · 2018-12-08T00:46:41Z

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

rockbmb

I took a careful look, but may have let some things slip. Will double-check later today.

test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

rockbmb · 2018-12-08T18:16:27Z

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

+{-# SPECIALISE inverseTotient :: Semiring b => (Integer -> b) -> Integer -> b #-}
+{-# SPECIALISE inverseTotient :: Semiring b => (Natural -> b) -> Natural -> b #-}
+
+-- | The inverse 'sigma' 1 function such that


Excellent comments.

rockbmb · 2018-12-08T18:17:49Z

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

+  zero = MinWord maxBound
+  one  = MinWord 1
+  plus  (MinWord a) (MinWord b) = MinWord (a `min` b)
+  times (MinWord a) (MinWord b) = MinWord (a * b)


TIL this how we define one valid semiring for the min operator, I wasn't quite sure (*) would work.

Bodigrim · 2018-12-08T19:36:41Z

@rockbmb thanks for valuable suggestions.

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

rockbmb · 2018-12-08T21:25:18Z

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.

rockbmb · 2018-12-19T19:11:07Z

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

-newtype MinNatural = MinNatural { unMinNatural :: Maybe Natural }
-  deriving (Show)
+data MinNatural
+  = MinNatural { unMinNatural :: !Natural }


Right, I didn't give much importance to this before, but why was this Maybe before, and why change it now?

It seems to me that switching to a newly defined sum type allows the non-infinite case to be more strict, since the indirection provided by Just is removed

Operations on DirichletSeries never force values of the underlying Map: all multiplications and filtering is based solely on keys. So values of the Map inevitably accumulate huge thunks, unless a) Data.Map.Strict is used, which triggers weak head normal form (WHNF) evaluation of values, b) datatype for values is strict, meaning that WHNF conicides with a normal form. This is crucial for a decent performance.

Maybe a type is lazy: forcing WHNF will force only evaluation of Just constructor, but not the value of a. On contrary MinNatural is strict, thanks to bang.

b-mehta

Excellent work! Apologies for the delay.

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs

benchmark/Math/NumberTheory/InverseBench.hs

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

b-mehta · 2018-12-21T01:45:01Z

Math/NumberTheory/ArithmeticFunctions/Inverse.hs

+-- > x `elem` inverseTotient (totient x)
+--
+-- The return value is parametrized by a semiring, which allows
+-- various applications. E. g., list all preimages:


It might be nice to include examples of calculating the sum or sum of squares of preimages:

>>> getSum $ inverseTotient Sum 120 4904 >>> getSum $ inverseTotient (Sum . (^2)) 120 1573974

Actually Sum is redundant: inverseTotient id 120 works as well.

True! In that case, it might be nice to add a note saying that the semiring on Word is the one given by addition and multiplication, to explain why that works.

We can do even something like inverseTotient (sigma 1) 120 %)

b-mehta · 2018-12-21T20:30:11Z

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 · 2018-12-21T22:33:08Z

@rockbmb @b-mehta thanks for reviewing!

Bodigrim · 2018-12-22T00:20:34Z

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

rockbmb · 2018-12-22T17:18:29Z

@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 · 2018-12-22T17:26:38Z

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

rockbmb reviewed Oct 24, 2018

View reviewed changes

Bodigrim added 17 commits December 4, 2018 22:27

Implement inverseTotient

e907fa8

Use semirings instead of Applicative+Monoid

a5c9027

Add MaxValue/MinValue semirings

88bde95

Unit tests for inverse totient

a2521f2

Instead of representing any series focus on Dirichlet series with imp…

8a6a56e

…licit unit leading coefficient

Specialize functions for Integer

a876370

Pass upper bound to timesAndCrop explicitly

04a4fba

Represent Dirichlet series with maps

752bb52

Fix unions

054d3f9

Use strict Map

bdef230

Introduce PrimePowers datatype

20eed48

Implement inverse sigma 1

e24b4e1

More elaborated invSigma

769605d

Test sigma function

b3b9622

Unit tests for sigma function

bcad9ed

Optimize main algorithm

07f9047

More haddock comments

4c0c305

Bodigrim changed the title ~~[WIP] Implement inverseTotient~~ Implement inverseTotient Dec 6, 2018

Bodigrim added 2 commits December 7, 2018 22:03

Add benchmarks

5386ad0

Final polishing

27eda69

This was referenced Dec 8, 2018

Stackage chessai/semirings#27

Closed

Inverse the Jordan and the sum-of-powers-of-divisors functions #148

Closed

rockbmb approved these changes Dec 8, 2018

View reviewed changes

rockbmb reviewed Dec 8, 2018

View reviewed changes

Math/NumberTheory/ArithmeticFunctions/Inverse.hs Outdated Show resolved Hide resolved

Math/NumberTheory/ArithmeticFunctions/Inverse.hs Outdated Show resolved Hide resolved

Bodigrim added 2 commits December 9, 2018 00:06

Review suggestions

c3383ec

Make MinNatural strict

ceb94f1

rockbmb reviewed Dec 19, 2018

View reviewed changes

b-mehta approved these changes Dec 21, 2018

View reviewed changes

Bodigrim added 2 commits December 21, 2018 19:27

Review suggestions

076dfa3

Write a helper function asSetOfPreimages

b67bb41

Bodigrim merged commit 34af34c into Bodigrim:master Dec 21, 2018

Bodigrim deleted the inverse-totient branch December 21, 2018 22:32

Bodigrim mentioned this pull request Dec 22, 2018

Summatory functions #60

Open

7 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement inverseTotient #142

Implement inverseTotient #142

Bodigrim commented Sep 28, 2018

rockbmb commented Sep 29, 2018

Bodigrim commented Sep 29, 2018

rockbmb left a comment

rockbmb commented Oct 24, 2018

Bodigrim commented Dec 4, 2018 •

edited

Loading

Bodigrim commented Dec 6, 2018

rockbmb commented Dec 7, 2018 •

edited

Loading

rockbmb commented Dec 7, 2018

Bodigrim commented Dec 7, 2018

b-mehta commented Dec 8, 2018

rockbmb left a comment

rockbmb Dec 8, 2018

rockbmb Dec 8, 2018

Bodigrim commented Dec 8, 2018

rockbmb commented Dec 8, 2018

rockbmb Dec 19, 2018

b-mehta Dec 20, 2018

Bodigrim Dec 21, 2018

b-mehta left a comment •

edited

Loading

b-mehta Dec 21, 2018

Bodigrim Dec 21, 2018

b-mehta Dec 21, 2018

Bodigrim Dec 21, 2018

b-mehta commented Dec 21, 2018

Bodigrim commented Dec 21, 2018

Bodigrim commented Dec 22, 2018

rockbmb commented Dec 22, 2018

b-mehta commented Dec 22, 2018

Implement inverseTotient #142

Implement inverseTotient #142

Conversation

Bodigrim commented Sep 28, 2018

rockbmb commented Sep 29, 2018

Bodigrim commented Sep 29, 2018

rockbmb left a comment

Choose a reason for hiding this comment

rockbmb commented Oct 24, 2018

Bodigrim commented Dec 4, 2018 • edited Loading

Bodigrim commented Dec 6, 2018

rockbmb commented Dec 7, 2018 • edited Loading

rockbmb commented Dec 7, 2018

Bodigrim commented Dec 7, 2018

b-mehta commented Dec 8, 2018

rockbmb left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Bodigrim commented Dec 8, 2018

rockbmb commented Dec 8, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

b-mehta left a comment • 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

b-mehta commented Dec 21, 2018

Bodigrim commented Dec 21, 2018

Bodigrim commented Dec 22, 2018

rockbmb commented Dec 22, 2018

b-mehta commented Dec 22, 2018

Bodigrim commented Dec 4, 2018 •

edited

Loading

rockbmb commented Dec 7, 2018 •

edited

Loading

b-mehta left a comment •

edited

Loading