-- | -- Module: Math.NumberTheory.Canon.Internals -- Copyright: (c) 2015-2019 Frederick Schneider -- Licence: MIT -- Maintainer: Frederick Schneider <fws.nyc@gmail.com> -- Stability: Provisional -- -- This module defines the internal canonical representation of numbers (CR_), a product of pairs (prime and exponent). -- PLEASE NOTE: It's not meant to be called directly. The Canon and SimpleCanon modules should be used instead. {-# LANGUAGE PatternSynonyms, ViewPatterns, ScopedTypeVariables, DataKinds, RankNTypes #-} module Math.NumberTheory.Canon.Internals ( CanonRep_, CR_, CanonElement_, crValidIntegral, crValidIntegralViaUserFunc, crValidRational, crValidRationalViaUserFunc, crExp, crRoot, crMaxRoot, crShow, ceShow, crFromInteger, crFromI, crToInteger, crToI, crCmp, crMult, crNegate, crAbs, crDivStrict, crSignum, crNumer, crDenom, crSplit, crDivRational, crIntegral, crShowRational, crToRational, crGCD, crLCM, crNegative, crPositive, crLog, crLogDouble, crDiv, crRadical, integerApply, crSimpleApply, crPrime, crHasSquare, crRecip, crMin, crMax, crValid, crMod, crModI, crSimplified, crDivisors, crNumDivisors, crWhichDivisor, crNthDivisor, crDivsPlus, crTau, crTotient, crPhi, crN1, cr0, cr1, creN1, crFactCutoff, crSmallFactCutoff, pattern PZero, pattern PZeroBad, pattern POne, pattern PNeg, pattern PNotPos, pattern PN1, pattern PIntNeg, pattern PIntNPos, pattern PIntPos, -- Functions deprecated from arithmoi that needed to be included here totient, pmI -- Short for powerModInteger ) where {- Canon or canon rep is short for canonical representation. In this context, it refers to the prime factorization of a number. Example: 250 = 2 * 5^3. It would be represented internally here as: [(2,1),(5,3)] So, this library along with Canon.hs can be used as shorthand for extremely large numbers. Multiplicative actions are cheap. But, addition and subtraction is potentially very expensive as any "raw" sums or differences are factorized. There are optimizations for sums/differences including those with special forms (algebraic, Aurifeuillean). Here are the possibilities: Zero: [(0,1)] One [] Other Positive Numbers: A set of (p,e) pairs where the p(rime)s are in ascending order. For "integers", the e(xponents) must be positive For "rationals", the e(xponents) must not be zero All "integers" are "rationals" but not vice-versa. Negative Numbers: (-1,1):P where P is a canon rep for any positive number Note: Much of the code is spent on handling these special cases (Zero, One, Negatives). Each integer and rational will have a unique "canon rep". Caveats: The behavior is undefined when directly working with "canon reps" which are not valid. The behavior of using rational CRs directly (when integral CRs are specified) is also undefined. The Canon library should be used as it hides these internal details. -} import Data.List (intersperse, partition) import Math.NumberTheory.Primes.Factorisation (smallFactors, stdGenFactorisation, factorise') import System.Random (mkStdGen) import Data.Bits (xor) import Data.List (sortBy) import Math.NumberTheory.Primes.Testing (isPrime) import GHC.Real (Ratio(..)) -- | Canon element: prime and exponent pair type CanonElement_ = (Integer,Integer) -- | Canonical representation: list of canon elements type CanonRep_ = [CanonElement_] -- | Shorthand for canonical representation type CR_ = CanonRep_ -- | Pattern to match the CR_ equivalent of 1 pattern POne :: forall t. [t] pattern POne = [] -- | Pattern to match the CR_ equivalent of zero pattern PZero :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)] pattern PZero = [(0,1)] -- | Pattern to match the CR_ equivalent of -1 pattern PN1 :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)] pattern PN1 = [(-1,1)] -- | Pattern to match a badly formed zero, meaning it's an invalid CR_ pattern PZeroBad :: forall t a. (Num a, Eq a) => [(a, t)] pattern PZeroBad <- ((0,_):_) -- MUST check after PZero -- No longer necessary -- pattern PNotIntegral :: forall a t. (Num a, Ord a) => [(t, a)] -- pattern PNotIntegral <- ( (_, compare 0 -> GT):_ ) -- negative exponent in the 2nd member of pair -- | Pattern to match a non-positive CR_ pattern PNotPos :: forall t a. (Num a, Ord a) => [(a, t)] pattern PNotPos <- ( (compare 1 -> GT, _):_ ) -- first term is 0, -1 and so not positive -- | Pattern to match a negative number pattern PIntNeg :: forall a. (Num a, Ord a) => a pattern PIntNeg <- (compare 0 -> GT) -- | Pattern to match a positive number pattern PIntPos :: forall a. (Num a, Ord a) => a pattern PIntPos <- (compare 0 -> LT) -- | Pattern to match a non-positive number pattern PIntNPos :: forall a. (Num a, Ord a) => a pattern PIntNPos <- (compare 1 -> GT) -- | Canonical values for a few special numbers creN1, cre0 :: CanonElement_ creN1 = (-1,1) cre0 = (0,1) crN1, cr0, cr1 :: CanonRep_ -- | Canon rep for -1 crN1 = [creN1] -- | Canon rep for 0 cr0 = [cre0] -- | Canon rep for 1 cr1 = [] -- Yes, a canonical "1" is just an empty list. -- | Pattern for a negative CR_ pattern PNeg :: forall a a1. (Num a, Num a1, Eq a, Eq a1) => [(a1, a)] pattern PNeg <- ((-1,1):_) crNegative, crPositive :: CR_ -> Bool -- | Check if a CR_ is negative. crNegative PNeg = True crNegative _ = False -- | Check if a CR_ is positive. crPositive PZero = False crPositive x = not $ crNegative x -- | Canon rep validity check: -- The 2nd param checks the validity of the base, the 3rd of the exponent. -- The base pred should be some kind of prime or psuedo-prime test unless you knew for -- certain the bases are prime. There are two choices for the exp pred: -- positiveOnly (True) or nonZero (False) (which allows for "rationals"). crValid :: CR_ -> (Integer -> Bool) -> Bool -> Bool crValid POne _ _ = True crValid PZero _ _ = True crValid PZeroBad _ _ = False crValid c bp ef | crNegative c = f (tail c) 1 | otherwise = f c 1 where f POne _ = True f ((cp,ce):cs) n | cp <= n || not (expPred ef ce) || not (bp cp) = False | otherwise = f cs cp f _ _ = error "Logic error in crValid'. Issue with pattern matching?" expPred b e = if b then (e > 0) else (e /= 0) crValidIntegral, crValidRational :: CR_ -> Bool crValidIntegralViaUserFunc, crValidRationalViaUserFunc :: CR_ -> (Integer -> Bool) -> Bool -- | Checks if a CR_ represents an integral number. crValidIntegral n = crValid n isPrime True -- | Checks if a CR_ is Integral and valid per user-supplied criterion. crValidIntegralViaUserFunc n f = crValid n f True -- | Checks if a CR_ is represents a rational number (inclusive of integral numbers). crValidRational n = crValid n isPrime False -- | Checks if a CR_ is Rational and valid per user-supplied criterion. crValidRationalViaUserFunc n f = crValid n f False crFromInteger, crFromI :: Integer -> (CR_, Bool) -- | Factor the number to convert it to a canonical rep. This is of course can be extremely expensive. -- (This now also returns a flag indicating True if the number was completely factored.) crFromInteger n | n == 0 = (cr0, True) | n == 1 = (cr1, True) | n == -1 = (crN1, True) | n < 0 = (creN1:cr, ff) -- in case arithmoi excludes -1 as a factor in the future | otherwise = (cr, ff) where cr = crSort cr' -- the prime factors must be in ascending order (cr', ff) = factorize (abs n) -- ToDo: pass to fcn crSort :: CR_ -> CR_ crSort = sortBy prd where prd (p1, _) (p2, _) | p1 < p2 = LT | otherwise = GT -- Note: if crFactCutoff is <= 0, complete factorization is attempted -- and all of the cutoff / spFactor logic is not used. crFactCutoff, crTrialDivCutoff, crSmallFactCutoff, crTrialDivCutoffSq :: Integer crFactCutoff = (10 :: Integer) ^ (80 :: Int) -- Note: if this is <= 0, complete factorization is attempted crTrialDivCutoff = 100000 crSmallFactCutoff = 10000000 -- use this higher cutoff if the number is beyond the factorization cutoff crTrialDivCutoffSq = crTrialDivCutoff * crTrialDivCutoff -- factorize and deftStdGenFact were adapted from arithmoi factorize :: Integer -> ([(Integer,Integer)], Bool) factorize n = (map (\(p,e) -> (p, toInteger e)) f, b) where (f, b) = if crFactCutoff > 0 then defStdGenFact (mkStdGen $ fromInteger n `xor` 0xdeadbeef) else (factorise' n, True) -- call underlying routine from arithmoi defStdGenFact sg = let (sfs,mb) = smallFactors (if n <= crFactCutoff then crTrialDivCutoff else crSmallFactCutoff) n in (sfs ++ case mb of Nothing -> [] Just m -> if m > crFactCutoff then [(m, 1)] else stdGenFactorisation (Just crTrialDivCutoffSq) sg Nothing m, case mb of Nothing -> True Just m -> m < crFactCutoff || isPrime m -- if less, do a complete factorization ) -- | Shorthand for crFromInteger function crFromI n = crFromInteger n crToInteger, crToI :: CR_ -> Integer -- | Convert a canon rep back to an Integer. crToInteger POne = 1 crToInteger PZero = 0 crToInteger c | (head c) == creN1 = -1 * (crToInteger $ tail c) -- negative number | otherwise = product $ map (\(x,y) -> x ^ y) c -- | Alias to crToInteger. crToI = crToInteger -- | Compute the modulus between a CR_ and Integer and return an Integer. crModI :: CR_ -> Integer -> Integer crModI _ 0 = error "Divide by zero error when computing n mod 0" crModI _ 1 = 0 crModI _ (-1) = 0 crModI POne PIntPos = 1 crModI PZero _ = 0 crModI c m | cn && mn = -1 * crModI (crAbs c) am | (cn && not mn) || (mn && not cn) = (signum m) * (am - f (crAbs c) am) | otherwise = f c m where cn = crNegative c mn = m < 0 am = abs m f c' m' = mod (product $ map (\(x,y) -> if (mod x m' == 0) then 0 else (pmI x (mmt y) m')) c') m' mmt e | e >= 1 = mod e $ totient m -- optimization | otherwise = error "Negative exponents are not allowed in crModI" -- | Compute modulus with all CR_ parameters. This wraps crModI. crMod :: CR_ -> CR_ -> (CR_, Bool) crMod c m = crFromI $ crModI c (crToI m) -- | Display a Canon Element (as either p^e or p). ceShow :: CanonElement_ -> String ceShow (p,e) = show p ++ if e == 1 then "" else "^" ++ (if e < 0 then "(" ++ se ++ ")" else se) where se = show e crShow, crShowRational :: CR_ -> String -- | Display a canonical representation. crShow POne = show (1 :: Integer) crShow x | null (tail x) = ceShow $ head x | otherwise = concat $ intersperse " * " $ map ceShow x -- | Display a CR_ rationally, as a quotient of its numerator and denominator. crShowRational c | d == cr1 = crShow n | otherwise = crShow n ++ "\n/\n" ++ crShow d where (n, d) = crSplit c crNegate, crAbs, crSignum :: CR_ -> CR_ -- | Negate a CR_. crNegate PZero = cr0 crNegate x | crNegative x = tail x | otherwise = creN1 : x -- | Take the absolute value of a CR_. crAbs x | crNegative x = tail x | otherwise = x -- | Compute the signum and return as CR_. crSignum PZero = cr0; crSignum x | crNegative x = crN1 | otherwise = cr1 -- | CR_ Compare Function crCmp :: CR_ -> CR_ -> Ordering crCmp POne y = crCmp1 y True crCmp x POne = crCmp1 x False crCmp x y | x == y = EQ | xN && yN = crCmp (tail y) (tail x) | xN = LT | yN = GT | eqZero x = LT | eqZero y = GT | otherwise = case compare (crLogDouble x) (crLogDouble y) of EQ -> compare (crLog x) (crLog y) -- We have to break out the big guns, both evaluate to infinity cmp -> cmp where xN = crNegative x yN = crNegative y eqZero PZero = True; eqZero _ = False -- Internal: Compare when either term is 1. crCmp1 :: CR_ -> Bool -> Ordering crCmp1 PNeg b = if b then GT else LT crCmp1 PZero b = if b then GT else LT crCmp1 _ b = if b then LT else GT crMin, crMax :: CR_ -> CR_ -> CR_ -- | Min function crMin x y = case crCmp x y of GT -> y _ -> x -- | Max functon crMax x y = case crCmp x y of LT -> y _ -> x divisionError, divByZeroError, zeroDivZeroError, negativeLogError :: String divisionError = "For this function, the dividend must be a multiple of the divisor." divByZeroError = "Division by zero error!" zeroDivZeroError = "Zero divided by zero is undefined!" negativeLogError = "The log of a negative number is undefined!" -- | Strict division: Generate error if exact division is not possible. crDivStrict :: CR_ -> CR_ -> CR_ crDivStrict x y = case crDiv x y of Left errorMsg -> error errorMsg Right quotient -> quotient -- | Attempt to take the quotient. crDiv :: CR_ -> CR_ -> Either String CR_ crDiv PZero PZero = Left zeroDivZeroError crDiv PZero _ = Right cr0 crDiv _ PZero = Left divByZeroError crDiv x' y' = f x' y' where -- call this after handling zeroes above, then division just occurs within here f x POne = Right x f POne _ = Left divisionError f x y | crNegative y = f (crNegate x) (crAbs y) | otherwise = case compare xp yp of LT -> case f xs y of Left _ -> Left divisionError Right r -> Right ((xp, xe):r) EQ| (xe > ye) -> case f xs ys of Left _ -> Left divisionError Right r -> Right ((xp,xe-ye):r) EQ| (xe == ye) -> f xs ys _ -> Left divisionError where ((xp,xe):xs) = x ((yp,ye):ys) = y crMult, crDivRational, crGCD, crLCM :: CR_ -> CR_ -> CR_ -- | Multiply two crs by summing the exponents for each prime. crMult PZero _ = cr0 crMult _ PZero = cr0 crMult POne y = y crMult x POne = x crMult x y = case compare xp yp of LT -> xh : crMult xs y -- cancel double negs or expnts adding to zero EQ -> if crNegative x || expSum == 0 then r else (xp, expSum) : r where r = crMult xs ys GT -> yh : crMult x ys where (xh@(xp,xe):xs) = x (yh@(yp,ye):ys) = y expSum = xe + ye -- | Division of rationals is equivalent to multiplying with negated exponents. crDivRational x y = crMult (crRecip y) x -- multiply by the reciprocal -- | For the GCD (Greatest Common Divisor), take the lesser of two exponents for each prime encountered. crGCD PZero y = y crGCD x PZero = x crGCD x y | crNegative x || crNegative y = f (crAbs x) (crAbs y) | crFactCutoff > 0 && ((spIncompFact x && spUnsmooth y) || -- either has an imcomplete factorization and the other (spIncompFact y && spUnsmooth x)) = f spx spy -- in case of unfactored bits | otherwise = f x y where f POne _ = cr1 f _ POne = cr1 f x' y' = case compare xp yp of LT -> f xs y' EQ -> (xp, min xe ye) : f xs ys GT -> f x' ys where ((xp,xe):xs) = x' ((yp,ye):ys) = y' (spx, spy) = spFactor x y -- | For the LCM (Least Common Multiple), take the max of two exponents for each prime encountered. crLCM PZero _ = cr0 crLCM _ PZero = cr0 crLCM x y | crNegative x || crNegative y = f (crAbs x) (crAbs y) | crFactCutoff > 0 && ((spIncompFact x && spUnsmooth y) || -- either has an imcomplete factorization and the other (spIncompFact y && spUnsmooth x)) = f spx spy -- in case of unfactored bits | otherwise = f x y where f POne y' = y' f x' POne = x' f x' y' = case compare xp yp of LT -> xh : f xs y' EQ -> (xp, max xe ye) : f xs ys GT -> yh : f x' ys where (xh@(xp,xe):xs) = x' (yh@(yp,ye):ys) = y' (spx, spy) = spFactor x y -- special factor: This covers the case where we have unfactored large components but on comparison with another -- cr we can see that the component can be reduced. We partition the cr into -- three pieces: small factor cutoff, prime, composite (implying it's > factor. cutoff) spFactor :: CR_ -> CR_ -> (CR_, CR_) spFactor x y = (sx ++ (grp $ crSort $ px ++ spF cx (py ++ cy)), sy ++ (grp $ crSort $ py ++ spF cy (px ++ cx)) ) where spl n = (s, p, c) where (s, r) = partition spSmooth n (p, c) = partition (\ce -> not $ spBigComposite ce) r (sx, px, cx) = spl x (sy, py, cy) = spl y grp (n:ns) = g' n ns grp _ = error "The list to be grouped in spFactor must have at least one element" g' (b,e) (r:rs) = if b == b' then g' (b, e + e') rs -- group by common base on sorted list else (b,e):g' (b', e') rs where (b', e') = r g' ce _ = [ce] -- take each entry in f and compute the gcd. ToDo: replace with fold spF n (f:fs) = spF (concat $ map (proc f) n) fs -- this is quadratic but likely with very short lists where proc (pf, _) (pn, en) = if g == 1 || g == pn then [(pn, en)] else [(g, en), (div pn g, en)] where g = gcd pn pf spF n _ = n -- Predicates used for special cases of GCD and LCM spUnsmooth, spIncompFact :: CR_ -> Bool spUnsmooth = any (\ce -> not $ spSmooth ce) spIncompFact = any spBigComposite spBigComposite, spSmooth :: CanonElement_ -> Bool spSmooth (b,_) = b <= crSmallFactCutoff spBigComposite (b,_) = b > crFactCutoff && (not $ isPrime b) -- | Take the reciprocal by raising a CR to the -1 power (equivalent to multiplying exponents by -1). crRecip :: CR_ -> CR_ crRecip x = crExp x (-1) True rootError :: CR_ -> Integer -> String rootError c r = "crRoot: All exponents must be multiples of " ++ (show r) ++ ". Not so with " ++ (show c) -- | Attempt to compute a particular root of a CR_. crRoot :: CR_ -> Integer -> CR_ crRoot _ PIntNeg = error "r-th roots are not allowed when r <= 0" crRoot POne _ = cr1 crRoot c r | crNegative c = if mod r 2 == 1 then creN1 : crRoot (crAbs c) r else error "Imaginary numbers not allowed: Even root of negative number requested" | otherwise = if mod ce r == 0 then (cp, div ce r) : crRoot cs r else error $ rootError c r where ((cp,ce):cs) = c -- | Takes the maximum root of the number. Generally, the abs value would be passed to the function. crMaxRoot :: CR_ -> Integer crMaxRoot c = foldr (\x -> flip gcd $ snd x) 0 c -- | Exponentiation. Note: this does allow for negative exponentiation if bool flag is True. crExp :: CR_ -> Integer -> Bool -> CR_ crExp _ PIntNeg False = error "Per param flag, negative exponentiation is not allowed here." crExp PZero PIntNPos _ = error "0^e where e <= 0 is either undefined or illegal" crExp PZero _ _ = cr0 crExp POne _ _ = cr1 crExp _ 0 _ = cr1 crExp c em _ = ce c where ce c' | crNegative c' = if mod em 2 == 1 then creN1 : absTail else absTail | otherwise = map (\(p,e) -> (p, e * em)) c' where absTail = ce $ crAbs c' -- | This log function is much more expensive but accurate. You have an "infinity" problem potentially with crLogDouble. crLog :: CR_ -> Rational crLog PNeg = error negativeLogError crLog c = sum $ map (\(p,e) -> (toRational $ logDouble $ fromIntegral p) * (fromIntegral e)) c where logDouble :: Double -> Double logDouble n = log n -- | Returns log of CR_ as a Double. crLogDouble :: CR_ -> Double crLogDouble PNeg = error negativeLogError crLogDouble c = sum $ map (\(x,y) -> log (fromIntegral x) * fromIntegral y) c crNumer, crDenom, crRadical :: CR_ -> CR_ -- | Compute numerator (by filtering on positive exponents). crNumer c = filter (\(_,e) -> e > 0) c -- | Compute denominator. (Grab the primes with negative exponents and then flip the sign of the exponents.) crDenom c = map (\(p,e) -> (p, (-1) * e)) $ filter (\(_,e) -> e < 0) c -- | Check if a CR_ represents an integer. crIntegral :: CR_ -> Bool crIntegral x = all (\(_,e) -> e > 0) x -- all exponents must be positive -- | Split a CR_ into its Numerator and Denominator. crSplit :: CR_ -> (CR_, CR_) crSplit c = (crNumer c, crDenom c) -- | Convert a CR_ to a Rational number. crToRational :: CR_ -> Rational crToRational c = (crToI $ crNumer c) :% (crToI $ crDenom c) -- | Compute the Radical of a CR_ (http://en.wikipedia.org/wiki/Radical_of_an_integer). -- Its the product of the unique primes in its factorization. crRadical n = map (\(p,_) -> (p, 1)) n -- | The Op(eration) is intended to be plus or minus. integerApply :: (Integer -> Integer -> Integer) -> CR_ -> CR_ -> Integer integerApply op x y = op (crToI x) (crToI y) -- | Calls integerApply and returns a CR_. crSimpleApply :: (Integer -> Integer -> Integer) -> CR_ -> CR_ -> (CR_, Bool) crSimpleApply op x y = crFromI $ integerApply op x y {- No longer needed. Different criteria used now pattern PPrime :: forall a a1. (Eq a, Num a, Num a1, Ord a1) => [(a1, a)] pattern PPrime <- [(compare 1 -> LT, 1)] -- of form x^1 where x > 1 -- prime (assumption p has already been verified to be prime) -} crPrime, crHasSquare :: CR_ -> Bool -- | Check if a number is a prime. crPrime cr | length cr > 1 || null cr = False | e == 1 && b > 1 && (crFactCutoff == 0 || b <= crFactCutoff || (b > crFactCutoff && isPrime b)) = True | otherwise = False where (b,e) = head cr -- | Check if a number has a squared (or higher) factor. crHasSquare = any (\(_,e) -> e > 1) -- | Divisor functions should be called with integral CRs (no negative exponents). crNumDivisors, crTau, crTotient, crPhi :: CR_ -> Integer crNumDivisors cr = product $ map (\(_,e) -> 1 + e) cr -- does return 1 for cr1 crTau = crNumDivisors crTotient cr = product $ map (\(p,e) -> (p-1) * p^(e-1)) cr crPhi = crTotient -- | Compute the nth divisor. This is zero based. -- Note: This is deterministic but it's not ordered by the value of the divisor. crNthDivisor :: Integer -> CR_ -> CR_ crNthDivisor 0 _ = cr1 crNthDivisor _ POne = error "Bad div num requested" crNthDivisor n c | m == 0 = r | otherwise = (cb,m) : r where (cb,ce):cs = c r = crNthDivisor (div n (ce + 1)) cs -- first param is the next n m = mod n (ce + 1) -- | Consider this to be the inverse of the crNthDivisor function. crWhichDivisor :: CR_ -> CR_ -> Integer crWhichDivisor d c | crPositive d == False || crPositive c == False = error "crWhichDivisor: Both params must be positive" | otherwise = f d c where f POne _ = 0 f _ POne = error "Not a valid divisor" f d' c' | dp < cp || (dp == cp && de > ce) = error "Not a valid divisor" | dp == cp = de + (ce + 1) * (f ds cs) | otherwise = (ce + 1) * (f d cs) where ((dp, de):ds) = d' ((cp, ce):cs) = c' -- | Efficiently compute all of the divisors based on the canonical representation. crDivisors :: CR_ -> [CR_] crDivisors c = foldr1 cartProd $ map pwrDivList c where cartProd xs ys = [x ++ y | y <- ys, x <- xs] pwrDivList (n,e) = [if y == 0 then cr1 else [(n,y)] | y <- [0..(fromInteger e)]] -- | Like the crDivisors function, except that it returns pairs of the CR_ and resp. numeric value, instead of just the CR_. crDivsPlus :: CR_ -> [(CR_, Integer)] crDivsPlus c = foldr1 cartProd (map pwrDivList c) where cartProd xs ys = [(xl ++ yl, xv * yv) | (yl, yv) <- ys, (xl, xv) <- xs] pwrDivList (e,n) = map tr $ pwrList e n powers x = 1 : map (* x) (powers x) pwrList n e = [(n,y) | y <- zip [0..e] (take (e'+1) $ powers n)] where e' = fromInteger e tr (a,b) = (if fb == 0 then cr1 else [(a, fb)], sb) -- this just transforms the data structure where (fb, sb) = b -- | Check if underlying rep is simplified crSimplified :: CR_ -> Bool crSimplified POne = True crSimplified PZero = True crSimplified PN1 = True crSimplified c = crPrime c -- | Compute totient: Logic from deprecated arithmoi function used here. totient :: Integer -> Integer totient n | n < 1 = error "Totient only defined for positive numbers" | n == 1 = 1 | ff == False = error "Totient not computable. The number could not be totally factored based on the factorizaton cutoff." | otherwise = product $ map (\(p,e) -> (p-1) * p ^ (e-1)) cr where (cr, ff) = factorize n -- ToDo: pass the param to function -- | powerModInteger adapted here from deprecated arithmoi function. pmI :: Integer -> Integer -> Integer -> Integer pmI x p m | x < 1 || p < 0 || m < 1 = error "pmI (powerModInteger) requires: x >= 1 &&, p >= 0, m >= 1" | otherwise = f p 1 (mod x m) -- last is the running exp of mod initially where f q w e | w == 0 || q == 0 = w | q == 1 = mod (w*e) m | otherwise = f (div q 2) nw (mod (e*e) m) where nw | mod q 2 == 1 = mod (w*e) m | otherwise = w