-- | -- Module: Math.NumberTheory.Canon.AurifCyclo -- Copyright: (c) 2015-2019 Frederick Schneider -- Licence: MIT -- Maintainer: Frederick Schneider -- Stability: Provisional -- -- Aurifeullian and Cyclotomic factorization method functions. {-# LANGUAGE PatternSynonyms, ViewPatterns #-} module Math.NumberTheory.Canon.AurifCyclo ( aurCandDec, aurDec, applyCycloPair, applyCycloPairWithMap, cyclo, cycloWithMap, cycloDivSet, cycloDivSetWithMap, chineseAurif, chineseAurifWithMap, crCycloAurifApply, applyCrCycloPair, divvy, CycloMap, getIntegerBasedCycloMap, showCyclo, crCycloInitMap, multPoly, divPoly, addPoly ) where import Math.NumberTheory.Canon.Internals import Math.NumberTheory.Moduli.Jacobi (JacobiSymbol(..), jacobi) import Data.Array (array, (!), Array(), elems) -- to do: convert to unboxed? https://wiki.haskell.org/Arrays import GHC.Real (numerator, denominator) import Data.List (sort, sortBy, (\\)) import qualified Data.Map as M -- CR_ Rep of 2 cr2 :: CR_ cr2 = fst $ crFromI 2 -- | This function checks if the inputs along with operator flag have a cyclotomic or Aurifeuillian form to greatly simplify factoring. -- If they do not, potentially much more expesive simple factorization is used via crSimpleApply. -- Note: The cyclotomic map is threaded into the functions crCycloAurifApply :: Bool -> CR_ -> CR_ -> CR_ -> Integer -> CycloMap -> (CR_, CycloMap) crCycloAurifApply b x y g gi m -- Optimization for prime g: If g is a prime (and exp not of from x^2 + y^2) but not Aurifeullian (verify | (crPrime g) && not (g == cr2 && b) = eA ([term1, termNM1], m) -- split into (x +/- y) and (x^(n-1) ... -/+ y^(n-1)) -- Factorize: grtx^g - grty^g via cyclotomic polynomials | not b = eA (cycA grtx grty g) -- Factorize x^n + y^n using cyclotomic polynomials (if n = 2^x*m where x >= 0 and m > 2) | b && not gpwr2 = eA (cycA (oddRoot x) (-1 * oddRoot y) odd') | otherwise = (fst $ crSimpleApply op x y, m) where op = if b then (+) else (-) ((gp, _):gs) = g gpwr2 = gp == 2 && gs == [] gth_root v = crToI $ crRoot v gi grtx = gth_root x grty = gth_root y -- used when factoring x^p +/- 1 where p is prime term1 = integerApply op (crRoot x gi) (crRoot y gi) -- a +/- b termNM1 = div (integerApply op x y) term1 -- divide a^g +/- b^g by the term above cycA x' y' n = (sort ia, m') -- sort the integers returned from low to high, should help if there are larger terms where (ia, m') = applyCrCycloPair x' y' n m -- eA stands for "enriched apply", v unlike to exceed factor cutoff eA (a,mp) = (foldr1 crMult $ map (\v' -> fst $ crFromI v') v, m') where (v, m') = case aurCandDecI x y gi g b of Nothing -> auL a mp -- can't do anything Brent Aurif-wise, try Chinese method Just (a1, a2) -> auL (divvy a a1 a2) mp -- meld in the 2 Aurif factors with input array auL al ma = case c of -- aL stands for "augmented list) Nothing -> (al, mp') -- just return what was passed in Just (a3, a4) -> (divvy al a3 a4, mp') -- additional "Chinese" factors where (c, mp') = chineseAurifCr x y b ma odd' | gp == 2 = tail g -- grabs number sans any power of 2 | otherwise = g oddRoot v = crToI $ crRoot v (crToI odd') {- The following functions implement Richard Brent's algorithm for computing Aurifeullian factors. His logic is used in conjuction with cyclotomic polynomial decomposition. http://maths-people.anu.edu.au/~brent/pd/rpb127.pdf http://maths-people.anu.edu.au/~brent/pd/rpb135.pdf -} -- | This function checks if the input is a candidate for Aurifeuillian decomposition. -- If so, split it into two and evaluate it. Otherwise, return nothing. -- The code will "prep" the input params for the internal function so they will be relatively prime. -- Possible solutions: Non-zero multiples of xi and yi (where xi and yi are relatively prime and of the form -- xi = a^(a*b), yi = 1 where a and b are positive integers and b is odd. (xi and y1 may be interchanged) -- If a is even, we will find factors of a^(a*b) + 1. The bool flag must be True -- If a is odd, we will find factors of a^(a*b) - 1. The bool flag must be False -- aurCandDec :: Integer -> Integer -> Bool -> Maybe (Integer, Integer) aurCandDec xi yi b = f (fst $ crFromI xi) (fst $ crFromI yi) where f xp yp = aurCandDecI x y n (fst $ crFromI n) b where n = gcd (crMaxRoot $ crAbs x) (crMaxRoot $ crAbs y) gxy = crGCD xp yp (x, y) = (crDivStrict xp gxy, crDivStrict yp gxy) -- this will fix the input to be relatively prime -- toDo: incorporate the factorization status when determining the correct gcd. -- Don't call this I(nternal) function directly. The function assumes that x and y are relatively prime. Currently uses Brent logic only. aurCandDecI :: CR_ -> CR_ -> Integer -> CR_ -> Bool -> Maybe (Integer, Integer) aurCandDecI x y n cr b| (nm4 == 1 && b) || (nm4 /= 1 && not b) || (xdg == x && ydg == y) || (m /= 0) = Nothing -- | otherwise = case aurDecI n' cr' of Nothing -> Nothing Just (gamma, delta) -> apply gamma delta where -- override of n, to attempt decomp for g = gcd when number of form: g^gd +/-1, -- this will only work when either x or y = 1 and not for any other divisor of g. -- If both terms are not 1, we just attempt an Aurif. decomp for n -- TODO: Need to integrate chineseAurif as it does something different (n', cr') | x /= cr1 && y /= cr1 = (n, cr) | otherwise = (gcd1i, gcd1) where x1 = if y /= cr1 then y else x gcd1 = crRadical $ crGCD x1 cr gcd1i = crToI gcd1 nm4 = mod n' 4 divTry a = case crDiv a (crExp cr' n' False) of -- check to divide by n^n, if not return original Left _ -> a Right quotient -> quotient xdg = divTry x ydg = divTry y mrGCD = gcd (crMaxRoot $ crAbs xdg) (crMaxRoot $ crAbs ydg) m = mod mrGCD (2*n') -- need to consider cyclotomic translations, order the terms (x', ml) | ydg /= y = ( (crDivRational ydg x), if (not b) then (-1) else 1) | otherwise = ( (crDivRational xdg y), 1); {- The more familar form of the below is (C(x)^2 - nxD(x)^2): gm(x)^2 -nx*dt(x)^2 => gamma +/- sqrt(nx) * delta -} xrtn = crMult cr' (crRoot x' n') xrtnr = crToRational xrtn sqrtnxr = crToRational $ crRoot (crMult cr' xrtn) 2 apply gm dt = Just (ml * numerator f1, numerator f2) where f1 = c - sqrtnxr * d f2 = c + sqrtnxr * d c = aA gm xrtnr d = aA dt xrtnr -- aA means applyArray. array/lists are treated like polynomials (zero-base assumed) aA a x'' = f (elems a) 1 0 where f [] _ a' = a' f (c:cs) m' a' = f cs (m'*x'') (a' + (toRational c)*m') -- Example Aurif. decomp: C5(x) = x^2 + 3x + 1, D5(x) = x + 1 => Cyclotomic5(x) = C5(x)^2 − 5x*D5(x)^2 -- | This function returns a pair of polynomials (in array form) or Nothing (if it's squareful). -- An illogical n (n <= 1) will generate an error. aurDec :: Integer -> Maybe (Array Integer Integer, Array Integer Integer) aurDec n | n <= 1 = error "aurifDecomp: n must be greater than 1" | otherwise = aurDecI n (fst $ crFromI n) -- no concern about factorization cutoff here -- | Internal Aurifeullian Decomposition Workhorse Function aurDecI :: Integer -> CR_ -> Maybe (Array Integer Integer, Array Integer Integer) aurDecI n cr | crHasSquare cr || n < 2 || n' < 2 = Nothing | otherwise = Just (gamma, delta) where nm4 = mod n 4 n' = if (nm4 == 1) then n else (2*n) d = div (totient n') 2 dm2 = mod d 2 -- max gamma and delta subscripts to explicitly compute (add'l terms come from symmetry) mg | dm2 == 1 = div (d-1) 2 | otherwise = div d 2 md | dm2 == 1 = div (d-1) 2 | otherwise = div (d-2) 2 -- create q array of size td2: q(2k+1) = jacobi n (2k+1), q(2k+2 = mn * totient (2k+2) q = array (1, d) ([(i, f i) | i <- [1..d]]) where f i | mod i 2 == 1 = convJacobi $ jacobi n i | otherwise = eQ where eQ | ff == True = moeb cr' * (totient g) * (cos' $ (n-1)*i) | otherwise = error "Moebius fcn can't be called if number not totally factored" (cr', ff)= crFromI $ div n' g g = gcd n' i -- moebius fcn: 0 if has square, otherwise based on number of distinct prime factors moeb crm | crHasSquare crm = 0 | mod (length crm) 2 == 1 = -1 | otherwise = 1 cos' c | m8 == 2 || m8 == 6 = 0 -- "cosine" function | m8 == 4 = -1 | m8 == 0 = 1 | otherwise = error "Logic error: bad/odd value passed to cos'" where m8 = mod c 8 -- These two arrarys have mutually recursive definitions gamma = array (0, d) ([(0,1)] ++ [(i, gf i) | i <- [1..d]]) where gf k | k > mg = gamma!(d-k) | otherwise = div gTerm (2*k) where gTerm = sum $ map f [0..k-1] where f j = n * q!(2*k-2*j-1) * delta!j - q!(2*k-2*j) * gamma!j delta = array (0, d-1) ([(0,1)] ++ [(i, df i) | i <- [1..d-1]]) where df k | k > md = delta!(d-k-1) | otherwise = div dTerm (2*k+1) where dTerm = gamma!k + sum (map f [0..k-1]) where f j = q!(2*k-2*j+1) * gamma!j - q!(2*k-2*j) * delta!j {- Pseudocode for computing gammas G and deltas D G(0) = 1 D(0) = 1 Evaluate G(k) for 1 .. floor(d/2) G(k) = (1/2k) * sum(n*q(2k-2j-1)*D(j) - q(2k-2j)G(j)) [for j= 0 to k-1) Evaluate D(k) for 1 .. floor(d-1/2) D(k) = (1/2k+1) * ( G(k) + sum(q(2k+1-2j)*D(j) - q(2k-2j)D(j)) ) Evaluate G(k) for (floor(d/2)+1) to d => G(k) = G(d-k) Evaluate D(k) for (floor(d+1/2)) to d-1 => D(k) = D(d-k) Cyc(n) = C(x)^2 -nxD(x)^2 where gamma and delta are the coeffs for C(x) and D(x) respectively -} -- | Internal function requires two integers (computed via Aurif. methods) along with a list of Integers. The product of -- the Integers must be a divisor of the list's product otherwise an error will be generated. -- It's called divvy because it splits the 2 integers across the array using the gcd. -- This will help factoring because the larger term(s) will be broken up into smaller pieces. divvy :: [Integer] -> Integer -> Integer -> [Integer] divvy a x y = d (sortBy rev a) (abs x) (abs y) where rev a' b' = if (a' > b') then LT else GT d [] x' y' | x' == 1 && y' == 1 = [] | abs x' == 1 && abs y' == 1 = [x' * y'] | otherwise = error "Empty list passed as first param but x' and y' weren't both 1" d (c:cs) x' y' | x' == 1 && y' == 1 = c:cs | otherwise = v ++ d cs (div x' gnx) (div y' gny) where v = filter (>1) $ [div q gny, gnx, gny] gnx = gcd c x' q = div c gnx gny = gcd q y' {- Test: Run divvy on this: let a = [4,7,8401,62324358089100907319521,682969,61,374857981681,547] let x = 50031545486420197 let y = 50031544711579219 -} {- Cyclotomic factorizations for numbers of the form x^n +/- y^n Example: x^15-y^15 = (x-y) (x^2+x y+y^2) (x^4+x^3 y+x^2 y^2+x y^3+y^4) (x^8-x^7*y +x^5*y^3 - x^4*y^4 +x^3*y^5 -x*y^7+y^8) C(15) is a form of the last term where y = 1 It's possible in some cases to do an additional Aurifeullian factorization (of the last term). -} type Poly = [Integer] -- look into optimizing at a later date -- | CycloPair: Pair of an Integer and its corresponding cyclotomic polynomial type CycloPair = (Integer, Poly) -- | CycloMapInternal: Map internal to CycloMap newtype type CycloMapInternal = M.Map CR_ CycloPair -- | CycloMap is a newtype hiding the details of a map of CR_ to pairs of integers and corresponding cyclotomic polynomials. newtype CycloMap = MakeCM CycloMapInternal deriving (Eq, Show) -- | Unwrap the CycloMap newtype. fromCM :: CycloMap -> CycloMapInternal fromCM (MakeCM cm) = cm -- | Unwrap the CycloMap and convert the internal canon rep keys to Integers, returning a "raw" map getIntegerBasedCycloMap :: CycloMap -> M.Map Integer CycloPair getIntegerBasedCycloMap cm = M.mapKeys crToI (fromCM cm) -- | This is an initial map with the cyclotomic polynomials for 1. crCycloInitMap :: CycloMap crCycloInitMap = MakeCM $ M.insert cr1 (1, gen_xNm1 1) M.empty -- | Wrapper function to query map internals cmLookup :: CR_ -> CycloMap -> Maybe CycloPair cmLookup c m = M.lookup c (fromCM m) -- Internal function for updating map internals cmInsert :: CR_ -> CycloPair -> CycloMap -> CycloMap cmInsert c p m = MakeCM $ M.insert c p (fromCM m) -- Computing the cyclotomic polynomials for the divisor set of a number: -- Begin with with a map of 2 elements to the cylotomic polynomial for 1 and 2 -- Check the radical and return the map: crCycloRad -- If the cr doesn't equal the radical, then open it up to the other factors and the square-free factors must be in the map -- Identity: x^n -1 is the product of cyclotomic polynomial for each d where d | n. -- | Integer wrapper for crCyclo with default CycloMap parameter cyclo :: Integer -> (CycloPair, CycloMap) cyclo n = crCyclo (fst $ crFromI n) crCycloInitMap -- no concern over factorization cutoff here -- | Integer wrapper for crCyclo cycloWithMap :: Integer -> CycloMap -> (CycloPair, CycloMap) cycloWithMap n m = crCyclo (fst $ crFromI n) m -- no concern over factorization cutoff here -- | Integer wrapper for crCycloDivSet with default CycloMap parameter cycloDivSet :: Integer -> CycloMap cycloDivSet n = fst $ crCycloDivSet (fst $ crFromI n) crCycloInitMap -- no concern over factorization cutoff here -- | Integer wrapper for crCycloDivSet cycloDivSetWithMap :: Integer -> CycloMap -> (CycloMap, CycloMap) cycloDivSetWithMap n m = crCycloDivSet (fst $ crFromI n) m -- no concern over factorization cutoff here -- | Return pair of expon. multiplier and radical's polynomial along with working cyclotomic map. crCyclo :: CR_ -> CycloMap -> (CycloPair, CycloMap) crCyclo cr m | crPositive cr = ((crToI $ crDivStrict cr r, p), m') | otherwise = error "crCyclo: Positive integer needed" where r = crRadical cr ((_,p), m') = crCycloRad r m -- | Return a pair of cyclo maps just for the divisors and then a master map. crCycloDivSet :: CR_ -> CycloMap -> (CycloMap, CycloMap) crCycloDivSet cr m | crPositive cr = m2 | otherwise = error "crCycloDivSet: Positive integer needed" where (_,m2) = c crd = crDivisors cr c = case cmLookup cr m of Nothing -> c' -- need to compute it Just p -> (p, (mf, m)) -- found it. add a filtered version. ToDo: optimize this where mf = MakeCM $ M.fromList $ filter (\(n,_) -> elem n crd) $ M.toList $ fromCM m {- Performance note: The filter version above tended to be somewhat faster than the lookup version below so I used that where mf = MakeCM $ M.fromList $ map l crd -- "lookup version" where l d = case cmLookup d m of Nothing -> error $ e d Just p -> (d, p) e d = error "crCycloDivSet: Logic error: Divisor = '" ++ (show d) ++ "' not found!" -} c' | r == cr = (pr, (pm, pm)) | otherwise = (cm, (mm, mm)) where (cm, mm) = mfn sqFulDivs pm r = crRadical cr (pr, pm) = crCycloRad r m sqFulDivs = crd \\ crDivisors r -- "squareful" divisors mfn [] _ = error "Logic Error in mfn: Empty list is forbidden" mfn (n:ns) mp | ns == [] = cp | otherwise = mfn ns mp' where cp@(_, mp') = crCycloAll n mp -- | Compute the "radical" divisors first and then the non-square free entries. crCycloRad :: CR_ -> CycloMap -> (CycloPair, CycloMap) crCycloRad cr m = case cmLookup cr m of Nothing -> c' -- need to compute it Just p -> (p, m) -- found it where c' | cs == [] = (cycpr, cmInsert cr cycpr m) | otherwise = (cyc_n, cmInsert cr cyc_n mp) where (_ : cs) = cr -- for primes, because the tail of the cr is [] meaning only one prime factor r = fromInteger $ crToI $ crRadical cr -- ToDo: Optimize cycpr to be quotient of (r^n -1)/(r-1) when r is a big prime cycpr = genPrimePoly r -- for composites -- Create polynomial of the form : x^n -1 xNm1 = gen_xNm1 r (cPrd, mp) = mf (init $ crDivisors cr) cyc_n = (1, divPoly xNm1 cPrd) -- mf (Memo Fold) takes a list of divisors and returns the pair: (cyclotomic product, memoized map) mf (n:ns) = f ns p m' where ((_,p), m') = crCycloRad n m f (n':ns') p' mp = f ns' (multPoly p' p'') m'' -- cycloMap-threaded mult. fold where ((_,p''), m'') = crCycloRad n' mp f _ p' mp = (p', mp) mf [] = error "Logic error: Blank list can't be passed to mf aka crCycloMemoFold" -- | Return a pair of Integer and its cyclotomic polynomial while efficiently building up a map crCycloAll :: CR_ -> CycloMap -> (CycloPair, CycloMap) crCycloAll cr m | p == cr1 = case cmLookup cr m of Nothing -> error "Logic error: Radical value not found for crCycloAll" Just cb -> (cb, m) -- found it | otherwise = (crp, cmInsert cr crp md) where (p, d) = crPullSq ((i,y), md) = case cmLookup d m of Nothing -> crCycloAll d m -- need to compute it Just c -> (c, m) -- found it -- Optimization: (1, x + 1) => (4, x + 1). Note: Cyc2(x) = x + 1 -- The first value is the exponential multiplier Cyc8(x) = C2(x^4) = (x^4) + 1 crp = ((fst $ head p) * i, y) crPullSq = f [] cr where f h [] = (cr1, h) f h (c@(cp,ce):cs) | ce > 1 = ([(cp, 1)], h ++ (cp, ce-1):cs) | otherwise = f (h ++ [c]) cs -- | These "apply cyclo" functions will use cyclotomic polynomial methods to factor x^e - b^e. applyCrCycloPair :: Integer -> Integer -> CR_ -> CycloMap -> ([Integer], CycloMap) applyCrCycloPair l r cr m = (applyCrCycloPairI l r cr (M.elems $ fromCM md), mn) where (md, mn) = crCycloDivSet cr m applyCrCycloPairI :: Integer -> Integer -> CR_ -> [CycloPair] -> [Integer] applyCrCycloPairI l r cr cds = map applyPoly cds where nd = crTotient cr pA v = a where a = array (0,nd) ([(0,1)] ++ [(i, v*a!(i-1)) | i <- [1..nd]]) -- array of powers lpa = pA l rpa = pA r applyPoly (m,p) = foldr1 (+) (map f $ zip p [0..]) where f (a, b) | a == 0 = 0 | otherwise = a * lpa!(m*b) * rpa!(m*(maxExp - b)) maxExp = toInteger $ length p - 1 -- | Wraps applyCycloPairWithMap with default CycloMap argument. applyCycloPair :: Integer -> Integer -> Integer -> [Integer] applyCycloPair x y e = fst $ applyCycloPairWithMap x y e crCycloInitMap -- | This will use cyclotomic polynomial methods to factor x^e - b^e. applyCycloPairWithMap :: Integer -> Integer -> Integer -> CycloMap -> ([Integer], CycloMap) applyCycloPairWithMap x y e m = applyCrCycloPair x y (fst $ crFromI e) m -- no concern over factorization cutoff here -- | This will display the cyclotomic polynomials for a CR. showCyclo :: CR_ -> CycloMap -> [Char] showCyclo n m = p $ snd $ fst $ crCyclo n m where p (c:cs) = show c ++ (p' cs (1 :: Int)) -- "LE" endianness is assumed here p _ = [] p' (c:cs) s | c == 0 = r | otherwise = (if c > 0 then " + " else " - ") ++ (if ac == 1 then "" else show ac) ++ "x" ++ (if s == 1 then "" else "^" ++ show s) ++ r where r = p' cs (s+1) ac = abs c p' _ _ = [] -- | All the exponents must be even to return True. crSquareFlag :: CR_ -> Bool crSquareFlag = all (\(_, ce) -> mod ce 2 == 0) -- | Wrapper for chineseAurifWithMap with default CycloMap parameter chineseAurif :: Integer -> Integer -> Bool -> Maybe (Integer, Integer) chineseAurif x y b = fst $ chineseAurifWithMap x y b crCycloInitMap -- | Integer wrapper for chineseAurifCr chineseAurifWithMap :: Integer -> Integer -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap) chineseAurifWithMap x y b m = chineseAurifCr (fst $ crFromI x) (fst $ crFromI y) b m -- practically speaking, factor cutoff is a non-issue -- The source for this algorithm is the paper by Sun Qi, Ren Debin, Hong Shaofang, Yuan Pingzhi and Han Qing -- http://www.jams.or.jp/scm/contents/Vol-2-3/2-3-16.pdf (The formula at 2.7 there is implemented below) -- This will handle a subset of the cases that the main Aurif. routines handle -- | This function reduces the two CR parameters by gcd before calling an internal function to find a "Chinese" Aurifeullian factorization. -- Solutions will be found for any non-zero multiple of xp yp b (where xp and yp are relatively prime) -- where xp is of the form: (q^2m * p) ^ (p * k) -- yp is of the form: (r^2n) ^ (p * k) -- (That said, the order of xp and yp can be switched and the same result or non-result would be obtained.) -- op = + if p == 3 mod 4 and op -- op = - if p == 1 mod 4 -- For xp and yp, all of variables (q, m, p, k, r, n) are (CRs equivalent to) positive integers. -- p and k are also odd. p must be square-free. q must be relatively prime to k. -- -- Integral Examples to Try: -- chineseAurif ((q^(2*m) * p) ^ (p * k)) ((r^(2*n))^(p * k)) True -- equivalent to 931^247 + 1 -- where (q, m, p, k, r, n) = (7, 1, 19, 13, 1, 1) -- -- chineseAurif ((q^(2*m) * p) ^ (p * k)) ((r^(2*n))^(p * k)) False -- equivalent to (5^16*29)^377 - 121^377 -- where (q, m, p, k, r, n) = (5, 8, 29, 13, 11, 1) -- chineseAurifCr :: CR_ -> CR_ -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap) chineseAurifCr xp yp b m = case c of Nothing -> chineseAurifI mbyx n myx (crToI myx) b m' -- if first try fails, try the reverse r -> (r, m') where (c, m') = chineseAurifI mbxy n mxy (crToI mxy) b m gcdxy = crGCD xp yp (x, y) = (crDivStrict xp gcdxy, crDivStrict yp gcdxy) -- strip out any commonality n = gcd (crMaxRoot $ crAbs x) (crMaxRoot $ crAbs y) ncr = fst $ crFromI n -- factorization cutoff not an issue here mbxy = crRoot (crDivRational x y) n mxy = crGCD (crNumer mbxy) ncr mbyx = crRecip mbxy myx = crGCD (crNumer mbyx) ncr -- | Internal function to find factor of mb^n +/- 1 (mb would be M from paper, mb meaning m "big". -- Solution forms: Any non-zero multiple of (q^2m * p) ^ (p * k) op (r^2n)^(p * k) where k is an odd, postive number, m, n > 0. -- This will work if the op is "+" when mod p 4 = 3 OR when op is "-" for when mod p 4 = 1. chineseAurifI :: CR_ -> Integer -> CR_ -> Integer -> Bool -> CycloMap -> (Maybe (Integer, Integer), CycloMap) chineseAurifI mbcr n mcr m b mp | mod n 2 == 0 || mod m 2 == 0 || -- n and m must both be odd m < 3 || km /= 0 || -- m must be odd and > 1 and m | n (mm4 == 1 && b) || -- sign and modulus (mm4 == 3 && not b) || -- must be in synch mbdm == cr0 || not (crSquareFlag mbdm) -- mb/m must be a square and integral = (Nothing, mp) | otherwise = case cv - (gd1 * gd2) of 0 -> (Just (gd1, gd2), mp') _ -> (Nothing, mp) -- addl check since paper doesn't indicate rationals are supported where mm4 = mod m 4 e = toRational $ if (mm4 == 3) then (-1) else mm4 (k, km) = quotRem n m mbdm = case crDiv mbcr mcr of Left _ -> cr0 -- error condition if not a multiple Right q -> q r = crToRational $ crRoot (crMult mbcr mcr) 2 -- sqrt (m*M) from paper mb = crToRational mbcr jR c = toRational $ convJacobi $ jacobi c m eM = e * mb v1 = (toRational m) * mb^(div (k * (m + 1)) 2) v2 = t * s where t = (jR 2) * r * (mb ^ (div (k-1) 2)) s = sum $ map (\c -> (jR c) * eM^(k*c)) $ filter (\c -> gcd c m == 1) [1..m] -- rel. prime ncr = fst $ crFromI n -- factor cutoff not an issue here -- get cyclotomic value cv = head $ applyCrCycloPairI (numerator eM) (denominator eM) ncr [cp] (cp, mp') = crCyclo ncr mp gd1 = gcd cv (numerator $ v1 - v2) -- delta1 gd2 = gcd cv (numerator $ v1 + v2) -- delta2 -- workaround after arithmoi changes convJacobi :: JacobiSymbol -> Integer convJacobi j = case j of MinusOne -> -1 Zero -> 0 One -> 1 -- Simple Polynomial functions just used internally -- Generate a polynomial of the form: x^n - 1 gen_xNm1 :: Int -> Poly gen_xNm1 r = -1 : (replicate (r-1) 0) ++ [1] genPrimePoly :: Int -> (Integer, Poly) genPrimePoly r = (1, replicate r 1) -- | Multiply two polynomials multPoly :: Num a => [a] -> [a] -> [a] multPoly [] _ = [] multPoly (p:p1) p2 = let pTimesP2 = multiplyBy p p2 xTimesP1Timesp2 = multiplyByX $ multPoly p1 p2 in addPoly pTimesP2 xTimesP1Timesp2 where multiplyBy a p' = map (a*) p' multiplyByX p' = 0:p' -- | Div function for polynomials. The assumption is that p1 is a multiple of p2 divPoly :: Integral a => [a] -> [a] -> [a] divPoly p1 p2 = go [] p1 (length p1 - length p2) where go q u n | n < 0 = q | otherwise = go (q0:q) u' (n-1) where q0 = div (head u) (head p2) u' = tail (addPoly u (map (\t -> -1 * t * q0) p2)) -- | Add two polynomials addPoly :: Num a => [a] -> [a] -> [a] addPoly p1 p2 = if (length p1 >= length p2) then (add' p1 p2) else (add' p2 p1) where add' p1' p2' = zipWith (+) p1' (p2' ++ repeat 0)