!      !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVW X Y Z [ \ ] ^ _ ` a b c d e f g h i j k l m n o p q r s tuvwxyz{|}~       !"#$%&'()*+,-./0123456789:;<=>?@ABC D E F G H I J K L M N O P QRSTUVWXYZ[\]^_`abcdefghijklmnopqr!s!t!u!v!w"x#y#z#{#|#}#~####$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%&&&&&B'Safe"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>None=?FMX# arithmoi/Check whether an element is a unit of the ring.arithmoidCalculate the greatest common divisor of two numbers and coefficients for the linear combination.For signed types satisfies: Tcase extendedGCD a b of (d, u, v) -> u*a + v*b == d && d == gcd a bCFor unsigned and bounded types the property above holds, but since u and vS must also be unsigned, the result may look weird. E. g., on 64-bit architecture 5extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1)'For unsigned and unbounded types (like ()) the result is undefined.For signed types we also have 8abs u < abs b || abs b <= 1 abs v < abs a || abs a <= 1(except if one of a and b is  of a signed type).      (c) 2017-2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NonePX5yarithmoi>A list of pairwise coprime numbers with their multiplicities.arithmoiUnwrap.arithmoi2Wrap a non-zero number with its multiplicity into .singleton 210 1!Coprimes {unCoprimes = [(210,1)]}arithmoi/Add a non-zero number with its multiplicity to .insert 360 1 (singleton 210 1)1Coprimes {unCoprimes = [(7,1),(5,2),(3,3),(2,4)]}4insert 2 4 (insert 7 1 (insert 5 2 (singleton 4 3))),Coprimes {unCoprimes = [(7,1),(5,2),(2,10)]}arithmoiThe input list is assumed to be a factorisation of some number into a list of powers of (possibly, composite) non-zero factors. The output list is a factorisation of the same number such that all factors are coprime. Such transformation is crucial to continue factorisation (lazily, in parallel or concurrent fashion) without having to merge multiplicities of primes, which occurs more than in one composite factor.&splitIntoCoprimes [(140, 1), (165, 1)]-Coprimes {unCoprimes = [(28,1),(33,1),(5,2)]}&splitIntoCoprimes [(360, 1), (210, 1)]1Coprimes {unCoprimes = [(7,1),(5,2),(3,3),(2,4)]}(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None&'.1FHMSX_farithmoiThis type represents residues with unknown modulo and rational numbers. One can freely combine them in arithmetic expressions, but each operation will spend time on modulo's recalculation:2 `modulo` 10 + 4 `modulo` 15(1 `modulo` 5)!(2 `modulo` 10) * (4 `modulo` 15)(3 `modulo` 5)$2 `modulo` 10 + fromRational (3 % 7)(1 `modulo` 10)$2 `modulo` 10 * fromRational (3 % 7)(8 `modulo` 10)8If performance is crucial, it is recommended to extract Mod m4 for further processing by pattern matching. E. g., fcase modulo n m of SomeMod k -> process k -- Here k has type Mod m InfMod{} -> error "impossible"arithmoirThis type represents elements of the multiplicative group mod m, i.e. those elements which are coprime to m. Use  toMultElement to construct.arithmoiUnwrap a residue.arithmoiWrapper for residues modulo m.Mod 3 :: Mod 10 stands for the class of integers, congruent to 3 modulo 10 ( &"17, "7, 3, 13, 23 &). The modulo is stored on type level, so it is impossible, for example, to add up by mistake residues with different moduli.:set -XDataKinds(3 :: Mod 10) + (4 :: Mod 12),error: Couldn't match type 12  with 10 ...(3 :: Mod 10) + 8(1 `modulo` 10)$Note that modulo cannot be negative.arithmoi.Linking type and value levels: extract modulo m as a value. arithmoi.Linking type and value levels: extract modulo m as a value.!arithmoiHThe canonical representative of the residue class, always between 0 and m-1 inclusively."arithmoiHThe canonical representative of the residue class, always between 0 and m-1 inclusively.#arithmoiHComputes the modular inverse, if the residue is coprime with the modulo.:set -XDataKindsinvertMod (3 :: Mod 10)3Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10invertMod (4 :: Mod 10)Nothing$arithmoiDrop-in replacement for , with much better performance.:set -XDataKindspowMod (3 :: Mod 10) 4(1 `modulo` 10)%arithmoiInfix synonym of $.&arithmoi4Attempt to construct a multiplicative group element.'arithmoitFor elements of the multiplicative group, we can safely perform the inverse without needing to worry about failure.(arithmoiCreate modular value by representative of residue class and modulo. One can use the result either directly (via functions from  and 5), or deconstruct it by pattern matching. Note that ( never returns .)arithmoi)Computes the inverse value, if it exists.invertSomeMod (3 `modulo` 10)3Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10invertSomeMod (4 `modulo` 10)Nothing$invertSomeMod (fromRational (2 % 5)) Just 5 % 2*arithmoiDrop-in replacement for `, with much better performance. When -O is enabled, there is a rewrite rule, which specialises  to *.powSomeMod (3 `modulo` 10) 4(1 `modulo` 10)+arithmoiuBeware that division by residue, which is not coprime with the modulo, will result in runtime error. Consider using # instead.2arithmoiuBeware that division by residue, which is not coprime with the modulo, will result in runtime error. Consider using ) instead. !"#$%&'()*!" #$%&'()*%8(7(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>Fz)=arithmoi.Calculate the integer cube root of an integer n!, that is the largest integer r such that r^3 <= n+. Note that this is not symmetric about 0, for example integerCubeRoot (-2) = (-2) while integerCubeRoot 2 = 1.>arithmoi9Calculate the integer cube root of a nonnegative integer n", that is, the largest integer r such that r^3 <= n. The precondition n >= 0 is not checked.?arithmoiReturns Nothing# if the argument is not a cube, Just r if n == r^3.@arithmoi"Test whether an integer is a cube.Aarithmoi8Test whether a nonnegative integer is a cube. Before = is calculated, a few tests of remainders modulo small primes weed out most non-cubes. For testing many numbers, most of which aren't cubes, this is much faster than  let r = cubeRoot n in r*r*r == n. The condition n >= 0 is not checked.Barithmoi}Test whether a nonnegative number is possibly a cube. Only about 0.08% of all numbers pass this test. The precondition n >= 0 is not checked.arithmoiHapproximate cube root, about 50 bits should be correct for large numbers=>?@AB=>?@AB(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>FOCarithmoi[Calculate the integer fourth root of a nonnegative number, that is, the largest integer r with r^4 <= n'. Throws an error on negaitve input.Darithmoi[Calculate the integer fourth root of a nonnegative number, that is, the largest integer r with r^4 <= n. The condition is not checked.EarithmoiReturns Nothing if n is not a fourth power, Just r if n == r^4 and r >= 0.FarithmoisTest whether an integer is a fourth power. First nonnegativity is checked, then the unchecked test is called.GarithmoiITest whether a nonnegative number is a fourth power. The condition is not$ checked. If a number passes the H0 test, its integer fourth root is calculated.HarithmoiRTest whether a nonnegative number is a possible fourth power. The condition is not6 checked. This eliminates about 99.958% of numbers.CDEFGHCDEFGH*(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None>F(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>Fk Iarithmoi:Calculate the integer square root of a nonnegative number n", that is, the largest integer r with r*r <= n'. Throws an error on negative input.Jarithmoi:Calculate the integer square root of a nonnegative number n", that is, the largest integer r with r*r <= n. The precondition n >= 0 is not checked.KarithmoiCalculate the integer square root of a nonnegative number as well as the difference of that number with the square of that root, that is if (s,r) = integerSquareRootRem n then s^2 <= n == s^2+r < (s+1)^2.LarithmoiCalculate the integer square root of a nonnegative number as well as the difference of that number with the square of that root, that is if (s,r) = integerSquareRootRem' n then s^2 <= n == s^2+r < (s+1)^2. The precondition n >= 0 is not checked.MarithmoiReturns % if the argument is not a square,  r if r*r == n and r >= 0<. Avoids the expensive calculation of the square root if n is recognized as a non-square before, prevents repeated calculation of the square root if only the roots of perfect squares are needed. Checks for negativity and P.NarithmoiXTest whether the argument is a square. After a number is found to be positive, first P@ is checked, if it is, the integer square root is calculated.Oarithmoi.Test whether the input (a nonnegative number) n is a square. The same as NW, but without the negativity test. Faster if many known positive numbers are tested.The precondition n >= 0J is not tested, passing negative arguments may cause any kind of havoc.ParithmoiTest whether a non-negative number may be a square. Non-negativity is not checked, passing negative arguments may cause any kind of havoc.First the remainder modulo 256 is checked (that can be calculated easily without division and eliminates about 82% of all numbers). After that, the remainders modulo 9, 25, 7, 11 and 13 are tested to eliminate altogether about 99.436% of all numbers.This is the test used by M@. For large numbers, the slower but more discriminating test Q is faster.QarithmoiTest whether a non-negative number may be a square. Non-negativity is not checked, passing negative arguments may cause any kind of havoc.First the remainder modulo 256 is checked (that can be calculated easily without division and eliminates about 82% of all numbers). After that, the remainders modulo several small primes are tested to eliminate altogether about 99.98954% of all numbers.JFor smallish to medium sized numbers, this hardly performs better than P, which uses smaller arrays, but for large numbers, where calculating the square root becomes more expensive, it is much faster (if the vast majority of tested numbers aren't squares). IJKLMNOPQ IJKLMNOPQ+(c) 2012 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>X^arithmoi totientSum n is, for n > 0 , the sum of otient k | k <- [1 .. n]]7, computed via generalised Mbius inversion. See  :http://mathworld.wolfram.com/TotientSummatoryFunction.html for the formula used for  totientSum.arithmoigeneralInversion g n/ evaluates the generalised Mbius inversion of g at the argument n.xThe generalised Mbius inversion implemented here allows an efficient calculation of isolated values of the function  f : N -> Z if the function g defined by , g n = sum [f (n `quot` k) | k <- [1 .. n]] Jcan be cheaply computed. By the generalised Mbius inversion formula, then 8 f n = sum [moebius k * g (n `quot` k) | k <- [1 .. n]]  which allows the computation in O(n) steps, if the values of the Mbius function are known. A slightly different formula, used here, does not need the values of the Mbius function and allows the computation in O(n^0.75) steps, using O(n^0.5) memory.}An example of a pair of such functions where the inversion allows a more efficient computation than the direct approach is w f n = number of reduced proper fractions with denominator <= n g n = number of proper fractions with denominator <= n (a proper fraction is a fraction  0 < n/d < 1). Then f n6 is the cardinality of the Farey sequence of order nr (minus 1 or 2 if 0 and/or 1 are included in the Farey sequence), or the sum of the totients of the numbers  2 <= k <= n. f n0 is not easily directly computable, but then g n = n*(n-1)/2\ is very easy to compute, and hence the inversion gives an efficient method of computing f n.For  valued functions, unboxed arrays can be used for greater efficiency. That bears the risk of overflow, however, so be sure to use it only when it's safe. The value f n is then computed by generalInversion g n#. Note that when many values of ft are needed, there are far more efficient methods, this method is only appropriate to compute isolated values of f.,(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>SafeIRarithmoiFollowing property holds: LapproxPrimeCount n >= primeCount n || n >= approxPrimeCountOverestimateLimitSarithmoiS nD gives an approximation of the number of primes not exceeding n'. The approximation is fairly good for n large enough.TarithmoiFollowing property holds: GnthPrimeApprox n <= nthPrime n || n >= nthPrimeApproxUnderestimateLimitUarithmoiU nV gives an approximation to the n-th prime. The approximation is fairly good for n large enough.RSTU(c) 2019 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NoneVV-(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>Safe7HVfWarithmoiWrapper for prime elements of a'. It is supposed to be constructed by  . /  /. and eliminated by X.One can leverage ? instance to generate lists of primes. Here are some examples.(Generate primes from the given interval: [nextPrime 101 .. precPrime 130]=[Prime 101,Prime 103,Prime 107,Prime 109,Prime 113,Prime 127]$Generate an infinite list of primes:[nextPrime 101 ..]?[Prime 101,Prime 103,Prime 107,Prime 109,Prime 113,Prime 127...9Generate primes from the given interval of form p = 6k+5:/[nextPrime 101, nextPrime 107 .. precPrime 150]=[Prime 101,Prime 107,Prime 113,Prime 131,Prime 137,Prime 149]Get next prime:succ (nextPrime 101) Prime 103Get previous prime:prec (nextPrime 101)Prime 97+Count primes less than a given number (cf. 0):fromEnum (precPrime 100)25Get 25-th prime number (cf. 1):toEnum 25 :: Prime IntPrime 97XarithmoiUnwrap prime element.WX (c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>SafeYarithmoi(Infinite zero-based table of factorials.take 5 factorial [1,1,2,6,24] The time-and-space behaviour of Y is similar to described in -Math.NumberTheory.Recurrences.Bilinear#memory.ZarithmoiZ k calculates the k-th Fibonacci number in O( log (abs k)1) steps. The index may be negative. This is efficient for calculating single Fibonacci numbers (with large index), but for computing many Fibonacci numbers in close proximity, it is better to use the simple addition formula starting from an appropriate pair of successive Fibonacci numbers.[arithmoi[ k returns the pair (F(k), F(k+1)) of the k-th Fibonacci number and its successor, thus it can be used to calculate the Fibonacci numbers from some index on without needing to compute the previous. The pair is efficiently calculated in O( log (abs k)#) steps. The index may be negative.\arithmoi\ k computes the k%-th Lucas number. Very similar to Z.]arithmoi] k computes the pair (L(k), L(k+1)) of the k7-th Lucas number and its successor. Very similar to [.^arithmoi^ p q k calculates the quadruple (U(k), U(k+1), V(k), V(k+1)) where U(i)- is the Lucas sequence of the first kind and V(i)= the Lucas sequence of the second kind for the parameters p and q, where  p^2-4q /= 04. Both sequences satisfy the recurrence relation A(j+2) = p*A(j+1) - q*A(j), the starting values are U(0) = 0, U(1) = 1 and V(0) = 2, V(1) = p[. The Fibonacci numbers form the Lucas sequence of the first kind for the parameters  p = 1, q = -1 and the Lucas numbers form the Lucas sequence of the second kind for these parameters. Here, the index must be non-negative, since the terms of the sequence for negative indices are in general not integers.YZ[\]^YZ[\]^ (c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>SafeXdK _arithmoiUInfinite zero-based table of binomial coefficients (also known as Pascal triangle): (binomial !! n !! k == n! / k! / (n - k)!.take 5 (map (take 5) binomial))[[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1]] Complexity: binomial !! n !! kI is O(n) bits long, its computation takes O(k n) time and forces thunks binomial !! n !! i for  0 <= i <= k'. Use the symmetry of Pascal triangle .binomial !! n !! k == binomial !! n !! (n - k) to speed up computations.One could also consider 23 to compute stand-alone values.`arithmoiInfinite zero-based table of  @https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind"Stirling numbers of the first kind.take 5 (map (take 5) stirling1)*[[1],[0,1],[0,1,1],[0,2,3,1],[0,6,11,6,1]] Complexity: stirling1 !! n !! kU is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks stirling1 !! i !! j for  0 <= i <= n and max(0, k - n + i) <= j <= k.One could also consider 24 to compute stand-alone values.aarithmoiInfinite zero-based table of  Ahttps://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind#Stirling numbers of the second kind.take 5 (map (take 5) stirling2))[[1],[0,1],[0,1,1],[0,1,3,1],[0,1,7,6,1]] Complexity: stirling2 !! n !! kU is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks stirling2 !! i !! j for  0 <= i <= n and max(0, k - n + i) <= j <= k.One could also consider 25 to compute stand-alone values.barithmoiInfinite one-based table of  (https://en.wikipedia.org/wiki/Lah_number Lah numbers.  lah !! n !! k equals to lah(n + 1, k + 1).take 5 (map (take 5) lah)3[[1],[2,1],[6,6,1],[24,36,12,1],[120,240,120,20,1]] Complexity:  lah !! n !! kS is O(n ln n) bits long, its computation takes O(k n ln n) time and forces thunks  lah !! n !! i for  0 <= i <= k.carithmoiInfinite zero-based table of  -https://en.wikipedia.org/wiki/Eulerian_number"Eulerian numbers of the first kind.take 5 (map (take 5) eulerian1)"[[],[1],[1,1],[1,4,1],[1,11,11,1]] Complexity: eulerian1 !! n !! kU is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks eulerian1 !! i !! j for  0 <= i <= n and max(0, k - n + i) <= j <= k.darithmoiInfinite zero-based table of  Qhttps://en.wikipedia.org/wiki/Eulerian_number#Eulerian_numbers_of_the_second_kind#Eulerian numbers of the second kind.take 5 (map (take 5) eulerian2)#[[],[1],[1,2],[1,8,6],[1,22,58,24]] Complexity: eulerian2 !! n !! kU is O(n ln n) bits long, its computation takes O(k n^2 ln n) time and forces thunks eulerian2 !! i !! j for  0 <= i <= n and max(0, k - n + i) <= j <= k.earithmoi Infinite zero-based sequence of  .https://en.wikipedia.org/wiki/Bernoulli_numberBernoulli numbers, computed via  bhttps://en.wikipedia.org/wiki/Bernoulli_number#Connection_with_Stirling_numbers_of_the_second_kind connection with a.take 5 bernoulli&[1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30] Complexity: bernoulli !! nS is O(n ln n) bits long, its computation takes O(n^3 ln n) time and forces thunks stirling2 !! i !! j for  0 <= i <= n and  0 <= j <= i.One could also consider 26 to compute stand-alone values.farithmoi 3https://en.wikipedia.org/wiki/Faulhaber%27s_formulaFaulhaber's formula.sum (map (^ 10) [0..100])9599241424342419242508sum $ zipWith (*) (faulhaberPoly 10) (iterate (* 100) 1)959924142434241924250 % 1arithmoiInfinite zero-based list of  *https://en.wikipedia.org/wiki/Euler_number Euler numbers'. The algorithm used was derived from  9http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf2Algorithms for Bernoulli numbers and Euler numbersI by Kwang-Wu Chen, second formula of the Corollary in page 7. Sequence  https://oeis.org/A122045A122045 in OEIS.take 10 euler' :: [Rational]G[1 % 1,0 % 1,(-1) % 1,0 % 1,5 % 1,0 % 1,(-61) % 1,0 % 1,1385 % 1,0 % 1]garithmoiThe same sequence as euler', but with type [a] instead of  [Ratio a] as the denominators in euler' are always 1.take 10 euler :: [Integer][1,0,-1,0,5,0,-61,0,1385,0]harithmoi Infinite zero-based list of the n)-th order Euler polynomials evaluated at 1'. The algorithm used was derived from  9http://www.emis.ams.org/journals/JIS/VOL4/CHEN/AlgBE2.pdf2Algorithms for Bernoulli numbers and Euler numbersh by Kwang-Wu Chen, third formula of the Corollary in page 7. Element-by-element division of sequences  https://oeis.org/A198631A1986631 and  https://oeis.org/A006519A006519 in OEIS."take 10 eulerPolyAt1 :: [Rational]E[1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2]arithmoiHelper for common code in =bernoulli, euler, eulerPolyAt1. All three sequences rely on  stirling2 and have the same general structure of zipping four lists together with multiplication, with one of those lists being the sublists in  stirling24, and two of them being the factorial sequence and  cycle [1, -1]#. The remaining list is passed to  helperForB_E_EP@ as an argument.Note: This function has a ([Ratio a] -> [Ratio a]) argument because bernoulli !! n will use, for all nonnegative n, every element in stirling2 !! n, while euler, eulerPolyAt1 only use tail $ stirling2 !! n(. As such, this argument serves to pass id in the former case, and tail in the latter. _`abcdefgh _`abcdeghf7"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>SafeSXqarithmoi:Infinite list of generalized pentagonal numbers. Example: take 10 pents[0,1,2,5,7,12,15,22,26,35]arithmoiWhen calculating the n-th partition number p(n) using the sum 8p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-11) + ...U, the signs of each term alternate every two terms, starting with a positive sign. pentagonalSignsN takes a list of numbers and produces such an alternated sequence. Examples:pentagonalSigns [1..5] [1,2,-3,-4,5]pentagonalSigns [1..6][1,2,-3,-4,5,6]iarithmoiInfinite zero-based table of  https://oeis.org/A000041partition numbers.take 10 partition[1,1,2,3,5,7,11,15,22,30]:set -XDataKinds%import Math.NumberTheory.Moduli.Classpartition !! 1000 :: Mod 1000(991 `modulo` 1000)i "(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>Safes^YZ[\]^_`abcdefghii8(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>Noneu 9(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>Nonev:(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NoneFarithmoiRemove factors of 2 and count them. If  n = 2^k*m with m odd, the result is (k, m) . Precondition: argument not 0 (not checked).arithmoiSpecialised version for <. Precondition: argument strictly positive (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiCount trailing zeros in a E. Precondition: argument nonzero (not checked, Integer invariant).arithmoiRemove factors of 2. If  n = 2^k*m with m odd, the result is m . Precondition: argument not 0 (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiSpecialised version for 2. Precondition: argument nonzero (not checked).arithmoiLShift argument right until the result is odd. Precondition: argument not 0, not checked.arithmoiLike ., but count the number of places to shift too.arithmoiNumber of 1-bits in a .arithmoiNumber of 1-bits in a .arithmoiNumber of 1-bits in an .arithmoi3It is difficult to convince GHC to unbox output of  and  'splitOff.go'K, so we fallback to a specialized unboxed version to minimize allocations.arithmoieMerges two ordered lists into an ordered list. Checks for neither its precondition or postcondition.arithmoi Work around -https://ghc.haskell.org/trac/ghc/ticket/14085 4(c) 2011 Daniel Fischer, 2017-2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None_jarithmoi%Represents three possible values of  +https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol.narithmoi +https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symboll of two arguments. The lower argument ("denominator") must be odd and positive, this condition is checked.2If arguments have a common factor, the result is l, otherwise it is k or m.jacobi 1001 9911)Zero -- arguments have a common factor 11jacobi 1001 9907MinusOnejmlknjmlkn;/(c) 2011 Daniel Fischer, 2017 Andrew LelechenkoMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NoneFtarithmoi isPrime n tests whether nb is a prime (negative or positive). It is a combination of trial division and Baillie-PSW test.If  isPrime n returns False then nE is definitely composite. There is a theoretical possibility that  isPrime n is True, but in fact nJ is not prime. However, no such numbers are known and none exist below 2^64N. If you have found one, please report it, because it is a major discovery.uarithmoiMiller-Rabin probabilistic primality test. It consists of the trial division test and several rounds of the strong Fermat test with different bases. The choice of trial divisors and bases are implementation details and may change in future silently.yFirst argument stands for the number of rounds of strong Fermat test. If it is 0, only trial division test is performed.If millerRabinV k n returns False then n% is definitely composite. Otherwise n' may appear composite with probability 1/4^k.varithmoiv n b tests whether non-negative n/ is a strong Fermat probable prime for base b.'Apart from primes, also some composite numbers have the tested property, but those are rare. Very rare are composite numbers having the property for many bases, so testing a large prime candidate with several bases can identify composite numbers with high probability. An odd number n > 3 is prime if and only if v n b holds for all b with 2 <= b <= (n-1)/2 , but of course checking all those bases would be less efficient than trial division, so one normally checks only a relatively small number of bases, depending on the desired degree of certainty. The probability that a randomly chosen base doesn't identify a composite number n is less than 1/4J, so five to ten tests give a reasonable level of certainty in general.Please consult  https://miller-rabin.appspot.com9Deterministic variants of the Miller-Rabin primality test! for the best choice of bases.warithmoiw n b tests whether n, is a Fermat probable prime for the base b, that is, whether b^(n-1)  n == 1. This is a weaker but simpler condition. However, more is lost in strength than is gained in simplicity, so for primality testing, the strong check should be used. The remarks about the choice of bases to test from v( apply with the modification that if a and b are Fermat bases for n, then a*b always is a Fermat base for n too. A Charmichael number is a composite number n3 which is a Fermat probable prime for all bases b coprime to n. By the above, only primes p <= n/2 not dividing nI need to be tested to identify Carmichael numbers (however, testing all those primes would be less efficient than determining Carmichaelness from the prime factorisation; but testing an appropriate number of prime bases is reasonable to find out whether it's worth the effort to undertake the prime factorisation).xarithmoiPrimality test after Baillie, Pomerance, Selfridge and Wagstaff. The Baillie-PSW test consists of a strong Fermat probable primality test followed by a (strong) Lucas primality test. This implementation assumes that the number n to test is odd and larger than 3. Even and small numbers have to be handled before. Also, before applying this test, trial division by small primes should be performed to identify many composites cheaply (although the Baillie-PSW test is rather fast, about the same speed as a strong Fermat test for four or five bases usually, it is, for large numbers, much more costly than trial division by small primes, the primes less than 1000X, say, so eliminating numbers with small prime factors beforehand is more efficient).The Baillie-PSW test is very reliable, so far no composite numbers passing it are known, and it is known (Gilchrist 2010) that no Baillie-PSW pseudoprimes exist below 2^64. However, a heuristic argument by Pomerance indicates that there are likely infinitely many Baillie-PSW pseudoprimes. On the other hand, according to  :http://mathworld.wolfram.com/Baillie-PSWPrimalityTest.html there is reason to believe that there are none with less than several thousand digits, so that for most use cases the test can be considered definitive.arithmoipThe Lucas-Selfridge test, including square-check, but without the Fermat test. For package-internal use only.tuvwx/(c) 2011 Daniel Fischer, 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None PSUVXk&yarithmoiy (n1, m1) (n2, m2) returns n such that n `mod` m1 == n1 and n `mod` m2 == n2 . Moduli m1 and m2 must be coprime, otherwise  is returned.&This function is slightly faster than z, but more restricted.chineseCoprime (1, 2) (2, 3)Just 5chineseCoprime (3, 4) (5, 6)!Nothing -- moduli must be coprimezarithmoiz (n1, m1) (n2, m2) returns n such that n `mod` m1 == n1 and n `mod` m2 == n2, if exists. Moduli m1 and m2$ are allowed to have common factors.chinese (1, 2) (2, 3)Just 5chinese (3, 4) (5, 6)Just 11chinese (3, 4) (2, 6)Nothing{arithmoiSame as y, but operates on residues.:set -XDataKinds%import Math.NumberTheory.Moduli.Class5(1 `modulo` 2) `chineseCoprimeSomeMod` (2 `modulo` 3)Just (5 `modulo` 6)5(3 `modulo` 4) `chineseCoprimeSomeMod` (5 `modulo` 6)Nothing|arithmoiSame as z, but operates on residues.:set -XDataKinds%import Math.NumberTheory.Moduli.Class.(1 `modulo` 2) `chineseSomeMod` (2 `modulo` 3)Just (5 `modulo` 6).(3 `modulo` 4) `chineseSomeMod` (5 `modulo` 6)Just (11 `modulo` 12).(3 `modulo` 4) `chineseSomeMod` (2 `modulo` 6)Nothing}arithmoi Given a list [(r_1,m_1), ..., (r_n,m_n)] of (residue,modulus) pairs, chineseRemainder; calculates the solution to the simultaneous congruences  r "a r_k (mod m_k) Lif all moduli are positive and pairwise coprime. Otherwise the result is Nothing, regardless of whether a solution exists.~arithmoi%chineseRemainder2 (r_1,m_1) (r_2,m_2) calculates the solution of  r "a r_k (mod m_k)if m_1 and m_2 are coprime.yz{|}~zy|{}~<(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None &'.FHX# arithmoiPoint on unknown curve.arithmoinWe use the Montgomery form of elliptic curve: b Y = X + a X + X (mod n). See Eq. (10.3.1.1) at p. 260 of  dhttp://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/S0025-5718-1987-0866113-7.pdf@Speeding the Pollard and Elliptic Curve Methods of Factorization by P. L. Montgomery.DSwitching to projective space by substitutions Y = y / z, X = x / z, we get b y z = x + a x z + x z (mod n). The point on projective elliptic curve is characterized by three coordinates, but it appears that only x- and z-components matter for computations. By the same reason there is no need to store coefficient b.That said, the chosen curve is represented by a24 = (a + 2) / 4 and modulo n at type level, making points on different curves incompatible.arithmoiExtract x-coordinate.arithmoiExtract z-coordinate.arithmoiBExtract (a + 2) / 4, where a is a coefficient in curve's equation.arithmoiExtract modulo of the curve.arithmoi s n- creates a point on an elliptic curve modulo n, uniquely determined by seed s. Some choices of s and nE produce ill-parametrized curves, which is reflected by return value .KWe choose a curve by Suyama's parametrization. See Eq. (3)-(4) at p. 4 of  9http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdfNImplementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware by K. Gaj, S. Kwon et al.arithmoiIf p0 + p1 = p2, then  p0 p1 p2 equals to p1 + p2-. It is also required that z-coordinates of p0, p1 and p2A are coprime with modulo of elliptic curve; and x-coordinate of p0G is non-zero. If preconditions do not hold, return value is undefined.*Remarkably such addition does not require  a24 constraint.,Computations follow Algorithm 3 at p. 4 of  9http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdfNImplementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware by K. Gaj, S. Kwon et al.arithmoiMultiply by 2.,Computations follow Algorithm 3 at p. 4 of  9http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdfNImplementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware by K. Gaj, S. Kwon et al.arithmoi1Multiply by given number, using binary algorithm.arithmoiFor debugging.arithmoiIn projective space Cs are equal, if they are both at infinity or if respective ratios  /  are equal. =(c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%Xk(arithmoiLet as = sum_i k_ia_i^s and bs = sum_i l_ib_i^s be Dirichlet series, and all a_i and b_i are divisors of n. Return Dirichlet series cs, which contains all terms as * bs = sum_i m_i/c_i^s such that c_i divides n.  >(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>Safe*    ?(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>X_I arithmoi!Compact store of primality flags.arithmoiSieve primes up to (and including) a bound (or 7, if bound is smaller). For small enough bounds, this is more efficient than using the segmented sieve.Since arrays are <-indexed, overflow occurs when the sieve size comes near  :: -, that corresponds to an upper bound near 15/8*. On 32Y-bit systems, that is often within memory limits, so don't give bounds larger than 8*10^9 there.arithmoi4Generate a list of primes for consumption from a .arithmoiAReturns true if integer is beyond representation range of type a.arithmoiQExtracts the longest strictly increasing prefix of the list (possibly infinite).arithmoiAscending list of primes.take 10 primesW[Prime 2,Prime 3,Prime 5,Prime 7,Prime 11,Prime 13,Prime 17,Prime 19,Prime 23,Prime 29] is a polymorphic list, so the results of computations are not retained in memory. Make it monomorphic to take advantages of memoization. Compare:set +sprimes !! 1000000 :: Prime IntPrime 15485867 (5.32 secs, 6,945,267,496 bytes)primes !! 1000000 :: Prime IntPrime 15485867 (5.19 secs, 6,945,267,496 bytes)against#let primes' = primes :: [Prime Int]primes' !! 1000000 :: Prime IntPrime 15485867 (5.29 secs, 6,945,269,856 bytes)primes' !! 1000000 :: Prime IntPrime 15485867(0.02 secs, 336,232 bytes)arithmoi(List of primes in the form of a list of s, more compact than , thus it may be better to use  >>= @ than keeping the list of primes alive during the entire run.arithmoiSieve up to bound in one go.arithmoi n* creates the list of primes not less than n.arithmoi n creates the list of s starting roughly at nW. Due to the organisation of the sieve, the list may contain a few primes less than n%. This form uses less memory than []B, hence it may be preferable to use this if it is to be reused. Deprecated(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NoneK@(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NoneYiarithmoi Factorise an  using a given list of numbers considered prime. If the list is not a list of primes containing all relevant primes, the result could be surprising.arithmoi bound n produces a factorisation of n using the primes <= bound. If n has prime divisors > bound@, the last entry in the list is the product of all these. If  n <= bound^24, this is a full factorisation, but very slow if n has large prime divisors.arithmoiCheck whether a number is coprime to all of the numbers in the list (assuming that list contains only numbers > 1 and is ascending).arithmoi bound n tests whether n is coprime to all primes <= bound. If  n <= bound^2., this is a full prime test, but very slow if n has no small prime divisors.A(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>XiarithmoiMaximal allowed argument of . Currently 8e18.arithmoi n == (n)2 is the number of (positive) primes not exceeding n.OFor efficiency, the calculations are done on 64-bit signed integers, therefore n must not exceed . Requires O(n^0.5)' space, the time complexity is roughly O(n^0.7). For small bounds,  n( simply counts the primes not exceeding n, for bounds from 30000T on, Meissel's algorithm is used in the improved form due to D.H. Lehmer, cf.  \http://en.wikipedia.org/wiki/Prime_counting_function#Algorithms_for_evaluating_.CF.80.28x.29.arithmoiMaximal allowed argument of . Currently 1.5e17.arithmoi n calculates the n%-th prime. Numbering of primes is 1 -based, so  1 == 2. Requires O((n*log n)^0.5)' space, the time complexity is roughly O((n*log n)^0.7A. The argument must be strictly positive, and must not exceed .(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NonejRSTUSRUT(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NoneF{ arithmoipowMod b e m computes (b^e) `mod` m in effective way. An exponent e must be non-negative, a modulo m/ must be positive. Otherwise the behaviour of powMod is undefined. powMod 2 3 53"powMod 3 12345678901234567890 10011 See also B and C.For finite numeric types (, , etc.) modulo m should be such that m^2 does not overflow, otherwise the behaviour is undefined. If you need both to fit into machine word and to handle large moduli, take a look at  and . powMod 3 101 (2^60-1 :: Integer)1018105167100379328 -- correctpowMod 3 101 (2^60-1 :: Int)01115647832265427613 -- incorrect due to overflowpowModInt 3 101 (2^60-1 :: Int)1018105167100379328 -- correctarithmoiSpecialised version of (, able to handle large moduli correctly.powModWord 3 101 (2^60-1)1018105167100379328arithmoiSpecialised version of (, able to handle large moduli correctly. > powModInt 3 101 (2^60-1)1018105167100379328(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>NoneFFarithmoiCalculate an integer root,  k n computes the (floor of) the k-th root of n, where k must be positive. r =  k n means r^k <= n < (r+1)^k4 if that is possible at all. It is impossible if k is even and n < 0 , since then r^k >= 0 for all r, then, and if k <= 0,  raises an error. For k < 5, a specialised version is called which should be more efficient than the general algorithm. However, it is not guaranteed that the rewrite rules for those fire, so if kO is known in advance, it is safer to directly call the specialised versions.arithmoi k n returns  if n is not a k-th power,  r if n == r^k. If k is divisible by 4, 3 or 2i, a residue test is performed to avoid the expensive calculation if it can thus be determined that n is not a k -th power.arithmoi k n checks whether n is a k -th power.arithmoi n checks whether n == r^k for some k > 1.arithmoi n produces the pair (b,k) with the largest exponent k such that n == b^k , except for  n <= 1V, in which case arbitrarily large exponents exist, and by an arbitrary decision (n,3) is returned.^First, by trial division with small primes, the range of possible exponents is reduced (if p^e exactly divides n, then k must be a divisor of e!, if several small primes divide n, kX must divide the greatest common divisor of their exponents, which mostly will be 1:, generally small; if none of the small primes divides n, the range of possible exponents is reduced since the base is necessarily large), if that has not yet determined the result, the remaining factor is examined by trying the divisors of the gcdg of the prime exponents if some have been found, otherwise by trying prime exponents recursively.arithmoi bd n produces the pair (b,k) with the largest exponent k such that n == b^k, where bd > 1 (it is expected that bd is much larger, at least 1000 or so), n > bd^2 and n has no prime factors p <= bd*, skipping the trial division phase of 0 when that is a priori known to be superfluous. It is only present to avoid duplication of work in factorisation and primality testing, it is not expected to be generally useful. The assumptions are not checked, if they are not satisfied, wrong results and wasted work may be the consequence.D(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None.HX_ arithmoi n% produces the prime factorisation of n.  0) is an error and the factorisation of 1 is empty. Uses a  < produced in an arbitrary manner from the bit-pattern of n.Warning:n there are no guarantees of any particular order of prime factors, do not expect them to be ascending. E. g.,factorise 10251562501[(101701,1),(100801,1)]arithmoiLike $, but without input checking, hence n > 1 is required.arithmoi is like , except that it doesn't use a pseudo random generator but steps through the curves in order. This strategy turns out to be surprisingly fast, on average it doesn't seem to be slower than the   based variant.arithmoi first strips off all small prime factors and then, if the factorisation is not complete, proceeds to curve factorisation. For negative numbers, a factor of -1# is included, the factorisation of 1 is empty. Since 0@ has no prime factorisation, a zero argument causes an error.arithmoiLike $, but without input checking, so n must be larger than 1.arithmoiA wrapper around = providing a few default arguments. The primality test is x, the prng function - naturally - !R. This function also requires small prime factors to have been stripped before.arithmoi is the driver for the factorisation. Its performance (and success) can be influenced by passing appropriate arguments. If you know that n has no prime divisors below b, any divisor found less than b*b must be prime, thus giving  Just (b*b)f as the first argument allows skipping the comparatively expensive primality test for those. If n is such that all prime divisors must have a specific easy to test for structure, a custom primality test can improve the performance (normally, it will make very little difference, since n has not many divisors, and many curves have to be tried to find one). More influence has the pseudo random generator (a function prng with 6 <= fst (prng k s) <= k-2 and an initial state for the PRNG) used to generate the curves to try. A lucky choice here can make a huge difference. So, if the default takes too long, try another one; or you can improve your chances for a quick result by running several instances in parallel. n0 requires that small (< 65536) prime factors of n^ have been stripped before. Otherwise it is likely to cycle forever. When in doubt, use . is unlikely to succeed if n/ has more than one (really) large prime factor.arithmoi n b1 b2 s tries to find a factor of n5 using the curve and point determined by the seed s ( 6 <= s < n-1H), multiplying the point by the least common multiple of all numbers <= b1 and all primes between b1 and b2. The idea is that there's a good chance that the order of the point in the curve over one prime factor divides the multiplier, but the order over another factor doesn't, if b1 and b2 are appropriately chosen. If they are too small, none of the orders will probably divide the multiplier, if they are too large, all probably will, so they should be chosen to fit the expected size of the smallest factor.It is assumed that n has no small prime factors.,The result is maybe a nontrivial divisor of n."arithmoi7The implementation follows the algorithm at p. 6-7 of  9http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdfNImplementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware by K. Gaj, S. Kwon et al.#arithmoi^Same as map (id *** flip multiply p) [from, thn .. to], but calculated in more efficient way.arithmoi n finds all prime divisors of n > 1 up to 2^16 by trial division and returns the list of these together with their multiplicities, and a possible remaining factor which may be composite.$arithmoiFor a given estimated decimal length of the smallest prime factor ("tier") return parameters B1, B2 and the number of curves to try before next "tier". Roughly based on `http://www.mersennewiki.org/index.php/Elliptic_Curve_Method#Choosing_the_best_parameters_for_ECMarithmoi"Lower bound for composite divisorsarithmoi Standard PRNGarithmoi3Estimated number of digits of smallest prime factorarithmoiThe number to factorisearithmoi#List of prime factors and exponentsarithmoi"Lower bound for composite divisorsarithmoiA primality testarithmoiA PRNGarithmoiInitial PRNG statearithmoi7Estimated number of digits of the smallest prime factorarithmoiThe number to factorisearithmoi#List of prime factors and exponents % Deprecated(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>Noneh (c) 2016-2018 Andrew.LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None=?_ arithmoi)A class for unique factorisation domains.arithmoiFactorise a number into a product of prime powers. Factorisation of 0 is an undefined behaviour. Otherwise following invariants hold: eabs n == abs (product (map (\(p, k) -> unPrime p ^ k) (factorise n))) all ((> 0) . snd) (factorise n)factorise (1 :: Integer)[]factorise (-1 :: Integer)[]factorise (6 :: Integer)[(Prime 2,1),(Prime 3,1)]factorise (-108 :: Integer)[(Prime 2,2),(Prime 3,3)]$This function is a replacement for E6. If you were looking for the latter, please import &Math.NumberTheory.Primes.Factorisation instead of this module.Warning:n there are no guarantees of any particular order of prime factors, do not expect them to be ascending. E. g.,factorise 10251562501#[(Prime 101701,1),(Prime 100801,1)]arithmoiNCheck whether an argument is prime. If it is then return an associated prime.isPrime (3 :: Integer)Just (Prime 3)isPrime (4 :: Integer)NothingisPrime (-5 :: Integer)Just (Prime 5)$This function is a replacement for F6. If you were looking for the latter, please import  Math.NumberTheory.Primes.Testing instead of this module.arithmoi-Smallest prime, greater or equal to argument. JnextPrime (-100) == 2 nextPrime 1000 == 1009 nextPrime 1009 == 1009arithmoiGLargest prime, less or equal to argument. Undefined, when argument < 2. 'precPrime 100 == 97 precPrime 97 == 97 WX WX(c) 2018 Frederick SchneiderMIT7Frederick Schneider <frederick.schneider2011@gmail.com>None>Xk* arithmoiRAn abstract representation of a smooth basis. It consists of a set of numbers "e2.arithmoiBuild a  from a set of numbers "e2. import qualified Data.Set as SetfromSet (Set.fromList [2, 3])*Just (SmoothBasis {unSmoothBasis = [2,3]})fromSet (Set.fromList [2, 4])*Just (SmoothBasis {unSmoothBasis = [2,4]})/fromSet (Set.fromList [1, 3]) -- should be >= 2NothingarithmoiBuild a  from a list of numbers "e2.fromList [2, 3]*Just (SmoothBasis {unSmoothBasis = [2,3]})fromList [2, 2](Just (SmoothBasis {unSmoothBasis = [2]})fromList [2, 4]*Just (SmoothBasis {unSmoothBasis = [2,4]})!fromList [1, 3] -- should be >= 2NothingarithmoiBuild a ) from a list of primes below given bound.fromSmoothUpperBound 10.Just (SmoothBasis {unSmoothBasis = [2,3,5,7]})fromSmoothUpperBound 1NothingarithmoiHelper used by  smoothOver (Integral constraint) and  smoothOver' ( Euclidean5 constraint) Since the typeclass constraint is just Num, it receives a norm comparison function for the generated smooth numbers. This function relies on the fact that for any element of a smooth basis p and any a it is true that norm (a * p) > norm a!. This condition is not checked.arithmoi(Generate an infinite ascending list of  +https://en.wikipedia.org/wiki/Smooth_numbersmooth numbers over a given smooth basis.import Data.Maybe1take 10 (smoothOver (fromJust (fromList [2, 5])))[1,2,4,5,8,10,16,20,25,32]arithmoiGenerate an ascending list of  +https://en.wikipedia.org/wiki/Smooth_numbersmooth numbers- over a given smooth basis in a given range.JIt may appear inefficient for short, but distant ranges; consider using  in such cases.import Data.Maybe6smoothOverInRange (fromJust (fromList [2, 5])) 100 200[100,125,128,160,200]arithmoiGenerate an ascending list of  +https://en.wikipedia.org/wiki/Smooth_numbersmooth numbers- over a given smooth basis in a given range.HIt is inefficient for large or starting near 0 ranges; consider using  in such cases.MSuffix BF stands for the brute force algorithm, involving a lot of divisions.import Data.Maybe8smoothOverInRangeBF (fromJust (fromList [2, 5])) 100 200[100,125,128,160,200]arithmoiisSmooth4 checks if a given number is smooth under a certain  SmoothBasis. Does not check if the  SmoothBasis is valid. G(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>Nonep&arithmoi5An argument for primality of a number (which must be > 1). s translate directly to s, correct arguments can be transformed into proofs. This type allows the manipulation of proofs while maintaining their correctness. The only way to access components of a * except the prime is through this type.arithmoi#A suggested Pocklington certificatearithmoi2Primality should be provable by trial division to alimitarithmoiaprime6 is said to be obviously prime, that holds for primes < 30arithmoiPrimality assumedarithmoidA proof of primality of a positive number. The type is abstract to ensure the validity of proofs.arithmoi%The number whose primality is proved.arithmoi9An argument for compositeness of a number (which must be > 1). s translate directly to s, correct arguments can be transformed into proofs. This type allows the manipulation of proofs while maintaining their correctness. The only way to access components of a . except the composite is through this type.arithmoicompo == firstDiv*secondDiv, where all are > 1arithmoicompo" fails the strong Fermat test for  fermatBasearithmoicompo fails the Lucas-Selfridge testarithmoiNo particular reason givenarithmoihA proof of compositeness of a positive number. The type is abstract to ensure the validity of proofs.arithmoi)The number whose compositeness is proved.arithmoi<A certificate of either compositeness or primality of an . Only numbers > 1W can be certified, trying to create a certificate for other numbers raises an error.arithmoi Eliminate .arithmoi@ transforms a proof of primality into an argument for primality.arithmoi checks the given argument and constructs a proof from it, if it is valid. For the explicit arguments, this is simple and resonably fast, for an , the verification uses  and hence may take a long time.&arithmoiKnot exported, this is the one place where invalid proofs can be constructedarithmoiK transforms a proof of compositeness into an argument for compositeness.arithmoi checks the given argument and constructs a proof from it, if it is valid. For the explicit arguments, this is simple and resonably fast, for a , the verification uses  and hence may take a long time.'arithmoiTnot exported, here is where invalid proofs can be constructed, they must not leakarithmoiCheck the validity of a z. Since it should be impossible to construct invalid certificates by the public interface, this should never return (.arithmoiCheck the validity of a q. Since it should be impossible to create invalid proofs by the public interface, this should never return (.arithmoiCheck the validity of a q. Since it should be impossible to create invalid proofs by the public interface, this should never return (.)arithmoi)\ records a trivially known prime. If the argument is not one of them, an error is raised.*arithmoi*\ finds out if its argument is a trivially known prime or not and returns the appropriate.+arithmoi+; checks whether its argument is a trivially known prime.,arithmoiList of trivially known primes.-arithmoiCertify a small number. This is not exposed and should only be used where correct. It is always checked after use, though, so it shouldn't be able to lie.arithmoi n constructs, for n > 15, a proof of either primality or compositeness of nR. This may take a very long time if the number has no small(ish) prime divisors.arithmoiCertify a number known to be not too small, having no small prime divisors and having passed the Baillie PSW test. So we assume it's prime, erroring if not. Since it's presumably a large number, we don't bother with trial division and construct a Pocklington certificate./arithmoimFind a decomposition of p-1 for the pocklington certificate. Usually bloody slow if p-1 has two (or more) large prime divisors.0arithmoiFind a factor of a known composite with approximately digits digits, starting with curve s. Actually, this may loop infinitely, but the loop should not be entered before the heat death of the universe.1arithmoi2Find a factor or say with which curve to continue.2arithmoi@Message in the unlikely case a Baillie PSW pseudoprime is found.3arithmoiFound a factor4arithmoiFermat failure356789:;<=>?@A)-.H(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>Nonexarithmoi n tests primality of n, first trial division by small primes is performed, then a Baillie PSW test and finally a prime certificate is constructed and verified, provided no step before found n@ to be composite. Constructing prime certificates can take a very" long time, so use this with care.(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>Noneyytuvwxtxuvw Deprecated(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None{@## Deprecated(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None> arithmoi n( produces the prime factorisation of nF, proving the primality of the factors, but doesn't report the proofs.arithmoi n produces a .arithmoi n) constructs a the prime factorisation of n (which must be positive) together with proofs of primality of the factors, using trial division up to 2^16 and elliptic curve factorisation for the remaining factors if necessary.,Construction of primality proofs can take a veryM long time, so this will usually be slow (but should be faster than using ; and proving the primality of the factors from scratch).Barithmoi4verify that we indeed have a correct primality proofCarithmoiproduce a proof of primality for primes Only called for (not too small) numbers known to have no small prime factors, so we can directly use BPSW without trial division.Darithmoiproduce a certified factorisation Assumes all small prime factors have been stripped before. Since it is not exported, that is known to hold. This is a near duplicate of &, I should some time clean this up.EarithmoiVmerge two lists of factors, so that the result is strictly increasing (wrt the primes)Farithmoi[merge a list of lists of factors so that the result is strictly increasing (wrt the primes)Garithmoi2message for an invalid proof, should never be usedDarithmoi"Lower bound for composite divisorsarithmoiA primality testarithmoiA PRNGarithmoiInitial PRNG statearithmoi7Estimated number of digits of the smallest prime factorarithmoiThe number to factorisearithmoi5List of prime factors, exponents and primality proofs(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NoneHVİarithmoiA container for a number and its pairwise coprime (but not neccessarily prime) factorisation. It is designed to preserve information about factors under multiplication. One can use this representation to speed up prime factorisation and computation of arithmetic functions.For instance, let p and q be big primes:2let p = 1000000000000000000000000000057 :: Integer2let q = 2000000000000000000000000000071 :: IntegerIt would be difficult to compute the totient function of their product as is, because once we multiplied them the information of factors is lost and I (p * qJ) would take ages. Things become different if we simply change types of p and q to prefactored ones:>let p = 1000000000000000000000000000057 :: Prefactored Integer>let q = 2000000000000000000000000000071 :: Prefactored IntegerNow the I% function can be computed instantly:,import Math.NumberTheory.ArithmeticFunctionsprefValue $ totient (p^2 * q^3)8000000000000000000000000001752000000000000000000000000151322000000000000000000000006445392000000000000000000000135513014000000000000000000001126361040)prefValue $ totient $ totient (p^2 * q^3)2133305798262843681544648472180210822742702690942899511234131900112583590230336435053688694839034890779375223070157301188739881477320529552945446912000Let us look under the hood:,import Math.NumberTheory.ArithmeticFunctions!prefFactors $ totient (p^2 * q^3)Coprimes {unCoprimes = [(1000000000000000000000000000057,1),(41666666666666666666666666669,1),(2000000000000000000000000000071,2),(111111111111111111111111111115,1),(2,4),(3,3)]}+prefFactors $ totient $ totient (p^2 * q^3)Coprimes {unCoprimes = [(39521,1),(6046667,1),(22222222222222222222222222223,1),(2000000000000000000000000000071,1),(361696272343,1),(85331809838489,1),(227098769,1),(199937,1),(5,3),(41666666666666666666666666669,1),(2,22),(3,8)]}Pairwise coprimality of factors is crucial, because it allows us to process them independently, possibly even in parallel or concurrent fashion.*Following invariant is guaranteed to hold: Eabs (prefValue x) = abs (product (map (uncurry (^)) (prefFactors x)))arithmoiNumber itself.arithmoicList of pairwise coprime (but not neccesarily prime) factors, accompanied by their multiplicities.arithmoiCreate  from a given number. fromValue 123NPrefactored {prefValue = 123, prefFactors = Coprimes {unCoprimes = [(123,1)]}}arithmoiCreate ` from a given list of pairwise coprime (but not neccesarily prime) factors with multiplicities.4fromFactors (splitIntoCoprimes [(140, 1), (165, 1)])\Prefactored {prefValue = 23100, prefFactors = Coprimes {unCoprimes = [(28,1),(33,1),(5,2)]}}4fromFactors (splitIntoCoprimes [(140, 2), (165, 3)])bPrefactored {prefValue = 88045650000, prefFactors = Coprimes {unCoprimes = [(28,2),(33,3),(5,5)]}}J(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None&'arithmoiA typical arithmetic function operates on the canonical factorisation of a number into prime's powers and consists of two rules. The first one determines the values of the function on the powers of primes. The second one determines how to combine these values into final result.In the following definition the first argument is the function on prime's powers, the monoid instance determines a rule of combination (typically KL or KME), and the second argument is convenient for unwrapping (typically, KN or KO).arithmoi3Convert to a function. The value on 0 is undefined.arithmoi-Convert to a function on prime factorisation.HarithmoiFactorisation is expensive, so it is better to avoid doing it twice. Write 'runFunction (f + g) n' instead of 'runFunction f n + runFunction g n'.P(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None|=?@CEFIMNINM=@?CFE(c) 2019 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%&'-.17=?HPSUVX_garithmoiHThis singleton data type establishes a correspondence between a modulo mC on type level and a cyclic group of the same order on term level.IarithmoiResidues modulo 2.JarithmoiResidues modulo 4.KarithmoiResidues modulo p^k for odd prime p.LarithmoiResidues modulo 2p^k for odd prime p.arithmoiHThis singleton data type establishes a correspondence between a modulo m4 on type level and its factorisation on term level.arithmoi Factors of m.MarithmoiFrom Data.Constraint.Nat.arithmoi.Wrapper to hide an unknown type-level natural.arithmoi,Unidirectional pattern for residues modulo 2p^k for odd prime p.arithmoi+Unidirectional pattern for residues modulo p^k for odd prime p.arithmoi-Unidirectional pattern for residues modulo 4.arithmoi-Unidirectional pattern for residues modulo 2.arithmoi5Create a singleton from a type-level positive modulo m, passed in a constraint.:set -XDataKindssfactors :: SFactors Integer 13+SFactors {sfactorsFactors = [(Prime 13,1)]}arithmoi#Create a singleton from factors of m-. Factors must be distinct, as in output of .import Math.NumberTheory.PrimessomeSFactors (factorise 98)6SFactors {sfactorsFactors = [(Prime 2,1),(Prime 7,2)]}arithmoi$Convert a singleton to a proof that m is known. Usage example: gtoModulo :: SFactors Integer m -> Natural toModulo t = case proofFromSFactors t of Sub Dict -> natVal tarithmoi5Create a singleton from a type-level positive modulo m, passed in a constraint.:set -XDataKindsimport Data.Maybe&cyclicGroup :: CyclicGroup Integer 169CGOddPrimePower' (Prime 13) 2$sfactorsToCyclicGroup (fromModulo 4) Just CG4'/sfactorsToCyclicGroup (fromModulo (2 * 13 ^ 3))*Just (CGDoubleOddPrimePower' (Prime 13) 3)+sfactorsToCyclicGroup (fromModulo (4 * 13))Nothingarithmoi Similar to  . e, but much faster, because it but performes only one primality test instead of full factorisation.arithmoi'Convert a cyclic group to a proof that m is known. Usage example: mtoModulo :: CyclicGroup Integer m -> Natural toModulo t = case proofFromCyclicGroup t of Sub Dict -> natVal tarithmoiwCheck whether a multiplicative group of residues, characterized by its modulo, is cyclic and, if yes, return its form.$sfactorsToCyclicGroup (fromModulo 4) Just CG4'/sfactorsToCyclicGroup (fromModulo (2 * 13 ^ 3))*Just (CGDoubleOddPrimePower' (Prime 13) 3)+sfactorsToCyclicGroup (fromModulo (4 * 13))NothingarithmoiInvert .import Data.MaybeGcyclicGroupToSFactors (fromJust (sfactorsToCyclicGroup (fromModulo 4)))>SFactors {sfactorsModulo = 4, sfactorsFactors = [(Prime 2,2)]}(c) 2012 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None>X"Parithmoi totientSum n is, for n > 0 , the sum of otient k | k <- [1 .. n]]7, computed via generalised Mbius inversion. See  :http://mathworld.wolfram.com/TotientSummatoryFunction.html for the formula used for  totientSum.import Data.ProxyAtotientSum (Proxy :: Proxy Data.Vector.Unboxed.Vector) 100 :: Int3044=totientSum (Proxy :: Proxy Data.Vector.Vector) 100 :: Integer3044arithmoigeneralInversion g n/ evaluates the generalised Mbius inversion of g at the argument n.xThe generalised Mbius inversion implemented here allows an efficient calculation of isolated values of the function  f : N -> Z if the function g defined by , g n = sum [f (n `quot` k) | k <- [1 .. n]] Jcan be cheaply computed. By the generalised Mbius inversion formula, then 8 f n = sum [moebius k * g (n `quot` k) | k <- [1 .. n]]  which allows the computation in O(n) steps, if the values of the Mbius function are known. A slightly different formula, used here, does not need the values of the Mbius function and allows the computation in O(n^0.75) steps, using O(n^0.5) memory.}An example of a pair of such functions where the inversion allows a more efficient computation than the direct approach is x f n = number of reduced proper fractions with denominator <= n g n = number of proper fractions with denominator <= n (a proper fraction is a fraction  0 < n/d < 1). Then f n6 is the cardinality of the Farey sequence of order nr (minus 1 or 2 if 0 and/or 1 are included in the Farey sequence), or the sum of the totients of the numbers  2 <= k <= n. f n0 is not easily directly computable, but then g n = n*(n-1)/2\ is very easy to compute, and hence the inversion gives an efficient method of computing f n.ISince the function arguments are used as array indices, the domain of f is restricted to . The value f n is then computed by generalInversion g n#. Note that when many values of ft are needed, there are far more efficient methods, this method is only appropriate to compute isolated values of f.(c) 2011 Daniel FischerMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%9arithmoiList all modular square roots.:set -XDataKindssqrtsMod sfactors (1 :: Mod 60)[(1 `modulo` 60),(49 `modulo` 60),(41 `modulo` 60),(29 `modulo` 60),(31 `modulo` 60),(19 `modulo` 60),(11 `modulo` 60),(59 `modulo` 60)]arithmoibList all square roots modulo a number, the factorisation of which is passed as a second argument.&sqrtsModFactorisation 1 (factorise 60)[1,49,41,29,31,19,11,59]arithmoi2List all square roots modulo the power of a prime.import Data.Maybeimport Math.NumberTheory.Primes-sqrtsModPrimePower 7 (fromJust (isPrime 3)) 2[4,5]-sqrtsModPrimePower 9 (fromJust (isPrime 3)) 3[3,12,21,24,6,15]arithmoi&List all square roots by prime modulo.import Data.Maybeimport Math.NumberTheory.Primes&sqrtsModPrime 1 (fromJust (isPrime 5))[1,4]&sqrtsModPrime 0 (fromJust (isPrime 5))[0]&sqrtsModPrime 2 (fromJust (isPrime 5))[]NarithmoisqrtModP' square prime finds a square root of square modulo prime. prime must be a (positive) prime, and square must+ be a positive quadratic residue modulo prime, i.e. 'jacobi square prime == 1.Oarithmoip must be of form 4k + 1ParithmoitonelliShanks square prime calculates a square root of square modulo prime, where prime is a prime of the form 4*k + 1 and square( is a positive quadratic residue modulo prime(, using the Tonelli-Shanks algorithm.Qarithmoi/prime must be odd, n must be coprime with prime'(c) 2016 Chris Fredrickson, Google Inc.MIT1Chris Fredrickson <chris.p.fredrickson@gmail.com>None7HVN, arithmoi<A Gaussian integer is a+bi, where a and b are both integers.arithmoiThe imaginary unit, where   .^ 2 == -1arithmoiConjugate a Gaussian integer.arithmoi2The square of the magnitude of a Gaussian integer.Rarithmoi2Compute whether a given Gaussian integer is prime.arithmoiAn infinite list of the Gaussian primes. Uses primes in Z to exhaustively generate all Gaussian primes (up to associates), in order of ascending magnitude.arithmoiTFind a Gaussian integer whose norm is the given prime number of form 4k + 1 using  dhttp://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdfHermite-Serret algorithm.SarithmoiTCompute the prime factorisation of a Gaussian integer. This is unique up to units (+- 1, +3- i). Unit factors are not included in the result.TarithmoiHRemove all (1:+1) factors from the argument, avoiding complex division.UarithmoiJRemove p and conj p factors from the argument, avoiding complex division.UarithmoiGaussian prime parithmoi%Precomputed norm of p, of form 4k + 1arithmoiHExpected number of factors (either p or conj p) in Gaussian integer zarithmoiGaussian integer z  6"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>None7HSVX(arithmoiAn Eisenstein integer is a + b, where a and b are both integers.,arithmoi1The imaginary unit for Eisenstein integers, where , == (-1/2) + ((sqrt 3)/2) == exp(2*pi*/3)and " is the usual imaginary unit with  == -1.Varithmoi Returns an EisensteinInteger1's sign, and its associate in the first sextant.-arithmoiSList of all Eisenstein units, counterclockwise across all sextants, starting with 1..arithmoiProduce a list of an EisensteinInteger's associates.Warithmoi3Function that does most of the underlying work for divMod and quotRemu, apart from choosing the specific integer division algorithm. This is instead done by the calling function (either divMod which uses div, or quotRem , which uses quot.)/arithmoiConjugate a Eisenstein integer.0arithmoi4The square of the magnitude of a Eisenstein integer.XarithmoiChecks if a given EisensteinInteger is prime. EisensteinInteger&s whose norm is a prime congruent to 0 or 1 modulo 3 are prime. See  "http://thekeep.eiu.edu/theses/2467BBandara, Sarada, "An Exposition of the Eisenstein Integers" (2016) , page 12.YarithmoiRemove 1 -  factors from an EisensteinIntegerI, and calculate that prime's multiplicity in the number's factorisation.1arithmoiMFind an Eisenstein integer whose norm is the given prime number in the form 3k + 1 using a modification of the  dhttp://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdfHermite-Serret algorithm.The maintainer  Dhttps://github.com/cartazio/arithmoi/pull/121#issuecomment-415010647Andrew Lelechenko derived the following:Each prime of the form 3n + 1 is actually of the form 6k + 1.One has X(z + 3k)^2 "a z^2 + 6kz + 9k^2 "a z^2 + (6k + 1)z - z + 9k^2 "a z^2 - z + 9k^2 (mod 6k + 1).The goal is to solve z^2 - z + 1 "a 0 (mod 6k + 1) . One has:  z^2 - z + 1 "a 0 (mod 6k + 1) z^2 - z "a -1 (mod 6k + 1) &z^2 - z + 9k^2 "a 9k^2 - 1 (mod 6k + 1) "(z + 3k)^2 "a 9k^2 - 1 (mod 6k + 1) -z + 3k = sqrtsModPrime(9k^2 - 1) (mod 6k + 1) /z = (sqrtsModPrime(9k^2 - 1) (mod 6k + 1)) - 3kFor example, let p = 7, then k = 1. Square root of 9*1^2-1 "a 1 (mod 7), and z = 1 - 3*1 = -2 "a 5 (mod 7).Truly, +norm (5 :+ 1) = 25 - 5 + 1 = 21 "a 0 (mod 7).2arithmoi6An infinite list of Eisenstein primes. Uses primes in ZL to exhaustively generate all Eisenstein primes in order of ascending norm.HEvery prime is in the first sextant, so the list contains no associates.MEisenstein primes from the whole complex plane can be generated by applying . to each prime in this list.Zarithmoi +Implementation notes for factorise function8Compute the prime factorisation of a Eisenstein integer. This function works by factorising the norm of an Eisenstein integer and then, for each prime factor, finding the Eisenstein prime whose norm is said prime factor with  findPrime.zThis is only possible because the norm function of the Euclidean Domain of Eisenstein integers is multiplicative: #norm (e1 * e2) == norm e1 * norm e2 for any two EisensteinIntegers e1, e2.!In the previously mentioned work  "http://thekeep.eiu.edu/theses/2467BBandara, Sarada, "An Exposition of the Eisenstein Integers" (2016)`, in Theorem 8.4 in Chapter 8, a way is given to express any Eisenstein integer  as :(-1)^a * ^b * (1 - )^c * product [_i^a_i | i <- [1..N]] where  a, b, c, a_i are nonnegative integers, N > 1 is an integer and _i are Eisenstein primes.Aplying norm0 to both sides of the equation from Theorem 8.4:  Pnorm  = norm ( (-1)^a * ^b * (1 - )^c * product [ _i^a_i | i <- [1..N]] ) == cnorm  = norm ((-1)^a) * norm (^b) * norm ((1 - )^c) * norm (product [ _i^a_i | i <- [1..N]]) == cnorm  = (norm (-1))^a * (norm )^b * (norm (1 - ))^c * product [ norm (_i^a_i) | i <- [1..N]] == dnorm  = (norm (-1))^a * (norm )^b * (norm (1 - ))^c * product [ (norm _i)^a_i) | i <- [1..N]] == Fnorm  = 1^a * 1^b * 3^c * product [ (norm _i)^a_i) | i <- [1..N]] == :norm  = 3^c * product [ (norm _i)^a_i) | i <- [1..N]] ==where  a, b, c, a_i are nonnegative integers, and N > 1 is an integer.nThe remainder of the Eisenstein integer factorisation problem is about finding appropriate Eisenstein primes [e_i | i <- [1..M]] such that <map norm [e_i | i <- [1..M]] == map norm [_i | i <- [1..N]] where  1 < N <= M are integers and ==5 is equality on sets (i.e.duplicates do not matter).NB: The reason M >= N is because the prime factors of an Eisenstein integer may include a prime factor and its conjugate (both have the same norm), meaning the number may have more Eisenstein prime factors than its norm has integer prime factors.[arithmoiRemove p and  conjugate p# factors from the argument, where p is an Eisenstein prime.\arithmoi0Divide an Eisenstein integer by an even integer.3arithmoi1See the source code and Haddock comments for the  factorise and isPrimeO functions in this module (they are not exported) for implementation details.[arithmoiEisenstein prime parithmoi Conjugate of parithmoiPrecomputed norm of p , of form 4k + 1arithmoi#Expected number of factors (either p or  conjugate p) in Eisenstein integer zarithmoiEisenstein integer z ()*+,-./012 ()*+,/0.-12)6(c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%X&>arithmoi)Find all solutions of ax + b "a 0 (mod m).:set -XDataKinds:solveLinear (6 :: Mod 10) 4 -- solving 6x + 4 "a 0 (mod 10)![(1 `modulo` 10),(6 `modulo` 10)]?arithmoi/Find all solutions of ax + bx + c "a 0 (mod m).:set -XDataKindsMsolveQuadratic sfactors (1 :: Mod 32) 0 (-17) -- solving x - 17 "a 0 (mod 32)C[(9 `modulo` 32),(25 `modulo` 32),(7 `modulo` 32),(23 `modulo` 32)]>arithmoiaarithmoibarithmoi list of x?arithmoiaarithmoibarithmoicarithmoi list of x>?>?"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>None>X:@arithmoi Evaluate the Q" function over a block. Value at 0), if zero falls into block, is undefined.This function should **not**_ be used with a negative lower bound. If it is, the result is undefined. Furthermore, do not: use a block length greater than maxBound :: Int.use a power that is either of 0, 1.VNone of these preconditions are checked, and if any occurs, the result is undefined, if the function terminates.sieveBlockNFree 2 1 106[True,True,True,False,True,True,True,False,False,True]Aarithmoi&For a given nonnegative integer power n, generate all n/-free numbers in ascending order, starting at 1.When n is 0 or 1, the resulting list is [1].Barithmoi Generate n-free numbers in a block starting at a certain value. The length of the list is determined by the value passed in as the third argument. It will be lesser than or equal to this value.aThis function should not be used with a negative lower bound. If it is, the result is undefined.The block length cannot exceed maxBound :: Int$, this precondition is not checked.As with nFrees , passing n = 0, 1 results in an empty list.@arithmoi Power whose n-freedom will be checked.arithmoiLower index of the block.arithmoiLength of the block.arithmoiVector of flags, where True at index i means the i-th element of the block is n-free.AarithmoiPower n to be used to generate n-free numbers.arithmoiGenerated infinite list of n-free numbers.BarithmoiPower n to be used to generate n-free numbers.arithmoiStarting number in the block.arithmoi,Maximum length of the block to be generated.arithmoiGenerated list of n-free numbers.@ABAB@ (c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None @AFHVCarithmoi$Represents three possible values of  -https://en.wikipedia.org/wiki/Mbius_functionMbius function.Darithmoi1Earithmoi0Farithmoi1GarithmoiConvert to any numeric type.Harithmoi5Evaluate the Mbius function over a block. Value of f. at 0, if zero falls into block, is undefined.,Based on the sieving algorithm from p. 3 of  $https://arxiv.org/pdf/1610.08551.pdfRComputations of the Mertens function and improved bounds on the Mertens conjecture1 by G. Hurst. It is approximately 5x faster than #R.sieveBlockMoebius 1 10[[MoebiusP,MoebiusN,MoebiusN,MoebiusZ,MoebiusN,MoebiusP,MoebiusN,MoebiusZ,MoebiusZ,MoebiusP]CDEFGHCDEFGHS(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NoneX!QarithmoiYCreate a multiplicative function from the function on prime's powers. See examples below.RarithmoiSee S.Sarithmoi2The set of all (positive) divisors of an argument.TarithmoiSee U.UarithmoiVThe unsorted list of all (positive) divisors of an argument, produced in lazy fashion.VarithmoiSee W.WarithmoiSame as R:, but with better performance on cost of type restriction.Xarithmoi Synonym for Y.map divisorCount [1..10][1,2,2,3,2,4,2,4,3,4]YarithmoiSee Z.Zarithmoi1The number of (positive) divisors of an argument. %tauA = multiplicative (\_ k -> k + 1)[arithmoiSee \.\arithmoiThe sum of the k1-th powers of (positive) divisors of an argument. HsigmaA = multiplicative (\p k -> sum $ map (p ^) [0..k]) sigmaA 0 = tauA]arithmoiSee ^.^arithmoi,Calculates the totient of a positive number n, i.e. the number of k with  1 <= k <= n and ] n k == 18, in other words, the order of the group of units in !$/(n)._arithmoiSee `.`arithmoi3Calculates the k-th Jordan function of an argument. jordanA 1 = totientAaarithmoiSee b.barithmoiCalculates the  4https://en.wikipedia.org/wiki/Ramanujan_tau_functionRamanujan tau function of a positive number n, using formulas given (http://www.numbertheory.org/php/tau.htmlherecarithmoiSee d.darithmoi.Calculates the Mbius function of an argument.earithmoiSee f.farithmoi1Calculates the Liouville function of an argument.garithmoiSee h.harithmoixCalculates the Carmichael function for a positive integer, that is, the (smallest) exponent of the group of units in !$/(n).iarithmoiTCreate an additive function from the function on prime's powers. See examples below.jarithmoiSee k.karithmoi!Number of distinct prime factors. "smallOmegaA = additive (\_ _ -> 1)larithmoiSee m.marithmoi3Number of prime factors, counted with multiplicity.  bigOmegaA = additive (\_ k -> k)narithmoiSee o.oarithmoi+The exponent of von Mangoldt function. Use log expMangoldtA) to recover von Mangoldt function itself.parithmoiSee q.qarithmoiCheck if an integer is n-free. An integer x is nf-free if in its factorisation into prime factors, no factor has an exponent larger than or equal to n.(ABCDEFGQRSTUVWXYZ[\]^_`abcdefghijklmnopq(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None,ABCDEFGQRSTUVWXYZ[\]^_`abcdefghijklmnopq,QRSTUVWXYZ[\]^_`abcdCDEFGefijklmghnopqAB!(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%17>P_rarithmoir) m is a type which is only inhabited by  5https://en.wikipedia.org/wiki/Primitive_root_modulo_nprimitive roots of m.sarithmoiExtract primitive root value.tarithmoi,Check whether a given modular residue is a  5https://en.wikipedia.org/wiki/Primitive_root_modulo_nprimitive root.:set -XDataKindsimport Data.Maybe4isPrimitiveRoot (fromJust cyclicGroup) (1 :: Mod 13)Nothing4isPrimitiveRoot (fromJust cyclicGroup) (2 :: Mod 13)PJust (PrimitiveRoot {unPrimitiveRoot = MultMod {multElement = (2 `modulo` 13)}})rstrst"(c) 2018 Bhavik MehtaMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None%X warithmoiComputes the discrete logarithm. Currently uses a combination of the baby-step giant-step method and Pollard's rho algorithm, with Bach reduction.:set -XDataKindsimport Data.Maybe7let cg = fromJust cyclicGroup :: CyclicGroup Integer 13(let rt = fromJust (isPrimitiveRoot cg 2)$let x = fromJust (isMultElement 11)discreteLogarithm cg rt x7wwT(c) 2011 Daniel FischerMIT1Daniel Fischer <daniel.is.fischer@googlemail.com>None d& !"#$%&'()*jmlknyz{|}~rstw#(c) 2017 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>NoneFX1 xarithmoi>A record, which specifies a function to evaluate over a block.>For example, here is a configuration for the totient function: SieveBlockConfig { sbcEmpty = 1 , sbcFunctionOnPrimePower = \p a -> (unPrime p - 1) * unPrime p ^ (a - 1) , sbcAppend = (*) }zarithmoivalue of a function on 1{arithmoi*how to evaluate a function on prime powers|arithmoi8how to combine values of a function on coprime arguments}arithmoiRCreate a config for a multiplicative function from its definition on prime powers.~arithmoiMCreate a config for an additive function from its definition on prime powers.arithmoi f x l8 evaluates an arithmetic function for integers between x and x+l-1 and returns a vector of length ld. It completely avoids factorisation, so it is asymptotically faster than pointwise evaluation of f. Value of f. at 0, if zero falls into block, is undefined.vBeware that for underlying non-commutative monoids the results may potentially differ from pointwise application via .This is a thin wrapper over , read more details there.,import Math.NumberTheory.ArithmeticFunctions%runFunctionOverBlock carmichaelA 1 10[1,1,2,2,4,2,6,2,6,4]arithmoiTEvaluate a function over a block in accordance to provided configuration. Value of f. at 0, if zero falls into block, is undefined.Based on Algorithm M of  #https://arxiv.org/pdf/1305.1639.pdf\Parity of the number of primes in a given interval and algorithms of the sublinear summation by A. V. Lelechenko. See Lemma 2 on p. 5 on its algorithmic complexity. For the majority of use-cases its time complexity is O(x^(1+)).9For example, following code lists smallest prime factors:CsieveBlock (SieveBlockConfig maxBound (\p _ -> unPrime p) min) 2 10[2,3,2,5,2,7,2,3,2,11]4And this is how to factorise all numbers in a block:EsieveBlock (SieveBlockConfig [] (\p k -> [(unPrime p, k)]) (++)) 2 10^[[(2,1)],[(3,1)],[(2,2)],[(5,1)],[(2,1),(3,1)],[(7,1)],[(2,3)],[(3,2)],[(2,1),(5,1)],[(11,1)]]arithmoiThis is  specialized to unboxed vectors.?sieveBlockUnboxed (SieveBlockConfig 1 (\_ a -> a + 1) (*)) 1 10[1,2,2,3,2,4,2,4,3,4] Hxyz{|}~ xyz{|}~H$(c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None_9arithmoiKCompute 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],The implementation follows Theorem 3.1 from  $https://arxiv.org/pdf/1610.08551.pdfRComputations of the Mertens function and improved bounds on the Mertens conjecture/ by G. Hurst, excluding segmentation of sieves.^arithmoi0Compute sum (map (x -> runMoebius (mu x) * f x))%(c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None>SXm arithmoi#Wrapper to use in conjunction with  and -. Extracts the minimal preimage of function.arithmoi#Wrapper to use in conjunction with  and -. Extracts the maximal preimage of function.arithmoi#Wrapper to use in conjunction with  and -. Extracts the minimal preimage of function.arithmoi#Wrapper to use in conjunction with  and -. Extracts the maximal preimage of function._arithmoiYConvert a list of powers of a prime into an atomic Dirichlet series (Section 4, Step 2).`arithmoiSee section 5.1 of the paper.aarithmoiSee section 5.2 of the paper.barithmoibInstead of multiplying all atomic series and filtering out everything, which is not divisible by np, we'd rather split all atomic series into a couple of batches, each of which corresponds to a prime factor of n:. This allows us to crop resulting Dirichlet series (see c calls in d9 below) at the end of each batch, saving time and memory.darithmoiMain workhorse.arithmoiThe inverse for ] function.'The return value is parameterized by a e, which allows various applications by providing different (multiplicative) embeddings. E. g., list all preimages (see a helper ):import qualified Data.Set as Simport Data.SemigroupFS.mapMonotonic getProduct (inverseTotient (S.singleton . Product) 120)NfromList [143,155,175,183,225,231,244,248,286,308,310,350,366,372,396,450,462]Count preimages:inverseTotient (const 1) 12017Sum preimages:inverseTotient id 1204904#Find minimal and maximal preimages:&unMinWord (inverseTotient MinWord 120)143&unMaxWord (inverseTotient MaxWord 120)462arithmoiThe inverse for [ 1 function.'The return value is parameterized by a e, which allows various applications by providing different (multiplicative) embeddings. E. g., list all preimages (see a helper ):import qualified Data.Set as Simport Data.SemigroupDS.mapMonotonic getProduct (inverseSigma (S.singleton . Product) 120)fromList [54,56,87,95]Count preimages:inverseSigma (const 1) 1204Sum preimages:inverseSigma id 120292#Find minimal and maximal preimages:$unMinWord (inverseSigma MinWord 120)54$unMaxWord (inverseSigma MaxWord 120)95arithmoi)Helper to extract a set of preimages for  or ._arithmoi&How to inject a number into a semiringarithmoi,Arithmetic function, which we aim to inversearithmoiList of powers of a primearithmoiAtomic Dirichlet series`arithmoi0Factorisation of a value of the totient functionarithmoi=Possible prime factors of an argument of the totient functionaarithmoi8Factorisation of a value of the sum-of-divisors functionarithmoiEPossible prime factors of an argument of the sum-of-divisors functionbarithmoi,Arithmetic function, which we aim to inversearithmoi3Factorisation of a value of the arithmetic functionarithmoi@Possible prime factors of an argument of the arithmetic functionarithmoiFBatches, corresponding to each element of the factorisation of a valuedarithmoi&How to inject a number into a semiringarithmoi,Arithmetic function, which we aim to inversearithmoi2How to find possible prime factors of the argumentarithmoi9Value of the arithmetic function, which we aim to inversearithmoi(Semiring element, representing preimagesU(c) 2018 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>None_z1farithmoi$Straightforward computation of [ n gx x | x <- [hi, hi - 1 .. lo] ]. Unfortunately, such list generator performs poor, so we fall back to manual recursion.harithmoi0bresenham(x+1) -> bresenham(x) for x >= (2n)^1/3iarithmoi"Division-free computation of [ n gv x | x <- [hi, hi - 1 .. lo] ]. In other words, we compute y-coordinates of highest integral points under hyperbola  x * y = n between x = lo and x = hi in reverse order.(The implementation follows section 5 of  #https://arxiv.org/pdf/1206.3369.pdfQA successive approximation algorithm for computing the divisor summatory functionD by R. Sladkey. It is 2x faster than a trivial implementation for .iV"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>SafeWjarithmoihJoins two lists element-by-element together into one, starting with the first one provided as argument.(take 10 $ intertwine [0, 2 ..] [1, 3 ..][0,1,2,3,4,5,6,7,8,9]karithmoiUSkips every odd-indexed element from an infinite list. Do NOT use with finite lists.take 10 (skipOdds [0, 1 ..])[0,2,4,6,8,10,12,14,16,18]larithmoiVSkips every even-indexed element from an infinite list. Do NOT use with finite lists.take 10 (skipEvens [0, 1 ..])[1,3,5,7,9,11,13,15,17,19]jklW"(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>SafeXarithmoi-Values of Hurwitz zeta function evaluated at (s, a) for  s " [0, 1 ..].RThe algorithm used was based on the Euler-Maclaurin formula and was derived from  %http://fredrikj.net/thesis/thesis.pdfDFast and Rigorous Computation of Special Functions to High Precision by F. Johansson, chapter 4.8, formula 4.8.5. The error for each value in this recurrence is given in formula 4.8.9 as an indefinite integral, and in formula 4.8.12 as a closed form formula. It is the user's responsibility: to provide an appropriate precision for the type chosen.For instance, when using Double/s, it does not make sense to provide a number   < 1e-53 as the desired precision. For Floats, providing an   < 1e-24@ also does not make sense. Example of how to call the function:zetaHurwitz 1e-15 0.25 !! 51024.3489745265808X(c) 2016 Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>SafeX:arithmoi\Infinite sequence of exact values of Riemann zeta-function at even arguments, starting with (0)3. Note that due to numerical errors conversion to m may return values below 1:,approximateValue (zetasEven !! 25) :: Double0.9999999999999996DUse your favorite type for long-precision arithmetic. For instance, YZ works fine:import Data.Number.Fixed2approximateValue (zetasEven !! 25) :: Fixed Prec5041.00000000000000088817842111574532859293035196051773narithmoibInfinite sequence of approximate values of Riemann zeta-function at odd arguments, starting with (1).arithmoi~Infinite sequence of approximate (up to given precision) values of Riemann zeta-function at integer arguments, starting with (0). take 5 (zetas 1e-14) :: [Double]H[-0.5,Infinity,1.6449340668482264,1.2020569031595942,1.0823232337111381]Beware to force evaluation of  zetas !! 1 if the type a2 does not support infinite values (for instance, YZ).n["(c) 2018 Alexandre Rodrigues BaldMIT4Alexandre Rodrigues Bald <alexandrer_b@outlook.com>SafeXarithmoi]Infinite sequence of exact values of Dirichlet beta-function at odd arguments, starting with (1).+approximateValue (betasOdd !! 25) :: Double0.9999999999999987import Data.Number.Fixed1approximateValue (betasOdd !! 25) :: Fixed Prec5040.99999999999999999999999960726927497384196726751694oarithmoi9Infinite sequence of approximate values of the Dirichlet = function at positive even integer arguments, starting with (0).arithmoiInfinite sequence of approximate (up to given precision) values of Dirichlet beta-function at integer arguments, starting with (0). take 5 (betas 1e-14) :: [Double]Q[0.5,0.7853981633974483,0.9159655941772189,0.9689461462593694,0.9889445517411051]o&5(c) 2018 Alexandre Rodrigues Bald, Andrew LelechenkoMIT/Andrew Lelechenko <andrew.lelechenko@gmail.com>SafeXp\]^_`a_`b_`c_`d_`e_`f_`g_`h_`i_`j_`k_`l_`mnopqrstuvwxyyz{|}~BC,,0,,1-- 3 6 7     ;F;;;;???????@@AAAABDEDDDDDDDD  E F . / G G G G G GGGGGGGGGGGGGGGGGGG G!GG"G#G$G%G&G'G(G)G*H+,-./0123456J7J7J8J9:;<==>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvdefwxyhijz{|}~              SSSSSSSSSSSSSISSSSSSSSSSSSSSSSSSSQS!!!!!"##########R$%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%WXX[[\]\]\]\]\]\]\]\]\]\]\]\]\]\]\]\\\\**\\+]+^\-  7 7                    999:::: :!\"):#$:%:&:':(:):*:+,-:.:/:0:1:2:3:4:5:6\7;8<9<:<;<<<=<><?<@<A<B<C<D<9=E=F=G=H=I==J=K=L>M>N>O>P>Q>R>S>T\U?V?W?X?Y?Z?[?\?]?^?_@`@a\bcdecdfDgDhDiDjGkGlmGnGoGpGqGrGsGtGuGvGwGxGyGzG{G|G}G~GGGGGGGGJFEFE\c$%%%%\I%_U\hUUVVVX[(arithmoi-0.10.0.0-ERYv2VuF0BeGX7OVoZ5V00Math.NumberTheory.Moduli.ClassMath.NumberTheory.Euclidean$Math.NumberTheory.Euclidean.CoprimesMath.NumberTheory.Powers.CubesMath.NumberTheory.Powers.Fourth Math.NumberTheory.Powers.Squares!Math.NumberTheory.Primes.CountingMath.NumberTheory.Primes.SmallMath.NumberTheory.Primes$Math.NumberTheory.Recurrences.Linear&Math.NumberTheory.Recurrences.BilinearMath.NumberTheory.RecurrencesMath.NumberTheory.Moduli.Jacobi Math.NumberTheory.Primes.Testing Math.NumberTheory.Moduli.ChineseMath.NumberTheory.Primes.Sieve&Math.NumberTheory.Primes.Factorisation Math.NumberTheory.Powers.Modular Math.NumberTheory.Powers.GeneralMath.NumberTheory.SmoothNumbers-Math.NumberTheory.Primes.Testing.Certificates0Math.NumberTheory.Primes.Factorisation.CertifiedMath.NumberTheory.Prefactored%Math.NumberTheory.ArithmeticFunctions"Math.NumberTheory.Moduli.Singleton"Math.NumberTheory.MoebiusInversionMath.NumberTheory.Moduli.Sqrt,Math.NumberTheory.Quadratic.GaussianIntegers.Math.NumberTheory.Quadratic.EisensteinIntegers"Math.NumberTheory.Moduli.Equations.Math.NumberTheory.ArithmeticFunctions.NFreedom-Math.NumberTheory.ArithmeticFunctions.Moebius&Math.NumberTheory.Moduli.PrimitiveRoot*Math.NumberTheory.Moduli.DiscreteLogarithm0Math.NumberTheory.ArithmeticFunctions.SieveBlock-Math.NumberTheory.ArithmeticFunctions.Mertens-Math.NumberTheory.ArithmeticFunctions.InverseMath.NumberTheory.ZetaGHC.TypeNats.CompatNumeric.NaturalNatural)Math.NumberTheory.Powers.Squares.Internal&Math.NumberTheory.MoebiusInversion.Int-Math.NumberTheory.Primes.Counting.ApproximateMath.NumberTheory.Primes.Types nextPrime precPrimeapproxPrimeCountnthPrimeApproxMath.Combinat.NumbersbinomialunsignedStirling1st stirling2nd bernoulli(Math.NumberTheory.Recurrences.PentagonalMath.NumberTheory.Unsafe'Math.NumberTheory.Primes.Sieve.IndexingMath.NumberTheory.Utils.Math.NumberTheory.Primes.Testing.Probabilistic#Math.NumberTheory.Curves.Montgomery'Math.NumberTheory.Utils.DirichletSeries$Math.NumberTheory.Utils.FromIntegral+Math.NumberTheory.Primes.Sieve.Eratosthenes4Math.NumberTheory.Primes.Factorisation.TrialDivision&Math.NumberTheory.Primes.Counting.ImplpowMod powSomeMod1Math.NumberTheory.Primes.Factorisation.Montgomery factoriseisPrime6Math.NumberTheory.Primes.Testing.Certificates.Internal*Math.NumberTheory.Primes.Testing.Certifiedtotient+Math.NumberTheory.ArithmeticFunctions.ClassData.SemigroupProductSum getProductgetSumMath.NumberTheory.PowersisNFreesieveBlockUnboxed.Math.NumberTheory.ArithmeticFunctions.StandardMath.NumberTheory.Moduli!Math.NumberTheory.Utils.HyperbolaMath.NumberTheory.Zeta.UtilsMath.NumberTheory.Zeta.HurwitzMath.NumberTheory.Zeta.RiemannData.Number.FixedFixed Math.NumberTheory.Zeta.Dirichletbase GHC.TypeNatsKnownNat&semirings-0.5.1-70vLMI5c3D2BCVQ0mWB8qzData.Euclideancoprimelcmgcddivide GcdDomaindegreeremquotquotRem EuclideanunwrapIntegral WrapIntegralWrappedIntegralisUnit extendedGCDCoprimes unCoprimes singletoninsertsplitIntoCoprimes$fMonoidCoprimes$fSemigroupCoprimes $fEqCoprimes$fShowCoprimesSomeModInfModMultMod multElementModgetMod getNatModgetVal getNatVal invertMod^% isMultElement invertGroupmodulo invertSomeMod$fFractionalMod$fNumMod $fBoundedMod $fShowMod$fBoundedMultMod$fMonoidMultMod$fSemigroupMultMod$fFractionalSomeMod $fNumSomeMod $fShowSomeMod $fOrdSomeMod $fEqSomeMod$fEqMod$fOrdMod $fEnumMod $fEqMultMod $fOrdMultMod $fShowMultModintegerCubeRootintegerCubeRoot' exactCubeRootisCubeisCube'isPossibleCubeintegerFourthRootintegerFourthRoot'exactFourthRoot isFourthPowerisFourthPower'isPossibleFourthPowerintegerSquareRootintegerSquareRoot'integerSquareRootRemintegerSquareRootRem'exactSquareRootisSquare isSquare'isPossibleSquareisPossibleSquare2!approxPrimeCountOverestimateLimit nthPrimeApproxUnderestimateLimit smallPrimesPrimeunPrime factorial fibonacci fibonacciPairlucas lucasPair generalLucas stirling1 stirling2lah eulerian1 eulerian2 faulhaberPolyeuler eulerPolyAt1 partition JacobiSymbolMinusOneZeroOnejacobi$fMonoidJacobiSymbol$fSemigroupJacobiSymbol$fEqJacobiSymbol$fOrdJacobiSymbol$fShowJacobiSymbol millerRabinVisStrongFermatPP isFermatPP bailliePSWchineseCoprimechinesechineseCoprimeSomeModchineseSomeModchineseRemainderchineseRemainder2 PrimeSieve primeSieve primeListprimes psieveList sieveFrom psieveFromtrialDivisionTotrialDivisionPrimeToprimeCountMaxArg primeCountnthPrimeMaxArgnthPrime powModWord powModInt integerRoot exactRoot isKthPowerisPerfectPower highestPower largePFPower factorise'stepFactorisationdefaultStdGenFactorisationdefaultStdGenFactorisation'stdGenFactorisationcurveFactorisationmontgomeryFactorisation smallFactorsUniqueFactorisation factorBack $fEnumPrime $fEnumPrime0 $fEnumPrime1 $fEnumPrime2$fUniqueFactorisationNatural$fUniqueFactorisationInteger$fUniqueFactorisationWord$fUniqueFactorisationInt SmoothBasisfromSetfromListfromSmoothUpperBound smoothOver' smoothOversmoothOverInRangesmoothOverInRangeBFisSmooth$fEqSmoothBasis$fShowSmoothBasisPrimalityArgumentPockDivisionObvious Assumptionaprime largeFactor smallFactor factorListalimitPrimalityProofcprimeCompositenessArgumentDivisorsFermatLucasBeliefcompo firstDivisor secondDivisor fermatBaseCompositenessProof composite Certificate CompositeargueCertificatearguePrimalityverifyPrimalityArgumentargueCompositenessverifyCompositenessArgumentcheckCertificatecheckCompositenessProofcheckPrimalityProofcertifyisCertifiedPrimecertifiedFactorisationcertificateFactorisationprovenFactorisation Prefactored prefValue prefFactors fromValue fromFactors $fUniqueFactorisationPrefactored$fNumPrefactored$fShowPrefactoredArithmeticFunction runFunctionrunFunctionOnFactors CyclicGroupSFactors unSFactorsSomeCGDoubleOddPrimePowerCGOddPrimePowerCG4CG2sfactors someSFactorsproofFromSFactors cyclicGroupcyclicGroupFromFactorscyclicGroupFromModuloproofFromCyclicGroupsfactorsToCyclicGroupcyclicGroupToSFactors $fNFDataSome $fShowSome $fOrdSome$fEqSome$fNFDataSFactors $fOrdSFactors $fEqSFactors $fNFDataSome0 $fShowSome0 $fOrdSome0 $fEqSome0$fNFDataCyclicGroup$fOrdCyclicGroup$fEqCyclicGroup$fShowSFactors$fGenericSFactors$fShowCyclicGroup$fGenericCyclicGroup totientSumgeneralInversionsqrtsModsqrtsModFactorisationsqrtsModPrimePower sqrtsModPrimeGaussianInteger:+realimagι conjugatenorm findPrime$$fUniqueFactorisationGaussianInteger$fEuclideanGaussianInteger$fGcdDomainGaussianInteger$fRingGaussianInteger$fSemiringGaussianInteger$fNumGaussianInteger$fShowGaussianInteger$fNFDataGaussianInteger$fEqGaussianInteger$fOrdGaussianInteger$fGenericGaussianIntegerEisensteinIntegerωids associates&$fUniqueFactorisationEisensteinInteger$fEuclideanEisensteinInteger$fGcdDomainEisensteinInteger$fRingEisensteinInteger$fSemiringEisensteinInteger$fNumEisensteinInteger$fShowEisensteinInteger$fNFDataEisensteinInteger$fEqEisensteinInteger$fOrdEisensteinInteger$fGenericEisensteinInteger solveLinearsolveQuadraticsieveBlockNFreenFrees nFreesBlockMoebiusMoebiusNMoebiusZMoebiusP runMoebiussieveBlockMoebius$fMonoidMoebius$fSemigroupMoebius$fVectorVectorMoebius$fMVectorMVectorMoebius$fUnboxMoebius $fEqMoebius $fOrdMoebius $fShowMoebiusmultiplicativedivisors divisorsA divisorsList divisorsListA divisorsSmalldivisorsSmallA divisorCounttautauAsigmasigmaAtotientAjordanjordanA ramanujan ramanujanAmoebiusmoebiusA liouville liouvilleA carmichael carmichaelAadditive smallOmega smallOmegaAbigOmega bigOmegaA expMangoldt expMangoldtAisNFreeA PrimitiveRootunPrimitiveRootisPrimitiveRoot$fEqPrimitiveRoot$fShowPrimitiveRootdiscreteLogarithmSieveBlockConfigsbcEmptysbcFunctionOnPrimePower sbcAppendmultiplicativeSieveBlockConfigadditiveSieveBlockConfigrunFunctionOverBlock sieveBlockmertens MinNaturalInfinity unMinNatural MaxNatural unMaxNaturalMinWord unMinWordMaxWord unMaxWordinverseTotient inverseSigmaasSetOfPreimages$fShowPrimePowers$fSemiringMaxWord$fSemiringMinWord$fSemiringMaxNatural$fSemiringMinNatural $fEqMaxWord $fOrdMaxWord $fShowMaxWord $fEqMinWord $fOrdMinWord $fShowMinWord$fEqMaxNatural$fOrdMaxNatural$fShowMaxNatural$fEqMinNatural$fOrdMinNatural$fShowMinNatural zetaHurwitz zetasEvenzetasbetasOddbetasghc-prim GHC.TypesNat+*^<=?-CmpNatDivLog2sameNat someNatValnatVal'natValSomeNat<=GHC.EnumminBoundGHC.RealGHC.NumNum FractionalappCuRtisqrtA karatsubaSqrt GHC.MaybeNothingJustIntEnumeuler'helperForB_E_EPpentspentagonalSigns array-0.5.3.0Data.Array.Base castSTUArray unsafeThaw unsafeFreezeunsafeAtboundsUArray unsafeWrite unsafeReadunsafeNewArray_idxPrtoPrimrhoshiftToOddCount shiftOCWordWord shiftOCIntshiftOCInteger integer-gmpGHC.Integer.TypeIntegershiftOCNatural GHC.NaturalbigNatZeroCountBigNat shiftToOdd shiftOInt shiftOWord shiftOInteger shiftToOdd#shiftToOddCount# bitCountWord#GHC.PrimWord# bitCountWord bitCountInt splitOff#splitOffmergeByrecipModuncheckedShiftR toWheel30 fromWheel30mod lucasTest SomePointPointpointXpointZpointA24pointNnewPointadddoublemultiply $fShowPoint $fEqPoint timesAndCropDirichletSeriesfromDistinctAscListlookupfilterunionsunionsize wordToInt wordToInteger intToWord intToIntegernaturalToIntegerintegerToNatural integerToWord integerToIntmaxBound doesNotFittakeWhileIncreasingsieveToPS sieveBits sieveRange countFromTo nthPrimeCt countToNthcountAlltrialDivisionWithtrialDivisionPrimeWithabs!random-1.1-3ypV4EIycgb35PKjTYYr5q System.RandomStdGenrandomRbigStepenumAndMultiplyFromThenTo testParms findParms primProof compProofFalsetrivial maybeTrivialisTrivialPrime trivialPrimes smallCert certifyBPSWfindDecomposition findFactorfindLoop bpswMessagefoundfermat TrialDivision PocklingtonTrivialfactorisedPartcofactor knownFactorstdLimitFactors StrongFermatLucasSelfridge firstFactor secondFactorwitnesstest primeCheckcertiFactorisationmergemergeAllinvalid$fNumArithmeticFunctionCG2'CG4'CGOddPrimePower'CGDoubleOddPrimePower'Magic sqrtModP'sqrtOfMinusOne tonelliShanks sqrtModPP' divideByTwo divideByPrime absSignum divHelper divideByThree quotEvenIsumMultMoebius atomicSeries invTotientinvSigmastrategyGHC.ListinvertFunction Data.SemiringSemiringpointsUnderHyperbola0stepBackpointsUnderHyperbola intertwineskipOdds skipEvensDoublezetasOdd betasEven