{-# LANGUAGE RebindableSyntax #-} ----------------------------------------------------------------------------- -- | -- Module : Data.YAP.DirichletSeries -- Copyright : (c) Ross Paterson 2021 -- License : BSD-style (see the file LICENSE) -- -- Maintainer : R.Paterson@city.ac.uk -- Stability : provisional -- Portability : portable -- -- An example instance of the algebraic classes: formal Dirichlet series. -- ----------------------------------------------------------------------------- module Data.YAP.DirichletSeries ( -- * Formal Dirichlet series DirichletSeries, constant, fromCoefficients, -- * Queries coefficients, approximations, -- * Division reciprocal, recipSimple, divSimple, -- * Series with integral coefficients mulInt, divSimpleInt, -- * Special compositions (.-), (.*), -- * Special series -- ** Generating multiplicative functions zeta, beta, -- ** Generating characteristic functions primeZeta, primePowerZeta, -- ** Polylogarithms polylogarithm, -- * Dirichlet characters Character(..), modulus, multCharacter, lSeries, -- ** Particular real Dirichlet characters principal, primePrincipal, legendreSymbol, primitive4, primitive8s, -- * Bell series and Euler product -- | A multiplicative arithmetic function (which can be generated by a -- Dirichlet series) is determined by its values on powers of primes. -- These values can be represented by a family of power series, each -- with constant term 1, for each prime (the /Bell series/ of the -- arithmetic function). The Dirichlet series of the multiplicative -- function can be conveniently constructed using the /Euler product/ -- of this family of power series. bellSeries, eulerProduct, -- * Examples -- $examples -- ** Sum-of-divisors transform #sum_of_divisors_transform# -- $sum_of_divisors_transform -- ** Möbius transform #moebius_transform# -- $moebius_transform -- ** More integer sequences -- $more_integer_sequences -- ** Enumerating real Dirichlet characters -- $dirichlet_characters -- ** Cyclotomic polynomials -- $cyclotomic_polynomials ) where import Prelude.YAP import Data.YAP.Algebra import Data.YAP.PowerSeries (PowerSeries) import qualified Data.YAP.PowerSeries as PS import Data.YAP.Utilities.List import Data.List (intercalate, sort) import Numeric.Natural infixl 9 .- infixl 9 .* infixl 7 `divSimple` -- | Formal Dirichlet series: -- \[ f(s) = \sum_{n=1}^\infty {a_n \over n^s} \] -- -- The series can be viewed as the Dirichlet generating function of the -- sequence \( \{a_n\} \) for some arithmetic function \(a\) of positive -- integers \(n\). -- -- * The arithmetic function \(a\) is /multiplicative/ -- if \( a_{mn} = a_m a_n \) when \( (m,n) = 1 \). -- This property is preserved by multiplication. -- -- * The arithmetic function \(a\) is /completely multiplicative/ -- if \( a_{mn} = a_m a_n \) for all \(m,n\). newtype DirichletSeries a = DS [a] -- ^ Coefficients, from 1 -- | Constant series: all coefficients after the first are zero. constant :: (AdditiveMonoid a) => a -> DirichletSeries a constant a = DS [a] -- | Dirichlet series formed from a list of coefficients. -- If the list is finite, the remaining coefficients are 'zero'. fromCoefficients :: (AdditiveMonoid a) => [a] -> DirichletSeries a fromCoefficients as = DS as -- | The infinite list of coefficients of the Dirichlet series. coefficients :: (AdditiveMonoid a) => DirichletSeries a -> [a] coefficients (DS as) = as ++ repeat zero -- | The infinite list of evaluations of truncations of the series. approximations :: (Floating a) => DirichletSeries a -> a -> [a] approximations (DS as) s = scanl1 (+) [a/fromInteger n**s | (a, n) <- zip as [1..]] -- | pointwise addition of coefficients instance (AdditiveMonoid a) => AdditiveMonoid (DirichletSeries a) where DS as + DS bs = DS (add as bs) zero = DS [] atimes n (DS as) | n == zero = zero | otherwise = DS (map (atimes n) as) add :: (AdditiveMonoid a) => [a] -> [a] -> [a] add = longZipWith (+) -- | pointwise subtraction of coefficients instance (AbelianGroup a) => AbelianGroup (DirichletSeries a) where negate (DS as) = DS (map negate as) gtimes n (DS as) | n == zero = zero | otherwise = DS (map (gtimes n) as) -- | Multiplication is Dirichlet convolution (commutative). instance (Semiring a) => Semiring (DirichletSeries a) where DS _ * DS [] = DS [] DS as * DS bs = DS (multiply 1 as bs) one = constant one fromNatural i = constant (fromNatural i) -- Dirichlet convolution multiply :: (Semiring a) => Int -> [a] -> [a] -> [a] multiply n0 as bs = sumTriangle [spread n (map (a*) bs) | (n, a) <- zip [n0..] as] -- | A special case of '(*)' where the second series has integral coefficients. mulInt :: (AbelianGroup a) => DirichletSeries a -> DirichletSeries Int -> DirichletSeries a mulInt (DS as) (DS bs) = DS (multiplyIntegral 1 bs as) multiplyIntegral :: (AbelianGroup a, ToInteger a, AbelianGroup b) => Int -> [a] -> [b] -> [b] multiplyIntegral n0 as bs = sumTriangle [spread n (map (gtimes a) bs) | (n, a) <- zip [n0..] as] recipSimpleIntegral :: (AbelianGroup a, ToInteger a) => DirichletSeries a -> DirichletSeries a recipSimpleIntegral (DS []) = error "reciprocal of zero" recipSimpleIntegral (DS (a:as)) = DS bs where bs = a : map negate (multiplyIntegral 2 as bs) -- | A special case of 'divSimple' where the second series has integral -- coefficients. divSimpleInt :: (AbelianGroup a) => DirichletSeries a -> DirichletSeries Int -> DirichletSeries a divSimpleInt f g = mulInt f (recipSimpleIntegral g) -- sum a list of lists, each offset by one from the previous list sumTriangle :: (AdditiveMonoid a) => [[a]] -> [a] sumTriangle [] = [] sumTriangle ([]:xss) = zero : sumTriangle xss sumTriangle ((x:xs):xss) = x : add xs (sumTriangle xss) -- spread coefficients to every nth position spread :: (AdditiveMonoid a) => Int -> [a] -> [a] spread n as = intercalate (replicate (n-1) zero) (map (:[]) as) -- separate elements of xs with corresponding elements of seps separate :: [a] -> [[a]] -> [a] separate [] _ = [] separate (x:xs) (sep:seps) = x : case xs of [] -> [] _ -> sep ++ separate xs seps separate _ _ = error "separate: fewer separators than elements" instance (Ring a) => Ring (DirichletSeries a) where fromInteger i = constant (fromInteger i) instance (FromRational a) => FromRational (DirichletSeries a) where fromRational x = constant (fromRational x) -- | multiplies the \(n\)th coefficient by \( - \log n \) instance (Floating a) => Differentiable (DirichletSeries a) where derivative = unsafeMapWithIndex $ \ n a -> - a * log (fromIntegral n) -- | divides the \(n\)th coefficient by \( - \log n \) instance (Floating a) => Integrable (DirichletSeries a) where integral = unsafeMapWithIndex $ \ n a -> - a / log (fromIntegral n) instance AdditiveFunctor DirichletSeries where mapAdditive f (DS as) = DS (fmap f as) -- | Reciprocal of a Dirichlet series whose first coefficient is non-zero. reciprocal :: (Field a) => DirichletSeries a -> DirichletSeries a reciprocal (DS []) = error "reciprocal of zero" reciprocal (DS (a:as)) = DS bs where bs = map (/a) (1 : map negate (multiply 2 as bs)) -- | Reciprocal of a Dirichlet series whose first coefficient is 1. -- This includes any series that generates a multiplicative function, -- in which case the reciprocal also generates a multiplicative function. recipSimple :: (Ring a) => DirichletSeries a -> DirichletSeries a recipSimple (DS []) = 1 recipSimple (DS (_:as)) = DS bs where bs = 1 : map negate (multiply 2 as bs) -- | Division by a Dirichlet series whose first coefficient is 1. divSimple :: (Ring a) => DirichletSeries a -> DirichletSeries a -> DirichletSeries a divSimple p q = p * recipSimple q -- | A variant of `fmap` that is well-defined only if @f n 'zero' = 'zero'@. unsafeMapWithIndex :: (Int -> a -> b) -> DirichletSeries a -> DirichletSeries b unsafeMapWithIndex f (DS as) = DS (zipWith f [1..] as) {- | Maps a function \(f(s)\) to \(f(s - k)\) for non-negative \(k\), by multiplying the \(n\)th coefficient by \(n^k\). If \(f(s)\) generates a (completely) multiplicative function, so does \(f(s - k)\). == __Example__ @'zeta' .- k@ produces the series whose coefficients are /k/th powers: >>> coefficients $ zeta .- 1 [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,... >>> coefficients $ zeta .- 2 [1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400,... -} (.-) :: (AdditiveMonoid a) => DirichletSeries a -> Int -> DirichletSeries a f .- k | k < 0 = error "negative argument to (.-)" | otherwise = unsafeMapWithIndex (\ n a -> atimes (toInteger n ^ k) a) f {- | Maps a function \(f(s)\) to \(f(k s)\) for positive \(k\), by moving the \(n\)th coefficient to position \(n^k\), with zeros in between. If \(f(s)\) generates a multiplicative function, so does \(f(k s)\). == __Example__ @'zeta' .* k@ produces the series whose coefficients are 1 if /n/ is an exact /k/th power and 0 otherwise: >>> coefficients $ zeta .* 2 [1,0,0,1,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,... >>> coefficients $ zeta .* 3 [1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,... -} (.*) :: (AdditiveMonoid a) => DirichletSeries a -> Int -> DirichletSeries a DS as .* k | k < 0 = error "non-positive argument to (.*)" | otherwise = DS (separate as zeros) where powers = [n^k | n <- [1..]]::[Int] zeros = [replicate (nnk - nk - 1) zero | (nnk, nk) <- zip (tail powers) powers] -- | Riemann zeta function: -- \[ -- \zeta(s) = \sum_{n=1}^\infty {1 \over n^s} -- \] -- All coefficients are 1, and thus trivially completely multiplicative. zeta :: (Semiring a) => DirichletSeries a zeta = DS (repeat one) -- | The Dirichlet \(\beta\) function, with coefficients -- \(1, 0, -1, 0, 1, 0, -1, 0, \ldots\), is completely multiplicative. -- It is equal to @'lSeries' 'primitive4'@. beta :: (Ring a) => DirichletSeries a beta = lSeries primitive4 -- | The polylogarithm function -- \[ -- Li_s(z) = \sum_{n=1}^\infty {z^n \over n^s} -- \] -- applied to a given value of \(z\). -- Some special cases: -- -- * \( Li_s(1) = \zeta(s) \) -- -- * \( Li_s(- 1) = - \eta(s) \) where \(\eta\) is the Dirichlet eta function -- polylogarithm :: (Semiring a) => a -> DirichletSeries a polylogarithm z = DS (iterate (*z) z) -- | The prime zeta function, whose the \(n\)th coefficient is 1 if \(n\) -- is prime and 0 otherwise (). -- The arithmetic function is not multiplicative. primeZeta :: (Semiring a) => DirichletSeries a -- zeta*zeta gives the number of factors of n, which is 2 for primes primeZeta = DS [fromBool (n == (2::Int)) | n <- coefficients (zeta*zeta)] -- | The prime power zeta function, whose the \(n\)th coefficient is 1 -- if \(n = p^k\) for prime \(p\) and positive \(k\) and 0 otherwise -- (). -- The arithmetic function is not multiplicative. primePowerZeta :: (Semiring a) => DirichletSeries a primePowerZeta = DS [fromBool (n == (1::Int)) | n <- coefficients (zeta*primeZeta)] fromBool :: (Semiring a) => Bool -> a fromBool False = zero fromBool True = one -- Dirichlet characters -- | A Dirichlet character of modulus @n@ is a completely multiplicative -- function from integers that repeats with period @n@ and is equal to -- 0 if and only if the argument is not coprime with @n@. Because it -- is periodic, we can represent a Dirichlet character as a list of the -- values for 1 to @n@. newtype Character a = Character [a] deriving (Eq, Ord, Show) -- | Pointwise multiplication of Dirichlet characters of moduli @m@ -- and @n@, yielding a character of modulus @m*n@. instance (Semiring a) => Semigroup (Character a) where x <> y = prod (multCharacter (modulus y) x) (multCharacter (modulus x) y) -- | The modulus of a Dirichlet character modulus :: Character a -> Int modulus (Character as) = length as -- | Make a new Dirichlet character with modulus @n@ times as large -- defining the same function of the integers. -- A non-principal Dirichlet character is called /primitive/ if it is -- not such a repetition of a Dirichlet character with a smaller modulus. multCharacter :: Int -> Character a -> Character a multCharacter n (Character as) = Character (concat (replicate n as)) -- characters assumed to be of the same length prod :: (Semiring a) => Character a -> Character a -> Character a prod (Character as) (Character bs) = Character (zipWith (*) as bs) -- | The identity is the unique Dirichlet character of modulus 1. instance (Semiring a) => Monoid (Character a) where mempty = Character [one] -- | The Dirichlet L-series for a Dirichlet character \(\chi\) is -- \[ -- L(s,\chi) = \sum_{n=1}^\infty \frac{\chi(n)}{n^s} -- \] lSeries :: (AdditiveMonoid a) => Character a -> DirichletSeries a lSeries (Character as) = fromCoefficients $ cycle as -- | The principal character of modulus @n@, taking the value @1@ for -- numbers coprime with @n@ and @0@ otherwise. principal :: (Semiring a) => Int -> Character a principal m = Character [if gcd m n == 1 then one else zero | n <- [1..m]] -- | A simpler implementation of the principal character, assuming that -- the parameter is prime. primePrincipal :: (Semiring a) => Int -> Character a primePrincipal p = Character (replicate (p-1) one ++ [zero]) -- | Given an odd prime \(p\), the Legendre symbol \((n|p)\) defines -- the unique quadratic Dirichlet character of modulus \(p\) (and also -- a primitive character): for non-zero \(n\), 1 if \(n\) is a quadratic -- residue of \(p\) and -1 otherwise. legendreSymbol :: (Ring a) => Int -> Character a legendreSymbol p = Character (symbols 1 (sort [n*n `mod` p | n <- [1..(p-1) `div` 2]])) where symbols n _ | n == p = [0] symbols n (sq:sqs) | n == sq = 1 : symbols (n+1) sqs symbols n sqs = -1 : symbols (n+1) sqs -- | The unique primitive character of modulus 4 primitive4 :: (Ring a) => Character a primitive4 = Character [1, 0, -1, 0] -- | The two primitive characters of modulus 8 primitive8s :: (Ring a) => [Character a] primitive8s = [ Character [1, 0, 1, 0, -1, 0, -1, 0], Character [1, 0, -1, 0, -1, 0, 1, 0]] {- | The /Bell series/ modulo \(p\) (usually assumed to be prime) of an arithmetic function \( a_n \) is the formal power series \[ g_p(x) = a_1 + a_p x + a_{p^2} x^2 + \cdots \] If the arithmetic function \(a\) is multiplicative, then \(a_1 = 1\). If it is completely multiplicative (and \(p\) is prime), then \(a_{p^k} = a_p^k\), so that the series simplifies to \[ g_p(x) = {1 \over 1 - {a_p x}} \] -} bellSeries :: (AdditiveMonoid a) => DirichletSeries a -> Int -> PowerSeries a bellSeries (DS as) p | p < 2 = error "Bell series is not defined for p < 2" | otherwise = PS.fromCoefficients (bell_series 1 as) where bell_series _ [] = [] bell_series n (c:cs) = c:bell_series np (drop (np - n - 1) cs) where np = n*p {- | If \(g_p\) is a family of power series for each prime \(p\), such that each constant term is 1, then the /Euler product/ \[ f(s) = \prod_{p\text{ prime}} g_p({1 \over p^s}) \] yields a Dirichlet series generating a multiplicative function. If \(g_p\) is the Bell series of a multiplicative function \(a\), then \(f\) is the Dirichlet generating function of \(a\): \[ f(s) = \prod_{p\text{ prime}} g_p({1 \over p^s}) = \prod_{p\text{ prime}} \left( 1 + {a_p \over p^s} + {a_{p^2} \over p^{2s}} + {a_{p^3} \over p^{3s}} + \cdots \right) = \sum_{n=1}^\infty {a_n \over n^s} \] Thus 'eulerProduct' and 'bellSeries' define an isomorphism between Dirichlet series generating multiplicative functions and prime-indexed families of power series with constant term 1. This isomorphism preserves multiplication (but not addition): prop> one = eulerProduct $ \ p -> one prop> eulerProduct g * eulerProduct h = eulerProduct $ \ p -> g p * h p The Euler product also satisfies: prop> eulerProduct g .* k = eulerProduct $ \ p -> g p .^ k prop> eulerProduct g .- k = eulerProduct $ \ p -> g p .* p^k -} eulerProduct :: (Eq a, Semiring a) => (Int -> PowerSeries a) -> DirichletSeries a eulerProduct f = fromCoefficients $ map snd $ foldr mult [] (map powers primes) where mult ps rest = take 2 ps ++ drop 2 (mergeProduct ps rest) powers p = zip (iterate (*p) one) (PS.coefficients (getSimple (f (fromNatural p)))) getSimple ps | PS.atZero ps == one = ps | otherwise = error "power series does not have constant term 1" -- indices are ordered within each list and coprime between lists mergeProduct :: (Ord a, Semiring a, Semiring b) => [(a, b)] -> [(a, b)] -> [(a, b)] mergeProduct ((k1, v1):kvs1') kvs2@((k2, v2):kvs2') = (k1*k2, v1*v2): merge [(k1*k2', v1*v2') | (k2', v2') <- kvs2'] (mergeProduct kvs1' kvs2) mergeProduct _ _ = [] -- indices are ordered within each list and non-overlapping between lists merge :: (Ord a) => [(a, b)] -> [(a, b)] -> [(a, b)] merge kvs1@((k1, v1):kvs1') kvs2@((k2, v2):kvs2') | k1 < k2 = (k1, v1):merge kvs1' kvs2 | otherwise = (k2, v2):merge kvs1 kvs2' merge [] kvs2 = kvs2 merge kvs1 [] = kvs1 -- all prime numbers primes :: [Natural] primes = two:filter isPrime (iterate (+two) three) where two = one+one three = one+two isPrime :: Natural -> Bool isPrime n = and [n `mod` p /= zero | p <- takeWhile small primes] where small p = p*p <= n {- $examples See also * "A catalog of interesting Dirichlet series", H. W. Gould and Temba Shonhiwa, /Missouri Journal of Mathematical Sciences/ 20:1 (2008) 2-18. * "Survey of Dirichlet Series of Multiplicative Arithmetic Functions", Richard J. Mathar, 2011-2012. Dirichlet series generating multiplicative functions can also be expressed as products of Bell series. For convenience, define: >>> import qualified Data.YAP.PowerSeries as PS >>> let x = PS.identity -} {- $sum_of_divisors_transform Multiplication of `zeta` with a Dirichlet series with coefficients \( \{b_n\} \) yields a series with coefficients \[ a_n = \sum_{d \mid n} b_d \] The same transformation on sequences from ordinary generating functions is performed by `Data.YAP.PowerSeries.lambertTransform`. As the zeta function is completely multiplicative, it can be expressed as the product of the Bell series \( g_p(x) = {1 \over 1 - x} \): >>> coefficients $ eulerProduct $ const $ PS.recipSimple (1-x) [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,... The number of divisors of \(n\), i.e. \(\tau(n)\), \(d(n)\) or \(\sigma_0(n)\) () is generated by \( \zeta^2(s) \), the Euler product of the Bell series \( g_p(x) = {1 \over (1-x)^2} \): >>> coefficients $ zeta^2 [1,2,2,3,2,4,2,4,3,4,2,6,2,4,4,5,2,6,2,6,... >>> coefficients $ eulerProduct $ const $ PS.recipSimple ((1-x)^2) [1,2,2,3,2,4,2,4,3,4,2,6,2,4,4,5,2,6,2,6,... In general, @zeta^k@ is the number of factorizations of \(n\) as \(k\) factors. The sum of the divisors of \(n\), i.e. \(\sigma(n)\) () is generated by \( \zeta(s)\zeta(s-1) \), the Euler product of the Bell series \( g_p(x) = {1 \over (1-x)(1-px)} \): >>> coefficients $ zeta*(zeta .- 1) [1,3,4,7,6,12,8,15,13,18,12,28,14,24,24,31,18,39,20,42,... >>> coefficients $ eulerProduct $ \ p -> PS.recipSimple ((1-x)*(1 - atimes p x)) [1,3,4,7,6,12,8,15,13,18,12,28,14,24,24,31,18,39,20,42,... The sum of squares of the divisors of \(n\), i.e. \(\sigma_2(n)\) () is generated by \( \zeta(s)\zeta(s-2) \), the Euler product of the Bell series \( g_p(x) = {1 \over (1-x)(1-p^2 x)} \): >>> coefficients $ zeta*(zeta .- 2) [1,5,10,21,26,50,50,85,91,130,122,210,170,250,260,341,290,455,362,546,... >>> coefficients $ eulerProduct $ \ p -> PS.recipSimple ((1-x)*(1 - atimes (p^2) x)) [1,5,10,21,26,50,50,85,91,130,122,210,170,250,260,341,290,455,362,546,... and so on for the divisor function \(\sigma_k(n) = \sum_{d \mid n} d^k\) for any power \(k\). The number of ways of writing \(n\) as the sum of at most two nonzero squares (): >>> coefficients $ zeta*beta [1,1,0,1,2,0,0,1,1,2,0,0,2,0,0,1,2,1,0,2,... In the special case where \(\{b_n\}\) is a characteristic sequence of some property (consisting of 0s and 1s), \(a_n\) is the number of divisors of \(n\) having that property. The number of distinct prime divisors of \(n\), i.e. \(\omega(n)\) (): >>> coefficients $ zeta*primeZeta [0,1,1,1,1,2,1,1,1,2,1,2,1,2,2,1,1,2,1,2,... The number of prime divisors of \(n\) counted with multiplicity, i.e. \(\Omega(n)\) (): >>> coefficients $ zeta*primePowerZeta [0,1,1,2,1,2,1,3,2,2,1,3,1,2,2,4,1,3,1,3,... -} {- $moebius_transform The 'zeta' function has a multiplicative inverse, whose coefficients are given by the Möbius function \( \mu(n) \) (): * 0 if \(n\) is divisible by the square of a prime, and otherwise * 1 if \(n\) has an even number of prime factors, or * \(-1\) if \(n\) has an odd number of prime factors. >>> coefficients $ recipSimple zeta [1,-1,-1,0,-1,1,-1,0,0,1,-1,0,-1,1,1,0,-1,0,-1,0,... Thus division by 'zeta' performs the inverse of the sum-of-divisors transform, i.e. the Möbius transform. The same transformation on sequences from ordinary generating functions is performed by `Data.YAP.PowerSeries.inverseLambertTransform`. The Möbius function is multiplicative (but not completely multiplicative), and expressed by the product of the reciprocal Bell series \( g_p(x) = 1 - x \): >>> coefficients $ eulerProduct $ const $ 1-x [1,-1,-1,0,-1,1,-1,0,0,1,-1,0,-1,1,1,0,-1,0,-1,0,... The number of numbers less than \(n\) and coprime with \(n\) () is the Euler totient function \( \varphi(n) \), which satisfies \[ \sum_{d \mid n}\varphi(d) = n \] i.e. @phi*'zeta' = 'zeta' .- 1@. Hence the generating function for @phi@ can be obtained by dividing the righthand side of this equation by 'zeta' (Möbius inversion), i.e. \( {\zeta(s-1) \over \zeta(s)} \), the Euler product of the Bell series \( g_p(x) = {1-x \over 1-px} \): >>> coefficients $ zeta .- 1 `divSimple` zeta [1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,8,... >>> coefficients $ eulerProduct $ \ p -> (1-x) `PS.divSimple` (1 - atimes p x) [1,1,2,2,4,2,6,4,6,4,10,4,12,6,8,8,16,6,18,8,... The generalization from @1@ to arbitrary \(k\) is the Jordan totient function \(J_k(n)\), which satisfies \[ \sum_{d \mid n} J_k(d) = n^k \] The Liouville function \( \lambda(n) = (-1)^{\Omega(n)} \) () has the property that \( \sum_{d|n}\lambda(d) \) is 1 if \(n\) is a perfect square and 0 otherwise, i.e. @liouville*'zeta' = 'zeta' .* 2@. Hence the Liouville function can be obtained by dividing the righthand side of this equation by 'zeta', i.e. \( {\zeta(2s) \over \zeta(s)} \), the Euler product of the Bell series \( g_p(x) = {1-x \over 1-x^2} = {1 \over 1+x} \): >>> coefficients $ zeta .* 2 `divSimple` zeta [1,-1,-1,1,-1,1,-1,-1,1,1,-1,-1,-1,1,1,1,-1,-1,-1,-1,... >>> coefficients $ eulerProduct $ const $ PS.recipSimple $ 1+x [1,-1,-1,1,-1,1,-1,-1,1,1,-1,-1,-1,1,1,1,-1,-1,-1,-1,... The von Mangoldt function \( \Lambda(n) \) (log of ): * \( \log p \) if \(n=p^k\) for prime \(p\) and positive integer \(k\) * 0 otherwise This function satisfies \[ \sum_{d \mid n} \Lambda(d) = \log~n \] so its Dirichlet generating function \( - {\zeta'(s) \over \zeta(s)} \) can be obtained by Möbius inversion: >>> coefficients $ - derivative zeta `divSimple` zeta [0.0,0.6931472,1.0986123,0.6931472,1.609438,-0.0,1.9459101,0.6931472,... If \(a(n)\) is the number of aperiodic binary strings of length \(n\) (), then \( 2^n = \sum_{d \mid n} a(n) \), so \(a(n)\) can be obtained by Möbius inversion of the sequence of powers of 2: >>> coefficients $ polylogarithm 2 `divSimple` zeta [2,2,6,12,30,54,126,240,504,990,2046,4020,8190,16254,... -} {- $more_integer_sequences The reciprocal of the series generating the Liouville function generates \(\mu(n)^2\), which is 1 if \(n\) is squarefree, and otherwise 0 (). Its generating function is \( {\zeta(s) \over \zeta(2s)} \), the Euler product of the Bell series \( g_p(x) = {1-x^2 \over 1-x} = 1+x \): >>> coefficients $ zeta `divSimple` zeta .* 2 [1,1,1,0,1,1,1,0,0,1,1,0,1,1,1,0,1,0,1,0,... >>> coefficients $ eulerProduct $ const $ 1+x [1,1,1,0,1,1,1,0,0,1,1,0,1,1,1,0,1,0,1,0,... Applying the sum-of-divisors transform to that, by multiplying by 'zeta', yields a series generating the number of squarefree divisors of \(n\), or equivalently the number of divisors of \(n\) that are coprime with the complementary divisor, i.e. \( \theta(n) = 2^{\omega(n)}\) (). Its generating function is \( {\zeta^2(s) \over \zeta(2s)} \), the Euler product of the Bell series \( g_p(x) = {1 - x^2 \over (1-x)^2} = {1+x \over 1-x} \): >>> coefficients $ zeta^2 `divSimple` zeta .* 2 [1,2,2,2,2,4,2,2,2,4,2,4,2,4,4,2,2,4,2,4,... >>> coefficients $ eulerProduct $ const $ (1+x) `PS.divSimple` (1-x) [1,2,2,2,2,4,2,2,2,4,2,4,2,4,4,2,2,4,2,4,... The Dedekind \(\psi\) function () has generating function \( {\zeta(s)\zeta(s-1) \over \zeta(2s)} \), the Euler product of the Bell series \( g_p(x) = {1-x^2 \over (1-x)(1-px)} = {1+x \over 1-px} \): >>> coefficients $ zeta*(zeta .- 1) `divSimple` zeta .* 2 [1,3,4,6,6,12,8,12,12,18,12,24,14,24,24,24,18,36,20,36,... >>> coefficients $ eulerProduct $ \ p -> (1+x) `PS.divSimple` (1 - atimes p x) [1,3,4,6,6,12,8,12,12,18,12,24,14,24,24,24,18,36,20,36,... The squarefree part of \(n\) () has generating function \( {\zeta(s-1)\zeta(2s) \over \zeta(2s - 2)} \), the Euler product of the Bell series \( g_p(x) = {1 - (p x)^2 \over (1-p x)(1-x^2)} \): >>> coefficients $ (zeta .- 1)*(zeta .* 2) `divSimple` zeta .- 2 .* 2 [1,2,3,1,5,6,7,2,1,10,11,3,13,14,15,1,17,2,19,5,... >>> coefficients $ eulerProduct $ \ p -> (1 - atimes p x^2) `PS.divSimple` ((1-x^2)*(1 - atimes p x)) [1,2,3,1,5,6,7,2,1,10,11,3,13,14,15,1,17,2,19,5,... The cubefree part of \(n\) () has generating function \( {\zeta(s-1)\zeta(3s) \over \zeta(3s - 3)} \), the Euler product of the Bell series \( g_p(x) = {1 - (p x)^3 \over (1-p x)(1-x^3)} \): >>> coefficients $ (zeta .- 1)*(zeta .* 3) `divSimple` zeta .- 3 .* 3 [1,2,3,4,5,6,7,1,9,10,11,12,13,14,15,2,17,18,19,20,... >>> coefficients $ eulerProduct $ \ p -> (1 - atimes p x^3) `PS.divSimple` ((1-x^3)*(1 - atimes p x)) [1,2,3,4,5,6,7,1,9,10,11,12,13,14,15,2,17,18,19,20,... The number of perfect partitions of \(n\) () has generating function \( {1 \over 2 -\zeta(s)} \): >>> coefficients $ recipSimple (2 - zeta) [1,1,1,2,1,3,1,4,2,3,1,8,1,3,3,8,1,8,1,8,3,3,1,20,2,3,4,8,1,13,... -} {- $dirichlet_characters The number of real Dirichlet characters modulo \(n\) is the number of square roots of unity modulo \(n\) (), whose series is generated by \((1 - {1 \over 2^s} + {2 \over 4^s}) {\zeta(s)^2 \over \zeta(2s)}\): >>> coefficients $ fromCoefficients [1, -1, 0, 2] * zeta^2 `divSimple` zeta .* 2 [1,1,2,2,2,2,2,4,2,2,2,4,2,2,4,4,2,2,2,4,4,2,2,8,2,2,2,4,2,4,... We can list these characters by defining series over the semiring @'Data.YAP.FiniteSet.FiniteSet' 'Character'@ using the Euler product. >>> import Data.YAP.FiniteSet First, we define a Bell series of principal characters for prime \(p\) as powers of the principal character modulo \(p\): >>> principals p = PS.recipOneMinus `PS.compose` PS.fromCoefficients [zero, singleton (primePrincipal p)] >>> PS.coefficients $ principals 3 [fromList [Character [1]], fromList [Character [1,1,0]], fromList [Character [1,1,0,1,1,0,1,1,0]], fromList [Character [1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0,1,1,0]],... The Euler product of this Bell series yields a series of singleton sets, each containing the principal character of each number \(n\): >>> coefficients $ eulerProduct principals [fromList [Character [1]], fromList [Character [1,0]], fromList [Character [1,1,0]], fromList [Character [1,0,1,0]], fromList [Character [1,1,1,1,0]], fromList [Character [1,0,0,0,1,0]],... The Bell series of (non-principal) primitive real Dirichlet characters for each prime \(p\) is finite: >>> primitives p = PS.fromCoefficients $ if p == 2 then [one, zero, singleton primitive4, fromList primitive8s] else [one, fromList [primePrincipal p, legendreSymbol p]] The Euler product of the product of these two series yields the set of real Dirichlet characters modulo each \(n\): >>> coefficients $ eulerProduct $ \ p -> principals p * primitives p [fromList [Character [1]], fromList [Character [1,0]], fromList [Character [1,-1,0],Character [1,1,0]], fromList [Character [1,0,-1,0],Character [1,0,1,0]], fromList [Character [1,-1,-1,1,0],Character [1,1,1,1,0]], fromList [Character [1,0,0,0,-1,0],Character [1,0,0,0,1,0]], fromList [Character [1,1,-1,1,-1,-1,0],Character [1,1,1,1,1,1,0]], fromList [Character [1,0,-1,0,-1,0,1,0],Character [1,0,-1,0,1,0,-1,0],Character [1,0,1,0,-1,0,-1,0],Character [1,0,1,0,1,0,1,0]], fromList [Character [1,-1,0,1,-1,0,1,-1,0],Character [1,1,0,1,1,0,1,1,0]], fromList [Character [1,0,-1,0,0,0,-1,0,1,0],Character [1,0,1,0,0,0,1,0,1,0]],... -} {- $cyclotomic_polynomials The sequence of cyclotomic polynomials \( \Phi_n(x) \) () satisfies \[ x^n - 1 = \prod_{d \mid n} \Phi_d(x) \] The list of left hand sides can be defined by >>> import qualified Data.YAP.Polynomial as Poly >>> let x = Poly.identity >>> let xnm1s = [xn - 1 | xn <- iterate (*x) x] Applying the 'Data.YAP.Logarithm.Logarithm' adaptor yields a summation: \[ \log(x^n - 1) = \sum_{d \mid n} \log(\Phi_d(x)) \] From this, the cyclotomic polynomials can be obtained using Möbius inversion in the field of rational polynomials. The denominators will all be 1, and can be discarded. >>> import Data.YAP.Logarithm >>> import Data.YAP.Ratio >>> map (mapAdditive numerator . numerator . exponential) $ coefficients $ fromCoefficients (map (logarithm . (% 1)) xnm1s) `divSimpleInt` zeta [fromCoefficients [-1,1], fromCoefficients [1,1], fromCoefficients [1,1,1], fromCoefficients [1,0,1], fromCoefficients [1,1,1,1,1], fromCoefficients [1,-1,1], fromCoefficients [1,1,1,1,1,1,1],... The numerators are always 0, 1 or -1, and the denominators are always 1. -}