bodigrim / arithmoi

Number theory: primes, arithmetic functions, modular computations, special sequences

arithmoi

arithmoi's Issues

New release of `arithmoi`

Can the current master be released to Hackage. It will help with Stackage's commercialhaskell/stackage#1691 going forward (we don't want to disable tests!)

Btw, I checked that master is supports all recent versions of other dependencies as well.

Implement Hurwitz Zeta function

After #44 and #120, both the Riemann zeta function and the Dirichlet beta function have been implemented.

The Hurwitz zeta function is another important function that can be used to define both the Riemann zeta and the Dirichlet beta functions, so this issue is about refactoring the two latter functions to use the former.

This should all go into the newly created Math.NumberTheory.Zeta namespace.

Ramanujan tau-function

Ramanujan tau function is probably the largest omission in Math.NumberTheory.ArithmeticFunctions.Standard. It is a multiplicative function, important in theory of modular forms.

A straightforward approach to its computation is to implement these formulas. This will be a nice first step.

More efficient algorithm is given by Charles 2006, but this refinement may be done latter, as a separate task.

Polymorphic fibonacci/lucas

Currently Math.NumberTheory.Recurrencies.Linear exposes monomorphic functions fibonacci, fibonacciPair, lucas, lucasPair, generalLucas, which return Integer values. This is not suitable for some scenarios.

For instance, recently I needed a Fibonacci number with index 10^15 modulo 123456789. The Fibonacci number itself is eminently huge, but, since modulo is small, the reduced value should be computable in a blink of an eye. If fibonacci function would be polymorphic by return value, allowing any Num, I could specify Mod 123456789 as a type and be done. But it is not flexible this way.

The proposal is to change type signatures to

fibonacci :: Num a => Int -> a
fibonacciPair :: Num a => Int -> (a, a)
lucas :: Num a => Int -> a
lucasPair :: Num a => Int -> (a, a)
generalLucas :: Num a => a -> a -> Int -> (a, a, a, a)

This change should be pretty mechanical; the only thing to deal with is shiftL, which should be replaced by addition. The main computation cost is caused by multiplications, any way.

This change is also more in line with Math.NumberTheory.Recurrencies.Bilinear, where all recurrent sequences are as polymorphic as possible.

Redesign API for modular square roots

Currently Math.NumberTheory.Moduli.Sqrt is a fractal of bad design. There are eight functions exposed and zero of them are actually suitable to find a modular square root x^2 ≡ n (mod m). User is obliged to factorise a modulo himself and satisfy a list of funny preconditions, stated in comments, but not enforced by types. And if some of preconditions are not satisfied, "the computation may not terminate".

The proposal is to redesign API from the scratch. There should be a simple function sqrtMod :: KnownNat m => Mod m -> Maybe (Mod m), which actually computes modular square root. Sure, since such function involves factorisation of the modulo at each call, more specific routines should be exposed as well. But:

• Most of preconditions must be encoded in types. If a function expects a prime argument, it must be Prime Integer instead of Integer.
• Other preconditions (e. g., primes of special form) must be checked in all publicly exposed functions.
• All computations must terminate, even if the input does not match expectations.

Do not report -1 as a factor

Current factorise routines (both Math.NumberTheory.Primes.Factorisation.Montgomery and Math.NumberTheory.GaussianIntegers) list -1 as a member of canonical factorisation of negative numbers. This is weird from mathematical perspective (canonical factorisation should consist of primes only) and will eventually be removed, when all primes throughout the library will be represented by newtypes from Math.NumberTheory.Primes.Types.

As a small step towards the goal, one can at least stop reporting -1 and Gaussian units in factorisations. E. g., factorise (-1) should be []. This is a small change, but take care to fix the test suite accordingly.

Not in scope: type constructor or class ‘Word64’

I had problems when compiling the package, getting this error for a few files.

• GHC version: 7.10.2
• cabal version:1.22.6.0, cabal library 1.22.4.0
• environment: (stackage-lts-3.17)[http://www.stackage.org/lts-3.17]
• System: 32-bit Debian Linux

The diff below fixed the problem for me, although it's probably not an universal solution:

diff --git a/Math/NumberTheory/Moduli.hs b/Math/NumberTheory/Moduli.hs
index 1daa9e0..29390d5 100644
--- a/Math/NumberTheory/Moduli.hs
+++ b/Math/NumberTheory/Moduli.hs
@@ -42,6 +42,7 @@ import Data.Array.Unboxed
 import Data.Array.Base (unsafeAt)
 import Data.Maybe (fromJust)
 import Data.List (nub)
+import Data.Word
 import Control.Monad (foldM, liftM2)

 import Math.NumberTheory.Utils (shiftToOddCount, splitOff)
diff --git a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
index f57f629..d094dad 100644
--- a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
+++ b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs
@@ -42,9 +42,9 @@ import Data.Array.ST
 #endif
 import Control.Monad (when)
 import Data.Bits
-#if __GLASGOW_HASKELL__ < 709
+-- #if __GLASGOW_HASKELL__ < 709
 import Data.Word
-#endif
+-- #endif

 import Math.NumberTheory.Powers.Squares (integerSquareRoot)
 import Math.NumberTheory.Utils

---

I also got a bunch of warnings like

Math/NumberTheory/Primes/Sieve/Misc.hs:153:93: Warning:
    Literal 3737841677 is out of the Int range -2147483648..2147483647

It's probably ok, but it'd be nice to avoid them so that people don't get scared something is broken :)

Check whether `evenI` is still faster than `even`

Module Math.NumberTheory.Moduli.Jacobi contains an interesting helper function:

-- For large Integers, going via Int is much faster than bit-fiddling
-- on the Integer, so we do that.
evenI :: Integral a => a -> Bool
evenI n = fromIntegral n .&. 1 == (0 :: Int)

Actually, it was definitely much faster than even from Prelude before GHC 8.4 due to trac 14437. We should check (looking at Core and CMM), whether evenI is still faster. And if it is, push its definition into base for the profit of the whole community.

Getting the smallest prime factors from a FactorSieve

It does not seem that there is an inexpensive way to look up the contents of a FactorSieve, which is unfortunate because a factor sieve can have other useful purposes beyond generating prime factorizations.

Due to memory optimizations I imagine that directly exposing the container is out of the question - but even just a lookup function like sieveSmallestFactor :: Sieve -> Integer -> Integer (or Maybe Integer) would be immensely useful.

Of course, there are the edge cases of how to handle inputs of 0 or 1.

Support Natural on the same level with Int, Word and Integer

There are three "first-class" integer (implementing Integral instances) types in arithmoi: Int, Word and Integer. It means that polymorphic functions Integral a => ... has been written, keeping in mind only these three possibilities for a. Usually such functions get {-# SPECIALISE #-} or/and {-# RULES #-} to improve performance. Here is an example from Math.NumberTheory.Powers.Squares:

integerSquareRoot :: Integral a => a -> a
integerSquareRoot n
  | n < 0       = error "integerSquareRoot: negative argument"
  | otherwise   = integerSquareRoot' n
{-# SPECIALISE integerSquareRoot :: Int -> Int,
                                    Word -> Word,
                                    Integer -> Integer
  #-}

integerSquareRoot' :: Integral a => a -> a
integerSquareRoot' = isqrtA
{-# RULES
"integerSquareRoot'/Int"     integerSquareRoot' = isqrtInt'
"integerSquareRoot'/Word"    integerSquareRoot' = isqrtWord
"integerSquareRoot'/Integer" integerSquareRoot' = isqrtInteger
  #-}

Further, such function are tested using Math.NumberTheory.TestUtils.testIntegralProperty to run each test thrice with different types.

Nowadays another integer type gains more traction: Natural, included in base since GHC 8.0. It represents non-negative integers of arbitrary size; Natural corresponds to Integer as Word corresponds to Int. We should support it through arithmoi on the same level with Int, Word and Integer.

There are three steps:

1. Check that polymorphic functions Integral a => ... are actually compatible with Natural arguments. This may be not the case, because (in contrast to Word) an attempt to create a negative Natural throws Underflow :: ArithException.

E. g., f x = (x - 1) + 2 is a perfectly valid function, when x is Int or Integer and even has an expected semantics for Word, but raises an exception on x = 0 :: Natural.

One can change Math.NumberTheory.TestUtils in the following way and adapt testIntegralProperty, testSameIntegralProperty and testIntegral2Property correspondingly. This should cover majority of cases, and other ones should be dealt manually.

 type TestableIntegral wrapper =
-  ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer]
-  , Matrix '[Arbitrary, Show] wrapper '[Large Int, Large Word, Huge Integer]
+  ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer, Natural]
+  , Matrix '[Arbitrary, Show] wrapper '[Large Int, Large Word, Huge Integer, Huge Natural]
   , Matrix '[Bounded, Integral] wrapper '[Int, Word]
   , Num (wrapper Integer)
+  , Num (wrapper Natural)
   , Functor wrapper
   )

2. Add {-# SPECIALISE #-} pragmas to cover all four "first-class" integer types. This is the easiest, almost mechanical part of the task.

3. Add rewrite rules {-# RULES #-} and provide corresponding functions, benefiting from Natural arguments, where appropriate. May be easy, may be hard, may be superfluous.

Lazy version of divisors

At the moment there are two functions to list divisors of an integer:

• divisors, which returns Set n,
• and its specialised version with better performance divisorsSmall, which returns IntSet.

However, there are use cases, when it is desirable to generate divisors lazily (I stumbled on this during solving of Problem 357). This is certainly possible, because underlying factorise function produces prime factors in lazy fashion.

sieveFactor (factorSieve 2048) (29*73) == [(2117,1)]

This may just be a documentation bug or perhaps I am misreading this:

sieveFactor fs n finds the prime factorisation of n using the FactorSieve fs. For negative n, a factor of -1 is included with multiplicity 1. After stripping any present factors 2, the remaining cofactor c (if larger than 1) is factorised with fs. This is most efficient of course if c does not exceed the bound with which fs was constructed. If it does, trial division is performed until either the cofactor falls below the bound or the sieve is exhausted. In the latter case, the elliptic curve method is used to finish the factorisation.

I understood this to mean that the FactorSieve is used to speed up the factorisation, but if an occasional argument to (sieveFactor fs) is higher than the argument to factorSieve, the results would still be correct (if less sped-up by the sieve).

However, this appears not to be the case. For example with arithmoi 0.4.1.2 on Haskell Platform 2014.2 (i.e., ghc 7.8.3), I get, for example: '''sieveFactor (factorSieve 2048) (29*73)==[(2117,1)]'''

A some checks seem to find this discrepancy only in cases where the smallest prime factor is 29.

Find Gaussian prime with a given norm in effective way

There are three kinds of Gaussian primes (see Math.NumberTheory.GaussianIntegers.primes):

• 1 + i,
• p + 0i, where p is an integer prime and p = 4k+3,
• a + bi, where a^2+b^2 is an integer prime and p = 4k+1.

In the third case it is non-trivial how to find a and b, satisfying a^2+b2 = p for a given p. Currently Math.NumberTheory.GaussianIntegers.findPrime' does it via brute force.

We should rather employ the Hermite-Serret algorithm. Resources:
https://math.stackexchange.com/a/5883
http://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdf

Dirichlet characters

https://en.wikipedia.org/wiki/Dirichlet_character

Despite of rather complicated applications, Dirichlet characters itself are a simple algebraic concept.
It would be nice to have a data type, representing them, annotated by modulo on type level:

data DirichletCharacter (mod :: Nat) = TBD

Since Dirichlet characters of given modulo form a group by multiplication, we should be able to define correspondent Semigroup and Monoid instances. Moreover, it is useful to be able to enumerate all Dirichlet characters, which can be reflected in Haskell by Enum and Bounded instances. Finally, each character is a multiplicative function, so there is a mapping between DirichletCharacter and ArithmeticFunction from Math.NumberTheory.ArithmeticFunctions.

Discrete logarithms

Once #86 lands into master, it would be nice to implement any non-trivial (or even trivial one!) algorithm for computation of discrete logarithms. I do not know much about this field, but hopefully wiki could serve a good starting point.

Migrate from array to vector

Older parts of arithmoi tend to use Data.Array instead of Data.Vector, as newer modules do. It would be nice to migrate to vectors completely. The migration should be pretty straightforward, so it is a nice first issue.

Partition function

The partition function p(n) represents the number of possible partitions of a natural number n, which is to say the number of distinct ways of representing n as a sum of natural numbers (with order irrelevant).

It would be nice to implement the partition function as a self-recurrent sequence partition :: Num a => [a] , using the pentagonal number theorem. The right place for it will be a new module Math.NumberTheory.Recurrencies.

sieve design and performance experiments

@ExpHP @Bodigrim lets talk out concrete ideas / types / experiments / data structures for fancy nice sieving here :)

Use Chudnovsky rational approximation of pi^2 for zetasEven

To evaluate Riemann zeta function on even arguments (zetasEven) we use formula:

The main issue here is that for reasonably big n we divide one huge number by another one, meaning a catastrophic loss of precision. For instance,

> approximateValue (zetasEven !! 25) :: Double
0.9999999999999996
> approximateValue (zetasEven !! 999) :: Double
NaN

which does not make any sense, because the zeta function is always finite and > 1. Currently we fix this behaviour by flooring values of zetasEven by 1.

This discrepancy hurts evaluation of Riemann zeta function on odd arguments and propagates further to numeric errors in Dirichlet beta function.

---

This is a rough idea, I have not given it much thought.

Chudnovsky algorithm gives us a rapidly converging series of rational approximations to pi^2. We can use it to generate a rational approximation of zetasEven !! n and consequently calculate not only zetaEven !! n itself (which is any way close to 1), but also zetasEven !! n - 1 using usual Double with a fairly good relative precision.

Rings of quadratic integers

It would be nice to implement a basic support for quadratic integers in arithmoi. There is a proof of concept in Math.NumberTheory.Quadratic.GaussianIntegers and Math.NumberTheory.Quadratic.EisensteinIntegers modules.

1. Implement quadratic integers as a data type with two Integer fields and discriminant of the field D on the type level. One can copy the kind of discriminants from Numeric.QuadraticIrrational.SquareFree into arithmoi tree.

{-# LANGUAGE DataKinds      #-}
{-# LANGUAGE KindSignatures #-}

import GHC.TypeNats.Compat
import Numeric.QuadraticIrrational.SquareFree

data QuadraticInteger (d :: SquareFree)
  = QuadraticInteger { a :: Integer, b :: Integer }

2. Define Num (QuadraticInteger d) instance. This will require definitions of conjugation and norm as well.

3. For principal rings (D = −1, −2, −3, −7, −11, −19, −43, −67, −163) define a function primes, generating an infinite sequence of prime elements, ordered by norm. Use solveQuadratic to solve norm x = p and norm x = p^2 as described in this comment.

4. Define UniqueFactorisation instance for principal rings.

5. Define Euclidean instance for Euclidean rings.

Segmentation fault with primeList and primeSieve

Reported by @treeowl (dafis/arithmoi#16).

To reproduce:

$ cabal sandbox init
$ cabal install arithmoi-0.4.1.1 hspec-2.1.2 --ghc-options=-fasm
$ cat >Spec.hs <<'EOF'

import Test.Hspec
import Test.QuickCheck
import Math.NumberTheory.Primes

primesSumWonk :: Int -> Int
primesSumWonk upto = sum . takeWhile (< upto) . map fromInteger . primeList $ primeSieve (toInteger upto)

primesSum :: Int -> Int
primesSum upto = sum . takeWhile (< upto) . map fromInteger $ primes

main = do
  hspec $
    describe "Primes.primesSumWonk" $
      it "works like primesSum" $
        property $ \n -> let upto = abs (n `rem` 45)
                         in primesSum upto == primesSumWonk upto

EOF
$ cabal exec runhaskell Spec.hs

Primes.primesSumWonk
Segmentation fault (core dumped)

@treeowl adds:

FYI, if you do decide to attempt to fix this, you might look at making it use the bitvec package for mutable bit vectors instead of weaving that junk into the prime sieve in a way that leaves me utterly baffled about just what it's supposed to be doing. Another very pleasant alternative might be to just steal sieves from the NumberSieves package. I have no idea about how their performance compares, theoretically, but anything is better than code that's blatantly broken.

GHC-8.0 compatible release

AFAICS the master branch builds ok with GHC-8.0.1. The Hackage version has restrictive bounds though (ghc-prim < 0.5).

Implement extendedGCD via integer-gmp

I tried to implement Math.NumberTheory.GCD.extendedGCD via GHC.Integer.GMP.Internals.extendedGCDInteger in a temporary branch.

extendedGCD :: Integral a => a -> a -> (a, a, a)
extendedGCD a b = (fromInteger d', fromInteger u', fromInteger v')
  where
    (d', u', v') = extendedGCDInteger (toInteger a) (toInteger b)

extendedGCDInteger :: Integer -> Integer -> (Integer, Integer, Integer)
extendedGCDInteger a b = (d, u, v)
  where
    (# d, u #) = gcdExtInteger a b
    v = if b == 0 then 0 else (d - u * a) `quot` b

But after this change our test-suite fails with a horrible error:

Assertion failed: (0 <= gn && gn <= gn0), function integer_gmp_gcdext, 
file libraries/integer-gmp/cbits/wrappers.c, line 309.

I suspect the issue is with integer-gmp (did anyone ever used integer_gmp_gcdext at all?) and not with arithmoi, but currently I do not have time to investigate further and fill a proper trac ticket.

Sort Gaussian primes by norm

It is preferable for applications (e. g., trial division) to list Gaussian primes sorted by ascending norm. Currently this is not so:

> take 10 Math.NumberTheory.GaussianIntegers.primes
[1+ι,3,2+ι,1+2*ι,7,11,3+2*ι,2+3*ι,4+ι,1+4*ι]
> map norm $ take 10 Math.NumberTheory.GaussianIntegers.primes
[2,9,5,5,49,121,13,13,17,17]

Generate a list of numbers with given prime divisors

It would be nice to have a function, which lists numbers with prime divisors in a given set. Here is a very inefficient implementation:

foo :: Set Integer -> Integer -> Integer -> [Integer]
foo ps from to = filter (all (flip S.member ps . fst) . factorise) [from..to]

Of course, real implementation should avoid factorisation completely.

Side note: a paper on smooth numbers, especially Ch. 5.7.

Some Useful Polynomials and Other Series

Here are some useful functions from my vault. Feel free to add or modify to your tastes any or all of them for future arithmoi releases. If there is interest I have some other number theoretic codes which may be of broader interest and which I'll be happy to share once I've cleaned them up for general use.

This is my Bernoulli series. I have checked that it produces the same results as the one in arithmoi, but not which is faster:

bernoullis =  [1,-1%2]++go (zipWith (\m pt -> map (% m) pt) [-3,-4..] $ map (init.init) $ drop 3 $ binomial)
  where
    go :: [[Rational]] -> [Rational]
    go (b:_:bs) = sum (zipWith (*) bernoullis b):0:go bs
    go _        = undefined 

The following code generates lists of polynomials, each of which in turn is represented as a simple list of coefficients starting with the constant term coefficient:

chebyshevPolysFirst :: forall a. (Num a) => [[a]]
chebyshevPolysFirst = [1]:[0,1]:(zipWith (\t1 t2 -> zipWith (-) (0:map (2*) t1) (t2++[0,0])) (tail chebyshevPolysFirst) chebyshevPolysFirst)

chebyshevPolysSecond :: forall a. (Num a) => [[a]]
chebyshevPolysSecond = [1]:[0,2]:(zipWith (\t1 t2 -> zipWith (-) (0:map (2*) t1) (t2++[0,0])) (tail chebyshevPolysSecond) chebyshevPolysSecond)

legendrePolys :: [[Rational]]
legendrePolys = [1]:[0,1]:(zipWith3 (\n p1 p2 -> zipWith (-) (0:map ((2*n-1)%n*) p1) ((map ((n-1)%n*) p2)++[0,0])) [2..] (tail legendrePolys) legendrePolys)

hermitePolysProb :: forall a. (Num a) => [[a]]
hermitePolysProb = [1]:[0,1]:(zipWith3 (\n h1 h2 -> zipWith (-) (0:h1) ((map (fromIntegral n*) h2)++[0,0])) [1..] (tail hermitePolysProb) hermitePolysProb)

hermitePolysPhysics :: forall a. (Num a) => [[a]]
hermitePolysPhysics = [1]:[0,2]:(zipWith3 (\n h1 h2 -> zipWith (-) (0:map (2*) h1) ((map ((2*fromIntegral n)*) h2)++[0,0])) [1..] (tail hermitePolysPhysics) hermitePolysPhysics)

besselPolys :: forall a. (Num a) => [[a]]
besselPolys = [1]:[1,1]:(zipWith3 (\n y1 y2 -> zipWith (+) (0:map ((2*fromIntegral n-1)*) y1) (y2++[0,0])) [2..] (tail besselPolys) besselPolys)

To evaluate any of these polynomials at a point, the following function is helpful:

polyEval p x = sum $ zipWith (*) p $ iterate (x*) 1

So, for example, the evaluate the 10th Chebyshev Polynomial of the First Kind at 5: polyEval (chevyshevPolysFirst!!10) 5 == 4517251249

The following offers the Faulhaber polynomial of degree p:

faulhaberPoly p = zipWith (*) ((0:) $ reverse $ take (p+1) $ bernoullis) $ map (% (fromIntegral p+1)) $ zipWith (*) (iterate negate (if odd p then 1 else -1)) $ binomials !! (fromIntegral p+1)

So, for example, to evaluate the sum of the 18th power of the first million integers:

sum [i^18|i<-[1..1000000]] == 52632078948868421052624778947368455052631578814768421052999912280701083128822055851844611528460355137844666500000

The advantage of this summation formula is that it takes time on the order of the exponent, rather than length of the original power sum. So even gigantic sums can be calculated instantly. For example, the sum of the 97th power of the first 10^18 integers (a 1,747 digit number) can be evaluated with numerator (polyEval (faulhaberPoly 97) 10^18) within milliseconds, even without compiling.

[A better implementation might generate an infinite list of polynomials (like the *Polys functions) and factor out the numerator in order to avoid needless ratio reduction, but this version is plenty fast for all the purposes I needed it.]

Finally, the following code evaluates the general hypergeometric function:

hypergeometric as bs x = sum $ zipWith4 (\ak bk xk fk -> ak/bk*xk*(1%fk)) aks bks xks factorial
  where
    xks = iterate (x*) 1
    aks = takeWhile (0/=) $ scanl (*) 1 $ map product $ iterate (map (1+)) as
    bks = scanl (*) 1 $ map product $ iterate (map (1+)) bs

The documentation for the latest version of this package is missing on Hackage.

The documentation for the latest version of this package is missing on Hackage. They've been pending on Hackage for some time, this usually indicated that Hackage wasn't able to build the package.

Is it possible for you to manually upload the package using cabal upload --doc, please?

Thanks!

Eisenstein integers

My feeling is that currently #73 is too broad and too vague. We should rather explore more special cases to hone API before taking a stab at the general problem.

We have already implemented the simplest quadratic ring, which is the ring of Gaussian integers. The next step is to handle the ring of Eisenstein integers.

The API may follow the outline of Math.NumberTheory.GaussianIntegers. Required mathematical background can be found here:

• https://en.wikipedia.org/wiki/Eisenstein_integer
• https://en.wikipedia.org/wiki/Eisenstein_prime

Change output of prime sieves from [Integer] to [Word]

Currently Eratosthenes sieves from Math.NumberTheory.Primes.Sieve.Eratosthenes output a list of Integer and operate internally on Integer and Int. This was reasonable historically, when x32 architecture was omnipresent. Nowadays it does not make much sense: Eratosthenes sieves consume far too much memory beyond 2^64 to be practical.

That said, it would be more memory efficient to use Int or Word to represent generated primes. Since primes are unsigned, Word looks good. It would be even better to use newtype wrapper Prm from Math.NumberTheory.Primes.Types, but it may be deferred to a separate patch.

Another idea to achieve better backward compatibility: sieves can be parametrised by unsigned type of offset (Word or Natural) to support rare cases, when one would like to use them beyond 2^64.

Benchmark Math.NumberTheory.GCD

There are two modules, implementing very low-level algorithms for GCD: Math.NumberTheory.GCD and Math.NumberTheory.GCD.LowLevel. I am wondering whether they are still faster than gcd in modern base (cf. #51 about Math.NumberTheory.Powers.Integer). It would be nice to benchmark them with criterion.

Namely:

1. Is binaryGCD faster than Prelude.gcd?
2. Is coprime faster than Prelude.gcd?
3. Is extendedGCD faster than GHC.Integer.GMP.Internals. gcdExtInteger?
4. Or maybe are their specialised and unboxed counterparts from Math.NumberTheory.GCD.LowLevel worth a trouble?

Dirichlet beta function

Dirichlet beta function is an important analogue of Riemann zeta function, implemented in #44. Beta function is especially handy when it comes to analytic number theory over Gaussian integers.

It would be nice to compute its values at positive integer points, similar to Riemann zeta function in Math.NumberTheory.Zeta.

Wrap output of splitIntoCoprimes into newtype

Currently splitIntoCoprimes guarantees certain properties of its output, but still returns a plain [(a, b)]. Such return type is easily corruptible by other functions. It would be much better to return something like

newtype Coprimes a b = Coprimes (Map a b)

The interface should not expose Coprimes constructor outside. Here is a quick draft of API:

instance Semigroup (Coprimes a b) where ...

instance Monoid (Coprimes a b) where ...

insert :: a -> b -> Coprimes a b -> Coprimes a b 

toList :: Coprimes a b -> [(a, b)]

---

Another vaguely related idea is to keep coprimes, which are known to be prime, separately. This can potentially speed up certain computations, especially with Prefactored. Like

data Coprimes a b = Coprimes (Map (Prime a) b) (Map a b)

test-suite/Math/NumberTheory/TestUtils.hs: Duplicate instance declarations

In Stackage build:

test-suite/Math/NumberTheory/TestUtils.hs:65:10: error:
    Duplicate instance declarations:
      instance Monad m => Serial m Word
        -- Defined at test-suite/Math/NumberTheory/TestUtils.hs:65:10
      instance [overlap ok] Monad m => Serial m Word
        -- Defined in ‘Test.SmallCheck.Series’
   |
65 | instance Monad m => Serial m Word where
   |          ^^^^^^^^^^^^^^^^^^^^^^^^

test-suite/Math/NumberTheory/TestUtils.hs:75:10: error:
    Duplicate instance declarations:
      instance Monad m => Serial m Natural
        -- Defined at test-suite/Math/NumberTheory/TestUtils.hs:75:10
      instance [overlap ok] Monad m => Serial m Natural
        -- Defined in ‘Test.SmallCheck.Series’
   |
75 | instance Monad m => Serial m Natural where
   |          ^^^^^^^^^^^^^^^^^^^^^^^^^^^

wordLog2 slower than naive approach

wordLog2:

import Math.NumberTheory.Logarithms

main = print . last . take 100000000 . map wordLog2 $ [1..]

$ time ./wordlog
26
./wordlog  2.87s user 0.00s system 99% cpu 2.870 total

my fastLog2:

import Data.Bits

fastLog2 :: Int -> Int
fastLog2 1 = 0
fastLog2 n = 1 + fastLog2 (shift n (-1))

main = print . last . take 100000000 . map fastLog2 $ [1..]

$ time ./fastlog
26
./fastlog  0.08s user 0.00s system 96% cpu 0.086 total

Is this a bug?

New Chinese Remainder Module?

Congratulations on 0.6.0.0! So much good stuff, so much stuff which I had already written and planned to contribute but never got around to!

Lest more of my efforts be wasted, let me start by offering an alternative implementation of the Chinese Remainder module. It differs from the current implementation in the following points:

1. It uses the M.NT.Moduli.Class types which makes the interface more consistent and generally neater. For example, the signature of chineseRemainder is [SomeMod] -> Maybe SomeMod. So, for example, chineseRemainder [3 \modulo` 4, 7 `modulo` 9, 0 `modulo` 5]' evaluates to Just (115 \modulo` 180)`. No longer mixing up which part of the pairs is the modulus, which the remainder, etc.

2. It correctly deals with InfMod constructor for SomeMod It is interpreted to represent the congruence class consisting of just a single rational (I am not sure I understand the reasoning for having a Rational, rather than Integer representation?). So, for example, chineseRemainder [3 \modulo` 4, 7 `modulo` 9, 0 `modulo` 5, 295] == Just 295 % 1andchineseRemainder [3 `modulo` 4, 7 `modulo` 9, 0 `modulo` 5, 296] == Nothing`

3. Most importantly for my usage cases, it correctly handles the cases where the moduli are not pairwise co-prime, rather than always returning Nothing, even if the congruences could be fulfilled: chineseRemainder [3 \modulo` 4, 7 `modulo` 6] = Just (7 `modulo` 12)` I am not positive that the algorithm is optimal, but after running many thousands of quickChecks, I'm pretty sure it is correct. As it only ever comes into play when the current implementation would just fail, I think even a suboptimal algorithm (if that is what it is), is an improvement.

import Control.Monad
import Data.Ratio
import Math.NumberTheory.Moduli.Class
import Math.NumberTheory.Primes.Factorisation

chineseRemainder :: [SomeMod] -> Maybe SomeMod
chineseRemainder (x:xs) = foldM chineseRemainder' x xs
chineseRemainder _ = Just (0 `modulo` 1)

chineseRemainder' :: SomeMod -> SomeMod -> Maybe SomeMod
chineseRemainder' x@(SomeMod x') y@(SomeMod y')
  | getMod x' == 1        = Just y
  | getMod y' == 1        = Just x
  | Just (SomeMod j) <- i = Just $ (xv+(yv-xv)*xm*getVal j) `modulo` (xn*yn)
  | (xv-yv) `mod` gm == 0 = chineseRemainder' (xv `modulo` xq) (yv `modulo` yq)
  | otherwise             = Nothing
  where
    xv = getVal x'
    xm = getMod x'
    xn = getNatMod x'

    yv = getVal y'
    ym = getMod y'
    yn = getNatMod y'

    i = invertSomeMod (xm `modulo` yn)

    gn = gcd xn yn
    gm = fromIntegral gn

    (xq, yq) = foldr enrich (xn `div` gn, yn `div` gn) $ factorise gm
    enrich (p',a) (xo, yo)
      | xo `mod` p == 0 = (xo*p^a, yo)
      | otherwise       = (xo, yo*p^a)
      where
        p = fromIntegral p'

chineseRemainder' x@(SomeMod x') y@(InfMod y')
  | denominator y' /= 1                       = Nothing
  | numerator y' `mod` getMod x' == getVal x' = Just y
  | otherwise                                 = Nothing

chineseRemainder' x@(InfMod _) y@(SomeMod _) = chineseRemainder' y x

chineseRemainder' x y
  | x == y    = Just x
  | otherwise = Nothing

prop_chineseRemainder' :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool
prop_chineseRemainder' xv (Positive xm) yv (Positive ym)
  | Just (SomeMod d) <- c = getMod d==l && [getVal d] == sols
  | otherwise             = null sols
  where
    l = lcm xm ym
    c = chineseRemainder' (xv `modulo` (fromIntegral xm)) (yv `modulo` (fromIntegral ym))
    sols = [i|i<-[0..l-1],(i-xv) `mod` xm==0,(i-yv) `mod` ym==0]

Euclidean rings

This is to avoid messy div{E,G}, mod{E,G}, gcd{E,G} (cf. the implementation of Gaussian integers and #121) and unify them with div, mod, gcd for usual Integral types.

class Euclidean a where
  quotRem :: a -> a -> (a, a)
  divMod  :: a -> a -> (a, a)

quot :: Euclidean a => a -> a -> a
quot = ...

rem :: Euclidean a => a -> a -> a
rem = ...

div :: Euclidean a => a -> a -> a
div = ...

mod :: Euclidean a => a -> a -> a
mod = ...

instance Integral a => Euclidean a where 
  ...

instance Euclidean GaussianInteger where
  ...

instance Euclidean EisensteinInteger where
  ...

Then we can replace Integral constraint with Euclidean in Math.NumberTheory.GCD, Math.NumberTheory.Prefactored and Math.NumberTheory.SmoothNumbers, making them both more general and applicable to quadratic fields.

https://en.wikipedia.org/wiki/Euclidean_domain

Spurious 1 factor reported for large repUnit

With version 0.6.0.1 of arithmoi (I'm using stack LTS), given this function:

repUnit :: Integer -> Integer
repUnit k = (10^k - 1) `div` 9

I get this result with factorise:

*Main Protolude> factorise (repUnit 80)
[(11,1),(17,1),(41,1),(73,1),(101,1),(137,1),(271,1),(3541,1),(9091,1),(27961,1),(1,1),(1676321,1),(5070721,1),(5882353,1),(5964848081,1),(19721061166646717498359681,1)]

You can see 1 is reported as a factor.

Sum up a function over primes

Math.NumberTheory.Primes.Counting.Impl.primeCount implements an algorithm to compute π(n), which is number of primes below n, in O(n^0.7) time and O(n^1/2) space. One can view π(n) as a sum of f p = 1 over prime p below n.

As it was pointed by @CarlEdman in #60 (primeLucy), this approach allows a generalization for any function f with the same time and space requirements. Moreover, there is a (difficult) way to improve its time complexity to O(n^2/3).

That said, it seems desirable to productionalize primeLucy, provide a nice API, write decent tests and benchmarks, and then express primeCount as a special case.

Additional resources:
https://en.wikipedia.org/wiki/Meissel–Lehmer_algorithm
http://sweet.ua.pt/tos/bib/5.4.pdf
https://projecteuler.net/thread=10;page=5#111677

Primitive roots

I discovered that arithmoi entirely lacks routines to deal with primitive roots modulo n. This gap should be filled, however I have not yet come up with a convenient and effective API.

The difficulty is that finding primitive roots depends not only on modulo n and its factorisation (this part can be dealt with layer by layer as in Math.NumberTheory.Moduli.Sqrt), but also on factorisation of totient n. There are use cases when one may have these factorisations beforehand and do not want them to be recalculated, and there are cases when one just wanna send n in and get an answer.

Further idea is to implement discrete logarithms, but this is a huge lump of work.

Split Math.NumberTheory.Logarithms into separate package

So people can depend on it. Unfortunately for some arithmoi seems to be "heavy" dependency.
E.g. scientific has an internal copy.

I could setup the package (and transfer to you / be co-maintainer, if you want), and make PR here to use it, if the proposal is ok.

Add support for `integer-gmp-1.0.0`

see also #3

Conditions on arguments of isStrongFermatPP

Hi,

isStrongFermatPP 4 9 returns True, meaning base 9 is a "strong liar" for 4
I discovered this accidentally today. I checked the definition and the algorithm seems to only work properly for odd integers > 3.

For n=1: It hung.
For n=2 and 3: It returned False for some bases

I suggest there be a hidden pre-condition function called within
isStrongFermatPP and millerRabinV (and perhaps the "new" function)
Then, an underlying function (isStrongFermatPP') would be called if it's a suitable candidate for Miller-Rabin. Also, isPrime works for negative or positive integers. I think isStrongFermatPP and related should as well.

-Fred

Summatory functions

Common difficulty of number theory is that arithmetic functions misbehave: their plots oscillate up and down in chaotic fashion. For instance, here is a plot of divisor function from wiki:

The pivot idea of analytic number theory is that instead of studying this jumble, we can study the summatory function, which is simply a sum (or an integral) of original function:

summatory :: (Int -> Int) -> Int -> Int 
summatory f n = sum $ map f [1..n]

Summatory functions behave well and a wide range of fruitful methods can be applied. E. g., divisor summatory function.

The interesting thing is that in many cases summatory function can be computed in sublinear time, faster than O(n). For example, let tau denote the divisor function. Then summatory tau n is equal to

2 * sum (map (\m -> x `quot` m) [1..sqrtx]) - sqrtx ^ 2

where sqrtx = integerSquareRoot x. This computation can be completed in O(x^{1/2}) time. There is a systematic way to derive such formulas, which I have explored in my old paper with a long title.

The proposal is to implement some of these algorithms, to provide users with a DSL to describe functions (in terms of Dirichlet convolution?) and to derive effective summatory formulas from DSL automatically.

Upd.: I think it is worse to create a separate type for multiplicative functions, accompanied by conversion multiplicative :: MultiplicativeFunction -> ArithmeticFunction, because there is so much more to do with multiplicative functions, which is undefined for general arithmetic ones. It would be nice to implement inverses of multiplicative functions as well, following the approach in Computing the inverses....

Implementation plan

• Function to evaluate arithmetic functions over intervals in quasilinear time (done in #77)
• Implement fast (O(x^2/3) time) summation of the Mobius function as per paper (done in #90)
• Implement fast (O(x^1/3) time) summation of the divisor function as per paper. Help wanted
• Grab Faulhaber polynomials from #70.
• Provide DSL to combine functions by Dirichlet convolution, deriving optimal summatory functions automatically as per paper.
• Review implementation of Math.NumberTheory.MoebiusInversion.totientSum.
• Compute inverses of multiplicative functions as per paper (done in #142).

Update license file / copyright field

It currently says (c) 2011 Daniel Fischer

Stackage

What about publishing this package on Stackage?

arithmoi-0.4.3.0 test suite does not compile on 32 bit platforms

Citing from http://hydra.nixos.org/build/49692523/nixlog/1/raw:

test-suite/Math/NumberTheory/Powers/GeneralTests.hs:98:9: warning: [-Woverflowed-literals]
    Literal 118372752099 is out of the Int range -2147483648..2147483647

test-suite/Math/NumberTheory/Powers/GeneralTests.hs:99:9: warning: [-Woverflowed-literals]
    Literal 118372752099 is out of the Int range -2147483648..2147483647
[12 of 23] Compiling Math.NumberTheory.Powers.FourthTests ( test-suite/Math/NumberTheory/Powers/FourthTests.hs, dist/build/spec/spec-tmp/Math/NumberTheory/Powers/FourthTests.dyn_o )

test-suite/Math/NumberTheory/Powers/FourthTests.hs:136:52: error:
    Variable not in scope:
      integerFourthRootSpecialCase1_Int :: Assertion

test-suite/Math/NumberTheory/Powers/FourthTests.hs:137:52: error:
    Variable not in scope:
      integerFourthRootSpecialCase1_Word :: Assertion

test-suite/Math/NumberTheory/Powers/FourthTests.hs:138:52: error:
    Variable not in scope: integerFourthRootSpecialCase2 :: Assertion
[13 of 23] Compiling Math.NumberTheory.Powers.CubesTests ( test-suite/Math/NumberTheory/Powers/CubesTests.hs, dist/build/spec/spec-tmp/Math/NumberTheory/Powers/CubesTests.dyn_o )

test-suite/Math/NumberTheory/Powers/CubesTests.hs:145:51: error:
    Variable not in scope: integerCubeRootSpecialCase1_Int :: Assertion

test-suite/Math/NumberTheory/Powers/CubesTests.hs:146:51: error:
    Variable not in scope:
      integerCubeRootSpecialCase1_Word :: Assertion

test-suite/Math/NumberTheory/Powers/CubesTests.hs:147:51: error:
    Variable not in scope: integerCubeRootSpecialCase2 :: Assertion
~~~

once random v2 is out, make sure arithmoi can use it/ benefits from it

millerRabinV performance enhancement?

millerRabinV : This routine in Math/NumberTheory/Primes/Testing/Probabilistic.hs
seems quite inefficient. According to the comments in the code, verifying a prime is actually "prime" with this routine would require testing 196 bases (primes through 1200).

In the last few years, there has been a project: https://miller-rabin.appspot.com/
to identify the maximum accurate solution for # of bases: 1 - 7.
A set of 7 bases exists (2, 325, 9375, 28178, 450775, 9780504, 1795265022)
that can verify the primality of a number <= 2^64.

Secondly, I wonder if this logic could be used to improve the performance of isPrime for numbers < 2^64. Not sure what the cutoff would be. http://www.trnicely.net/misc/bpsw.html

[Edit: Fixed a few typos]

curveFactorisaton : quirk with selecting "parms"

Hi,

In the internal sub-function "fact" in curveFactorisation in Math/NumberTheory/Primes/Factorisation/Montgomery.hs, there's this line:

mult j <$> fact k (if null pfs then digs+4 else digs)

digs = # of digits and is passed to findParms to return B1, B2 and curve count.
Looks a little awkward, I think it should be incrementing it so that it grabs the next tier. (For instance, jumping from 15 to 20).

When factoring a number, it looks like you would see this pattern due to the varying differences between tiers:
8, 12, 16, 20, 24(), 28, 32, 36, 40, 44 () -- num digits checked in findParms
12, 15, 20, 25, 25, 30, 35, 40, 45, 45 -- "tier" found in testParams

• means that that digit-tier would be executed twice which is clearly incorrect.

For a quick fix, we could change the hard-coded 4
to a conditional: (if digs <= 16 then 4 else 5)

This "jump of 4" logic is also present in the certiFactorisation function.

Elliptic curve testing (using GMP-ECM's B2 and curve counts)

Hi,

I noticed that GMP-ECM (and other code bases apparently) uses different B2 and curve count values:
http://www.mersennewiki.org/index.php/Elliptic_Curve_Method

Given it's GMP-ECM, I would assume these figures are the best known.

(The B2 and curve counts are the 3rd and 4th members of each tuplet in testParms list in Math/NumberTheory/Primes/Factorisation/Montgomery.hs)

The arithmoi code has an entry for 12 (digit factors). The GMP-ECM tables do not. Is the 12 entry really needed? Oddly, there is "no" optimal number of curves for 15 in the GMP-ECM. Not sure how to interpret that.

Do you think it's worth changing the arithmoi code?