-- |
-- Module:      Math.NumberTheory.ArithmeticFunctions.Inverse
-- Copyright:   (c) 2018 Andrew Lelechenko
-- Licence:     MIT
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--
-- Computing inverses of multiplicative functions.
-- The implementation is based on
-- <https://www.emis.de/journals/JIS/VOL19/Alekseyev/alek5.pdf Computing the Inverses, their Power Sums, and Extrema for Euler’s Totient and Other Multiplicative Functions>
-- by M. A. Alekseyev.

{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Math.NumberTheory.ArithmeticFunctions.Inverse
  ( inverseTotient
  , inverseJordan
  , inverseSigma
  , inverseSigmaK
  , -- * Wrappers
    MinWord(..)
  , MaxWord(..)
  , MinNatural(..)
  , MaxNatural(..)
  , -- * Utils
    asSetOfPreimages
  ) where

import Prelude hiding (rem, quot)
import Data.Bits (Bits)
import Data.Euclidean
import Data.Foldable
import Data.List (partition, sortOn)
import Data.Map (Map)
import qualified Data.Map as M
import Data.Maybe
import Data.Ord (Down(..))
import Data.Semiring (Semiring(..), Mul(..))
import Data.Set (Set)
import qualified Data.Set as S
import Data.Traversable
import Numeric.Natural

import Math.NumberTheory.ArithmeticFunctions
import Math.NumberTheory.Logarithms
import Math.NumberTheory.Roots (exactRoot, integerRoot)
import Math.NumberTheory.Primes
import Math.NumberTheory.Utils.DirichletSeries (DirichletSeries)
import qualified Math.NumberTheory.Utils.DirichletSeries as DS
import Math.NumberTheory.Utils.FromIntegral

data PrimePowers a = PrimePowers
  { forall a. PrimePowers a -> Prime a
_ppPrime  :: Prime a
  , forall a. PrimePowers a -> [Word]
_ppPowers :: [Word] -- sorted list
  }

instance Show a => Show (PrimePowers a) where
  show :: PrimePowers a -> String
show (PrimePowers Prime a
p [Word]
xs) = String
"PP " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show (forall a. Prime a -> a
unPrime Prime a
p) forall a. [a] -> [a] -> [a]
++ String
" " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show [Word]
xs

-- | Convert a list of powers of a prime into an atomic Dirichlet series
-- (Section 4, Step 2).
atomicSeries
  :: Num a
  => (a -> b)               -- ^ How to inject a number into a semiring
  -> ArithmeticFunction a c -- ^ Arithmetic function, which we aim to inverse
  -> PrimePowers a          -- ^ List of powers of a prime
  -> DirichletSeries c b    -- ^ Atomic Dirichlet series
atomicSeries :: forall a b c.
Num a =>
(a -> b)
-> ArithmeticFunction a c -> PrimePowers a -> DirichletSeries c b
atomicSeries a -> b
point (ArithmeticFunction Prime a -> Word -> m
f m -> c
g) (PrimePowers Prime a
p [Word]
ks) =
  forall a b. [(a, b)] -> DirichletSeries a b
DS.fromDistinctAscList (forall a b. (a -> b) -> [a] -> [b]
map (\Word
k -> (m -> c
g (Prime a -> Word -> m
f Prime a
p Word
k), a -> b
point (forall a. Prime a -> a
unPrime Prime a
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
k))) [Word]
ks)

-- | See section 5.1 of the paper.
invJordan
  :: forall a. (Integral a, UniqueFactorisation a, Eq a)
  => Word
  -- ^ Value of k in 'jordan' k
  -> [(Prime a, Word)]
  -- ^ Factorisation of a value of the totient function
  -> [PrimePowers a]
  -- ^ Possible prime factors of an argument of the totient function
invJordan :: forall a.
(Integral a, UniqueFactorisation a, Eq a) =>
Word -> [(Prime a, Word)] -> [PrimePowers a]
invJordan Word
k [(Prime a, Word)]
fs = forall a b. (a -> b) -> [a] -> [b]
map (\Prime a
p -> forall a. Prime a -> [Word] -> PrimePowers a
PrimePowers Prime a
p (Prime a -> [Word]
doPrime Prime a
p)) [Prime a]
ps
  where
    divs :: [a]
    divs :: [a]
divs = forall n a. ArithmeticFunction n a -> [(Prime n, Word)] -> a
runFunctionOnFactors forall n. Num n => ArithmeticFunction n [n]
divisorsListA [(Prime a, Word)]
fs

    ps :: [Prime a]
    ps :: [Prime a]
ps = forall a b. (a -> Maybe b) -> [a] -> [b]
mapMaybe (\a
d -> forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Word
k (a
d forall a. Num a => a -> a -> a
+ a
1) forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= forall a. UniqueFactorisation a => a -> Maybe (Prime a)
isPrime) [a]
divs

    doPrime :: Prime a -> [Word]
    doPrime :: Prime a -> [Word]
doPrime Prime a
p = case forall a b. Eq a => a -> [(a, b)] -> Maybe b
lookup Prime a
p [(Prime a, Word)]
fs of
      Maybe Word
Nothing -> [Word
1]
      Just Word
w  -> [Word
1 .. Word
wforall a. Num a => a -> a -> a
+Word
1]

-- | See section 5.2 of the paper.
invSigma
  :: forall a. (Euclidean a, Integral a, UniqueFactorisation a, Enum (Prime a), Bits a)
  => Word
  -- ^ Value of k in 'sigma' k
  -> [(Prime a, Word)]
  -- ^ Factorisation of a value of the sum-of-divisors function
  -> [PrimePowers a]
  -- ^ Possible prime factors of an argument of the sum-of-divisors function
invSigma :: forall a.
(Euclidean a, Integral a, UniqueFactorisation a, Enum (Prime a),
 Bits a) =>
Word -> [(Prime a, Word)] -> [PrimePowers a]
invSigma Word
k [(Prime a, Word)]
fs
  = forall a b. (a -> b) -> [a] -> [b]
map (\(Prime a
x, Set Word
ys) -> forall a. Prime a -> [Word] -> PrimePowers a
PrimePowers Prime a
x (forall a. Set a -> [a]
S.toList Set Word
ys))
  forall a b. (a -> b) -> a -> b
$ forall k a. Map k a -> [(k, a)]
M.assocs
  forall a b. (a -> b) -> a -> b
$ forall k a. Ord k => (a -> a -> a) -> Map k a -> Map k a -> Map k a
M.unionWith forall a. Semigroup a => a -> a -> a
(<>) Map (Prime a) (Set Word)
pksSmall Map (Prime a) (Set Word)
pksLarge
  where
    numDivs :: a
    numDivs :: a
numDivs = forall n a. ArithmeticFunction n a -> [(Prime n, Word)] -> a
runFunctionOnFactors forall a n. Num a => ArithmeticFunction n a
tauA [(Prime a, Word)]
fs

    divs :: [a]
    divs :: [a]
divs = forall n a. ArithmeticFunction n a -> [(Prime n, Word)] -> a
runFunctionOnFactors forall n. Num n => ArithmeticFunction n [n]
divisorsListA [(Prime a, Word)]
fs

    n :: a
    n :: a
n = forall a. Num a => [(Prime a, Word)] -> a
factorBack [(Prime a, Word)]
fs

    -- There are two possible strategies to find possible prime factors
    -- of an argument of the sum-of-divisors function.
    -- 1. Take each prime p and each power e such that p^e <= n,
    -- compute sigma(p^e) and check whether it is a divisor of n.
    -- (corresponds to 'pksSmall' below)
    -- 2. Take each divisor d of n and each power e such that e <= log_2 d,
    -- compute p = floor(e-th root of (d - 1)) and check whether sigma(p^e) = d
    -- and p is actually prime (correposnds to 'pksLarge' below).
    --
    -- Asymptotically the second strategy is beneficial, but computing
    -- exact e-th roots of huge integers (especially when they exceed MAX_DOUBLE)
    -- is very costly. That is why we employ the first strategy for primes
    -- below limit 'lim' and the second one for larger ones. This allows us
    -- to loop over e <= log_lim d which is much smaller than log_2 d.
    --
    -- The value of 'lim' below was chosen heuristically;
    -- it may be tuned in future in accordance with new experimental data.
    lim :: a
    lim :: a
lim = a
numDivs forall a. Ord a => a -> a -> a
`max` a
2

    pksSmall :: Map (Prime a) (Set Word)
    pksSmall :: Map (Prime a) (Set Word)
pksSmall = forall k a. [(k, a)] -> Map k a
M.fromDistinctAscList
      [ (Prime a
p, Set Word
pows)
      | Prime a
p <- [forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
nextPrime a
2 .. forall a.
(Bits a, Integral a, UniqueFactorisation a) =>
a -> Prime a
precPrime a
lim]
      , let pows :: Set Word
pows = Prime a -> Set Word
doPrime Prime a
p
      , Bool -> Bool
not (forall (t :: * -> *) a. Foldable t => t a -> Bool
null Set Word
pows)
      ]

    doPrime :: Prime a -> Set Word
    doPrime :: Prime a -> Set Word
doPrime Prime a
p' = let p :: a
p = forall a. Prime a -> a
unPrime Prime a
p' in forall a. [a] -> Set a
S.fromDistinctAscList
      [ Word
e
      | Word
e <- [Word
1 .. Int -> Word
intToWord (Integer -> Integer -> Int
integerLogBase (forall a. Integral a => a -> Integer
toInteger (a
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
k)) (forall a. Integral a => a -> Integer
toInteger a
n))]
      , a
n forall a. Euclidean a => a -> a -> a
`rem` ((a
p forall a b. (Num a, Integral b) => a -> b -> a
^ (Word
k forall a. Num a => a -> a -> a
* (Word
e forall a. Num a => a -> a -> a
+ Word
1)) forall a. Num a => a -> a -> a
- a
1) forall a. Euclidean a => a -> a -> a
`quot` (a
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
k forall a. Num a => a -> a -> a
- a
1)) forall a. Eq a => a -> a -> Bool
== a
0
      ]

    pksLarge :: Map (Prime a) (Set Word)
    pksLarge :: Map (Prime a) (Set Word)
pksLarge = forall (f :: * -> *) k a.
(Foldable f, Ord k) =>
(a -> a -> a) -> f (Map k a) -> Map k a
M.unionsWith forall a. Semigroup a => a -> a -> a
(<>)
      [ forall b a. b -> (a -> b) -> Maybe a -> b
maybe forall a. Monoid a => a
mempty (forall k a. k -> a -> Map k a
`M.singleton` forall a. a -> Set a
S.singleton Word
e) (forall a. UniqueFactorisation a => a -> Maybe (Prime a)
isPrime a
p)
      | a
d <- [a]
divs
      , Word
e <- [Word
1 .. Int -> Word
intToWord (forall a. Euclidean a => a -> a -> a
quot (Integer -> Integer -> Int
integerLogBase (forall a. Integral a => a -> Integer
toInteger a
lim) (forall a. Integral a => a -> Integer
toInteger a
d)) (Word -> Int
wordToInt Word
k)) ]
      , let p :: a
p = forall a b. (Integral a, Integral b) => b -> a -> a
integerRoot (Word
e forall a. Num a => a -> a -> a
* Word
k) (a
d forall a. Num a => a -> a -> a
- a
1)
      , a
p forall a b. (Num a, Integral b) => a -> b -> a
^ (Word
k forall a. Num a => a -> a -> a
* (Word
e forall a. Num a => a -> a -> a
+ Word
1)) forall a. Num a => a -> a -> a
- a
1 forall a. Eq a => a -> a -> Bool
== a
d forall a. Num a => a -> a -> a
* (a
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
k forall a. Num a => a -> a -> a
- a
1)
      ]

-- | Instead of multiplying all atomic series and filtering out everything,
-- which is not divisible by @n@, we'd rather split all atomic series into
-- a couple of batches, each of which corresponds to a prime factor of @n@.
-- This allows us to crop resulting Dirichlet series (see 'filter' calls
-- in @invertFunction@ below) at the end of each batch, saving time and memory.
strategy
  :: forall a c. (GcdDomain c, Ord c)
  => ArithmeticFunction a c
  -- ^ Arithmetic function, which we aim to inverse
  -> [(Prime c, Word)]
  -- ^ Factorisation of a value of the arithmetic function
  -> [PrimePowers a]
  -- ^ Possible prime factors of an argument of the arithmetic function
  -> [(Maybe (Prime c, Word), [PrimePowers a])]
  -- ^ Batches, corresponding to each element of the factorisation of a value
strategy :: forall a c.
(GcdDomain c, Ord c) =>
ArithmeticFunction a c
-> [(Prime c, Word)]
-> [PrimePowers a]
-> [(Maybe (Prime c, Word), [PrimePowers a])]
strategy (ArithmeticFunction Prime a -> Word -> m
f m -> c
g) [(Prime c, Word)]
factors [PrimePowers a]
args = (forall a. Maybe a
Nothing, [PrimePowers a]
ret) forall a. a -> [a] -> [a]
: [(Maybe (Prime c, Word), [PrimePowers a])]
rets
  where
    ([PrimePowers a]
ret, [(Maybe (Prime c, Word), [PrimePowers a])]
rets)
      = forall (t :: * -> *) s a b.
Traversable t =>
(s -> a -> (s, b)) -> s -> t a -> (s, t b)
mapAccumL [PrimePowers a]
-> (Prime c, Word)
-> ([PrimePowers a], (Maybe (Prime c, Word), [PrimePowers a]))
go [PrimePowers a]
args
      forall a b. (a -> b) -> a -> b
$ forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (forall a. a -> Down a
Down forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fst) [(Prime c, Word)]
factors

    go
      :: [PrimePowers a]
      -> (Prime c, Word)
      -> ([PrimePowers a], (Maybe (Prime c, Word), [PrimePowers a]))
    go :: [PrimePowers a]
-> (Prime c, Word)
-> ([PrimePowers a], (Maybe (Prime c, Word), [PrimePowers a]))
go [PrimePowers a]
ts (Prime c
p, Word
k) = ([PrimePowers a]
rs, (forall a. a -> Maybe a
Just (Prime c
p, Word
k), [PrimePowers a]
qs))
      where
        predicate :: PrimePowers a -> Bool
predicate (PrimePowers Prime a
q [Word]
ls) = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Word
l -> forall a. Maybe a -> Bool
isJust forall a b. (a -> b) -> a -> b
$ m -> c
g (Prime a -> Word -> m
f Prime a
q Word
l) forall a. GcdDomain a => a -> a -> Maybe a
`divide` forall a. Prime a -> a
unPrime Prime c
p) [Word]
ls
        ([PrimePowers a]
qs, [PrimePowers a]
rs) = forall a. (a -> Bool) -> [a] -> ([a], [a])
partition PrimePowers a -> Bool
predicate [PrimePowers a]
ts

-- | Main workhorse.
invertFunction
  :: forall a b c.
     (Num a, Semiring b, Euclidean c, UniqueFactorisation c, Ord c)
  => (a -> b)
  -- ^ How to inject a number into a semiring
  -> ArithmeticFunction a c
  -- ^ Arithmetic function, which we aim to inverse
  -> ([(Prime c, Word)] -> [PrimePowers a])
  -- ^ How to find possible prime factors of the argument
  -> c
  -- ^ Value of the arithmetic function, which we aim to inverse
  -> b
  -- ^ Semiring element, representing preimages
invertFunction :: forall a b c.
(Num a, Semiring b, Euclidean c, UniqueFactorisation c, Ord c) =>
(a -> b)
-> ArithmeticFunction a c
-> ([(Prime c, Word)] -> [PrimePowers a])
-> c
-> b
invertFunction a -> b
point ArithmeticFunction a c
f [(Prime c, Word)] -> [PrimePowers a]
invF c
n
  = forall a b.
(Ord a, Num a, Semiring b) =>
a -> DirichletSeries a b -> b
DS.lookup c
n
  forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (forall a b c. (a -> b -> c) -> b -> a -> c
flip (forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Maybe (Prime c, Word)
-> [PrimePowers a] -> DirichletSeries c b -> DirichletSeries c b
processBatch)) (forall a b. [(a, b)] -> DirichletSeries a b
DS.fromDistinctAscList []) [(Maybe (Prime c, Word), [PrimePowers a])]
batches
  where
    factors :: [(Prime c, Word)]
factors = forall a. UniqueFactorisation a => a -> [(Prime a, Word)]
factorise c
n
    batches :: [(Maybe (Prime c, Word), [PrimePowers a])]
batches = forall a c.
(GcdDomain c, Ord c) =>
ArithmeticFunction a c
-> [(Prime c, Word)]
-> [PrimePowers a]
-> [(Maybe (Prime c, Word), [PrimePowers a])]
strategy ArithmeticFunction a c
f [(Prime c, Word)]
factors forall a b. (a -> b) -> a -> b
$ [(Prime c, Word)] -> [PrimePowers a]
invF [(Prime c, Word)]
factors

    processBatch
      :: Maybe (Prime c, Word)
      -> [PrimePowers a]
      -> DirichletSeries c b
      -> DirichletSeries c b
    processBatch :: Maybe (Prime c, Word)
-> [PrimePowers a] -> DirichletSeries c b -> DirichletSeries c b
processBatch Maybe (Prime c, Word)
Nothing [PrimePowers a]
xs DirichletSeries c b
acc
      = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (forall a b.
(Num a, Euclidean a, Ord a, Semiring b) =>
a
-> DirichletSeries a b
-> DirichletSeries a b
-> DirichletSeries a b
DS.timesAndCrop c
n) DirichletSeries c b
acc
      forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c.
Num a =>
(a -> b)
-> ArithmeticFunction a c -> PrimePowers a -> DirichletSeries c b
atomicSeries a -> b
point ArithmeticFunction a c
f) [PrimePowers a]
xs

    -- This is equivalent to the next, general case, but is faster,
    -- since it avoids construction of many intermediate series.
    processBatch (Just (Prime c
p, Word
1)) [PrimePowers a]
xs DirichletSeries c b
acc
      = forall a b.
(a -> Bool) -> DirichletSeries a b -> DirichletSeries a b
DS.filter (\c
a -> c
a forall a. Euclidean a => a -> a -> a
`rem` forall a. Prime a -> a
unPrime Prime c
p forall a. Eq a => a -> a -> Bool
== c
0)
      forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (forall a b.
(Num a, Euclidean a, Ord a, Semiring b) =>
a
-> DirichletSeries a b
-> DirichletSeries a b
-> DirichletSeries a b
DS.timesAndCrop c
n) DirichletSeries c b
acc'
      forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c.
Num a =>
(a -> b)
-> ArithmeticFunction a c -> PrimePowers a -> DirichletSeries c b
atomicSeries a -> b
point ArithmeticFunction a c
f) [PrimePowers a]
xs2
      where
        ([PrimePowers a]
xs1, [PrimePowers a]
xs2) = forall a. (a -> Bool) -> [a] -> ([a], [a])
partition (\(PrimePowers Prime a
_ [Word]
ts) -> forall (t :: * -> *) a. Foldable t => t a -> Int
length [Word]
ts forall a. Eq a => a -> a -> Bool
== Int
1) [PrimePowers a]
xs
        vs :: DirichletSeries c b
vs = forall a b.
(Ord a, Semiring b) =>
[DirichletSeries a b] -> DirichletSeries a b
DS.unions forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c.
Num a =>
(a -> b)
-> ArithmeticFunction a c -> PrimePowers a -> DirichletSeries c b
atomicSeries a -> b
point ArithmeticFunction a c
f) [PrimePowers a]
xs1
        (DirichletSeries c b
ys, DirichletSeries c b
zs) = forall a b.
(a -> Bool)
-> DirichletSeries a b
-> (DirichletSeries a b, DirichletSeries a b)
DS.partition (\c
a -> c
a forall a. Euclidean a => a -> a -> a
`rem` forall a. Prime a -> a
unPrime Prime c
p forall a. Eq a => a -> a -> Bool
== c
0) DirichletSeries c b
acc
        acc' :: DirichletSeries c b
acc' = DirichletSeries c b
ys forall a b.
(Ord a, Semiring b) =>
DirichletSeries a b -> DirichletSeries a b -> DirichletSeries a b
`DS.union` forall a b.
(Num a, Euclidean a, Ord a, Semiring b) =>
a
-> DirichletSeries a b
-> DirichletSeries a b
-> DirichletSeries a b
DS.timesAndCrop c
n DirichletSeries c b
zs DirichletSeries c b
vs

    processBatch (Just (Prime c, Word)
pk) [PrimePowers a]
xs DirichletSeries c b
acc
      = (\(Prime c
p, Word
k) -> forall a b.
(a -> Bool) -> DirichletSeries a b -> DirichletSeries a b
DS.filter (\c
a -> c
a forall a. Euclidean a => a -> a -> a
`rem` (forall a. Prime a -> a
unPrime Prime c
p forall a b. (Num a, Integral b) => a -> b -> a
^ Word
k) forall a. Eq a => a -> a -> Bool
== c
0)) (Prime c, Word)
pk
      forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (forall a b.
(Num a, Euclidean a, Ord a, Semiring b) =>
a
-> DirichletSeries a b
-> DirichletSeries a b
-> DirichletSeries a b
DS.timesAndCrop c
n) DirichletSeries c b
acc
      forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall a b c.
Num a =>
(a -> b)
-> ArithmeticFunction a c -> PrimePowers a -> DirichletSeries c b
atomicSeries a -> b
point ArithmeticFunction a c
f) [PrimePowers a]
xs

{-# SPECIALISE invertFunction :: Semiring b => (Int -> b) -> ArithmeticFunction Int Int -> ([(Prime Int, Word)] -> [PrimePowers Int]) -> Int -> b #-}
{-# SPECIALISE invertFunction :: Semiring b => (Word -> b) -> ArithmeticFunction Word Word -> ([(Prime Word, Word)] -> [PrimePowers Word]) -> Word -> b #-}
{-# SPECIALISE invertFunction :: Semiring b => (Integer -> b) -> ArithmeticFunction Integer Integer -> ([(Prime Integer, Word)] -> [PrimePowers Integer]) -> Integer -> b #-}
{-# SPECIALISE invertFunction :: Semiring b => (Natural -> b) -> ArithmeticFunction Natural Natural -> ([(Prime Natural, Word)] -> [PrimePowers Natural]) -> Natural -> b #-}

-- | The inverse for 'totient' function.
--
-- The return value is parameterized by a 'Semiring', which allows
-- various applications by providing different (multiplicative) embeddings.
-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
--
-- >>> import qualified Data.Set as S
-- >>> import Data.Semigroup
-- >>> S.mapMonotonic getProduct (inverseTotient (S.singleton . Product) 120)
-- fromList [143,155,175,183,225,231,244,248,286,308,310,350,366,372,396,450,462]
--
-- Count preimages:
--
-- >>> inverseTotient (const 1) 120
-- 17
--
-- Sum preimages:
--
-- >>> inverseTotient id 120
-- 4904
--
-- Find minimal and maximal preimages:
--
-- >>> unMinWord (inverseTotient MinWord 120)
-- 143
-- >>> unMaxWord (inverseTotient MaxWord 120)
-- 462
inverseTotient
  :: (Semiring b, Integral a, Euclidean a, UniqueFactorisation a)
  => (a -> b)
  -> a
  -> b
inverseTotient :: forall b a.
(Semiring b, Integral a, Euclidean a, UniqueFactorisation a) =>
(a -> b) -> a -> b
inverseTotient = forall b a.
(Semiring b, Integral a, Euclidean a, UniqueFactorisation a) =>
Word -> (a -> b) -> a -> b
inverseJordan Word
1
{-# SPECIALISE inverseTotient :: Semiring b => (Int -> b) -> Int -> b #-}
{-# SPECIALISE inverseTotient :: Semiring b => (Word -> b) -> Word -> b #-}
{-# SPECIALISE inverseTotient :: Semiring b => (Integer -> b) -> Integer -> b #-}
{-# SPECIALISE inverseTotient :: Semiring b => (Natural -> b) -> Natural -> b #-}

-- | The inverse for 'jordan' function.
--
-- Generalizes the 'inverseTotient' function, which is 'inverseJordan' 1.
--
-- The return value is parameterized by a 'Semiring', which allows
-- various applications by providing different (multiplicative) embeddings.
-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
--
-- >>> import qualified Data.Set as S
-- >>> import Data.Semigroup
-- >>> S.mapMonotonic getProduct (inverseJordan 2 (S.singleton . Product) 192)
-- fromList [15,16]
--
-- Similarly to 'inverseTotient', it is possible to count and sum preimages, or
-- get the maximum/minimum preimage.
--
-- Note: it is the __user's responsibility__ to use an appropriate type for
-- 'inverseSigmaK'. Even low values of k with 'Int' or 'Word' will return
-- invalid results due to over/underflow, or throw exceptions (i.e. division by
-- zero).
--
-- >>> asSetOfPreimages (inverseJordan 15) (jordan 15 19) :: S.Set Int
-- fromList []
--
-- >>> asSetOfPreimages (inverseJordan 15) (jordan 15 19) :: S.Set Integer
-- fromList [19]
inverseJordan
  :: (Semiring b, Integral a, Euclidean a, UniqueFactorisation a)
  => Word
  -> (a -> b)
  -> a
  -> b
inverseJordan :: forall b a.
(Semiring b, Integral a, Euclidean a, UniqueFactorisation a) =>
Word -> (a -> b) -> a -> b
inverseJordan Word
k a -> b
point = forall a b c.
(Num a, Semiring b, Euclidean c, UniqueFactorisation c, Ord c) =>
(a -> b)
-> ArithmeticFunction a c
-> ([(Prime c, Word)] -> [PrimePowers a])
-> c
-> b
invertFunction a -> b
point (forall n. Num n => Word -> ArithmeticFunction n n
jordanA Word
k) (forall a.
(Integral a, UniqueFactorisation a, Eq a) =>
Word -> [(Prime a, Word)] -> [PrimePowers a]
invJordan Word
k)
{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Int -> b) -> Int -> b #-}
{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Word -> b) -> Word -> b #-}
{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Integer -> b) -> Integer -> b #-}
{-# SPECIALISE inverseJordan :: Semiring b => Word -> (Natural -> b) -> Natural -> b #-}

-- | The inverse for 'sigma' 1 function.
--
-- The return value is parameterized by a 'Semiring', which allows
-- various applications by providing different (multiplicative) embeddings.
-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
--
-- >>> import qualified Data.Set as S
-- >>> import Data.Semigroup
-- >>> :set -XFlexibleContexts
-- >>> S.mapMonotonic getProduct (inverseSigma (S.singleton . Product) 120)
-- fromList [54,56,87,95]
--
-- Count preimages:
--
-- >>> inverseSigma (const 1) 120
-- 4
--
-- Sum preimages:
--
-- >>> inverseSigma id 120
-- 292
--
-- Find minimal and maximal preimages:
--
-- >>> unMinWord (inverseSigma MinWord 120)
-- 54
-- >>> unMaxWord (inverseSigma MaxWord 120)
-- 95
inverseSigma
  :: (Semiring b, Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a)
  => (a -> b)
  -> a
  -> b
inverseSigma :: forall b a.
(Semiring b, Euclidean a, UniqueFactorisation a, Integral a,
 Enum (Prime a), Bits a) =>
(a -> b) -> a -> b
inverseSigma = forall b a.
(Semiring b, Euclidean a, UniqueFactorisation a, Integral a,
 Enum (Prime a), Bits a) =>
Word -> (a -> b) -> a -> b
inverseSigmaK Word
1
{-# SPECIALISE inverseSigma :: Semiring b => (Int -> b) -> Int -> b #-}
{-# SPECIALISE inverseSigma :: Semiring b => (Word -> b) -> Word -> b #-}
{-# SPECIALISE inverseSigma :: Semiring b => (Integer -> b) -> Integer -> b #-}
{-# SPECIALISE inverseSigma :: Semiring b => (Natural -> b) -> Natural -> b #-}

-- | The inverse for 'sigma' function.
--
-- Generalizes the 'inverseSigma' function, which is 'inverseSigmaK' 1.
--
-- The return value is parameterized by a 'Semiring', which allows
-- various applications by providing different (multiplicative) embeddings.
-- E. g., list all preimages (see a helper 'asSetOfPreimages'):
--
-- >>> import qualified Data.Set as S
-- >>> import Data.Semigroup
-- >>> :set -XFlexibleContexts
-- >>> S.mapMonotonic getProduct (inverseSigmaK 2 (S.singleton . Product) 850)
-- fromList [24,26]
--
-- Similarly to 'inverseSigma', it is possible to count and sum preimages, or
-- get the maximum/minimum preimage.
--
-- Note: it is the __user's responsibility__ to use an appropriate type for
-- 'inverseSigmaK'. Even low values of k with 'Int' or 'Word' will return
-- invalid results due to over/underflow, or throw exceptions (i.e. division by
-- zero).
--
-- >>> asSetOfPreimages (inverseSigmaK 17) (sigma 17 13) :: S.Set Int
-- fromList *** Exception: divide by zero
inverseSigmaK
  :: (Semiring b, Euclidean a, UniqueFactorisation a, Integral a, Enum (Prime a), Bits a)
  => Word
  -> (a -> b)
  -> a
  -> b
inverseSigmaK :: forall b a.
(Semiring b, Euclidean a, UniqueFactorisation a, Integral a,
 Enum (Prime a), Bits a) =>
Word -> (a -> b) -> a -> b
inverseSigmaK Word
k a -> b
point = forall a b c.
(Num a, Semiring b, Euclidean c, UniqueFactorisation c, Ord c) =>
(a -> b)
-> ArithmeticFunction a c
-> ([(Prime c, Word)] -> [PrimePowers a])
-> c
-> b
invertFunction a -> b
point (forall n a.
(Integral n, Num a, GcdDomain a) =>
Word -> ArithmeticFunction n a
sigmaA Word
k) (forall a.
(Euclidean a, Integral a, UniqueFactorisation a, Enum (Prime a),
 Bits a) =>
Word -> [(Prime a, Word)] -> [PrimePowers a]
invSigma Word
k)
{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Int -> b) -> Int -> b #-}
{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Word -> b) -> Word -> b #-}
{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Integer -> b) -> Integer -> b #-}
{-# SPECIALISE inverseSigmaK :: Semiring b => Word -> (Natural -> b) -> Natural -> b #-}

--------------------------------------------------------------------------------
-- Wrappers

-- | Wrapper to use in conjunction with 'inverseTotient' and 'inverseSigma'.
-- Extracts the maximal preimage of function.
newtype MaxWord = MaxWord { MaxWord -> Word
unMaxWord :: Word }
  deriving (MaxWord -> MaxWord -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MaxWord -> MaxWord -> Bool
$c/= :: MaxWord -> MaxWord -> Bool
== :: MaxWord -> MaxWord -> Bool
$c== :: MaxWord -> MaxWord -> Bool
Eq, Eq MaxWord
MaxWord -> MaxWord -> Bool
MaxWord -> MaxWord -> Ordering
MaxWord -> MaxWord -> MaxWord
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: MaxWord -> MaxWord -> MaxWord
$cmin :: MaxWord -> MaxWord -> MaxWord
max :: MaxWord -> MaxWord -> MaxWord
$cmax :: MaxWord -> MaxWord -> MaxWord
>= :: MaxWord -> MaxWord -> Bool
$c>= :: MaxWord -> MaxWord -> Bool
> :: MaxWord -> MaxWord -> Bool
$c> :: MaxWord -> MaxWord -> Bool
<= :: MaxWord -> MaxWord -> Bool
$c<= :: MaxWord -> MaxWord -> Bool
< :: MaxWord -> MaxWord -> Bool
$c< :: MaxWord -> MaxWord -> Bool
compare :: MaxWord -> MaxWord -> Ordering
$ccompare :: MaxWord -> MaxWord -> Ordering
Ord, Int -> MaxWord -> ShowS
[MaxWord] -> ShowS
MaxWord -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MaxWord] -> ShowS
$cshowList :: [MaxWord] -> ShowS
show :: MaxWord -> String
$cshow :: MaxWord -> String
showsPrec :: Int -> MaxWord -> ShowS
$cshowsPrec :: Int -> MaxWord -> ShowS
Show)

instance Semiring MaxWord where
  zero :: MaxWord
zero = Word -> MaxWord
MaxWord forall a. Bounded a => a
minBound
  one :: MaxWord
one  = Word -> MaxWord
MaxWord Word
1
  plus :: MaxWord -> MaxWord -> MaxWord
plus  (MaxWord Word
a) (MaxWord Word
b) = Word -> MaxWord
MaxWord (Word
a forall a. Ord a => a -> a -> a
`max` Word
b)
  times :: MaxWord -> MaxWord -> MaxWord
times (MaxWord Word
a) (MaxWord Word
b) = Word -> MaxWord
MaxWord (Word
a forall a. Num a => a -> a -> a
* Word
b)

-- | Wrapper to use in conjunction with 'inverseTotient' and 'inverseSigma'.
-- Extracts the minimal preimage of function.
newtype MinWord = MinWord { MinWord -> Word
unMinWord :: Word }
  deriving (MinWord -> MinWord -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MinWord -> MinWord -> Bool
$c/= :: MinWord -> MinWord -> Bool
== :: MinWord -> MinWord -> Bool
$c== :: MinWord -> MinWord -> Bool
Eq, Eq MinWord
MinWord -> MinWord -> Bool
MinWord -> MinWord -> Ordering
MinWord -> MinWord -> MinWord
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: MinWord -> MinWord -> MinWord
$cmin :: MinWord -> MinWord -> MinWord
max :: MinWord -> MinWord -> MinWord
$cmax :: MinWord -> MinWord -> MinWord
>= :: MinWord -> MinWord -> Bool
$c>= :: MinWord -> MinWord -> Bool
> :: MinWord -> MinWord -> Bool
$c> :: MinWord -> MinWord -> Bool
<= :: MinWord -> MinWord -> Bool
$c<= :: MinWord -> MinWord -> Bool
< :: MinWord -> MinWord -> Bool
$c< :: MinWord -> MinWord -> Bool
compare :: MinWord -> MinWord -> Ordering
$ccompare :: MinWord -> MinWord -> Ordering
Ord, Int -> MinWord -> ShowS
[MinWord] -> ShowS
MinWord -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MinWord] -> ShowS
$cshowList :: [MinWord] -> ShowS
show :: MinWord -> String
$cshow :: MinWord -> String
showsPrec :: Int -> MinWord -> ShowS
$cshowsPrec :: Int -> MinWord -> ShowS
Show)

instance Semiring MinWord where
  zero :: MinWord
zero = Word -> MinWord
MinWord forall a. Bounded a => a
maxBound
  one :: MinWord
one  = Word -> MinWord
MinWord Word
1
  plus :: MinWord -> MinWord -> MinWord
plus  (MinWord Word
a) (MinWord Word
b) = Word -> MinWord
MinWord (Word
a forall a. Ord a => a -> a -> a
`min` Word
b)
  times :: MinWord -> MinWord -> MinWord
times (MinWord Word
a) (MinWord Word
b) = Word -> MinWord
MinWord (Word
a forall a. Num a => a -> a -> a
* Word
b)

-- | Wrapper to use in conjunction with 'inverseTotient' and 'inverseSigma'.
-- Extracts the maximal preimage of function.
newtype MaxNatural = MaxNatural { MaxNatural -> Natural
unMaxNatural :: Natural }
  deriving (MaxNatural -> MaxNatural -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MaxNatural -> MaxNatural -> Bool
$c/= :: MaxNatural -> MaxNatural -> Bool
== :: MaxNatural -> MaxNatural -> Bool
$c== :: MaxNatural -> MaxNatural -> Bool
Eq, Eq MaxNatural
MaxNatural -> MaxNatural -> Bool
MaxNatural -> MaxNatural -> Ordering
MaxNatural -> MaxNatural -> MaxNatural
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: MaxNatural -> MaxNatural -> MaxNatural
$cmin :: MaxNatural -> MaxNatural -> MaxNatural
max :: MaxNatural -> MaxNatural -> MaxNatural
$cmax :: MaxNatural -> MaxNatural -> MaxNatural
>= :: MaxNatural -> MaxNatural -> Bool
$c>= :: MaxNatural -> MaxNatural -> Bool
> :: MaxNatural -> MaxNatural -> Bool
$c> :: MaxNatural -> MaxNatural -> Bool
<= :: MaxNatural -> MaxNatural -> Bool
$c<= :: MaxNatural -> MaxNatural -> Bool
< :: MaxNatural -> MaxNatural -> Bool
$c< :: MaxNatural -> MaxNatural -> Bool
compare :: MaxNatural -> MaxNatural -> Ordering
$ccompare :: MaxNatural -> MaxNatural -> Ordering
Ord, Int -> MaxNatural -> ShowS
[MaxNatural] -> ShowS
MaxNatural -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MaxNatural] -> ShowS
$cshowList :: [MaxNatural] -> ShowS
show :: MaxNatural -> String
$cshow :: MaxNatural -> String
showsPrec :: Int -> MaxNatural -> ShowS
$cshowsPrec :: Int -> MaxNatural -> ShowS
Show)

instance Semiring MaxNatural where
  zero :: MaxNatural
zero = Natural -> MaxNatural
MaxNatural Natural
0
  one :: MaxNatural
one  = Natural -> MaxNatural
MaxNatural Natural
1
  plus :: MaxNatural -> MaxNatural -> MaxNatural
plus  (MaxNatural Natural
a) (MaxNatural Natural
b) = Natural -> MaxNatural
MaxNatural (Natural
a forall a. Ord a => a -> a -> a
`max` Natural
b)
  times :: MaxNatural -> MaxNatural -> MaxNatural
times (MaxNatural Natural
a) (MaxNatural Natural
b) = Natural -> MaxNatural
MaxNatural (Natural
a forall a. Num a => a -> a -> a
* Natural
b)

-- | Wrapper to use in conjunction with 'inverseTotient' and 'inverseSigma'.
-- Extracts the minimal preimage of function.
data MinNatural
  = MinNatural { MinNatural -> Natural
unMinNatural :: !Natural }
  | Infinity
  deriving (MinNatural -> MinNatural -> Bool
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
/= :: MinNatural -> MinNatural -> Bool
$c/= :: MinNatural -> MinNatural -> Bool
== :: MinNatural -> MinNatural -> Bool
$c== :: MinNatural -> MinNatural -> Bool
Eq, Eq MinNatural
MinNatural -> MinNatural -> Bool
MinNatural -> MinNatural -> Ordering
MinNatural -> MinNatural -> MinNatural
forall a.
Eq a
-> (a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
min :: MinNatural -> MinNatural -> MinNatural
$cmin :: MinNatural -> MinNatural -> MinNatural
max :: MinNatural -> MinNatural -> MinNatural
$cmax :: MinNatural -> MinNatural -> MinNatural
>= :: MinNatural -> MinNatural -> Bool
$c>= :: MinNatural -> MinNatural -> Bool
> :: MinNatural -> MinNatural -> Bool
$c> :: MinNatural -> MinNatural -> Bool
<= :: MinNatural -> MinNatural -> Bool
$c<= :: MinNatural -> MinNatural -> Bool
< :: MinNatural -> MinNatural -> Bool
$c< :: MinNatural -> MinNatural -> Bool
compare :: MinNatural -> MinNatural -> Ordering
$ccompare :: MinNatural -> MinNatural -> Ordering
Ord, Int -> MinNatural -> ShowS
[MinNatural] -> ShowS
MinNatural -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [MinNatural] -> ShowS
$cshowList :: [MinNatural] -> ShowS
show :: MinNatural -> String
$cshow :: MinNatural -> String
showsPrec :: Int -> MinNatural -> ShowS
$cshowsPrec :: Int -> MinNatural -> ShowS
Show)

instance Semiring MinNatural where
  zero :: MinNatural
zero = MinNatural
Infinity
  one :: MinNatural
one  = Natural -> MinNatural
MinNatural Natural
1

  plus :: MinNatural -> MinNatural -> MinNatural
plus MinNatural
a MinNatural
b = MinNatural
a forall a. Ord a => a -> a -> a
`min` MinNatural
b

  times :: MinNatural -> MinNatural -> MinNatural
times MinNatural
Infinity MinNatural
_ = MinNatural
Infinity
  times MinNatural
_ MinNatural
Infinity = MinNatural
Infinity
  times (MinNatural Natural
a) (MinNatural Natural
b) = Natural -> MinNatural
MinNatural (Natural
a forall a. Num a => a -> a -> a
* Natural
b)

-- | Helper to extract a set of preimages for 'inverseTotient' or 'inverseSigma'.
asSetOfPreimages
  :: (Ord a, Semiring a)
  => (forall b. Semiring b => (a -> b) -> a -> b)
  -> a
  -> S.Set a
asSetOfPreimages :: forall a.
(Ord a, Semiring a) =>
(forall b. Semiring b => (a -> b) -> a -> b) -> a -> Set a
asSetOfPreimages forall b. Semiring b => (a -> b) -> a -> b
f = forall a b. (a -> b) -> Set a -> Set b
S.mapMonotonic forall a. Mul a -> a
getMul forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall b. Semiring b => (a -> b) -> a -> b
f (forall a. a -> Set a
S.singleton forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. a -> Mul a
Mul)