bodigrim / arithmoi Goto Github PK
View Code? Open in Web Editor NEWNumber theory: primes, arithmetic functions, modular computations, special sequences
Home Page: http://hackage.haskell.org/package/arithmoi
License: MIT License
Number theory: primes, arithmetic functions, modular computations, special sequences
Home Page: http://hackage.haskell.org/package/arithmoi
License: MIT License
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.
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 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.
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.
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:
Prime Integer
instead of Integer
.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.
I had problems when compiling the package, getting this error for a few files.
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 :)
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.
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
.
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:
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
)
Add {-# SPECIALISE #-}
pragmas to cover all four "first-class" integer types. This is the easiest, almost mechanical part of the task.
Add rewrite rules {-# RULES #-}
and provide corresponding functions, benefiting from Natural
arguments, where appropriate. May be easy, may be hard, may be superfluous.
At the moment there are two functions to list divisors of an integer:
Set n
,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.
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.
There are three kinds of Gaussian primes (see Math.NumberTheory.GaussianIntegers.primes):
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
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
.
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.
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
.
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.
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.
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 }
Define Num (QuadraticInteger d)
instance. This will require definitions of conjugation and norm as well.
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.
Define UniqueFactorisation
instance for principal rings.
Define Euclidean
instance for Euclidean rings.
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 theNumberSieves
package. I have no idea about how their performance compares, theoretically, but anything is better than code that's blatantly broken.
AFAICS the master
branch builds ok with GHC-8.0.1. The Hackage version has restrictive bounds though (ghc-prim < 0.5
).
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.
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]
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.
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. 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!
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:
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.
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:
binaryGCD
faster than Prelude.gcd
?coprime
faster than Prelude.gcd
?extendedGCD
faster than GHC.Integer.GMP.Internals. gcdExtInteger
?Math.NumberTheory.GCD.LowLevel
worth a trouble?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.
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)
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:
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?
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:
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.
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 % 1and
chineseRemainder [3 `modulo` 4, 7 `modulo` 9, 0 `modulo` 5, 296] == Nothing`
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]
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.
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.
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
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.
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.
see also #3
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
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
Math.NumberTheory.MoebiusInversion.totientSum
.It currently says (c) 2011 Daniel Fischer
What about publishing this package on Stackage?
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
~~~
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]
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
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.
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?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.