Math.NumberTheory.Powers.Squares: exported symbols usage examples

Symbols

  • exactSquareRoot See 1 Occurences [+] Collapse [-]
    Found in Math.NumberTheory.Powers from the package arithmoi
    -- Stability:   Provisional
    -- Portability: Non-portable (GHC extensions)
    --
    -- Calculating integer roots, modular powers and related things.
    -- This module reexports the most needed functions from the implementation
    -- modules. The implementation modules provide some additional functions,
    -- in particular some unsafe functions which omit some tests for performance
    -- reasons.
    --
    module Math.NumberTheory.Powers
      ( -- *  Integer Roots
        -- ** Square roots
        integerSquareRoot
      , isSquare
      , exactSquareRoot
        -- ** Cube roots
      , integerCubeRoot
      , isCube
      , exactCubeRoot
        -- ** Fourth roots
    

    integerSquareRoot See 12 Occurences [+] Collapse [-]
    Found in Math.NumberTheory.ArithmeticFunctions.Mertens from the package arithmoi
        u = (integerSquareRoot x + 1) `max` ((integerCubeRoot x) ^ (2 :: Word) `quot` 2)
    
        sfunc :: Word -> Int
        sfunc y
          = 1
          - sum [ U.unsafeIndex mes (fromIntegral $ y `quot` n) |  n <- [y `quot` u + 1 .. kappa] ]
          + fromIntegral kappa * U.unsafeIndex mes (fromIntegral nu)
          - sumMultMoebius lookupMus (\n -> fromIntegral $ y `quot` n) [1 .. nu]
          where
            nu = integerSquareRoot y
            kappa = y `quot` (nu + 1)
    
        -- cacheSize ~ u
        cacheSize :: Word
        cacheSize = u `max` (x `quot` u) `max` integerSquareRoot x
    
        -- 1-based index
        mus :: U.Vector Moebius
        mus = sieveBlockMoebius 1 cacheSize
    

    Found in Math.NumberTheory.ArithmeticFunctions.Mertens from the package arithmoi
    mertens :: Word -> Int
    mertens 0 = 0
    mertens 1 = 1
    mertens x = sumMultMoebius lookupMus (\n -> sfunc (x `quot` n)) [1 .. x `quot` u]
      where
        u = (integerSquareRoot x + 1) `max` ((integerCubeRoot x) ^ (2 :: Word) `quot` 2)
    
        sfunc :: Word -> Int
        sfunc y
          = 1
          - sum [ U.unsafeIndex mes (fromIntegral $ y `quot` n) |  n <- [y `quot` u + 1 .. kappa] ]
          + fromIntegral kappa * U.unsafeIndex mes (fromIntegral nu)
          - sumMultMoebius lookupMus (\n -> fromIntegral $ y `quot` n) [1 .. nu]
          where
            nu = integerSquareRoot y
            kappa = y `quot` (nu + 1)
    
        -- cacheSize ~ u
        cacheSize :: Word
        cacheSize = u `max` (x `quot` u) `max` integerSquareRoot x
    

    Found in Math.NumberTheory.ArithmeticFunctions.Mertens from the package arithmoi
    import Math.NumberTheory.Powers.Squares
    import Math.NumberTheory.ArithmeticFunctions.Moebius
    
    -- | Compute individual values of Mertens function in O(n^(2/3)) time and space.
    --
    -- > > map (mertens . (10 ^)) [0..9]
    -- > [1,-1,1,2,-23,-48,212,1037,1928,-222]
    --
    -- The implementation follows Theorem 3.1 from <https://arxiv.org/pdf/1610.08551.pdf Computations of the Mertens function and improved bounds on the Mertens conjecture> by G. Hurst, excluding segmentation of sieves.
    mertens :: Word -> Int
    mertens 0 = 0
    mertens 1 = 1
    mertens x = sumMultMoebius lookupMus (\n -> sfunc (x `quot` n)) [1 .. x `quot` u]
      where
        u = (integerSquareRoot x + 1) `max` ((integerCubeRoot x) ^ (2 :: Word) `quot` 2)
    
        sfunc :: Word -> Int
        sfunc y
          = 1
          - sum [ U.unsafeIndex mes (fromIntegral $ y `quot` n) |  n <- [y `quot` u + 1 .. kappa] ]
    

    Found in Math.NumberTheory.ArithmeticFunctions.Moebius from the package arithmoi
      where
        lowIndex :: Int
        lowIndex = wordToInt lowIndex'
    
        len :: Int
        len = wordToInt len'
    
        highIndex :: Int
        highIndex = lowIndex + len - 1
    
        -- Bit fiddling in 'mapper' is correct only
        -- if all sufficiently small (<= 191) primes has been sieved out.
        ps :: [Int]
        ps = takeWhile (<= (191 `max` integerSquareRoot highIndex)) $ map fromInteger primes
    
        mapper :: Int -> Word8 -> Word8
        mapper ix val
          | val .&. 0x80 == 0x00
          = 1
    

    Found in Math.NumberTheory.ArithmeticFunctions.Moebius from the package arithmoi
    import Data.Word
    #if __GLASGOW_HASKELL__ < 803
    import Data.Semigroup
    #endif
    import qualified Data.Vector.Generic         as G
    import qualified Data.Vector.Generic.Mutable as M
    import qualified Data.Vector.Primitive as P
    import qualified Data.Vector.Unboxed         as U
    import qualified Data.Vector.Unboxed.Mutable as MU
    import GHC.Exts
    import GHC.Integer.GMP.Internals
    import Unsafe.Coerce
    
    import Math.NumberTheory.Primes (primes)
    import Math.NumberTheory.Powers.Squares (integerSquareRoot)
    import Math.NumberTheory.Utils.FromIntegral (wordToInt)
    
    import Math.NumberTheory.Logarithms
    
    -- | Represents three possible values of <https://en.wikipedia.org/wiki/Möbius_function Möbius function>.
    

    Found in Math.NumberTheory.ArithmeticFunctions.SieveBlock from the package arithmoi
        let lowIndex :: Int
            lowIndex = wordToInt lowIndex'
    
            len :: Int
            len = wordToInt len'
    
        as <- V.unsafeThaw $ V.enumFromN lowIndex' len
        bs <- MV.replicate len empty
    
        let highIndex :: Int
            highIndex = lowIndex + len - 1
    
            ps :: [Int]
            ps = takeWhile (<= integerSquareRoot highIndex) $ map fromInteger primes
    
        forM_ ps $ \p -> do
    
          let p# :: Word#
              !p'@(W# p#) = intToWord p
    

    Found in Math.NumberTheory.ArithmeticFunctions.SieveBlock from the package arithmoi
    import Data.Coerce
    import qualified Data.Vector as V
    import qualified Data.Vector.Mutable as MV
    import GHC.Exts
    #if __GLASGOW_HASKELL__ < 709
    import Data.Monoid
    #endif
    
    import Math.NumberTheory.ArithmeticFunctions.Class
    import Math.NumberTheory.ArithmeticFunctions.Moebius (sieveBlockMoebius)
    import Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed
    import Math.NumberTheory.Logarithms (integerLogBase')
    import Math.NumberTheory.Primes (primes)
    import Math.NumberTheory.Primes.Types
    import Math.NumberTheory.Powers.Squares (integerSquareRoot)
    import Math.NumberTheory.Utils (splitOff#)
    import Math.NumberTheory.Utils.FromIntegral (wordToInt, intToWord)
    
    -- | 'runFunctionOverBlock' @f@ @x@ @l@ evaluates an arithmetic function
    -- for integers between @x@ and @x+l-1@ and returns a vector of length @l@.
    

    Found in Math.NumberTheory.MoebiusInversion from the package arithmoi
    --   The value @f n@ is then computed by @generalInversion g n@. Note that when
    --   many values of @f@ are needed, there are far more efficient methods, this
    --   method is only appropriate to compute isolated values of @f@.
    generalInversion :: (Int -> Integer) -> Int -> Integer
    generalInversion fun n
        | n < 1     = error "Möbius inversion only defined on positive domain"
        | n == 1    = fun 1
        | n == 2    = fun 2 - fun 1
        | n == 3    = fun 3 - 2*fun 1
        | otherwise = fastInvert fun n
    
    fastInvert :: (Int -> Integer) -> Int -> Integer
    fastInvert fun n = big `unsafeAt` 0
      where
        !k0 = integerSquareRoot (n `quot` 2)
        !mk0 = n `quot` (2*k0+1)
        kmax a m = (a `quot` m - 1) `quot` 2
        big = runSTArray $ do
            small <- newArray_ (0,mk0) :: ST s (STArray s Int Integer)
            unsafeWrite small 0 0
    

    Found in Math.NumberTheory.MoebiusInversion.Int from the package arithmoi
    --   The value @f n@ is then computed by @generalInversion g n@. Note that when
    --   many values of @f@ are needed, there are far more efficient methods, this
    --   method is only appropriate to compute isolated values of @f@.
    generalInversion :: (Int -> Int) -> Int -> Int
    generalInversion fun n
        | n < 1     = error "Möbius inversion only defined on positive domain"
        | n == 1    = fun 1
        | n == 2    = fun 2 - fun 1
        | n == 3    = fun 3 - 2*fun 1
        | otherwise = fastInvert fun n
    
    fastInvert :: (Int -> Int) -> Int -> Int
    fastInvert fun n = big `unsafeAt` 0
      where
        !k0 = integerSquareRoot (n `quot` 2)
        !mk0 = n `quot` (2*k0+1)
        kmax a m = (a `quot` m - 1) `quot` 2
        big = runSTUArray $ do
            small <- newArray_ (0,mk0) :: ST s (STUArray s Int Int)
            unsafeWrite small 0 0
    

    Found in Math.NumberTheory.Powers from the package arithmoi
    -- Licence:     MIT
    -- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
    -- Stability:   Provisional
    -- Portability: Non-portable (GHC extensions)
    --
    -- Calculating integer roots, modular powers and related things.
    -- This module reexports the most needed functions from the implementation
    -- modules. The implementation modules provide some additional functions,
    -- in particular some unsafe functions which omit some tests for performance
    -- reasons.
    --
    module Math.NumberTheory.Powers
      ( -- *  Integer Roots
        -- ** Square roots
        integerSquareRoot
      , isSquare
      , exactSquareRoot
        -- ** Cube roots
      , integerCubeRoot
      , isCube
    

    Found in Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed from the package arithmoi
        let lowIndex :: Int
            lowIndex = wordToInt lowIndex'
    
            len :: Int
            len = wordToInt len'
    
        as <- V.unsafeThaw $ V.enumFromN lowIndex' len
        bs <- MV.replicate len empty
    
        let highIndex :: Int
            highIndex = lowIndex + len - 1
    
            ps :: [Int]
            ps = takeWhile (<= integerSquareRoot highIndex) $ map fromInteger primes
    
        forM_ ps $ \p -> do
    
          let p# :: Word#
              !p'@(W# p#) = intToWord p
    

    Found in Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed from the package arithmoi
      , multiplicativeSieveBlockConfig
      , additiveSieveBlockConfig
      , sieveBlockUnboxed
      ) where
    
    import Control.Monad (forM_)
    import Control.Monad.ST (runST)
    import qualified Data.Vector.Unboxed as V
    import qualified Data.Vector.Unboxed.Mutable as MV
    import GHC.Exts
    
    import Math.NumberTheory.ArithmeticFunctions.Moebius (Moebius)
    import Math.NumberTheory.Logarithms (integerLogBase')
    import Math.NumberTheory.Primes (primes)
    import Math.NumberTheory.Powers.Squares (integerSquareRoot)
    import Math.NumberTheory.Utils (splitOff#)
    import Math.NumberTheory.Utils.FromIntegral (wordToInt, intToWord)
    
    -- | A record, which specifies a function to evaluate over a block.
    --
    

    integerSquareRoot' See 17 Occurences [+] Collapse [-]
    Found in Math.NumberTheory.Primes.Factorisation.Montgomery from the package arithmoi
        | n == 1    = []
        | ptest n   = [(n, 1)]
        | otherwise = evalState (fact n digits) seed
          where
            digits :: Int
            digits = fromMaybe 8 mbdigs
    
            ptest :: Integer -> Bool
            ptest = maybe primeTest (\bd k -> k <= bd || primeTest k) primeBound
    
            rndR :: Integer -> State g Integer
            rndR k = state (prng k)
    
            perfPw :: Integer -> (Integer, Int)
            perfPw = maybe highestPower (largePFPower . integerSquareRoot') primeBound
    
            fact :: Integer -> Int -> State g [(Integer, Int)]
            fact 1 _ = return mempty
            fact m digs = do
              let (b1, b2, ct) = findParms digs
    

    Found in Math.NumberTheory.Primes.Factorisation.Montgomery from the package arithmoi
    import Data.IntMap (IntMap)
    import qualified Data.IntMap as IM
    import Data.List (foldl')
    import Data.Maybe
    #if __GLASGOW_HASKELL__ < 803
    import Data.Semigroup
    #endif
    
    import GHC.TypeNats.Compat
    
    import Math.NumberTheory.Curves.Montgomery
    import Math.NumberTheory.GCD (splitIntoCoprimes)
    import Math.NumberTheory.Moduli.Class
    import Math.NumberTheory.Powers.General     (highestPower, largePFPower)
    import Math.NumberTheory.Powers.Squares     (integerSquareRoot')
    import Math.NumberTheory.Primes.Sieve.Eratosthenes
    import Math.NumberTheory.Primes.Sieve.Indexing
    import Math.NumberTheory.Primes.Testing.Probabilistic
    import Math.NumberTheory.Unsafe
    import Math.NumberTheory.Utils
    

    Found in Math.NumberTheory.Primes.Factorisation.TrialDivision from the package arithmoi
    --   in the list is the product of all these. If @n <= bound^2@, this is a
    --   full factorisation, but very slow if @n@ has large prime divisors.
    trialDivisionTo :: Integer -> Integer -> [(Integer,Int)]
    trialDivisionTo bd
        | bd < 100      = trialDivisionTo 100
        | bd < 10000000 = trialDivisionWith (primeList $ primeSieve bd)
        | otherwise     = trialDivisionWith (takeWhile (<= bd) $ (psieveList >>= primeList))
    
    -- | Check whether a number is coprime to all of the numbers in the list
    --   (assuming that list contains only numbers > 1 and is ascending).
    trialDivisionPrimeWith :: [Integer] -> Integer -> Bool
    trialDivisionPrimeWith prs n
        | n < 0     = trialDivisionPrimeWith prs (-n)
        | n < 2     = False
        | otherwise = go n (integerSquareRoot' n) prs
          where
            go !m !r (p:ps) = r < p || m `rem` p /= 0 && go m r ps
            go _ _ _ = True
    
    -- | @'trialDivisionPrimeTo' bound n@ tests whether @n@ is coprime to all primes @<= bound@.
    

    Found in Math.NumberTheory.Primes.Factorisation.TrialDivision from the package arithmoi
    trialDivisionWith :: [Integer] -> Integer -> [(Integer,Int)]
    trialDivisionWith prs n
        | n < 0     = trialDivisionWith prs (-n)
        | n == 0    = error "trialDivision of 0"
        | n == 1    = []
        | otherwise = go n (integerSquareRoot' n) prs
          where
            go !m !r (p:ps)
                | r < p     = [(m,1)]
                | otherwise =
                    case splitOff p m of
                      (0,_) -> go m r ps
                      (k,q) -> (p,k) : if q == 1
                                         then []
                                         else go q (integerSquareRoot' q) ps
            go m _ _    = [(m,1)]
    
    -- | @'trialDivisionTo' 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^2@, this is a
    

    Found in Math.NumberTheory.Primes.Factorisation.TrialDivision from the package arithmoi
        ) where
    
    import Math.NumberTheory.Primes.Sieve.Eratosthenes
    import Math.NumberTheory.Powers.Squares
    import Math.NumberTheory.Utils
    
    -- | Factorise an 'Integer' 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.
    trialDivisionWith :: [Integer] -> Integer -> [(Integer,Int)]
    trialDivisionWith prs n
        | n < 0     = trialDivisionWith prs (-n)
        | n == 0    = error "trialDivision of 0"
        | n == 1    = []
        | otherwise = go n (integerSquareRoot' n) prs
          where
            go !m !r (p:ps)
                | r < p     = [(m,1)]
                | otherwise =
                    case splitOff p m of
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
        go2 ps = go 1 ps
        go !acc ((p, 1) : pps) = go (lcm acc (p - 1)) pps
        go acc ((p, k) : pps)  = go ((lcm acc (p - 1)) * p ^ (k - 1)) pps
        go acc []              = acc
    
    -- NOTE: This is a legacy implementation of FactorSieve which uses the
    --       same (2,3,5) wheel optimization as the other sieves.
    --       It is still used to generate the other sieves.
    spfSieve :: Word -> ST s (STUArray s Int Word)
    spfSieve bound = do
      let (octs,lidx) = idxPr bound
          !mxidx = 8*octs+lidx
          mxval :: Word
          mxval = 30*fromIntegral octs + fromIntegral (rho lidx)
          !mxsve = integerSquareRoot' mxval
          (kr,r) = idxPr mxsve
          !svbd = 8*kr+r
      ar <- newArray (0,mxidx) 0
      let start k i = 8*(k*(30*k+2*rho i) + byte i) + idx i
          tick p stp off j ix
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
            | bound < n = tdLoop tt n (integerSquareRoot' n) 0
            | otherwise = case unsafeAt sve (toIdx n) of
                            nt -> tt `lcm` fromIntegral nt
        lstIdx = snd (bounds sve)
        tdLoop !tt n sr ix
            | lstIdx < ix   = curve tt n
            | sr < p'       = tt `lcm` (n-1)      -- n is a prime
            | pix /= p-1    = tdLoop tt n sr (ix+1)    -- not a prime, next
            | otherwise     = case splitOff p' n of
                                (0,_) -> tdLoop tt n sr (ix+1)
                                (k,m) -> let nt = (lcm tt (p'-1))*p'^(k-1)
                                         in case m of
                                              1 -> nt
                                              j | j <= bound -> nt `lcm` fromIntegral (unsafeAt sve (toIdx j))
                                                | otherwise  -> tdLoop nt j (integerSquareRoot' j) (ix+1)
              where
                p = toPrim ix
                p' = fromIntegral p
                pix = unsafeAt sve ix
        curve tt n = tt `lcm` carmichaelFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n)
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
                      (0,_) -> go5 tt n
                      (k,1) | tt == 1   -> 2*3^(k-1)
                            | otherwise -> tt*3^(k-1)
                      (k,m) | tt == 1   -> go5 (2*3^(k-1)) m
                            | otherwise -> go5 (tt*3^(k-1)) m
        go5 !tt n = case splitOff 5 n of
                      (0,_) -> sieveLoop tt n
                      (k,m) -> let tt' = case fromInteger tt .&. (3 :: Int) of
                                           0 -> tt
                                           2 -> 2*tt
                                           _ -> 4*tt
                                   nt = tt'*5^(k-1)
                               in if m == 1 then nt else sieveLoop nt m
        sieveLoop !tt n
            | bound < n = tdLoop tt n (integerSquareRoot' n) 0
            | otherwise = case unsafeAt sve (toIdx n) of
                            nt -> tt `lcm` fromIntegral nt
        lstIdx = snd (bounds sve)
        tdLoop !tt n sr ix
            | lstIdx < ix   = curve tt n
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
            | bound < n = tdLoop tt n (integerSquareRoot' n) 0
            | otherwise = case unsafeAt sve (toIdx n) of
                            nt -> tt*fromIntegral nt
        lstIdx = snd (bounds sve)
        tdLoop !tt n sr ix
            | lstIdx < ix   = curve tt n
            | sr < p'       = tt*(n-1)      -- n is a prime
            | pix /= p-1    = tdLoop tt n sr (ix+1)    -- not a prime, next
            | otherwise     = case splitOff p' n of
                                (0,_) -> tdLoop tt n sr (ix+1)
                                (k,m) -> let nt = tt*ppTotient (p',k)
                                         in case m of
                                              1 -> nt
                                              j | j <= bound -> nt*fromIntegral (unsafeAt sve (toIdx j))
                                                | otherwise  -> tdLoop nt j (integerSquareRoot' j) (ix+1)
              where
                p = toPrim ix
                p' = fromIntegral p
                pix = unsafeAt sve ix
        curve tt n = tt * totientFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n)
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
        check n
            | n < 1     = error "Totient only defined for positive numbers"
            | n == 1    = 1
            | otherwise = go2 n
        go2 n = case shiftToOddCount n of
                  (0,_) -> go3 1 n
                  (k,m) -> let tt = (shiftL 1 (k-1)) in if m == 1 then tt else go3 tt m
        go3 !tt n = case splitOff 3 n of
                      (0,_) -> go5 tt n
                      (k,m) -> let nt = tt*(2*3^(k-1)) in if m == 1 then nt else go5 nt m
        go5 !tt n = case splitOff 5 n of
                      (0,_) -> sieveLoop tt n
                      (k,m) -> let nt = tt*(4*5^(k-1)) in if m == 1 then nt else sieveLoop nt m
        sieveLoop !tt n
            | bound < n = tdLoop tt n (integerSquareRoot' n) 0
            | otherwise = case unsafeAt sve (toIdx n) of
                            nt -> tt*fromIntegral nt
        lstIdx = snd (bounds sve)
        tdLoop !tt n sr ix
            | lstIdx < ix   = curve tt n
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
                0 | p-3 == 2*i -> [(fromIntegral p,c+1)]
                  | otherwise  -> (fromIntegral p,c) : (2*fromIntegral i+3,1) : []
                q | fromIntegral q == p -> countLoop p (i `quot` p - 1) (c+1)
                  | otherwise -> (fromIntegral p, c) : intLoop i
        lstIdx = snd (bounds sve)
        tdLoop n sr ix
            | lstIdx < ix   = curve n
            | sr < p        = [(n,1)]
            | pix /= 0      = tdLoop n sr (ix+1)    -- not a prime
            | otherwise     = case splitOff p n of
                                (0,_) -> tdLoop n sr (ix+1)
                                (k,m) -> (p,k) : case m of
                                                   1 -> []
                                                   j | j <= bound -> intLoop (fromIntegral (j `shiftR` 1) - 1)
                                                     | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)
              where
                p = fromIntegral $ 2 * ix + 3
                pix = unsafeAt sve ix
        curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
        check 0 = error "0 has no prime factorisation"
        check 1 = []
        check n
          | n < 0       = (-1,1) : check (-n)
          | n <= bound  = go2w (fromIntegral n)     -- avoid expensive Integer ops if possible
          | fromInteger n .&. (1 :: Int) == 1 = sieveLoop n
          | otherwise   = go2 n
        go2w n
            | n .&. 1 == 1 = intLoop ((n-3) `shiftR` 1)
            | otherwise = case shiftToOddCount n of
                            (k,m) -> (2,k) : if m == 1 then [] else intLoop ((m-3) `shiftR` 1)
        go2 n = case shiftToOddCount n of
                  (k,m) -> (2,k) : if m == 1 then [] else sieveLoop m
        sieveLoop n
            | bound < n  = tdLoop n (integerSquareRoot' n) 0
            | otherwise = intLoop (fromIntegral (n `shiftR` 1)-1)
        intLoop :: Word -> [(Integer,Int)]
        intLoop !n = case unsafeAt sve (fromIntegral n) of
                      0 -> [(2*fromIntegral n+3,1)]
                      p -> let p' = fromIntegral p in countLoop p' (n `quot` p' - 1) 1
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
    -- | @'factorSieve' n@ creates a store of smallest prime factors of the numbers not exceeding @n@.
    --   If you need to factorise many smallish numbers, this can give a big speedup since it avoids
    --   many superfluous divisions. However, a too large sieve leads to a slowdown due to cache misses.
    --   The prime factors are stored as 'Word16' for compactness, so @n@ must be
    --   smaller than @2^32@.
    factorSieve :: Integer -> FactorSieve
    factorSieve bound
      | 4294967295 < bound  = error "factorSieve: overflow"
      | bound < 8   = FS 7 (array (0,2) [(0,0),(1,0),(2,0)])
      | otherwise   = FS bnd fSieve
        where
          bnd = fromInteger bound
          ibnd = fromInteger ((bound - 3) `quot` 2)
          svbd = (fromInteger (integerSquareRoot' bound) - 1) `quot` 2
          fSieve = runSTUArray $ do
              sieve <- newArray (0,ibnd) 0 :: ST s (STUArray s Int Word16)
              let sift i
                      | i < svbd = do
                          sp <- unsafeRead sieve i
    

    Found in Math.NumberTheory.Primes.Sieve.Misc from the package arithmoi
          -- ** Carmichael
        , carmichaelSieve
        , sieveCarmichael
        ) where
    
    import Control.Monad.ST
    import Data.Array.ST
    import Data.Array.Unboxed
    import Control.Monad (when)
    import Data.Bits
    import GHC.Word
    
    import System.Random
    
    import Math.NumberTheory.Powers.Squares (integerSquareRoot')
    import Math.NumberTheory.Primes.Sieve.Indexing
    import Math.NumberTheory.Primes.Factorisation.Montgomery
    import Math.NumberTheory.Unsafe
    import Math.NumberTheory.Utils
    

    Found in Math.NumberTheory.Primes.Testing.Certificates.Internal from the package arithmoi
                                                            ++ "Please report this to the maintainer.\n\n"
                                                            ++ show cpr)
                                Prime ppr ->(p,e,bs,ppr)
                  where
                    q = nm1 `quot` p
                    x = powModInteger bs q n
                    y = powModInteger x p n
                    g = gcd n (x-1)
    
    -- | Find a decomposition of p-1 for the pocklington certificate.
    --   Usually bloody slow if p-1 has two (or more) /large/ prime divisors.
    findDecomposition :: Integer -> (Integer, [(Integer,Int,Bool)], Integer)
    findDecomposition n = go 1 n [] prms
      where
        sr = integerSquareRoot' n
        pbd = min 1000000 (sr+20)
        prms = primeList (primeSieve $ pbd)
        go a b afs (p:ps)
            | a > b     = (a,afs,b)
            | otherwise = case splitOff p b of
    

    Found in Math.NumberTheory.Primes.Testing.Certificates.Internal from the package arithmoi
    smallCert n
        | n < 30    = Trivial n
        | otherwise = TrialDivision n (integerSquareRoot' n + 1)
    
    -- | @'certify' n@ constructs, for @n > 1@, a proof of either
    --   primality or compositeness of @n@. This may take a very long
    --   time if the number has no small(ish) prime divisors
    certify :: Integer -> Certificate
    certify n
        | n < 2     = error "Only numbers larger than 1 can be certified"
        | n < 31    = case trialDivisionWith trivialPrimes n of
                        ((p,_):_) | p < n     -> Composite (Factors n p (n `quot` p))
                                  | otherwise -> Prime (Trivial n)
                        _ -> error "Impossible"
        | n < billi = let r2 = integerSquareRoot' n + 2 in
                      case trialDivisionTo r2 n of
                        ((p,_):_) | p < n       -> Composite (Factors n p (n `quot` p))
                                  | otherwise   -> Prime (TrialDivision n r2)
                        _ -> error "Impossible"
        | otherwise = case smallFactors 100000 n of
    

    Found in Math.NumberTheory.Primes.Testing.Certificates.Internal from the package arithmoi
    --   known prime.
    isTrivialPrime :: Integer -> Bool
    isTrivialPrime n = n `elem` trivialPrimes
    
    -- | List of trivially known primes.
    trivialPrimes :: [Integer]
    trivialPrimes = [2,3,5,7,11,13,17,19,23,29]
    
    -- | Certify 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.
    smallCert :: Integer -> PrimalityProof
    smallCert n
        | n < 30    = Trivial n
        | otherwise = TrialDivision n (integerSquareRoot' n + 1)
    
    -- | @'certify' n@ constructs, for @n > 1@, a proof of either
    --   primality or compositeness of @n@. This may take a very long
    --   time if the number has no small(ish) prime divisors
    certify :: Integer -> Certificate
    

    integerSquareRootRem No usage example found for this symbol :( Collapse [-]
    integerSquareRootRem' No usage example found for this symbol :( Collapse [-]
    isPossibleSquare No usage example found for this symbol :( Collapse [-]
    isPossibleSquare2 No usage example found for this symbol :( Collapse [-]
    isSquare See 1 Occurences [+] Collapse [-]
    Found in Math.NumberTheory.Powers from the package arithmoi
    -- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
    -- Stability:   Provisional
    -- Portability: Non-portable (GHC extensions)
    --
    -- Calculating integer roots, modular powers and related things.
    -- This module reexports the most needed functions from the implementation
    -- modules. The implementation modules provide some additional functions,
    -- in particular some unsafe functions which omit some tests for performance
    -- reasons.
    --
    module Math.NumberTheory.Powers
      ( -- *  Integer Roots
        -- ** Square roots
        integerSquareRoot
      , isSquare
      , exactSquareRoot
        -- ** Cube roots
      , integerCubeRoot
      , isCube
      , exactCubeRoot
    

    isSquare' No usage example found for this symbol :( Collapse [-]