arithmoi 0.13.3.0 → 0.13.4.0

raw patch · 16 files changed

+87/−28 lines, 16 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Math/NumberTheory/ArithmeticFunctions/Class.hs +2/−0
• Math/NumberTheory/ArithmeticFunctions/Mertens.hs +1/−1
• Math/NumberTheory/ArithmeticFunctions/Moebius.hs +1/−1
• Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs +3/−1
• Math/NumberTheory/Diophantine.hs +19/−10
• Math/NumberTheory/DirichletCharacters.hs +1/−1
• Math/NumberTheory/Moduli/Cbrt.hs +1/−1
• Math/NumberTheory/Moduli/Singleton.hs +1/−1
• Math/NumberTheory/Primes/Counting/HowPrimeCountingWorks.md +8/−8
• Math/NumberTheory/Primes/Testing/Certified.hs +1/−1
• Math/NumberTheory/Quadratic/EisensteinIntegers.hs +1/−1
• Math/NumberTheory/Quadratic/GaussianIntegers.hs +15/−0
• arithmoi.cabal +1/−1
• changelog.md +15/−0
• test-suite/Math/NumberTheory/ArithmeticFunctions/InverseTests.hs +1/−1
• test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs +16/−0

@@ -56,6 +56,8 @@     = ArithmeticFunction (\_ _ -> ()) (const x)   (ArithmeticFunction f1 g1) <*> (ArithmeticFunction f2 g2)     = ArithmeticFunction (\p k -> (f1 p k, f2 p k)) (\(a1, a2) -> g1 a1 (g2 a2))+  liftA2 h (ArithmeticFunction f1 g1) (ArithmeticFunction f2 g2)+    = ArithmeticFunction (\p k -> (f1 p k, f2 p k)) (\(a1, a2) -> h (g1 a1) (g2 a2))  instance Semigroup a => Semigroup (ArithmeticFunction n a) where   (<>) = liftA2 (<>)

@@ -19,7 +19,7 @@ import Math.NumberTheory.ArithmeticFunctions.Moebius import Math.NumberTheory.Utils.FromIntegral --- | Compute individual values of Mertens function in O(n^(2/3)) time and space.+-- | Compute individual values of Mertens function in \( O(n^{2/3}) \) time and space. -- -- >>> map (mertens . (10 ^)) [0..9] -- [1,-1,1,2,-23,-48,212,1037,1928,-222]

@@ -40,7 +40,7 @@  -- | Represents three possible values of <https://en.wikipedia.org/wiki/Möbius_function Möbius function>. data Moebius-  = MoebiusN -- ^ -1+  = MoebiusN -- ^ \-1   | MoebiusZ -- ^  0   | MoebiusP -- ^  1   deriving (Eq, Ord, Show)

@@ -150,7 +150,9 @@               npHigh = highIndex'      `shiftR` 1           forM_ [npLow .. npHigh] $ \np@(W# np#) -> do             let ix = wordToInt (np `shiftL` 1) - lowIndex :: Int-                tz = I# (word2Int# (ctz# np#))+                -- Calling ctz# assumes that np /= 0, otherwise you get 64+                -- and can have an out-of-bounds read from 'fs'.+                tz = if np == 0 then 0 else I# (word2Int# (ctz# np#))             MU.unsafeModify as (\x -> x `shiftL` (tz + 1)) ix             MG.unsafeModify bs (\y -> y `append` V.unsafeIndex fs tz) ix 

@@ -18,7 +18,7 @@ import           Math.NumberTheory.Utils.FromIntegral  -- | See `cornacchiaPrimitive`, this is the internal algorithm implementation--- | as described at https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm +-- as described at https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm cornacchiaPrimitive' :: Integer -> Integer -> [(Integer, Integer)] cornacchiaPrimitive' d m = concatMap   (findSolution . Inf.head . Inf.dropWhile (\r -> r * r >= m) . gcdSeq m)@@ -37,16 +37,21 @@     (s2, rem1) = divMod (m - r * r) d     s          = integerSquareRoot s2 --- | Finds all primitive solutions (x,y) to the diophantine equation --- |    x^2 + d*y^2 = m--- | when 1 <= d < m and gcd(d,m)=1--- | Given m is square free these are all the positive integer solutions+-- | @cornacchiaPrimitive d m@ finds all primitive solutions \((x, y)\) (i.e. with \(\gcd(x, y) = 1\))+-- to the Diophantine equation+-- \[ x^2 + d \cdot y^2 = m \]+-- Preconditions: \(1 \le d < m\) and \(\gcd(d, m) = 1\).+--+-- Throws error if the preconditions are not met.+--+-- When \(m\) is square-free these are all the positive integer solutions;+-- use 'cornacchia' to find all solutions for arbitrary \(m\). cornacchiaPrimitive :: Integer -> Integer -> [(Integer, Integer)] cornacchiaPrimitive d m   | not (1 <= d && d < m) = error "precondition failed: 1 <= d < m"   | gcd d m /= 1          = error "precondition failed: d and m coprime"   |-  -- If d=1 then the algorithm doesn't generate symmetrical pairs +  -- If d=1 then the algorithm doesn't generate symmetrical pairs     d == 1                = concatMap genPairs solutions   | otherwise             = solutions  where@@ -60,10 +65,14 @@   squareProducts acc f = [ a * b | a <- acc, b <- squarePowers f ]   squarePowers (p, a) = map (unPrime p ^) [0 .. wordToInt a `div` 2] --- | Finds all positive integer solutions (x,y) to the--- | diophantine equation:--- |    x^2 + d*y^2 = m--- | when 1 <= d < m and gcd(d,m)=1+-- | @cornacchia d m@ finds all positive integer solutions \((x, y)\) to the Diophantine equation+-- \[ x^2 + d \cdot y^2 = m \]+-- Preconditions: \(1 \le d < m\) and \(\gcd(d, m) = 1\).+--+-- Throws error if the preconditions are not met.+--+-- Unlike 'cornacchiaPrimitive', this also finds non-primitive solutions (where \(\gcd(x, y) > 1\))+-- by solving the equation for each square divisor of \(m\) and scaling the results. cornacchia :: Integer -> Integer -> [(Integer, Integer)] cornacchia d m   | not (1 <= d && d < m) = error "precondition failed: 1 <= d < m"

@@ -307,7 +307,7 @@           where p' = unPrime p                 m = p'^(k-1)*(p'-1) --- the validity of the producted dirichletfactor list here requires the template to be valid+-- the validity of the produced dirichletfactor list here requires the template to be valid unroll :: [Template] -> Natural -> [DirichletFactor] unroll t m = snd (mapAccumL func m t)   where func :: Natural -> Template -> (Natural, DirichletFactor)

@@ -121,7 +121,7 @@ -- Otherwise, cubic reciprocity is called. cubicReciprocity alpha beta = cubicSymbolHelper beta alpha --- | This function takes two Eisenstein intgers @alpha@ and @beta@ and returns+-- | This function takes two Eisenstein integers @alpha@ and @beta@ and returns -- three arguments @(gamma, delta, newSymbol)@. @gamma@ and @delta@ are the -- associated primary numbers of alpha and beta respectively. @newSymbol@ -- is the cubic symbol measuring the discrepancy between the cubic residue

@@ -227,7 +227,7 @@  -- | Similar to 'cyclicGroupFromFactors' . 'factorise', -- but much faster, because it--- but performes only one primality test instead of full+-- but performs only one primality test instead of full -- factorisation. cyclicGroupFromModulo   :: (Integral a, UniqueFactorisation a)

@@ -45,7 +45,7 @@ In general, the method is as follows:  1. All arrays are of size the square root of the counting range limit, as follows:-   +     a. A `roughs` array to process the remaining rough numbers in culling waves, initialized to all odds starting from 1 (ie. 1, 3, ...),      b. A `phis` array which contains counts of remaining values after culling up to some level of base prime but is used for emulating the stack during processing initialized with the values as culled by two as in ((`limit`/roughs[`i`] + 1) `div` 2) where `i` is the zero-base index, and@@ -77,8 +77,8 @@                     (acc + phip2(limit / npmult) - phi npmult pi) -- recursion    ``` -   Phi(limit (Pi(sqrt(limit))) is the above function called with:   -   phi 1 (lenth paseprms).+   Phi(limit (Pi(sqrt(limit))) is the above function called with:+   phi 1 (length paseprms).     Using the recursive function is not efficient because it doesn't have the asymptotic improvements to performance due to partial sieving. @@ -89,11 +89,11 @@      b. this phase processes each of the remaining `roughs` values in turn for all remaining after being marked and splits the processing depending on whether the base prime value multiplied by the remaining rough element is:      - less than or equal to the square root limit, in which case it compresses the `roughs` and `pis` array by subtracting the next `pis` element from the current one and writing it into an index that doesn't included the marked/skipped elements, or-  -    - if greater than it subtracts a value looked up by rough value index in the `pisndxs` array and again writes the result into the compression index as per above.  -  ++    - if greater than it subtracts a value looked up by rough value index in the `pisndxs` array and again writes the result into the compression index as per above.+     - in all cases it moves the `roughs` to synchronizes with the pi's indices to move values to the left to fill the holes left by eliminating the marked rough values.-  +      c. the last phase is to adjust pisndxs values so they address the correct rough/pis values after the above compression of these arrays has taken place.  3. The partial sieving loop only runs the base prime values up to the square root of the square root of the limit; after that point all of the roughs are prime (other than one) and therefore, the `limit` divided by each of the unique products of pairs of these roughs are added to the final answer from the loop with the usual procedure of looking up the values to be added from the `pisndx` array by using as an index: toIndex (limit / p1 / p2) but with the limitation that p1 can not be larger than (limit / p1 / p1) which satisfies the terminating condition that `p1` must be less than or equal to the quotient:@@ -106,7 +106,7 @@  In more detail, the "bottom up" partial sieving method is as follows; this implementeds the odds-only optimization so the arrays represent only the odd numbers: -1. Keep track of the number of base primes that have been processed yet and the current effective number of rough numbers; INITIALIAZE `y` as the square root of the limit and three arrays each of the size `y`, with the first small counts array keeping track of the current number of possible primes up to the number represented by the index value so initialized with each element containing the index value, the second the roughs array containing the current level zero rough numbers as per the current culling passes (starts with nothing culled, so all odd values could be prime), and finally, the third containing the current Phi's/ (Pi's in this case) for each rough number in the rough numbers array (initially limit divided by the rough value corresponding to the index).+1. Keep track of the number of base primes that have been processed yet and the current effective number of rough numbers; INITIALIZE `y` as the square root of the limit and three arrays each of the size `y`, with the first small counts array keeping track of the current number of possible primes up to the number represented by the index value so initialized with each element containing the index value, the second the roughs array containing the current level zero rough numbers as per the current culling passes (starts with nothing culled, so all odd values could be prime), and finally, the third containing the current Phi's/ (Pi's in this case) for each rough number in the rough numbers array (initially limit divided by the rough value corresponding to the index).  2. For each base prime in order from the lowest to the limit^(1/4) (same as the square root of `y`) do the following: 

@@ -206,7 +206,7 @@                     , "The number in question is:\n"                     , show n                     , "\nOther parties like wikipedia might also be interested."-                    , "\nSorry for aborting your programm, but this is a major discovery."+                    , "\nSorry for aborting your program, but this is a major discovery."                     ]  -- | Found a factor

@@ -228,7 +228,7 @@ --        where @a, b, c, a_i@ are nonnegative integers, @N > 1@ is an integer and --        @π_i@ are Eisenstein primes. ----- Aplying @norm@ to both sides of the equation from Theorem 8.4:+-- Applying @norm@ to both sides of the equation from Theorem 8.4: -- -- 1. @norm μ = norm ( (-1)^a * ω^b * (1 - ω)^c * product [ π_i^a_i | i <- [1..N]] ) ==@ -- 2. @norm μ = norm ((-1)^a) * norm (ω^b) * norm ((1 - ω)^c) * norm (product [ π_i^a_i | i <- [1..N]]) ==@

@@ -17,6 +17,8 @@     ι,     conjugate,     norm,+    ids,+    associates,     primes,     findPrime, ) where@@ -93,6 +95,19 @@   | a <  0 && b <= 0 = (-z, -1)   -- fourth quadrant: [0, inf) x (-inf, 0)i   | otherwise        = ((-b) :+ a, -ι)++-- | List of all Gaussian units, counterclockwise across all quadrants,+-- starting with @1@.+--+-- @since 0.13.4.0+ids :: [GaussianInteger]+ids = [1, ι, -1, -ι]++-- | Produce a list of a @GaussianInteger@'s associates.+--+-- @since 0.13.4.0+associates :: GaussianInteger -> [GaussianInteger]+associates (a :+ b) = [a :+ b, (-b) :+ a, (-a) :+ (-b), b :+ (-a)]  instance GcdDomain GaussianInteger 

@@ -1,5 +1,5 @@ name:          arithmoi-version:       0.13.3.0+version:       0.13.4.0 cabal-version: 2.0 build-type:    Simple license:       MIT

@@ -1,5 +1,20 @@ # Changelog +## 0.13.4.0++### Fixed++* Fix `sieveBlockUnboxed` segfaulting when a block starts from 0.+* Fix typos.++### Added++* Add `ids` and `associates` to `Math.NumberTheory.Quadratic.GaussianIntegers`.++### Changed++* Define `liftA2` for `ArithmeticFunction` explicitly.+ ## 0.13.3.0  ### Changed

@@ -363,7 +363,7 @@ -- TestTree  -- Tests for 'Int', 'Word' are omitted because 'inverseSigmaK/inverseJordan'--- tests would quickly oveflow in these types.+-- tests would quickly overflow in these types. testIntegralPropertyNoLargeInverse   :: forall bool. (SC.Testable IO bool, QC.Testable bool)   => String -> (forall a. (Euclidean a, Semiring a, Integral a, Bits a, UniqueFactorisation a, Show a, Enum (Prime a)) => Power Word -> Positive a -> bool) -> TestTree

@@ -19,9 +19,11 @@ import Test.Tasty.HUnit  import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U  import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.ArithmeticFunctions.SieveBlock+import Math.NumberTheory.Primes (Prime(..))  pointwiseTest :: (Eq a, Show a) => ArithmeticFunction Word a -> Word -> Word -> IO () pointwiseTest f lowIndex len = assertEqual "pointwise"@@ -67,6 +69,19 @@       _ -> MoebiusZ   } +doesNotSegfaultOnZero :: IO ()+doesNotSegfaultOnZero = assertBool "should not segfault" $ xs == xs+  where+    xs = U.sum $ sieveBlockUnboxed sigmaConfig 0 5++sigmaConfig :: SieveBlockConfig Int+sigmaConfig = SieveBlockConfig+  { sbcEmpty = 1+  , sbcAppend = (*)+  , sbcFunctionOnPrimePower =+      \p n -> fromIntegral $ (unPrime p ^ (n+1) - 1) `quot` (unPrime p - 1)+  }+ testSuite :: TestTree testSuite = testGroup "SieveBlock"   [ testGroup "pointwise"@@ -79,4 +94,5 @@     , testCase "carmichael" $ pointwiseTest carmichaelA 1 1000     ]   , testGroup "special moebius" moebiusSpecialCases+  , testCase "does not segfault at 0" doesNotSegfaultOnZero   ]