-- | -- Module: Math.NumberTheory.Canon -- Copyright: (c) 2015-2019 Frederick Schneider -- Licence: MIT -- Maintainer: Frederick Schneider -- Stability: Provisional -- -- A Canon is exponentation-based representation for arbitrarily massive numbers, including prime towers and hyper-expressions. {-# LANGUAGE PatternSynonyms, ViewPatterns, RankNTypes #-} module Math.NumberTheory.Canon ( Canon, makeCanon, makeCanon', BareStatus(..), CanonValueType(..), cShowFull, cShowFullAsCode, cShowAsCode, cShowAsCodeUnf, cShowUnf, cMult, cDiv, cAdd, cSubtract, cExp, cReciprocal, (>^), (<^), cGCD, cLCM, cMod, cOdd, cEven, cTotient, cPhi, cNegative, cPositive, cIntegral, cRational, cIrrational, cPrime, cSimplified, cSplit, cNumerator, cDenominator, cCanonical, cBare, cBareStatus, cValueType, cDelve, cIsPrimeTower, cPrimeTowerLevel, cSuperLog, cSuperLogCmp, -- Hyper levels 4 through 9 for these 4 lines cTetration, cPentation, cHexation, cHeptation, cOctation, cNonation, cTetrationL, cPentationL, cHexationL, cHeptationL, cOctationL, cNonationL, (<^>), (<<^>>), (<<<^>>>), (<<<<^>>>>), (<<<<<^>>>>>), (|<^>|), (<^^>), (<<^^>>), (<<<^^>>>), (<<<<^^>>>>), (<<<<<^^>>>>>), (|<^^>|), -- Operators for hyper levels 10-50 (~^~), (~<^>~), (~<<^>>~), (~<<<^>>>~), (~<<<<^>>>>~), -- 10-14 (~|^|~), (~|<^>|~), (~|<<^>>|~), (~|<<<^>>>|~), (~|<<<<^>>>>|~), -- 15-19 (~~^~~), (~~<^>~~), (~~<<^>>~~), (~~<<<^>>>~~), (~~<<<<^>>>>~~), -- 20-24 (~~|^|~~), (~~|<^>|~~), (~~|<<^>>|~~), (~~|<<<^>>>|~~), (~~|<<<<^>>>>|~~), -- 25-29 (~~~^~~~), (~~~<^>~~~), (~~~<<^>>~~~), (~~~<<<^>>>~~~), (~~~<<<<^>>>>~~~), -- 30-34 (~~~|^|~~~), (~~~|<^>|~~~), (~~~|<<^>>|~~~), (~~~|<<<^>>>|~~~), (~~~|<<<<^>>>>|~~~), -- 35-39 (~~~~^~~~~), (~~~~<^>~~~~), (~~~~<<^>>~~~~), (~~~~<<<^>>>~~~~), (~~~~<<<<^>>>>~~~~), -- 40-44 (~~~~|^|~~~~), (~~~~|<^>|~~~~), (~~~~|<<^>>|~~~~), (~~~~|<<<^>>>|~~~~), (~~~~|<<<<^>>>>|~~~~), -- 45-49 (~~~~~^~~~~~), -- FIFTY cAddOpLevel, cMultOpLevel, cExpOpLevel, cTetrOpLevel, -- Hyper levels 1-4 cPentOpLevel, cHexOpLevel, cHeptOpLevel, cOctOpLevel, cNonOpLevel, -- Hyper levels 5-9 cGetHyperList, cGetHyperOp, maxHyperOpDispLevel, maxHyperOpDelveLevel, cFactorSum, cConvertToSum, cMaxExpoToExpand, cFactorHorizon, cApplyHy, cHyperOp, cHyperExpr, cHyperExprAny, cMaxHyperOp, cMinHyperOp, cHyperSum, cHyperProd, cHyperExpo, cHyperSumAny, cHyperize, cQuasiCanonize, cQuasiCanonized, cCleanup, cGetAddends, cGetFactors, cCleanupAsNumDenPair, CanonElement, cGetBase, cGetExponent, cGetBases, cGetBasesDeep, cGetExponents, cGetElements, cNumDivisors, cTau, cDivisors, cNthDivisor, cWhichDivisor, cRelativelyPrime, cGetFirstNDivisors, cN1, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9, CycloMap, getIntegerBasedCycloMap, showCyclo, crCycloInitMap -- Exposes cyclotomic map-related functionality from AurifCyclo ) where import Math.NumberTheory.Primes (primes) import Math.NumberTheory.Primes.Testing (isPrime) import Data.List import Data.Maybe (fromMaybe) import GHC.Real (Ratio(..)) import Math.NumberTheory.Canon.Internals import Math.NumberTheory.Canon.Additive import Math.NumberTheory.Canon.AurifCyclo import Math.NumberTheory.Canon.Simple (CanonConv(..)) -- import Debug.Trace (trace) -- | CanonValueType: 3 possibilities for this GADT (integral, non-integral rational, irrational). -- Imaginary/complex numbers are not supported data CanonValueType = IntC | NirC | IrrC deriving (Eq, Ord, Show) -- | This element is a base, exponent pair. The base is an integer and is generally prime or 0, -1. -- The exponent is also a Canon (allowing for arbitrary nesting) -- A Canon conceptually consists of a list of these elements. The first member of the pair will -- be a Canon raised to the first power. By doing this, we're allow for further generality -- in the definition of a Canon. type CanonElement = (Canon, Canon) -- | GCR_ stands for Generalized Canonical Representation. This is internal to Canon. type GCR_ = [GCRE_] type GCRE_ = (Integer, Canon) -- | Canon: GADT for either Bare (Integer) or some variation of a Can(onical) form (see CanonValueType). data Canon = Bare Integer BareStatus | Can GCR_ CanonValueType | HX Canon [Canon] CanonValueType -- | BareStatus: A "Bare Simp" number means a prime number, +/-1 or 0. The code must set the flag properly -- A "Bare NSim" number is an Integer that has not been checked (to see if it can be factored). data BareStatus = Simp | NSim deriving (Eq, Ord, Show) -- | Create a Canon from an Integer. This may involve expensive factorization. makeCanon :: Integer -> Canon makeCanon n = fst $ makeCanon' n -- | Create a Canon from an Integer. Also return True if the number is fully factored makeCanon' :: Integer -> (Canon, Bool) makeCanon' n = (f cr, ff) where f POne = Bare 1 Simp f c | null cs && eh == 1 = if superLogI bh > superLogICutoff -- because we assume bare < hyper expr then error "Lib limitation: Can't handle massive bare numbers > cutoff" else Bare bh (if ff then Simp else NSim) | otherwise = Can g (gcrCVT g) where (bh,eh):cs = c g = map (\(p,e) -> (p, makeCanon e)) c -- (can't be reduced any further) (cr, ff) = crFromI n -- 2nd param, the totally factored flag not used at this time -- | Convert from underlying canonical rep. to Canon. The 2nd param indicates whether or not to force factorization/simplification. crToC :: CR_ -> Bool -> Canon crToC POne _ = Bare 1 Simp crToC c _ | crSimplified c = Bare (fst $ head c) Simp -- a little ugly | otherwise = Can g (gcrCVT g) where g = map (\(p,e) -> (p, makeCanon e)) c -- | Instances for Canon instance Eq Canon where x == y = cEq x y -- | Internal value that corresponds with ~~~~~^~~~~~ (level 50 hyperoperation) maxHyperOpDispLevel :: Integer maxHyperOpDispLevel = 50; -- | Max hyper operaton level when converting to canonical form (for the sake of combining and reducing terms) maxHyperOpDelveLevel :: Canon maxHyperOpDelveLevel = makeCanon 100; -- These must correspond with the built-in and defined operators (from addition through hexation), except for ^ hyperOpStrings :: [String] -- ensure this is consistent with small canons / maxHyperOpDisplayLevel hyperOpStrings = [ "", "+", "*", "^", "<^>", "<<^>>", "<<<^>>>", "<<<<^>>>>", "<<<<<^>>>>>", "|<^>|", -- 0-9 "~^~", "~<^>~", "~<<^>>~", "~<<<^>>>~", "~<<<<^>>>>~", -- 10-14 "~|^|~", "~|<^>|~", "~|<<^>>|~", "~|<<<^>>>|~", "~|<<<<^>>>>|~", -- 15-19 "~~^~~", "~~<^>~~", "~~<<^>>~~", "~~<<<^>>>~~", "~~<<<<^>>>>~~", -- 20-24 "~~|^|~~", "~~|<^>|~~", "~~|<<^>>|~~", "~~|<<<^>>>|~~", "~~|<<<<^>>>>|~~", -- 25-29 "~~~^~~~", "~~~<^>~~~", "~~~<<^>>~~~", "~~~<<<^>>>~~~", "~~~<<<<^>>>>~~~", -- 30-34 "~~~|^|~~~", "~~~|<^>|~~~", "~~~|<<^>>|~~~", "~~~|<<<^>>>|~~~", "~~~|<<<<^>>>>|~~~", -- 35-39 "~~~~^~~~~", "~~~~<^>~~~~", "~~~~<<^>>~~~~", "~~~~<<<^>>>~~~~", "~~~~<<<<^>>>>~~~~", -- 40-44 "~~~~|^|~~~~", "~~~~|<^>|~~~~", "~~~~|<<^>>|~~~~", "~~~~|<<<^>>>|~~~~", "~~~~|<<<<^>>>>|~~~~", -- 45-49 "~~~~~^~~~~~"] -- FIFTY smallCanons :: [Canon] smallCanons = map (\n -> makeCanon n) [0..maxHyperOpDispLevel] -- | Levels starting with 1 in the hyperoperation hierarchy cAddOpLevel, cMultOpLevel, cExpOpLevel, cTetrOpLevel, cPentOpLevel, cHexOpLevel, cHeptOpLevel, cOctOpLevel, cNonOpLevel :: Canon (_: cAddOpLevel : cMultOpLevel : cExpOpLevel : cTetrOpLevel : cPentOpLevel : cHexOpLevel : cHeptOpLevel : cOctOpLevel : cNonOpLevel : _) = smallCanons -- | Various show functions: cShowFull - fully expand large primes and composites in Canon expression. -- "Unf" in name means don't factor unless it's too big too display -- "AsCode" in name means you can copy and paste the results and execute them. cShowFull, cShowFullAsCode, cShowAsCode, cShowAsCodeUnf, cShowUnf, cShowForEqChk :: Canon -> String cShowFull = cShow True False False False cShowFullAsCode = cShow True True False False cShowAsCode = cShow False True False False -- displays hyperexpr wrapped in parens cShowAsCodeUnf = cShow False True True False cShowUnf = cShow False False True False cShowForEqChk = cShow False False False True instance Show Canon where -- If debugging ... show = cShowAsCode -- so can it be pasted back in and run. Leave this way? Maybe default should not use { } and use use <^ show = cShow False False False False -- 1st bool = b: if True, full display of all integers, -- 2nd bool = p: if True, all parens for most hyperexprs, -- 3rd bool = i: if True, display unfactored integers where possible -- 4th bool = m: if True, when a hyOp is a sum or product, sort it, when to check for equality) -- Note: If parens flag is true and as long as the hyperOp doesn't exceed the max display level, -- you can copy and paste the expression back in as input cShow :: Bool -> Bool -> Bool -> Bool -> Canon -> String cShow b _ i _ (Bare n NSim) = showI b n False i -- False means composite cShow b _ i _ (Bare n Simp) = showI b n True i -- True means prime (or -1, 0, 1) cShow b p i m (HX h l _) | p && (cHyperExprAny h || h > maxSmallC) = "cApplyHy " ++ showH h ++ " [" ++ (concat $ intersperse ", " $ map (cShow b p i m) cl) ++ "] True" -- fmt as fcn call! | otherwise = fmt1 (head cl) ++ s' (tail cl) where cl | h == cAddOpLevel && any cNegative l = pR ++ nR -- put the negatives in back | otherwise = l where (pR, nR) = partition cPositive l -- there should always be at least one of each. The sum must be positive fmt1 hD | not (cHyperExpr hD) && h == cMultOpLevel = cShow b p i m hD | otherwise = showH hD showH c | h == cAddOpLevel || cBare c || (i && canConvToI c) = rep -- showH small helper function for clarity of expression | (not p) && cHyperExpr c = "{" ++ rep ++ "}" | otherwise = "(" ++ rep ++ ")" where rep = cShow b p i m c fmtHy f | f = "-" -- indicates a negative sign, flip a plus to minus | cHyperExprAny h || h > maxSmallC = fmt' h | p && h == cExpOpLevel = "<^" | otherwise = hyperOpStrings !! (fromInteger $ cToI h) -- and write cApplyHy for exp where fmt' c | (not p) && cHyperExpr c = "" | otherwise = "" where rep = cShow b p i m c rep' | not m = rep | cGetHyperOp c == cAddOpLevel = cShow b p i m (fst $ cConvertToSum c) | cGetHyperOp c == cMultOpLevel = cShow b p i m c | otherwise = rep s' (x:xs) = " " ++ fmtHy f' ++ " " ++ showH (if f' then (negate x) else x) ++ s' xs where f' = h == cAddOpLevel && cNegative x s' _ = "" maxSmallC = smallCanons !! (fromInteger $ maxHyperOpDispLevel) cShow b p i m c | denom == c1 = s numer False | otherwise = s numer True ++ " / " ++ s denom True where (numer, denom) = cSplit c s (Bare n f) _ = cShow b p i m (Bare n f) s v w | i && canConvToI v = show $ cToI v -- if the Canonical is not too big, convert it back to integer (when i flag is true) | w = "(" ++ catList ++ ")" | otherwise = catList -- w = with(out) parens where catList = concat $ intersperse " * " $ map sE $ gcr' $ cToGCR v gcr' g@(x:y:gs)| fst x == -1 && snd y == c1 = ((fst x * fst y, snd y) : gs) | otherwise = g -- above: display (-1,1), (2,1) as -2 gcr' g = g sE (p', e) | ptLevel > 2 = sp ++ " <^> " ++ s ptLevel True -- sE means show element | otherwise = case e of Bare 1 _ -> sp Bare _ _ -> sp ++ expOp ++ se _ -> sp ++ " " ++ expOp ++ " (" ++ se ++ ")" where ptLevel = cPrimeTowerLevelI e p' 1 sp = showI b p' (isPrime p' || p' == -1) i se = cShow b p i m e expOp = if p then "<^" else "^" canConvToI :: Canon -> Bool canConvToI c = not $ cSuperLogGT (fst $ cSuperLog c) cSuperLogCutoff -- Allow via first parameter to suppress full printing of massive integers and just indicate an "x-digit number" showI :: Bool -> Integer -> Bool -> Bool -> String showI b n pOrC i | i = show n -- just as is | not pOrC = "[" ++ txt ++ "]" -- composites go in brackets (either number or digit count) | truncFlag = "(" ++ txt ++ ")" -- if prime but tooBig, put in parens | otherwise = txt -- just the number where txt = if truncFlag then (show $ nd n) ++ "-digit " ++ (if pOrC then "prime" else "composite") else show n truncFlag = (not b) && n > integerShowCutoff nd :: Integer -> Integer nd n' = nd' n' 1 -- Count digits 1 by 1. ToDo: Optimize where nd' n'' ct = if (n'' >= 10) then nd' (div n'' 10) (ct + 1) else ct instance Enum Canon where toEnum n = makeCanon $ fromIntegral n fromEnum c = fromIntegral $ cToI c instance Ord Canon where compare x y = cCmp x y instance Real Canon where toRational c | cIrrational c = toRational $ cToD c | otherwise = (cToI $ cNumerator c) :% (cToI $ cDenominator c) instance Integral Canon where toInteger c | cIntegral c = cToI c | otherwise = floor $ cToD c quotRem n m = fst $ cQuotRem n m crCycloInitMap -- tries to use map but ultimately throws it away mod n m = cMod n m instance Fractional Canon where fromRational (n :% d) = makeCanon n / makeCanon d (/) x y = fst $ cDiv x y crCycloInitMap -- tries to use map but ultimately throws it away instance Num Canon where -- tries to use the map but ultimately throws it away when using +, - and * operators fromInteger n = makeCanon n x + y = fst $ cAdd x y crCycloInitMap x - y = fst $ cSubtract x y crCycloInitMap x * y = fst $ cMult x y crCycloInitMap negate x = cNegate x abs x = cAbs x signum x = cSignum x -- | Checks if the Canon is Canonical, a more complex expression. cCanonical :: Canon -> Bool cCanonical (Can _ _) = True cCanonical _ = False -- | Checks if the Canon just a "Bare" Integer. cBare :: Canon -> Bool cBare (Bare _ _) = True cBare _ = False -- | Returns the status for "Bare" numbers. cBareStatus :: Canon -> BareStatus cBareStatus (Bare _ b) = b cBareStatus _ = error "cBareStatus: Can only checked for 'Bare' Canons" -- | Return the CanonValueType (Integral, etc). cValueType :: Canon -> CanonValueType cValueType (Bare _ _) = IntC cValueType (Can _ v) = v cValueType (HX _ _ v) = v -- | Split a Canon into the numerator and denominator. cSplit :: Canon -> (Canon, Canon) cSplit c = (cNumerator c, cDenominator c) -- | Check for equality. cEq :: Canon -> Canon -> Bool -- cEq a b | trace ("cEq: (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined cEq (Bare x _ ) (Bare y _ ) = x == y cEq (Bare _ Simp) (Can _ _ ) = False cEq (Can _ _ ) (Bare _ Simp) = False cEq a@(HX _ _ _) b@(HX _ _ _) | signum a /= signum b = False | h1 /= h2 = False -- Confirm: Always true? | h1 < cExpOpLevel = cmpHyLists cL1 cL2 | otherwise = cShowForEqChk a' == cShowForEqChk b' -- Note: Kludge for deeply nested numbers where (cL1:cL2:_) = map cGetHyperList [a',b'] -- Quadratic compare: necessary to avoid sorting -- which can be problematic and expensive for hyperoperations cmpHyLists x y | length x /= length y = False | otherwise = c' x y [] c' l@(x:xs) (y:ys) bN | x == y = c' xs (bN++ys) [] | otherwise = c' l ys (bN ++ [y]) c' (_:_) _ _ = False c' _ y _ = null y -- "Endless" looping! (a', b') = (cQuasiCanonize $ fst $ cConvertToSum a, cQuasiCanonize $ fst $ cConvertToSum b) (a', b') = (fst $ cConvertToSum a, fst $ cConvertToSum b) -- ToDo: make this more robust? (h1, h2) = (cGetHyperOp a', cGetHyperOp b') cEq hx@(HX _ _ _) b | cBare b || cMaxHyperOp hx >= cPentOpLevel = False | otherwise = cValueType hx == cValueType b && cGetBases hx == cGetBases b && (cSuperLogCmp (fst $ cSuperLog hx) (fst $ cSuperLog b) == EQ) -- ToDo: Verify this is robust cEq a hx@(HX _ _ _) = cEq hx a cEq (Bare x NSim) y | cValueType y /= IntC = False | otherwise = x == cToI y cEq x y@(Bare _ NSim) = cEq y x cEq (Can x a ) (Can y b) = if a /= b then False else gcrEqCheck x y -- | Check if a Canon is integral and odd/even, respectively. Note: Return False for both if the Canon is not integral. -- See CanonValueType for possible cases. cOdd, cEven :: Canon -> Bool cOdd = cMod2Check 1 gcrOdd cEven = cMod2Check 0 gcrEven cMod2Check :: Int -> (GCR_ -> Bool) -> Canon -> Bool cMod2Check m _ (Bare x _) = mod x 2 == toInteger m cMod2Check _ f (Can c IntC) = f c cMod2Check m _ (HX PoA cL IntC) = mod (sum $ map (\c -> mod c c2) cL) 2 == smallCanons !! m -- match add "operator" cMod2Check m _ (HX PoM cL IntC) = mod (product $ map (\c -> mod c c2) cL) 2 == smallCanons !! m -- match on mult "operator" cMod2Check m _ (HX _ cL IntC) = mod (head cL) 2 == smallCanons !! m cMod2Check _ _ _ = False -- | GCD and LCM functions for Canon cGCD, cLCM :: Canon -> Canon -> Canon cGCD x y | cHyperExprAny x || cHyperExprAny y = head $ cMultiplicative x y Gcd | otherwise = cLGApply gcrGCD x y cLCM x y | cHyperExprAny x || cHyperExprAny y = head $ cMultiplicative x y Lcm | otherwise = cLGApply gcrLCM x y -- | Compare Function (cHyperCmp is internal) cCmp, cCmpH, cCmp' :: Canon -> Canon -> Ordering -- cCmp a b | trace ("cCmp: (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined cCmp (Bare x _) (Bare y _) = compare x y cCmp x@(Can _ _) y@(Bare _ _) = gcrCmp (cToGCR x) (cToGCR y) cCmp x@(Bare _ _) y@(Can _ _) = gcrCmp (cToGCR x) (cToGCR y) cCmp x@(Can _ _) y@(Can _ _) = gcrCmp (cToGCR x) (cToGCR y) cCmp x@(HX _ _ _) (Bare _ _) = if signum x == c1 then GT else LT -- Hyperexpr always has greater magnitude cCmp (Bare _ _) y@(HX _ _ _) = if signum y == c1 then LT else GT -- Inverse of above cCmp x y | signum y == c1 && signum x /= c1 = LT | signum x == c1 && signum y /= c1 = GT | signum x == cN1 && signum y == cN1 = cCmp (abs y) (abs x) | otherwise = cCmpH x y -- At this point, we are comparing positive hyper expressions. Should not be called directly. -- cCmpH a b | trace ("cCmpH: (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined -- Interferes with show cCmpH x@(Can _ _) y@(HX _ _ _) | not (cSuperLogGT (fst $ cSuperLog x) cSuperLogCutoff) = LT | otherwise = cCmpH (cConvertToHyperExpr x) y cCmpH x@(HX _ _ _) y@(Can _ _) | not (cSuperLogGT (fst $ cSuperLog y) cSuperLogCutoff) = GT | otherwise = cCmpH x (cConvertToHyperExpr y) cCmpH a@(HX h1 cL1 _) b@(HX h2 cL2 _) | a == b = EQ | (h1 == cAddOpLevel || h2 == cAddOpLevel) && aS /= a = cCmp aS bS | (h1 == cMultOpLevel || h2 == cMultOpLevel) && aR /= a -- we don't always take this. Otherwise, we can have an endless loop = cCmp aR bR | mOa > mOb + 1 = GT | mOb > mOa + 1 = LT | mOa > cHexOpLevel && mOb > cHexOpLevel && candPred a && candPred b -- To Verify: do the bases have to match = compare (tryLiftTail a) (tryLiftTail b) | flag1Less = cba | flag1More = cab | mP && exprDomination a b = GT | mP && exprDomination b a = LT | bP && lA' > lB' = GT -- ToDo: Further investigate if there are any exceptions to this? | bP && lA' < lB' = LT | h1 > cMultOpLevel && h1 == h2 && dominates cL1 cL2 True = GT | h2 > cMultOpLevel && h1 > h2 && dominates cL1 cL2 False = GT | h1 > cMultOpLevel && h2 == h1 && dominates cL2 cL1 True = LT | h1 > cMultOpLevel && h2 > h1 && dominates cL2 cL1 False = LT | bP = case compare (last cL1) (last cL2) of --For large enough hyOps, last entry says which is > EQ -> compare (reduce a) (reduce b) -- If equal try compare lists with all but last members cmp -> cmp | otherwise = cCmp' a b where (mOa, mOb, bP) = (cMaxHyperOp a, cMaxHyperOp b, h1 == h2 && h1 > cPentOpLevel) mP = mOa >= hyperOpCutoff && mOb >= hyperOpCutoff reduce c | len == 0 = error "Logic error in comparison: rFmt must have a hyper list with at least one entry" | len == 1 = head l | otherwise = simpleHX (cGetHyperOp c) (init l) -- create new hyper expr with all but the last entry where (l, len) = (cGetHyperList c, length l) (flag1Less, cba) = (h1 >= cPentOpLevel && h2 == h1 + 1, comp1Diff b a False) -- ToDo: Verify it handles embedded HEs. (flag1More, cab) = (h2 >= cPentOpLevel && h1 == h2 + 1, comp1Diff a b True) (hLa, hLb, lHa, lHb) = (cGetHyperList a, cGetHyperList b, length hLa, length hLb) (lA', lB') | lHa > lHb = (cApplyHy h1 (drop lD hLa) True, last hLb) | lHa < lHb = (last hLb, cApplyHy h2 (drop lD hLb) True) | otherwise = (last hLa, last hLb) where lD = abs (lHa - lHb) -- Modify this to see if there are any terms in common -- ((aS, bS), (aR, bR)) = (reduceSums a b, reduceProds a b) -- can cause endless looping ((aS, bS), (_,aR,bR)) = (reduceSums a b, simpleReduce a b False) cCmpH x y = error $ errorStrg -- We should never get to this spot in the code where errorStrg = "Logic error in cCmpH in program flow: " ++ show x ++ " vs. " ++ show y ++ "." {- Two known cases that will cause loops. compare (3 * ((7 <^> 5) * (5 <<^>> 8 <<^>> 6)) + 17 <<^>> 5 + 4) (3 * ((7 <^> 4) * (5 <<^>> 6 <<^>> 8)) + 2) compare (3 * 5 <<^>> 8 <<^>> 6 + 2) ( 3 * 5 <<^>> 6 <<^>> 8 + 4) -} reduceSums :: Canon -> Canon -> (Canon, Canon) reduceSums a b = (sum aS', negate $ sum bS') where (aS', bS') = partition cPositive (cGetAddends diff) -- low level diff, no infernal looping! diff = combineSum $ simpleHX cAddOpLevel (cGetAddends a ++ (map cNegate $ cGetAddends b)) reduceProds :: Canon -> Canon -> (Canon, Canon) -- reduceProds a b | trace ("reduceProds: (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined reduceProds a b = (aR', bR') where (_:aR':bR':_) = cMultiplicative a b Gcd dominates :: [Canon] -> [Canon] -> Bool -> Bool -- dominates a b _ | trace ("dominates: (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined dominates a' b' gtf = d' a' b' (0 :: Integer) -- gtf indicates the underlying hyper operation level was greater in a than b where d' (x:xs) (y:ys) pc | x < y = False | otherwise = d' xs ys (pc + if x > y then 1 else 0) d' _ (_:_) _ = False d' (_:_) _ _ = True d' _ _ pc = gtf || pc > 0 -- if flag set or positive ct, it dominates -- a has a hyper operation in its base one more than b's base. We are dealing with positive hyper expressions here -- ToDo: what if there are hyper expressions embedded comp1Diff :: Canon -> Canon -> Bool -> Ordering comp1Diff a' b' cF = if cF then r else flp r -- EQ in this context means inconclusive where hLA@(aB:aE:_) = cGetHyperList a' (lA2, lB2) = (length hLA, makeCanon $ toInteger $ length $ cGetHyperList b') flp r' = case r' of GT -> LT LT -> GT c -> c r | lA2 < ml6 = LT -- larger embedded hexation in pentated b | lA2 > 2 = GT -- tower for the larger hyperoperation. e.g. 6 <<<^>>> 7 <<<^>>> 3 -- would be larger than any pentation tower -- The above is equivalent to: 6 <<^>> (6 <<<^>>> 7 <<<^>>> 3 - 1) | aE > lB2 = GT -- The "exponent" for the larger is greater than the length of the smaller -- For instance: 5 <<^>> 7 <<^>> 8 would be less than 6 <<<^>>> 4. | aE < lB2 = LT -- e.g. Downgrade a = 6 [6,3] to 5 [6,6,6] and compare it to b (aE == lB) | otherwise = compare (simpleHX (cGetHyperOp b') (replicate (fromInteger $ cToI aE) aB)) b' ml6 = maxHypLen cHexOpLevel b' -- ToDo:adapt this so it finds the maximum chunk? -- maximum length of list based on hyper operation. Assumed to be the maximum in the expression maxHypLen :: Canon -> Canon -> Int maxHypLen h c = mhl c 0 where mhl c' mx | cHyperExprAny c' = if cGetHyperOp c' == h then max mx (length cL) else (foldl1 max $ map (maxHypLen h) cL) | otherwise = 0 where cL = cGetHyperList c' hyperOpCutoff :: Canon hyperOpCutoff = cTetrOpLevel -- unsigned values are assumed. This checks if s is less than d or less than a subexpression of d exprDomination :: Canon -> Canon -> Bool -- exprDomination d s | trace ("exprDomination: (s=" ++ show s ++ ") and (d=" ++ show d ++ ")") False = undefined exprDomination d s = eD d s False -- The flag indicates what whether we are already embedded or not in the structure eD :: Canon -> Canon -> Bool -> Bool -- eD d' s' b' | trace ("eD: (d' = " ++ show d' ++ ", s' = " ++ show s' ++ ", b' = " ++ show b' ++ ")") False = undefined eD d' s' b' | notBoth s' d' && not b' = s' < d' -- first level check | notBoth s' d' && b' = s' <= d' | s' == d' = b' -- equality shows domination if at an inner level | s' /= sRp = eD dRp sRp b' | rC = rC | b' && (compare d' s' /= LT) = True -- (if inside the nested expression). Could be expensive. | otherwise = False -- at last check if individual items in list dominate where notBoth x y = not (cMaxHyperOp x >= hyperOpCutoff && cMaxHyperOp y >= hyperOpCutoff) (sRs, dRs) = if b' && (cHyperSum s' || cHyperSum d') then reduceSums s' d' else (s', d') (sRp, dRp) = if b' && (cHyperProd sRs || cHyperProd dRs) then reduceProds sRs dRs else (sRs, dRs) rC = any (\e -> eD e s' False || eD e s' True) $ cGetHyperList d' -- Fall back comparison function. If the numbers are small enough and sufficiently close, -- they will be converted back to integers and compared. We are dealing with positive hyper expressions here. -- cCmp' a b | trace ("cCmp': (a=" ++ show a ++ ") and (b=" ++ show b ++ ")") False = undefined cCmp' a b | aH == cPentOpLevel && bH == cTetrOpLevel && any cHyperExprAny (cGetHyperList b) && cSuperLogGT slb sla = LT -- pentation vs. nested tetration | bH == cPentOpLevel && aH == cTetrOpLevel && any cHyperExprAny (cGetHyperList a) && cSuperLogGT sla slb = GT -- nested tetration vs. pentation | aH >= cPentOpLevel && aH > bH -- Note: comp1Diff handles the case where aH = bH + 1 = GT | bH >= cPentOpLevel && bH > aH = LT | aH <= cTetrOpLevel && bH <= cTetrOpLevel && cSuperLogGT sla slb = GT | aH <= cTetrOpLevel && bH <= cTetrOpLevel && cSuperLogGT slb sla = LT | bBh == cPentOpLevel && aBh == bBh = pCmp a b | aBh == cAddOpLevel || bBh == cAddOpLevel || cmpAddends /= EQ = cmpAddends | aBh == cMultOpLevel || bBh == cMultOpLevel || cmpFactors /= EQ = cmpFactors | aBh == bBh && aBh > cMultOpLevel && cmpHyperList /= EQ = cmpHyperList -- ToDo: -- Note: super log is only practical <= level 10 | cSuperLogGT sla slb = GT -- These two checks will handle cases like: compare (5 <^> 8 <<^>> 6) (17 <<^>> 5) | cSuperLogGT slb sla = LT | otherwise = error $ "Unable to accurately compare a = " ++ show a ++ " and b = " ++ show b where (aH, bH) = (cMaxHyperOp a, cMaxHyperOp b) (aBh, bBh) = (cGetHyperOp a, cGetHyperOp b) (sla, slb) = (fst $ cSuperLog a, fst $ cSuperLog b) cmpList f = compare (sort $ f a) (sort $ f b) cmpAddends = cmpList cGetAddends cmpFactors = cmpList cGetFactors cmpHyperList = cmpList cGetHyperList -- Only for pentation check pCmp :: Canon -> Canon -> Ordering pCmp a b | pA > pB = GT | pA < pB = LT | otherwise = cSuperLogCmp sla' slb' where pTail x = cApplyHy aBh (tail $ cGetHyperList x) False (pA, pB) = (pTail a, pTail b) sl x t = fst $ cSuperLog $ simpleHX aBh (x:[t - m + 2]) sla' = sl (head $ cGetHyperList a) pA slb' = sl (head $ cGetHyperList b) pB m = min pA pB aBh = cGetHyperOp a -- | wrapper to create apply a hyperoperation to a list cApplyHy :: Canon -> [Canon] -> Bool -> Canon -- the Bool says whether to raise an error for a null list cApplyHy ho a b = if length a == 0 && b then error "cApplyHy: Null list passed. Specified as fatal condition by calling fcn" else fst (cHyperOp ho a crCycloInitMap) -- This function will do any simplifications -- | Find the maximum hyperoperation embedded in a Canon cMaxHyperOp :: Canon -> Canon cMaxHyperOp = findSigHyOp max -- | Find the minimum hyperoperation embedded in a Canon. (If not at all, return zer0 cMinHyperOp :: Canon -> Canon cMinHyperOp = findSigHyOp mHo where mHo a b | a == b = a | a == c0 = b | b == c0 = a | otherwise = min a b -- Can be called with f = max or mHo findSigHyOp :: (Canon -> Canon -> Canon) -> Canon -> Canon findSigHyOp _ (Bare _ _) = c0 findSigHyOp f (Can g _) = foldl1 f $ map runningSig g where runningSig (_, e) | e == c1 = cMultOpLevel | otherwise = f cExpOpLevel (findSigHyOp f e) -- at least exp findSigHyOp f (HX h hl _) = f h (foldl1 f $ map (findSigHyOp f) hl) -- | QuotRem Function cQuotRem :: Canon -> Canon -> CycloMap -> ((Canon, Canon), CycloMap) cQuotRem x y m | cHyperExprAny x || cHyperExprAny y = ((hQ, c0), mR) -- ToDo: Handle non-zero modulus, say if x is a sum. | cIntegral x && cIntegral y = ((gcrToC q', md'), m'') | otherwise = error "cQuotRem: Must both parameters must be integral." where (q', md', m'') = case gcrDiv (cToGCR x) gy of Left _ -> (q, md, m') Right quotient -> (quotient, c0, m) where gy = cToGCR y md = cMod x y q = gcrDivStrict (cToGCR d) gy -- equivalent to: (x - mod x y) / y. (d, m') = cSubtract x md m (hQ, mR) = cDiv x y m -- | Mod function cMod :: Canon -> Canon -> Canon -- cMod c m | trace ("cMod: (c=" ++ show c ++ "), m=" ++ show m ++ ")") False = undefined cMod c m | not (cIntegral c) || not (cIntegral m) = error "cMod: Must both parameters must be integral" | c < m = c | m == cGCD c m = c0 -- c is a multiple of m. If m is a hyper expr, this all we can do for now | otherwise = makeCanon $ cModI c (cToI m) -- | Mod function between a Canon and integer. This is usually called by cMod cModI :: Canon -> Integer -> Integer -- cModI c m | trace ("cModI: (c=" ++ show c ++ "), m=" ++ show m ++ ")") False = undefined cModI _ 0 = error "cModI: Divide by zero error when computing n mod 0" cModI _ 1 = 0 cModI _ (-1) = 0 cModI Pc1 PIntPos = 1 cModI Pc0 _ = 0 cModI c m | cn && mn = -1 * cma | (cn && not mn) || (mn && not cn) = signum m * ((abs m) - cma) | otherwise = cModI' c m where (cn, mn, cma) = (cNegative c, m < 0, cModI' (abs c) (abs m)) -- cModI' b m' | trace ("cModI' (b=" ++ show b ++ "), m'=" ++ show m' ++ ")") False = undefined cModI' (Bare n _) m' = mod n m' cModI' (Can c' _) m' = if c == makeCanon m' then 0 else mod (product $ map pM c') m' where pM (b,e) = if (mod b m' == 0) then 0 else (pmI b (mmt e) m') mmt e = cModI e (totient m') -- optimization cModI' (HX PoA cL _) m' = mod (sum $ map (\ce -> cModI ce m') cL) m' cModI' h@(HX PoM cL _) m' = if (cModI (product $ cGetBases h) m' == 0) then 0 -- simple check if the bases are a multiple of the modulus else mod (product $ map (\ce -> cModI ce m') cL) m' cModI' (HX PoE cL _) m' = cModI (foldr1 (<^) (b':tail cL)) m' -- convert it a power tower where b' = makeCanon $ cModI (head cL) m' cModI' (HX h cL _) m' | h == cTetrOpLevel && length cL == 2 && not (cHyperExprAny (cL !! 1)) && totient twrHeight > m' = cModHyTwr (head cL) m' twrHeight -- exp | otherwise = cModHyper (head cL) m' -- to infinity and beyond :) where twrHeight = cToI $ cL !! 1 -- https://www.quora.com/What-would-be-the-remainder-if-Grahams-number-were-divided-by-2-4-5-6-7-8-9-or-10 -- ToDo : Optimize this for larger m especially powers of 10 to show trailing digits -- cModHyper b m' | trace ("cModHyper: (b=" ++ show b ++ "), m'=" ++ show m' ++ ")") False = undefined cModHyper b m' | (all (\e -> elem e bB) $ cGetBases $ mC') = 0 -- the base of a hyper expression must be multiple of m' | otherwise = cToI $ f mC' where (bB, mC') = (cGetBases b, makeCanon m') -- f mC | trace ("f: (mC=" ++ show mC ++ ")") False = undefined f mC | mC == c2 = if cEven b then c0 else c1 | otherwise = cMod (b <^ f phi) mC where phi = fst $ cPhi mC crCycloInitMap -- ToDo: Optimize: This runs in linear time. It could leverage the information above if "f phi" is low enough -- cModHyTwr b m' s | trace ("cModHyTwr: (b=" ++ show b ++ "), m'=" ++ show m' ++ ", s = " ++ show s ++ ")") False = undefined cModHyTwr b m' s | r == 0 = 0 | otherwise = cToI $ cm' (s-1) r where r = makeCanon $ cModI b m' cm' y lv | y == 0 = nv -- at end | otherwise = cm' (y-1) nv where nv = makeCanon $ cModI (r <^ lv) m' -- | Totient functions cTotient, cPhi :: Canon -> CycloMap -> (Canon, CycloMap) cTotient c m | (not $ cIntegral c) || cNegative c = error "Not defined for non-integral or negative numbers" | not $ cSimplified c = error "cTotient Can't compute if number not completely factored" | c == c0 = (c0, m) | otherwise = f (cToGCR c) c1 m where f [] prd m' = (prd, m') f ((p,e):gs) prd m' = f gs wp mw -- f is equivalent to the crTotient function but with threading of CycloMap -- => product $ map (\(p,e) -> (p-1) * p^(e-1)) cr where cp = makeCanon p -- "Canon-ize" cp above. Generally, this should be a prime already (pM1, mp) = cSubtract cp c1 m' (eM1, me) = cSubtract e c1 mp (pxeM1, mpm) = cExp cp eM1 False me (nprd, mprd) = cMult pM1 pxeM1 mpm (wp, mw) = cMult prd nprd mprd cPhi = cTotient -- | The thinking around the hyperoperators is that they should look progressively scarier :) -- | They range from level 4 / tetration (<^>) to level 50 (~~~~~^~~~~~). Please read .odp file for the naming convention. infixr <^>, <<^>>, <<<^>>>, <<<<^>>>>, <<<<<^>>>>>, |<^>|, -- 4-9 ~^~, ~<^>~, ~<<^>>~, ~<<<^>>>~, ~<<<<^>>>>~, -- 10-14 ~|^|~, ~|<^>|~, ~|<<^>>|~, ~|<<<^>>>|~, ~|<<<<^>>>>|~, -- 15-19 ~~^~~, ~~<^>~~, ~~<<^>>~~, ~~<<<^>>>~~, ~~<<<<^>>>>~~, -- 20-24 ~~|^|~~, ~~|<^>|~~, ~~|<<^>>|~~, ~~|<<<^>>>|~~, ~~|<<<<^>>>>|~~, -- 25-29 ~~~^~~~, ~~~<^>~~~, ~~~<<^>>~~~, ~~~<<<^>>>~~~, ~~~<<<<^>>>>~~~, -- 30-34 ~~~|^|~~~, ~~~|<^>|~~~, ~~~|<<^>>|~~~, ~~~|<<<^>>>|~~~, ~~~|<<<<^>>>>|~~~, -- 35-39 ~~~~^~~~~, ~~~~<^>~~~~, ~~~~<<^>>~~~~, ~~~~<<<^>>>~~~~, ~~~~<<<<^>>>>~~~~, -- 40-44 ~~~~|^|~~~~, ~~~~|<^>|~~~~, ~~~~|<<^>>|~~~~, ~~~~|<<<^>>>|~~~~, ~~~~|<<<<^>>>>|~~~~, -- 45-49 ~~~~~^~~~~~ -- FIFTY (<^>), (<<^>>), (<<<^>>>), (<<<<^>>>>), (<<<<<^>>>>>), (|<^>|), -- 4-9 (~^~), (~<^>~), (~<<^>>~), (~<<<^>>>~), (~<<<<^>>>>~), -- 10-14 (~|^|~), (~|<^>|~), (~|<<^>>|~), (~|<<<^>>>|~), (~|<<<<^>>>>|~), -- 15-19 (~~^~~), (~~<^>~~), (~~<<^>>~~), (~~<<<^>>>~~), (~~<<<<^>>>>~~), -- 20-24 (~~|^|~~), (~~|<^>|~~), (~~|<<^>>|~~), (~~|<<<^>>>|~~), (~~|<<<<^>>>>|~~), -- 25-29 (~~~^~~~), (~~~<^>~~~), (~~~<<^>>~~~), (~~~<<<^>>>~~~), (~~~<<<<^>>>>~~~), -- 30-34 (~~~|^|~~~), (~~~|<^>|~~~), (~~~|<<^>>|~~~), (~~~|<<<^>>>|~~~), (~~~|<<<<^>>>>|~~~), -- 35-39 (~~~~^~~~~), (~~~~<^>~~~~), (~~~~<<^>>~~~~), (~~~~<<<^>>>~~~~), (~~~~<<<<^>>>>~~~~), -- 40-44 (~~~~|^|~~~~), (~~~~|<^>|~~~~), (~~~~|<<^>>|~~~~), (~~~~|<<<^>>>|~~~~), (~~~~|<<<<^>>>>|~~~~), -- 45-49 (~~~~~^~~~~~) -- FIFTY :: Canon -> Canon -> Canon a <^> b = cTetration a b a <<^>> b = cPentation a b a <<<^>>> b = cHexation a b a <<<<^>>>> b = cHeptation a b a <<<<<^>>>>> b = cOctation a b a |<^>| b = cNonation a b a ~^~ b = cApplyHy (makeCanon 10) [a,b] True a ~<^>~ b = cApplyHy (makeCanon 11) [a,b] True a ~<<^>>~ b = cApplyHy (makeCanon 12) [a,b] True a ~<<<^>>>~ b = cApplyHy (makeCanon 13) [a,b] True a ~<<<<^>>>>~ b = cApplyHy (makeCanon 14) [a,b] True a ~|^|~ b = cApplyHy (makeCanon 15) [a,b] True a ~|<^>|~ b = cApplyHy (makeCanon 16) [a,b] True a ~|<<^>>|~ b = cApplyHy (makeCanon 17) [a,b] True a ~|<<<^>>>|~ b = cApplyHy (makeCanon 18) [a,b] True a ~|<<<<^>>>>|~ b = cApplyHy (makeCanon 19) [a,b] True a ~~^~~ b = cApplyHy (makeCanon 20) [a,b] True a ~~<^>~~ b = cApplyHy (makeCanon 21) [a,b] True a ~~<<^>>~~ b = cApplyHy (makeCanon 22) [a,b] True a ~~<<<^>>>~~ b = cApplyHy (makeCanon 23) [a,b] True a ~~<<<<^>>>>~~ b = cApplyHy (makeCanon 24) [a,b] True a ~~|^|~~ b = cApplyHy (makeCanon 25) [a,b] True a ~~|<^>|~~ b = cApplyHy (makeCanon 26) [a,b] True a ~~|<<^>>|~~ b = cApplyHy (makeCanon 27) [a,b] True a ~~|<<<^>>>|~~ b = cApplyHy (makeCanon 28) [a,b] True a ~~|<<<<^>>>>|~~ b = cApplyHy (makeCanon 29) [a,b] True a ~~~^~~~ b = cApplyHy (makeCanon 30) [a,b] True a ~~~<^>~~~ b = cApplyHy (makeCanon 31) [a,b] True a ~~~<<^>>~~~ b = cApplyHy (makeCanon 32) [a,b] True a ~~~<<<^>>>~~~ b = cApplyHy (makeCanon 33) [a,b] True a ~~~<<<<^>>>>~~~ b = cApplyHy (makeCanon 34) [a,b] True a ~~~|^|~~~ b = cApplyHy (makeCanon 35) [a,b] True a ~~~|<^>|~~~ b = cApplyHy (makeCanon 36) [a,b] True a ~~~|<<^>>|~~~ b = cApplyHy (makeCanon 37) [a,b] True a ~~~|<<<^>>>|~~~ b = cApplyHy (makeCanon 38) [a,b] True a ~~~|<<<<^>>>>|~~~ b = cApplyHy (makeCanon 39) [a,b] True a ~~~~^~~~~ b = cApplyHy (makeCanon 40) [a,b] True a ~~~~<^>~~~~ b = cApplyHy (makeCanon 41) [a,b] True a ~~~~<<^>>~~~~ b = cApplyHy (makeCanon 42) [a,b] True a ~~~~<<<^>>>~~~~ b = cApplyHy (makeCanon 43) [a,b] True a ~~~~<<<<^>>>>~~~~ b = cApplyHy (makeCanon 44) [a,b] True a ~~~~|^|~~~~ b = cApplyHy (makeCanon 45) [a,b] True a ~~~~|<^>|~~~~ b = cApplyHy (makeCanon 46) [a,b] True a ~~~~|<<^>>|~~~~ b = cApplyHy (makeCanon 47) [a,b] True a ~~~~|<<<^>>>|~~~~ b = cApplyHy (makeCanon 48) [a,b] True a ~~~~|<<<<^>>>>|~~~~ b = cApplyHy (makeCanon 49) [a,b] True a ~~~~~^~~~~~ b = cApplyHy (makeCanon 50) [a,b] True cTetration, cPentation, cHexation, cHeptation, cOctation, cNonation :: Canon -> Canon -> Canon -- | Tetration Function - Level 4 cTetration a b = cApplyHy cTetrOpLevel [a,b] True -- | Pentation Function - Level 5 cPentation a b = cApplyHy cPentOpLevel [a,b] True -- | Hexation Function - Level 6 cHexation a b = cApplyHy cHexOpLevel [a,b] True -- | Heptation Function - Level 7 cHeptation a b = cApplyHy cHeptOpLevel [a,b] True -- | Octation Function -- Level 8 cOctation a b = cApplyHy cOctOpLevel [a,b] True -- | Nonation Function -- Level 9 cNonation a b = cApplyHy cNonOpLevel [a,b] True -- | Hyperoperation List Operators. On display, the towers will have single caret operators interspersed. infixr 9 <^^>, <<^^>>, <<<^^>>>, <<<<^^>>>>, <<<<<^^>>>>>, |<^^>| (<^^>), (<<^^>>), (<<<^^>>>), (<<<<^^>>>>), (<<<<<^^>>>>>), (|<^^>|) :: Canon -> [Canon] -> Canon a <^^> b = fst $ cTetrationL a b crCycloInitMap a <<^^>> b = fst $ cPentationL a b crCycloInitMap a <<<^^>>> b = fst $ cHexationL a b crCycloInitMap a <<<<^^>>>> b = fst $ cHeptationL a b crCycloInitMap a <<<<<^^>>>>> b = fst $ cOctationL a b crCycloInitMap a |<^^>| b = fst $ cNonationL a b crCycloInitMap cTetrationL, cPentationL, cHexationL, cHeptationL, cOctationL, cNonationL :: Canon -> [Canon] -> CycloMap -> (Canon, CycloMap) -- | Tetration List Function cTetrationL a b m = cHyperOp cTetrOpLevel (a:b) m -- | Pentation List Function cPentationL a b m = cHyperOp cPentOpLevel (a:b) m -- | Hexation List Function cHexationL a b m = cHyperOp cHexOpLevel (a:b) m -- | Heptation List Function cHeptationL a b m = cHyperOp cHeptOpLevel (a:b) m -- | Octation List Function cOctationL a b m = cHyperOp cOctOpLevel (a:b) m -- | Nonation List Function cNonationL a b m = cHyperOp cNonOpLevel (a:b) m -- | Generalized Hyperoperation Function (https://en.wikipedia.org/wiki/Hyperoperation) cHyperOp :: Canon -> [Canon] -> CycloMap -> (Canon, CycloMap) -- cHyperOp n l _ | trace ("cHyperOp: (ho=" ++ show n ++ "), l=" ++ show l ++ ")") False = undefined cHyperOp n l@(a:b:cs) m | any (not . cIntegral) (n:l) = error "cHyperOp requires the 1st 2 parameters to be integral at this time." | b < cN1 && n > cExpOpLevel = error $ hyperLvlError b n | n > c2 && any cNegative (b:cs) = error "cHyperOp: At this time, all trailing entries must be >= 0 when using exponentiation or greater." | cNegative a && n > c3 = error "cHyperOp: At this time, the base must be >= 0 when using tetration or greater." | cNegative a && n == c3 = (if oddPwr then negate absHe else absHe, m) | n < c0 = error "cHyperOp: Requires the level n >= 0" | any (== c0) l = if n == cAddOpLevel then filterV c0 else (if n == cMultOpLevel then (c0, m) else stripVs c0) | any (== c1) l && n > cAddOpLevel = if n == cMultOpLevel then filterV c1 else stripVs c1 | (a /= c0 && a /= c1 && b > c1 && not (a == c2 && b == c2)) || n == c1 || n == c2 = tryToDemoteOrPromote | null cs' = cHyperOpSpecial (toInteger n) a b m | otherwise = error "Can not handle special cases with more than 2 params at this time" where -- ToDo: Weave in the cycloMap -- Note: This tetration demotion logic is closely tied to the cSuperLogCutoff -- The idea that anything internally considered as a hyperexpression must be greater than -- the cutoff which is currently 10^10^5. Even 22934 ^ 22934 absHe = fst $ cHyperOp n ((abs a):b:cs) m oddPwr = cOdd $ fst $ cHyperOp n (b:cs) m tryToDemoteOrPromote | n == cAddOpLevel = (sum l, m) | n == cMultOpLevel = (product l, m) | hyperFree && n == cExpOpLevel = (foldr1 (<^) l, m) -- Note: The underlying function calls cHyperOp for hyper expressions | n == cPentOpLevel && l == [c3, c3] -- expand to 3 <^> 3 <^> 3 so it can be reduced = (c3 <^> c3 <^> c3, m) | null cs && n == cTetrOpLevel && b == 2 = (a <^ a, m) -- tetration to exp | null cs && n > cTetrOpLevel && b == 2 = if cGetHyperOp a == nM1 then (cApplyHy nM1 (a:(cGetHyperList a)) True, m) -- e.g. (3<^>4) <<^>> 2 => (3<^>4)<^>3<^>4 else (cApplyHy nM1 [a,a] True, m) -- e.g. can't append: (3<<^>>4) <<^>> 2 = (3<<^>>4)<<^>>(3<<^>>4) | a == 2 && ( ((null cs && ((n == 5 && b == 3) || (n == 4 && b == 4)) ) || (length cs == 1 && head cs == 2 && n == 4 && b == 2) ) ) = (makeCanon (65536 :: Integer), m) -- 2^2^2^2 = 2 <^> 4 = 2 <^> 2 <^> 2 = 2 <<^>> 3 | a == 2 && b == 3 && null cs -- 2 <<<^>>> 3 = 2 <<^>> (2 <<^>> 2) = 2 <<^>> 4 -- Special demotion case for 2 = (cApplyHy nM1 [a, c4] True, m) | a == 2 && b == 4 && null cs && n == 5 = (cApplyHy nM1 [a, cApplyHy n [a, b - 1] True] True, m) -- another demotion: 2 <<^>> 4 = 2 <^> 65536. Both help with comparisons | null cs && ((a <= 6 && n == 4 && b == 3) || (a == 3 && n == 5 && b == 2)) = (a <^ a <^ a, m) | ((lL > 2 && n > cMultOpLevel) || (lL >= 2 && n <= cMultOpLevel)) && sameVal l -- e.g. (1 + 5<^>7) ^ (1 + 5<^>7) ^ (1 + 5<^>7) = (1+ 5<^>7) <^> 3 = (promotedC, m) | otherwise = (cleanup b, m) where nM1 = n - c1 hyperFree = not $ any cHyperExprAny l sameVal (x:xs) = s' xs where s' (v:vs) | v == x = s' vs | otherwise = False s' _ = True sameVal _ = error "cHyperOp: List with at least two items expected" (lL, lenC) = (length l, makeCanon $ toInteger lL) promotedC = case n of 1 -> (head l) * lenC 2 -> (head l) <^ lenC _ -> cApplyHy (n+1) [head l, lenC] True cs' = if b == c1 then [] else cs -- blank out cs if b == 1 -- ToDo : always correct? defHypEx = HX n l IntC -- this just takes the input and creates a HyExp. Might not be what's returned later filterV v = (cApplyHy n (filter (/= v) l) False, m) stripVs v = (cApplyHy n nl False, m) where nl = if v == c0 then s l [] else (fst $ span (/= v) l) -- e.g. [2,3,1,4,5] -> [2,3] s (x:xs) wl = if x == c0 then (if (ct0 xs 0) == 1 then wl -- two trailing zeros evaulate to 1 else (if length wl > 0 then init wl else []) ) else s xs (wl ++ [x]) where ct0 (y:ys) ct = if y == c0 then ct0 ys (ct+1) else ct ct0 _ ct = ct :: Integer s _ _ = error "Logic error in strip0: should not get to the end" -- Examples [2,3,4,0,5] => 0^5 = 0 so [2,3,4,0] -> [2,3,1] -> [2,3] -- for 0 [2,3,4,0,0] => 0^0 = 1 so [2,3,4,1] -> [2,3,4] -- [2,3,0,0,0] => [2,3,0,1] -> [2,3,0] -> [2,1] -> [2] -- Upgrade Chain Example Below n = 4, a = 7, b = HX (n+1) [a, 13] IntC => HX 5 [7,13] IntC -- Then, the answer should be HX (n+1) [a,t+1] IntC = HX 5 [7,14] IntC -- 7 <^> (7 <<^>> 13) simplifies to 7 <<^>> 14 -- Function for regular cases -- cleanup n' | trace ("cleanup: Processing: (ho=" ++ show n' ++ ")") False = undefined cleanup (HX h cL@(a':e:xs) _) | a == a' && h == n + c1 && null xs && null cs = cApplyHy h [a, e + c1] False -- Upgrade Chain | h == n && null cs = cApplyHy h (a:cL) False -- combine into longer chain: 5<^>(7<^>7) = 5<^>7<^>7 | n == eL && cGetHyperOp a == mL = distProdForExpo | otherwise = defHypEx cleanup _ | n == eL && cGetHyperOp a == mL = distProdForExpo | otherwise = defHypEx distProdForExpo = computeExpr mL $ map (\p -> f (p:es)) $ cGetHyperList b' -- dist expo if it can't be upgraded where (b':es) = l f l'@(x:xs) | cGetHyperOp x == eL = computeExpr eL [bX, eXEval * computeExpr eL xs] -- (x^a)^b = x^(a*b) | cHyperExprAny x = computeExpr eL l' | otherwise = foldr1 (<^) l' where (bX:eX) = cGetHyperList x eXEval = computeExpr eL eX f _ = error "Logic Error: Empty list found in cleanup" (mL, eL) = (cMultOpLevel, cExpOpLevel) cHyperOp h (a:_) m | h < c0 || not (cIntegral h) = error "cHyperOp: Hyper operator must be >= 0 and integral" | otherwise = (a, m) cHyperOp h l m | h < cAddOpLevel || not (cIntegral h) = error "cHyperOp: Nullary value not defined if hyper operator is lower than addition" | h == cAddOpLevel = (sum l, m) | otherwise = (product l, m) hyperLvlError :: Canon -> Canon -> String hyperLvlError b n = "cHyperOp: Hyperexpr not defined when b < -1 and n is beyond exponentiation. b = " ++ show b ++ ", n = " ++ show n ++ "." -- go through the map and flatten any sums/products in the list -- take the list. partition it by having "hyper" expressions or not. Collapse the non-hyper entries cFlattenAndGroup :: Canon -> [Canon] --cFlattenAndGroup c | trace ("cFlattenAndGroup: Processing: (c=" ++ show c ++ ")") False = undefined cFlattenAndGroup c = cFlattenAndGroup' c cMaxExpoToExpand cFlattenAndGroup' :: Canon -> Canon -> [Canon] cFlattenAndGroup' c mx | n1 == cAddOpLevel || n1 == cMultOpLevel = fAndG | otherwise = [c] where n1 = cGetHyperOp c -- hyper op from input (cA,n,cL) = (abs c, cGetHyperOp cA, cGetHyperList cA) fAndG = fmt (gF nH) (sF h) where (gF, tF, sF) | n == cAddOpLevel = (sum, tryFlatSum, sortByHpo) | otherwise = (product, tryFlatProd, id) (h, nH) = partition cHyperExprAny $ concat $ map tF cL fmt nonHC hyL | n == cAddOpLevel && nonHC == c0 = hyL | n == cMultOpLevel && nonHC == c1 = hyL | otherwise = (nonHC:hyL) tryFlatSum v | cGetHyperOp v' == cAddOpLevel = cGetHyperList v' | otherwise = [v'] where v' = fst $ cConvertToSum' v mx tryFlatProd v | cGetHyperOp v == cMultOpLevel = cGetHyperList v | otherwise = [v] -- When operating on a sum, we can flatten some products and distribute them -- ToDo: Factor algebraic expressions with hyperoperations. -- Elements with more hyper expressions in base are sorted first. The lists of hyper ops are sorted in descending -- order. Products of equal "hyper length" will then be compared by the lists. sortByHpo :: [Canon] -> [Canon] -- sortByHpo v | trace ("sortByHpo: Processing v = " ++ show v ++ ".") False = undefined sortByHpo l' | length l' == 1 = l' | otherwise = filter (/= c0) $ map collHy $ groupBy (\x y -> snd x == snd y) $ sortBy sHPO $ map hpo $ map combineProd l' -- This groups by hyOps, bases pair where collHy z = combineSum $ computeExpr cAddOpLevel (map fst z) -- allows for crude sorting without doing any heavy lifting type CanonInfo = ([Canon], ([[Canon]], [Canon])) -- ([bases], ([["exponents"]],[hyper ops'])) combineProd, combineSum :: Canon -> Canon -- combineProd c | trace ("combineProd: Processing c = " ++ show c ++ ".") False = undefined combineProd c | cGetHyperOp c' == cMultOpLevel = simpleHX c2 (map fst $ reverse $ sortBy sHPO $ map hpo $ cGetHyperList c') | otherwise = c' where c' = if cNegative c then negate aCm else aCm where aCm = combine cMultOpLevel $ abs c combineSum c = combine cAddOpLevel c combine :: Canon -> Canon -> Canon -- combine h c | trace ("combine: Processing h = " ++ show h ++ ", c = " ++ show c ++ ".") False = undefined combine h c | length (cGetHyperList c) < 2 || cGetHyperOp c /= h = c | nH == nullary = computeExpr h cLc | h == cAddOpLevel = computeExpr h (cLc ++ [nH]) -- leave non-hyper expressions at the end | otherwise = computeExpr h (nH:cLc) where (nH, cL') = (if h == cMultOpLevel then prepM else prepA) c cLc = if null cL' then [] else (combine' (tail cL') [head cL'] []) -- can be null for sums -- use the quadratic check logic and group them together -- factors (multiplicands) -- prepM c | trace ("prepM: Processing c=" ++ show c ++ ".") False = undefined prepM (HX PoE (b:es) _) = (c1, [(b, computeExpr nxtOp es)]) prepM c'@(HX PoM _ _) = (fld nHe, concat $ map (snd . prepM) hE) where (hE, nHe) = partition cHyperExpr $ cFlattenAndGroup c' prepM c'@(HX _ _ _) = (c1, [(c', c1)]) prepM c' = (c', []) -- addends prepA (HX PoM l _) = (c0, [(computeExpr nxtOp hE, product nHe)]) where (hE, nHe) = partition cHyperExpr l prepA (HX PoA l _) = (fld nHe, concat $ map (snd . prepA) hE) where (hE, nHe) = partition cHyperExpr l prepA c'@(HX _ _ _) = (c0, [(c', c1)]) prepA c' = (c', []) -- combine' c l wL | trace ("combine': Processing c=" ++ show c ++ ", l=" ++ show l ++", wL=" ++ show wL ++ ".") False = undefined combine' l@((xB,xE):xs) (yP@(yB,yE):ys) wL | xB == yB = combine' xs (wL ++ combinedTerm ++ ys) [] | otherwise = combine' l ys (yP:wL) where combinedTerm = if (xE + yE == c0) then [] else [(xB, xE + yE)] combine' (xP:xs) _ wL = combine' xs (xP:wL) [] combine' _ y _ = map f y -- this is the simplified list where f (b,e) | e == c1 = b | cGetHyperOp e == nxtOp = computeExpr nxtOp (b:cGetHyperList e) | otherwise = computeExpr nxtOp [b, e] (fld, nullary,nxtOp) | h == cAddOpLevel = (sum, c0, cMultOpLevel) | otherwise = (product, c1, cExpOpLevel) -- ToDo: Modify so that there are pairs of numbers when there are repeated exponents? Or is this close enough? hpo :: Canon -> (Canon, CanonInfo) hpo c' = (c', h' (abs c')) where h' c@(HX PoM l2 _) = (getHyperBases c, (reverse $ sort $ concat e2, reverse $ sort $ filter (\h -> h /= c0) $ concat h2)) where (e2, h2) = unzip $ map spHyOp $ filter cHyperExpr l2 h' c@(HX _ _ _) = (getHyperBases c, spHyOp c) -- Use the back hyperOp for now h' _ = ([], ([[]], [])) -- nothing to consider in sorting getHyperBases :: Canon -> [Canon] getHyperBases c = cGetBases' False False True c sHPO :: (Canon, CanonInfo) -> (Canon, CanonInfo) -> Ordering sHPO (_,(b1,(e1,hl1))) (_,(b2,(e2,hl2))) | hl1 > hl2 = LT | hl1 < hl2 = GT | otherwise = case compare hl2 hl1 of -- rev sort the hyper ops EQ -> case compare b1 b2 of -- and then the bases if needed EQ -> case compare (length e1) (length e2) of EQ -> compare e1 e2 cmpE -> cmpE cmpB -> cmpB cmp -> cmp -- give greater weight to hyper expressions raised to an exponent. ToDo: verify soundness spHyOp :: Canon -> ([[Canon]], [Canon]) -- spHyOp c | trace ("spHyOp: Processing: (" ++ show c ++ ")") False = undefined spHyOp c | h == cExpOpLevel && (cHyperExpr $ head hL) = (replicate nR $ tHl, replicate nR hH) | h == cMultOpLevel = (filter (/= []) $ sort $ concat e', sort $ concat h') -- handle product | otherwise = ([tHl], [h]) where (h, hL) = (cGetHyperOp c, cGetHyperList c) hH = cGetHyperOp $ head hL nR = if l > 2 || e > 1000 then 1000 else (fromInteger $ cToI e) -- ToDo: handle edge case when grouping where l = length $ cGetHyperList $ head hL e = head $ tail hL (e', h') = unzip $ map spHyOp $ filter cHyperExpr hL tHl | h == cExpOpLevel && length bHl > 1 = tail bHl | otherwise = if length hL > 1 then tail hL else [] where bHl = cGetHyperList $ head hL -- Function for special cases: -- Note: When n (first param) is zero, that is known as "succession" -- Cases when a is zero ... cHyperOpSpecial :: Integer -> Canon -> Canon -> CycloMap -> (Canon, CycloMap) cHyperOpSpecial 0 Pc0 b' m = cAdd b' c1 m cHyperOpSpecial 1 Pc0 b' m = (b', m) cHyperOpSpecial 2 Pc0 _ m = (c0, m) cHyperOpSpecial 3 Pc0 b' m = (if b' == c0 then c1 else c0, m) cHyperOpSpecial _ Pc0 b' m = (if cOdd b' then c0 else c1, m) -- Cases when b is zero ... cHyperOpSpecial 0 _ Pc0 m = (c1, m) cHyperOpSpecial 1 a' Pc0 m = (a', m) cHyperOpSpecial 2 _ Pc0 m = (c0, m) cHyperOpSpecial _ _ Pc0 m = (c1, m) -- Cases when b is -1 ... cHyperOpSpecial 0 _ PcN1 m = (c0, m) cHyperOpSpecial 1 a' PcN1 m = cSubtract a' c1 m cHyperOpSpecial 2 a' PcN1 m = (cNegate a', m) cHyperOpSpecial 3 a' PcN1 m = (cReciprocal a', m) cHyperOpSpecial _ _ PcN1 m = (c0, m) -- Other Cases ... cHyperOpSpecial h Pc2 Pc2 m | h == 0 = (smallCanons !! 3, m) | otherwise = (smallCanons !! 4, m) -- recursive identity cHyperOpSpecial _ Pc1 _ m = (c1, m) cHyperOpSpecial _ a' Pc1 m = (a', m) cHyperOpSpecial _ _ _ _ = error "Can't compute this hyperoperation. b must be >= -1. Not a 'special' case" -- | Return the list of canons from a hyper expression cGetHyperList :: Canon -> [Canon] cGetHyperList (HX _ cL _) = cL cGetHyperList _ = [] -- | Return the level of hyperoperation from a hyper expression. cGetHyperOp :: Canon -> Canon cGetHyperOp (HX h _ _) = h cGetHyperOp _ = c0 -- | Exponentiation and root operator declarations infixr 9 <^, >^ (<^), (>^) :: Canon -> Canon -> Canon a <^ b = fst $ cExp a b True crCycloInitMap r >^ n = cRoot n r -- | Convert a Canon back to its underlying rep (if possible). cToCR :: Canon -> CR_ cToCR (Can c v) | v /= IrrC = gcrToCR c | otherwise = error "cToCR: Cannot convert irrational canons to underlying data structure" cToCR (Bare 1 _ ) = cr1 cToCR (Bare n NSim) = fst $ crFromI n cToCR (Bare n Simp) = [(n,1)] -- not ideal cToCR c = gcrToCR $ cToGCR c -- this could be EXTREMELY expensive for hyper-expressions. -- | Convert generalized canon rep to Canon. gcrToC :: GCR_ -> Canon gcrToC g | gcrBare g = Bare (gcrToI g) Simp | otherwise = Can g (gcrCVT g) -- | For generalized canon rep, determine the CanonValueType. gcrCVT :: GCR_ -> CanonValueType gcrCVT POne = IntC gcrCVT g = g' g IntC -- start Integral, can only get "worse" where g' _ IrrC = IrrC -- short-circuits once irrational canon is found g' POne v = v g' ((_,ce):cs) v = g' cs (dcv v ce) -- check the exponents for expr's value type g' _ _ = error "gcrCVT : Logic error. Patterns should have been exhaustive" -- checking exponents dcv IrrC _ = IrrC dcv _ (Can _ IrrC) = IrrC dcv _ (Can _ NirC) = IrrC dcv IntC (Bare n _) = if n < 0 then NirC else IntC dcv v (Bare _ _) = v dcv v c = if cNegative c then NirC else v -- | Define some small canons for convenience cN1, c0, c1, c2, c3, c4, c5, c6, c7, c8, c9 :: Canon (cN1:c0:c1:c2:c3:c4:c5:c6:c7:c8:c9:_) = map makeCanon [-1..maxHyperOpDispLevel] impossibleHyperValue :: Canon impossibleHyperValue = cN1 -- used internally as a sentinel -- | Convert Canon to Integer if possible. cToI :: Canon -> Integer cToI (Bare i _ ) = i cToI c@(Can g v) | v == IntC && (cHyperExpr c || cSuperLogGT (fst $ cSuperLog c) cSuperLogCutoff) = error $ tooBigError c | v == IntC = gcrToI g | otherwise = error $ nonIntError c cToI h@(HX _ _ v) | v == IntC = error $ tooBigError h -- always too big. cHyperOp is in sync with cSuperLogCutoff | otherwise = error $ nonIntError h tooBigError, nonIntError :: Canon -> String tooBigError c = "This expression is too large to be converted to an integer: " ++ show c nonIntError c = "Can't convert a non-integral canon to an integer: " ++ show c -- | Convert Canon To Double. cToD :: Canon -> Double cToD (Bare i _) = fromIntegral i cToD (Can c _) = gcrToD c cToD (HX _ _ _) = error "This hyper-expression is too large to convert to a double" -- | Multiply Function: Generally speaking, this will be much cheaper than addition/subtraction which requires factoring. -- You are usually just merging lists of prime, exponent pairs and adding exponents where common primes are found. -- This notion is the crux of the library. -- -- Note: This can be used instead of the '*' operator if you want to maintain a CycloMap for performance -- reasons. cMult :: Canon -> Canon -> CycloMap -> (Canon, CycloMap) cMult Pc0 _ m = (c0, m) cMult _ Pc0 m = (c0, m) cMult Pc1 y m = (y, m) cMult x Pc1 m = (x, m) cMult x y m | not (cHyperExprAny x) && not (cHyperExprAny y) = (gcrToC g, m') | otherwise = (multH x y, m) -- This attempts to do deeper combining and can be problematic (head $ cMultiplicative x y Mult, m) where (g, m') = gcrMult (cToGCR x) (cToGCR y) m -- | Addition and subtraction is generally much more expensive because it requires refactorization. -- There is logic to look for algebraic forms which can greatly reduce simplify factorization. -- Note: This can be used instead of the +/- operators if you want to maintain a CycloMap for performance -- reasons. cAdd, cSubtract :: Canon -> Canon -> CycloMap -> (Canon, CycloMap) cAdd = cApplyAdtvOp True cSubtract = cApplyAdtvOp False -- | Internal Function to compute sum or difference based on first param. Much heavy lifting under the hood here. cApplyAdtvOp :: Bool -> Canon -> Canon -> CycloMap -> (Canon, CycloMap) -- cApplyAdtvOp _ x y _ | trace ("cApplyAdtvOp: Processing: (" ++ show x ++ ") and (" ++ show y ++ ")") False = undefined cApplyAdtvOp _ x Pc0 m = (x, m) cApplyAdtvOp True Pc0 y m = (y, m) -- True -> (+) cApplyAdtvOp False Pc0 y m = (negate y, m) -- False -> (-) cApplyAdtvOp b x y m | not b && x == y = (c0, m) | b && x == (negate y) = (c0, m) | not hax && not hay = (r, m') | otherwise = (addH x (if b then y else negate y), m) where (hax, hay) = (cHyperExprAny x, cHyperExprAny y) gcd' = cGCD x y -- non-hyper (x', y') = (x / gcd', y / gcd') r | tooBigToAdd x' || tooBigToAdd y' = simpleHX cAddOpLevel [x, if b then y else (negate y)] | otherwise = gcd' * (crToC c False) (c, m') = crApplyAdtvOptConv b (cToCR x') (cToCR y') m -- costly bit tooBigToAdd :: Canon -> Bool tooBigToAdd c@(Can _ _) | cHyperExprAny c = True | otherwise = cSuperLogGT (fst $ cSuperLog c) cSuperLogCutoff tooBigToAdd (HX _ _ _) = True tooBigToAdd (Bare _ _) = False -- | Exponentiation: This does allow for negative exponentiation if the Bool flag is True. -- Note: This can be used instead of the exponentiation operator if you want to maintain a CycloMap for performance -- reasons. cExp :: Canon -> Canon -> Bool -> CycloMap -> (Canon, CycloMap) -- cExp c e _ _ | trace ("cExp: Processing: " ++ show c ++ " <^ " ++ show e ++ ".") False = undefined cExp c e b m | cNegative e && (not b) = error "Per param flag, negative exponentiation is not allowed here." | cIrrational c && cIrrational e = error "cExp: Raising an irrational number to an irrational power is not currently supported." | otherwise = cExp' c e where cExp' _ Pc0 = (c1, m) cExp' Pc1 _ = (c1, m) cExp' Pc0 _ | cNegative e = error "0^e where e < 0 gives a div by zero error" | otherwise = (c0, m) cExp' _ Pc1 = (c, m) -- just return the value cExp' _ _ | cHyperExprAny c || cHyperExprAny e = (cApplyHy cExpOpLevel [c,e] True, m) | otherwise = (gcrToC g, mg) where (g, mg) = gE (cToGCR c) e m gE g' e' m' | gcrNegative g' = case cValueType e' of -- gcr exponentiation IntC -> if cOdd e' then (gcreN1:absTail, m'') else (absTail, m'') NirC -> if cOdd d then (gcreN1:absTail, m'') else error "gE: Imaginary numbers not supported" IrrC -> error "gE: Raising neg numbers to irr. powers not supported" | otherwise = f g' m' -- equivalent to multiplying each exp by e' (with CycloMap threaded) where (absTail, m'') = gE (gcrAbs g') e' m' (_, d) = cSplit e' -- even denom generates an imag. number f [] mf = ([], mf) f ((p,x):gs) mf = (fp, mf') where (prd, mx) = cMult e' x mf (t, mn) = f gs mx (fp, mf') = gcrMult [(p, prd)] t mn -- | Functions to check if a canon is negative/positive cNegative, cPositive :: Canon -> Bool -- cNegative c | trace ("cNegative: (l=" ++ show c ++ "))") False = undefined cNegative (Bare n _ ) = n < 0 cNegative (Can c _ ) = gcrNegative c cNegative (HX PoA cL _) | lp == 0 = True | ln == 0 = False | otherwise = (cCmp pH nH == LT) where (posCL, negCL') = partition cPositive cL negCL = map negate negCL' lp = length posCL ln = length negCL pH = cApplyHy cAddOpLevel posCL True nH = cApplyHy cAddOpLevel negCL True -- ToDo: Are there cases where combineSum could be used. Convert To Sum caused loops because it calls cNegative cNegative (HX PoM cL _) = cNegative $ head cL cNegative (HX _ _ _) = False -- tetration and beyond can only result in positive numbers cPositive (Bare n _ ) = n > 0 cPositive (Can c _ ) = gcrPositive c cPositive h@(HX PoA _ _) = not $ cNegative h -- zero is not possible in a hyper-expression. cPositive h@(HX PoM _ _) = not $ cNegative h -- same for products cPositive (HX _ _ _) = True -- tetration and beyond can only result in positive numbers -- | Functions for negation, absolute value and signum cNegate, cAbs, cSignum :: Canon -> Canon -- cNegate c | trace ("cNegate: Processing: v=("++show c ++ ")") False = undefined cNegate (Bare 0 _) = c0 cNegate (Bare 1 _) = cN1 cNegate (Bare (-1) _) = c1 cNegate (Bare x Simp) = Can (gcreN1 : [(x, c1)]) IntC -- prepend a "-1", not ideal cNegate (Bare x NSim) = Bare (-1 * x) NSim cNegate (Can x v) = gcrNegateCanonical x v -- HyperOp case: Product of canons. cNegate h@(HX PoA cL _) | cNegative h = simpleHX cAddOpLevel (reverse $ map negate cL) -- only should happen internally | otherwise = simpleHX cMultOpLevel [cN1,h] cNegate (HX PoM cL _) | hD == cN1 = cApplyHy cMultOpLevel (tail cL) True | nhH && cNegative hD = simpleHX cMultOpLevel ((abs hD):(tail cL)) -- change the leading term which should not | nhH && cPositive hD = simpleHX cMultOpLevel ((negate hD):(tail cL)) -- be hyper (if exists) in a product | otherwise = simpleHX cMultOpLevel (cN1:cL) -- prepend to existing list/product where (hD, nhH) = (head cL, not $ cHyperExpr hD) cNegate he@(HX _ _ _) = simpleHX cMultOpLevel [cN1, he] -- prepend a negative one to existing expression, making a new 2-element expr. -- this applies to hyper sums which internally will always be kept positive cAbs x | cNegative x = cNegate x | otherwise = x --cSignum c | trace ("cSignum: (c = " ++ show c ++ ")") False = undefined cSignum (Bare 0 _) = c0 cSignum g | cNegative g = cN1 | otherwise = c1 -- This internal function works for either gcrGCD or gcrLCM. cLGApply :: (GCR_ -> GCR_ -> GCR_) -> Canon -> Canon -> Canon cLGApply f x y | cNegative x || cNegative y = gcrToC $ f (cToGCR $ cAbs x) (cToGCR $ cAbs y) | otherwise = gcrToC $ f (cToGCR x) (cToGCR y) -- | This function tries to convert a hyper expression to "canonical" form. It is rather limited -- due to the way power towers branch for composite numbers in canonical form. Conversions can be used for non-integral division. tryToCanonizeHyperExpr :: Canon -> Maybe Canon tryToCanonizeHyperExpr c@(HX _ _ _) | cHyperSumAny c || cMaxHyperOp c > cTetrOpLevel || cMaxTetrLevel > 10 = Nothing | otherwise = Just $ conv c where cMaxTetrLevel = mtl c0 c mtl wM (HX h l@(x:xs) _) | h == cTetrOpLevel = foldl1 max [wM, mtl wM x, cApplyHy h xs True] | otherwise = foldl1 max (wM:(map (mtl wM) l)) mtl wM _ = wM conv (HX PoM l _) = product $ map conv l conv (HX PoE l _) = foldr1 (<^) $ map conv l conv (HX h l@(b:x:_) _) | h /= c4 = error "Logic error: Only tetration allowed here" | length l /= 2 = error "Logic error: Tetration list must only have two elements" | hB < c4 = convToTwr l | hB == c4 && simpleHyperExpr b && cMaxTetrBase <= 10 -- Note: Quite limited = nestedTetr x | otherwise = foldr1 (<^) $ replicate (fromInteger $ cToI $ l !! 1) (conv $ l !! 0) where (hB, bHl, bT) = (cGetHyperOp b, cGetHyperList b, convToTwr bHl) cMaxTetrBase = mtl c0 b convToTwr l' = foldr1 (<^) $ replicate (fromInteger $ cToI $ l' !! 1) (conv $ l' !! 0) nestedTetr oe | oe == 2 = b <^ b | otherwise = bT <^ nestedTetr (oe - 1) -- only relevant to tetration and above simpleHyperExpr c'@(HX h' l' _) | h' < cTetrOpLevel || not (cIntegral c') || length l' /= 2 || not (any cHyperExprAny l') = True | otherwise = False simpleHyperExpr _ = False conv c' = c' -- non-hyper expA tryToCanonizeHyperExpr c = Just c -- | Div function : Multiply by the reciprocal. cDiv :: Canon -> Canon -> CycloMap -> (Canon, CycloMap) cDiv _ Pc0 _ = error "cDiv: Division by zero error" cDiv Pc0 _ m = (c0, m) cDiv x y m | not (cHyperExprAny x) && not (cHyperExprAny y) = cMult (cReciprocal y) x m -- multiply by the reciprocal | y' == c1 = (x', m) -- x is a multiple of y (One or both is a hyper expr) | otherwise = case tryHyperDiv x y m of Right r -> r Left s -> error s where (x', y') = reduceProds x y -- do not call this directly. It assumes hyper operations tryHyperDiv :: Canon -> Canon -> CycloMap -> Either String (Canon, CycloMap) tryHyperDiv x y m | fmx /= hyDef && fmy /= hyDef = if (cCanonical fQ && cIntegral fQ && cHyperExprAny fQ) --canonical yet has "hyper exponents". Convert quotient to hyper expression. then Right (cConvertToHyperExpr fQ, m') else Right (fQ, m') | otherwise = Left ("At this time, one can only divide hyper expressions when x is a multiple of y, non-sums or limited tetrations: x = " ++ show x ++ ", y = " ++ show y) where fmch v = fromMaybe hyDef (tryToCanonizeHyperExpr v) (hyDef, fmx, fmy) = (impossibleHyperValue, fmch x, fmch y) (fQ, m') = cDiv fmx fmy m -- feed the canonical reps back into the function -- Converts an integral "Canonical" canon to a hyper product. Error if the canon is not integral Otherwise, it leaves the canon as is. -- ToDo: What if the result is not a hyper expr after going through the function? cConvertToHyperExpr :: Canon -> Canon cConvertToHyperExpr c | not (cIntegral c) = error "Cannot convert a non-integral canon to a hyper expression" | cCanonical c && cHyperExprAny c = cApplyHy cMultOpLevel ((product nHe):hE) False | otherwise = c where (hE, nHe) = partition cHyperExpr $ map hF $ cToGCR c hF (p, e) = if e == c1 then pC else (cApplyHy cExpOpLevel [pC, e] True) where pC = makeCanon p -- | Compute reciprocal (by negating exponents or equivalent). cReciprocal :: Canon -> Canon cReciprocal x | not (cHyperExprAny x) = fst $ cExp x cN1 True crCycloInitMap -- raise number to (-1)st power | fmx /= hyDef = cReciprocal fmx | otherwise = error $ "At this time, one can only take reciprocals of hyper expressions which are " ++ "non-sums and limited tetrations." where fmch v = fromMaybe hyDef (tryToCanonizeHyperExpr v) (hyDef, fmx) = (impossibleHyperValue, fmch x) -- | Functions to check if a Canon is Integral, (Ir)Rational, "Simplified", a prime or a prime tower cIntegral, cIrrational, cRational, cSimplified, cPrime, cIsPrimeTower :: Canon -> Bool cIntegral c = cValueType c == IntC cIrrational c = cValueType c == IrrC cRational c = not $ cIrrational c cSimplified (Bare _ Simp) = True cSimplified (Bare _ NSim) = False cSimplified (Can c _) = gcrSimplified c cSimplified c@(HX h l _) = h /= cAddOpLevel && ((cHyperProd c && all cSimplified l) || (cSimplified $ head l)) cPrime c = cSimplified c && c > c1 -- Simp includes 0, -1 cIsPrimeTower c = cPrimeTowerLevel c > 0 -- x^x would not be, but x^x^x would be -- | Utility functions regarding hyperoperations. "Any" functions search the entire expression cHyperExpr, cHyperExprAny, cHyperSum, cHyperSumAny, cHyperProd, cHyperExpo :: Canon -> Bool cHyperExpr = cHyperPredCheck (>= cAddOpLevel) False cHyperExprAny = cHyperPredCheck (>= cAddOpLevel) True cHyperSum (HX h (j:k:cs) _) = h == cAddOpLevel || (h == cMultOpLevel && j == cN1 && cGetHyperOp k == cAddOpLevel && null cs) cHyperSum _ = False cHyperSumAny = cHyperPredCheck (== cAddOpLevel) True -- when looking any we can just go by hyper op cHyperProd c@(HX PoM _ _) = not $ cHyperSum c -- Note: a negative sum is not considered a product cHyperProd _ = False; cHyperExpo = cHyperPredCheck (== cExpOpLevel) False -- checks if this is an exponential expression -- | Takes a predicate related to the hyper operation. It will search recursively if the 2nd flag is set. cHyperPredCheck :: (Canon -> Bool) -> Bool -> Canon -> Bool cHyperPredCheck f b c | f (cGetHyperOp c) = True | not b = False -- don't do the any check | otherwise = cHP' c where cHP' (HX _ l _) = any (cHyperPredCheck f b) l cHP' (Can g _) = any (cHyperPredCheck f b) $ map snd g cHP' _ = False -- | cNumerator and cDenominator are for processing "rational" canon reps. cNumerator, cDenominator :: Canon -> Canon cNumerator (Can c _ ) = gcrToC $ filter (\x -> cPositive $ snd x) c -- filter positive exponents cNumerator b = b cDenominator (Can c _ ) = gcrToC $ map (\(p,e) -> (p, cN1*e)) $ filter (\(_,e) -> cNegative e) c -- negate neg expnts cDenominator _ = c1 -- ToDo: For now, hyper expressions are always integral -- ToDo : Tweak cQuasiCanonize to make this function obsolete. The 2nd param isn't part of the QC function. -- cNestExpTail can be used whether or not the base is a hyper expression. It unravels tetration and beyond. -- For example: 7 <<<^>>> 8 = 7 <^ (7 <^> (-1 + 7 <<^>> (-1 + 7 <<<^>>> 7)). The expr after <^ would be the exp. tail cNestExpTail :: Canon -> Bool -> Canon cNestExpTail c'@(HX h (b:xs) IntC) bF | h == cAddOpLevel= c1 | h < cExpOpLevel = error errorMsg | otherwise = if bF then baseTail b * eH else eH -- If the flag is set, process the base as well. where expTail = cApplyHy h xs True eH = expRec h expTail expRec h' e' | h' == cExpOpLevel = e' -- otherwise, recursively demote down | otherwise = expRec (h' - c1) newE -- e.g. x<^>y<^>z = x^(x<^>((y<^>z)-1)) where newE = cApplyHy h' [b, e' - c1] True baseTail (HX PoA _ _) = c1 baseTail c@(HX PoM _ _) = error $ nonPrimePowerError c -- Limited but cQuasiCanonize kind of supplants this fcn baseTail c@(HX _ _ _) = cNestExpTail c bF baseTail c@(Can g _) | length g > 1 = error $ nonPrimePowerError c | otherwise = snd $ head g -- the exponent baseTail (Bare _ Simp) = c1 baseTail c = error $ nonPrimePowerError c errorMsg = "nestedExpTail: requires a hyper expression at level >= exponentiation: " ++ show c' cNestExpTail _ _ = c1 -- | Break code into a canonized cCleanup :: Canon -> Canon cCleanup = cHyperize . cQuasiCanonize -- | Split the hyperoperation into a cleaned-up numerator and denominator pair (if denom is 1). This still represents an integral value. e.g. 3 <^> 7 / 3 <^> 4 cCleanupAsNumDenPair :: Canon -> (Canon,Canon) cCleanupAsNumDenPair c = (n,d) where (n, d) = (cHyperize $ pr nL, cHyperize $ pr dL) qc = cQuasiCanonize c pr cL = simpleHX cMultOpLevel $ map expDemote $ filter (\(_, e) -> e /= c0) cL (nL, dL) = unzip $ map (\(p,(eP,eN)) -> ((p, eP), (p, eN)) ) $ map (\(p,e) -> (p, spl e)) $ map (\c' -> expPromote c') $ cGetFactors qc spl c' = (simpleHX cAddOpLevel pos, simpleHX cAddOpLevel (map negate neg)) where (pos, neg) = partition cPositive $ cGetAddends c' -- positive and negative entries in exponent sum expression -- | Hyperize will take a Canon in quasi-canonized form and try to clean it up in a tidier expression -- Example: 7 ^ ( 1 + 2 * (49 <^> 7) = 7 * 49 <^> 8. ToDo: Enhancement: Partial hyperizing? cHyperize :: Canon -> Canon cHyperize c | not (cQuasiCanonized c) || (h /= cExpOpLevel && h /= cMultOpLevel) || null iM = c | any cNegative $ concat $ map (\(_,e) -> cGetAddends e) $ map expPromote $ cGetFactors c -- = c -- For example, we can't cleanup 3 <^> 5 / 3 <^> 4 = 3 ^ (3<^>4 - 3<^>3) into a simple expression | not (null $ cGetBases' False True False $ simpleHX cMultOpLevel iM) -- in-scope bases are non-unique so not valid = c | not (foldl1 (&&) $ map snd process) -- not all "tail-convertible") = c | not (foldl1 (&&) $ map (\(_,l) -> allTheSame $ map snd l) grp) -- not all multipliers are the same = c | null grp' || not (foldl1 (&&) $ map snd grp') -- not all elements of each base accounted for = c | otherwise = iSp * oSp where h = cGetHyperOp c (iM, oM) = partition (\m -> cGetHyperOp m == cExpOpLevel) $ cGetFactors c process = map (\l -> hypMap (l !! 0, l !! 1)) $ map cGetHyperList iM oSp = product (oM ++ (map (snd . snd . fst) process)) -- everything that could not be rolled up grp = grpExpr $ concat $ map (\((p,(eMap,_)),_) -> map (\((_,t),m)-> (t, (p, m))) eMap) process grp' = map (\(e,l) -> ((e, snd $ head l), cBaseRadical e == product (map fst l))) grp iSp = product $ map (\((e,m),_) -> cApplyHy cExpOpLevel [e, m] True) grp' -- in scope product allTheSame l@(x:_:_) = and $ map (== x) (tail l) allTheSame _ = True grpExpr l@(_:_:_) = gE' l [] grpExpr ((e,p):_) = [(e, [p])] grpExpr _ = error $ "Blank list passed to grpExpr when processing c = " ++ show c gE' l@((xf,_):_) wL = gE' nM ((xf, map snd m):wL) -- all the add'l base info for that expression where (m,nM) = partition (\e -> xf == fst e) l gE' _ wL = wL -- Called by hyperize at this point, the constants should have been removed hypMap :: (Canon, Canon) -> ((Canon, ([((Canon, Canon), Canon)], Canon)), Bool) -- ToDo: Better to change this to a Maybe hypMap (p, e) = ((p, (mV', osProd)), not $ any (\((_,t),_) -> t == impossibleHyperValue) mV') where (iS, oS) = partition candPred $ cGetAddends e -- only process osProd = computeExpr cMultOpLevel (map (\x -> p <^ x) oS) mV' = mV p iS -- mV is short for mapped values mV :: Canon -> [Canon] -> [((Canon, Canon), Canon)] mV p iS | null iS = [] | otherwise = map (hypCheck p) $ cGetAddends $ grpAndSrtList p iS hypCheck :: Canon -> Canon -> ((Canon, Canon), Canon) hypCheck p c = ((c, liftedTail), p') -- e.g. p == 13 for (3 <<^>> 4) <^ 13 where (fs, base) = (cGetFactors c, head $ cGetHyperList $ head $ fs) liftedTail | b' /= c1 = impossibleHyperValue -- p must be a "clean" multiple of b | otherwise = tryLiftTail $ head fs (_, p', b') = simpleReduce (product $ tail fs) (qcBase p base) False -- -- group and sort list of canons grpAndSrtList :: Canon -> [Canon] -> Canon -- grpAndSrtList p iS | trace ("grpAndSrtList: (p = " ++ show p ++ ", iS = " ++ show iS ++ ")") False = undefined grpAndSrtList p iS = simpleHX cAddOpLevel $ map s $ cGetAddends $ factorSumIter True $ simpleHX cAddOpLevel $ map m' iS where m' a = applyFcnForHy a cMultOpLevel cFlattenAndGroup s a = applyFcnForHy a cMultOpLevel (srt p) -- sort the portions of the product so that the first item's base "derivative" will equal the product of the tail srt :: Canon -> Canon -> [Canon] srt p a' = ((reverse $ sortOn (nestLevel p) cs) ++ ncs) where (cs, ncs) = partition (\e -> candPred e && elem p (cGetBases e)) $ cGetFactors a' candPred :: Canon -> Bool candPred c'@(HX PoM _ _) = any candPred $ cGetHyperList c' candPred c' = powSq c' || cGetHyperOp c' == cTetrOpLevel where powSq (HX PoE (b:e:xs) _) = null xs && b == e -- some item raised to itself powSq _ = False -- 7 ^ {7 <^> {7 <<^>> (2^2) - 1}} would be lifted to Just 7 <<^>> 5. Useful when hyperizing tryLiftTail :: Canon -> Canon -- tryLiftTail c | trace ("tryLiftTail: (c = " ++ show c ++ ")") False = undefined tryLiftTail c | cGetHyperOp c < cExpOpLevel || length l < 2 || cGetHyperOp cLift == cExpOpLevel = impossibleHyperValue | otherwise = cLift where (l, b, cLift) = (cGetHyperList c, head l, cApplyHy cExpOpLevel [b, c] True) qcBase :: Canon -> Canon -> Canon -- qcBase p c | trace ("qcBase: (p = " ++ show p ++ ", c = " ++ show c ++ ")") False = undefined qcBase p c@(Bare n _) = if (n == cToI p) then c1 else error (errMsgQCB p c []) -- ToDo: Handle unfactored numbers qcBase p c@(Can g _) | cHyperExpr p || length pe == 0 = error (errMsgQCB p c []) | otherwise = snd $ head pe where pe = filter (\(p',_) -> p' == pI) g -- ToDo: Unfactored edge cases pI = cToI p qcBase p c | cGetHyperOp p == cAddOpLevel && p == c = c1 | length qce == 1 = snd $ head qce -- 3 <^ (3 <^> 4) -> 3 <^> 4 | otherwise = error $ errMsgQCB p c qce where qce = filter (\(b,_) -> p == b) $ map expPromote $ cGetFactors $ cQuasiCanonize c errMsgQCB :: Canon -> Canon -> [(Canon, Canon)] -> String errMsgQCB p c qce = "Logic error in qcBase: Canon: " ++ show c ++ " did not contain p = " ++ show p ++ ". Length = " ++ show (length qce) -- Checks how embedded a prime (or sum) is in a hyper expression nestLevel :: Canon -> Canon -> Int nestLevel p c | cBare p || cGetHyperOp p == cAddOpLevel = nL c 0 | otherwise = error $ "nestLevel: Only for a prime or sum: " ++ show p ++ " when checking: " ++ show c where nL c'@(HX PoA l _) lvl = if p == c' then (lvl + 1) else (getMax l lvl) nL (HX PoM l _) lvl = getMax l lvl nL (HX _ (b:_) _) lvl = nL b (lvl+1) nL (Bare n _) lvl = if cHyperExpr p || n /= pI then badLevel else lvl nL (Can g _) lvl | cHyperExpr p || null matchP = badLevel | otherwise = if x == 1 then lvl else (lvl + 1) where matchP = filter (\(p',_) -> p' == pI) g x = snd $ head matchP nL _ _ = error "Should not reach this point in the code" pI = cToI p badLevel = -1 -- no match getMax l lvl = if null l then 0 else (foldl1 max $ map (\m -> nL m lvl) l) -- simple utility function to create hyper expression if list specified simpleHX :: Canon -> [Canon] -> Canon simpleHX h c | ln == 1 = head c | ln > 1 = HX h c IntC | otherwise = if h == cAddOpLevel then c0 else c1 -- nullary values where ln = length c -- | If the Canon is a product, return the factors. Otherwise, return the Canon itself. cGetAddends :: Canon -> [Canon] cGetAddends c = if cGetHyperOp c == cAddOpLevel then cGetHyperList c else [c] -- | If the Canon is a sum, return the addends. Otherwise, return the Canon itself. cGetFactors :: Canon -> [Canon] cGetFactors c = if cGetHyperOp c == cMultOpLevel then cGetHyperList c else [c] -- | Take a canon and a list of indexes and delve into the canon. This operates on the internal hyper lists cDelve :: Canon -> [Int] -> Canon cDelve c xL | x < length hL = if null xs then e else cDelve e xs | otherwise = error $ "Can not find index " ++ show x ++ " in canon: " ++ show c where (hL, e) = (cGetHyperList c, hL !! x) (x:xs) = xL applyFcnForHy :: Canon -> Canon -> (Canon -> [Canon]) -> Canon applyFcnForHy c h f = if cGetHyperOp c == h then simpleHX h (f c) else c -- FactorSumIter: Performs steps like: a*b + b + c => b * (a+1)+c. factorSumIter :: Bool -> Canon -> Canon factorSumIter hF (HX PoA hL _) = fs (head hL) (tail hL) [] where fs a (x:xs) wL | gcd' == c1 || (hF && not (cHyperExpr gcd')) = fs a xs (wL ++ [x]) | otherwise = fs nv xs wL where (nv, gcd') = simpleFactorPair a x hF fs a _ wL = computeExpr cAddOpLevel (wL ++ [a]) factorSumIter _ c = c -- | This checks if the (hyper)expression is in quasi-canonical form cQuasiCanonized :: Canon -> Bool cQuasiCanonized (HX PoA _ _) = True cQuasiCanonized c@(HX PoM l _) = all cQuasiCanonized l && null (cGetBases' True True False c) cQuasiCanonized (HX PoE (b:_:xs) _) = (cBare b || cGetHyperOp b == cAddOpLevel) && null xs -- only b ^ e not b ^ e ^ x cQuasiCanonized (HX h _ _) = h > maxHyperOpDelveLevel -- anything else like tetration has not been simplified cQuasiCanonized _ = True -- | This is akin to canonical form except you may have sums in the bases. It converts expression up to a hyperoperational cutoff cQuasiCanonize :: Canon -> Canon -- cQuasiCanonize c | trace ("cQuasiCanonize: (c = " ++ show c ++ ")") False = undefined cQuasiCanonize c | cGetHyperOp c > maxHyperOpDelveLevel || (pF && null sM) -- nothing below the the hyper limit = c -- don't attempt to canonize | pF && not (null bM) -- there are entries beyond the hyper limit. = computeExpr cMultOpLevel ((cL sMp) ++ bM) | otherwise = computeExpr cMultOpLevel (cL c) -- all below the hyper limit where (bM, sM) = partition (\m -> cGetHyperOp m > maxHyperOpDelveLevel) $ cGetHyperList c -- partition product (sMp, pF) = (computeExpr cMultOpLevel sM, cGetHyperOp c == cMultOpLevel) -- "Endless" looping! cL c' = map (\l -> promote (fst $ head l, fst $ cConvertToSum $ sum $ map snd l)) $ cL c' = map (\l -> promote (fst $ head l, sum $ map snd l)) $ -- ToDo: Make this more robust? groupBy (\x y -> fst x == fst y) $ sortOn fst $ can' c' promote (b',e') | e' == c1 = b' | cHyperExprAny e' = (computeExpr cExpOpLevel [b',e']) | otherwise = b' <^ e' can' c'@(HX h l'@(b:xs) IntC) | h == cAddOpLevel = [(c', c1)] -- you need the base | h == cMultOpLevel = concat $ map can' l' | otherwise = map (\(b', e') -> (b', mul e')) $ can' b where expTail = cApplyHy h xs True eH = expRec h expTail mul :: Canon -> Canon mul e' | e' == c1 = eH | eH == c1 = e' | otherwise = applyFcnForHy (simpleHX cMultOpLevel [e', eH]) cMultOpLevel cFlattenAndGroup expRec h' e' | h' == cExpOpLevel = e' -- otherwise, recursively demote down | otherwise = expRec (h' - c1) newE -- e.g. x<^>y<^>z = x^(x<^>((y<^>z)-1)) where newE = cApplyHy h' [b, e' - c1] True can' (Can g _) = map (\(b', e') -> (makeCanon b', e')) g -- essentially the same thing can' c' = [(c', c1)] nonPrimePowerError :: Canon -> String -- collect all of the errors together nonPrimePowerError c = "cNestExpTail: Can't compute if base is a product or unfactored: " ++ show c cBaseRadical :: Canon -> Canon -- need a cutoff for what to process cBaseRadical = cBaseRadical' False False False cBaseRadical' :: Bool -> Bool -> Bool -> Canon -> Canon cBaseRadical' f d h c = product $ cGetBases' f d h c -- For cToGCR combine the bases or partial bases that match on a gcm (say with the other term) in cTryToCanonize cHyOpLvlOutOfRange :: Canon -> Bool cHyOpLvlOutOfRange h | h < cMultOpLevel || h > maxHyperOpDelveLevel = True | otherwise = False data FuncType = Mult | Gcd | Lcm deriving (Show, Eq) -- Several choices from multiplicative functions. "Let results escape" to cHyperOp only at end cMultiplicative :: Canon -> Canon -> FuncType -> [Canon] -- cMultiplicative v w t | trace ("cMultiplicative: Processing: v=("++show v++", w= "++show w++", t = "++show t++")") False = undefined cMultiplicative v w t | not (cHyperExpr v) && not (cHyperExpr w) = case t of Gcd -> [gvw, div v gvw, div w gvw] Lcm -> [cLCM v w] _ -> [v * w] | t == Mult = [v * w] -- No longer does anything distinct from multiplication | t == Lcm && relPrime = if v' == c1 then [abs w] else (if w' == c1 then [abs v] else [hyperize' $ cCleanup $ head $ cMultiplicative' vA wA t]) | t == Lcm = [hyperize' $ cCleanup $ head $ cMultiplicative' (cQuasiCanonize vA) (cQuasiCanonize wA) t] | t == Gcd && relPrime = [gHvw, f v' v, f w' w] | otherwise = [gHvw', f v2 v, f w2 w] where gvw = cGCD v w -- non-hyper (vA, wA) = (abs v, abs w) hyperize' c = simpleHX cMultOpLevel (concat $ map (\e -> cGetFactors $ if cQuasiCanonized e then cHyperize e else e) $ cGetFactors c) (gHvw:v':w':_) = map (hyperize' . cCleanup) $ cMultiplicative' vA wA Gcd -- first try (gHvw2:v2:w2:_) = map (hyperize' . cCleanup) $ cMultiplicative' (cQuasiCanonize v') (cQuasiCanonize w') Gcd gHvw' = gHvw * gHvw2 relPrime = null $ intersect (cGetBases v') (cGetBases w') f a' a = if signum a == cN1 then negate a' else a' -- efficient way to adjust by sign -- do not call directly. assumes unsigned "hyper" input cMultiplicative' :: Canon -> Canon -> FuncType -> [Canon] -- cMultiplicative' v w t | trace ("cMultiplicative': Processing: v=("++show v++", w= "++show w++", t = "++show t++")") False = undefined cMultiplicative' v w t = [apply r, apply xN, apply yN] where -- Internally we just manipulate lists rather than continually passing interim results to cApplyHy apply l = product $ (ccProd:hEs) where ccProd = product $ map conv cc where conv c | cHyperExpr c = (head hL) <^ (head $ tail hL) -- always of this form | otherwise = c where hL = cGetHyperList c (cc, hEs) = partition canConv l canConv c = not (cHyperExprAny c) || (cGetHyperOp c == cExpOpLevel && length hL == 2 && not (any cHyperExprAny hL)) where hL = cGetHyperList c (xN, yN, r) = p' (allFactors v) (allFactors w) -- this mode creates a hyper list p' x y -- We ignore the first two return values when not running for Gcd. ToDo: simplify the code below | eitherNull x1 y1 = if t == Gcd then (x1, y1, r1) else ([], [], r1 ++ x1 ++ y1) | eitherNull x2 y2 = if t == Gcd then (x2, y2, r2) else ([], [], r2 ++ x2 ++ y2) | eitherNull x3 y3 = if t == Gcd then (x3, y3, r3) else ([], [], r3 ++ x3 ++ y3) | otherwise = if t == Gcd then (x4, y4, r4) else ([], [], r4 ++ x4 ++ y4) where (x1, y1, r1) = rW x y [] False False (x2, y2, r2) = rW x1 y1 r1 True False (x3, y3, r3) = rW x2 y2 r2 False True (x4, y4, r4) = rW x3 y3 r3 True True -- rW x' y' g' bx by | trace ("rW: (x = "++show x'++",y="++show y'++", g' = "++show g'++",bx="++show bx++",by="++show by++")") False = undefined rW x' y' g' bx by = promote (fmt x') (fmt y') [] [] g' bx by t where fmt m = allFactors (simpleHX cMultOpLevel m) -- ToDo: change wrap kludge eitherNull j k = null j || null k -- promote a' b' aW' bW' g' aF bF _ | trace ("promote: (a' = "++show a'++", b' = "++show b'++", aW' = "++show aW'++", bW' = "++show bW'++", g' = "++show g'++", aF = "++show aF++", bF = "++show bF++")") False = undefined promote a' b' aW' bW' g' aF bF t' = if (aF && bF) then r' (sortOn getBase' a') (sortOn getBase' b') aW' bW' g' [] [] -- Loopage worries seem unfounded. else r' a' b' aW' bW' g' [] [] -- Doesn't seem necessary for both where -- r' aC bC aW bW g _ bN | trace ("r': (aC = "++show aC++", bC = "++show bC++", aW = "++show aW++", bW = "++show bW++", g = "++show g++", bN = "++show bN++")") False = undefined r' (a:as) (b:bs) aW bW g aN bN | aB == bB && (aT /= impossibleHyperValue && bT /= impossibleHyperValue) = r' as bs nAw nBw (g ++ [aB <^ m]) aN bN | otherwise = r' (a:as) bs aW bW g aN (bN ++ [b]) where (aB, bB) = (f' a aF, f' b bF) (aT, bT) = (e' a aF, e' b bF) m | t' == Gcd && aF && bF = if a < b then aT else bT -- leverage the comparision instead of using tail | t' == Gcd = min aT bT | t' == Lcm && aF && bF = if a > b then aT else bT -- leverage, pt. 2 | t' == Lcm = max aT bT | otherwise = aT + bT -- multiply, so add the exponents nGcdE w' | w' /= m && cMaxHyperOp w' > cTetrOpLevel && cMaxHyperOp m > cTetrOpLevel = simpleHX cAddOpLevel [w', negate m] -- to avoid infinite loops | otherwise = w' - m (nAw, nBw) | t' == Gcd = (aW ++ [aB <^ nGcdE aT], bW ++ [bB <^ nGcdE bT]) | otherwise = (aW, bW) r' (a:as) _ aW bW g aN bN = r' as bN newAw bW g newAn [] -- add to aW, use tail for #1 + feed in bN list for #2 where (newAw, newAn) | t' == Gcd = (aW ++ [a], aN) | otherwise = (aW, aN ++ [a]) -- not found r' _ b aW bW g aN bN -- at the end | t' == Gcd = (aW, bW ++ b ++ bN, g) -- add unprocessed members of b lists | otherwise = (aN, b ++ bN, g) -- feed the lists back in f' j fB = if fB then getBase' j else j e' j fB = if fB then (if cGetHyperOp j > maxHyperOpDelveLevel then impossibleHyperValue else cNestExpTail j False) -- ToDo: replace with cQuasiCanonize else c1 -- if whole expression, exp is just 1 getBase' s@(HX PoA _ _) = s getBase' c@(HX PoM _ _) = error $ getBaseHypError c getBase' (HX _ (b:_) _) = b getBase' (Bare b _) = makeCanon b getBase' c@(Can g _) | length g == 1 = gcrToC g --makeCanon $ fst $ head g | otherwise = error $ getBaseCanError c getBase' _ = error "getBase' unknown canon type encountered." getBaseCanError c = "getBase': canonical param with > 1 base not allowed: " ++ show c getBaseHypError c = "getBase': hyper prod param not allowed: " ++ show c -- Note: This does not delve into hyperoperations and try separate or group -- It does split canonical elements unlike cGetFactors allFactors :: Canon -> [Canon] -- allFactors c | trace ("allFactors: (c = " ++ show c ++ ")") False = undefined allFactors (HX PoM l _) = concat $ map allFactors l allFactors c@(HX _ _ _) = [c] allFactors c | cIntegral c = filter (/= c1) $ map expLift $ cToGCR c | otherwise = error "allFactors only takes integral canons" -- Note: We call the constructor directly so there are no demotions applied by cApplyHy where expLift (b', e') = if e' == c1 then bC else (simpleHX cExpOpLevel [bC,e']) -- avoids "2" bug where bC = makeCanon b' -- | Compute the rth-root of a Canon. cRoot :: Canon -> Canon -> Canon cRoot c r | not (cPositive r) = error "r-th roots are not allowed when r <= 0 or not integral" | r == c1 || c == c0 = c | cNegative c && cEven r = error "cRoot does not support imaginary numbers (even roots of negative numbers)." | all (\(_,e) -> cMod e r == 0) cL' = if cNegative c then negate root else root | cMaxHyperOp c > maxHyperOpDelveLevel = error $ "Root could not be found but that may be due to the level of hyper operation being beyond the cutoff: " ++ show c | otherwise = error $ "The root requested was not a multiple of all the exponents in the expansion of " ++ show c where cL' = map expPromote $ allFactors $ cQuasiCanonize $ abs c root = simpleHX cMultOpLevel $ map (\(p,e) -> expDemote (p, e / r)) cL' -- | This is used for tetration, etc. It defaults to zero for non-integral reps. cPrimeTowerLevel :: Canon -> Canon cPrimeTowerLevel (Bare _ Simp) = c1 cPrimeTowerLevel (Can g IntC) | gcrPrimePower g = cPrimeTowerLevelI (snd $ head g) (fst $ head g) (1 :: Integer) | otherwise = c0 cPrimeTowerLevel c@(HX h l@(b:xl) _) | h < cExpOpLevel || any cHyperExprAny l || not (cPrime b) = c0 -- ToDo: handle nested hyper expression cases properly | h == cExpOpLevel = if cQuasiCanonized c && cMaxHyperOp c > cExpOpLevel then (cPrimeTowerLevel $ cHyperize c) else (makeCanon $ toInteger $ length l) | h == cTetrOpLevel = simpleHX h xl | h <= maxHyperOpDelveLevel = cDelve (cQuasiCanonize c) [1,1] -- gets the tetration expression | otherwise = c -- it's so massive just return the number itself. Not that critical. cPrimeTowerLevel _ = c0 -- | Internal workhorse function to compute the height of a prime tower (e.g. 5^(5^7) => 3) cPrimeTowerLevelI :: Canon -> Integer -> Integer -> Canon cPrimeTowerLevelI (Bare b _) n l | b == n = makeCanon $ l + 1 | otherwise = c0 cPrimeTowerLevelI (Can g IntC) n l | gcrPrimePower g == False = c0 | n /= (fst $ head g) = c0 | otherwise = cPrimeTowerLevelI (snd $ head g) n (l+1) cPrimeTowerLevelI _ _ _ = 0 -- This is only for internal display. Not needed for hyper-expressions. -- | Function to convert Canon to generalized canon rep cToGCR :: Canon -> GCR_ --cToGCR c | trace ("cToGCR: (c = " ++ show c ++ ")") False = undefined -- Tracing here may cause stack overflow! cToGCR c = case gAtt of Just g -> g _ -> error $ noConvError where gAtt = cToGCR' c noConvError = "Could not to convert unwieldy hyper expression to canonical rep: " ++ show c ++ "." cToGCR' :: Canon -> Maybe GCR_ cToGCR' (Can x _) = Just x -- don't attempt to factor a large composite at this time. To Do: may not handle trial division. cToGCR' (Bare x s) | s == Simp || gcrBigComposite (x, c1) = Just $ if x == 1 then gcr1 else [(x, c1)] | otherwise = cToGCR' (makeCanon x) -- ToDo:Thread in CycloMap? cToGCR' c | any (\c' -> cHyperExpr c' && (cHyOpLvlOutOfRange $ cGetHyperOp c')) $ cGetFactors c = Nothing | any (not . cSimplified) f = Nothing -- if there are any unsimplified bases, fail | otherwise = Just g -- product of base / exponent pairs where f = cGetFactors $ cQuasiCanonize c g = sortOn fst $ concat $ map conv f conv (Bare n _) = [(n, c1)] conv (Can g' _) = g' conv (HX _ l _) = [(cToI (l !! 0), l !! 1)] -- already in canonized mode so this is safe -- Warning: Don't call this for 0 or +/- 1. The value type will not change by negating the value gcrNegateCanonical :: GCR_ -> CanonValueType -> Canon gcrNegateCanonical g v | gcrNegative g = case gcrPrime (tail g) of True -> Bare (fst $ head $ tail g) Simp False -> Can (tail g) v | otherwise = Can (gcreN1 : g) v -- just prepend gcrNegate :: GCR_ -> GCR_ gcrNegate Pg0 = gcr0 gcrNegate x | gcrNegative x = tail x | otherwise = gcreN1 : x gcrNegative :: GCR_ -> Bool gcrNegative PgNeg = True gcrNegative _ = False gcrPositive :: GCR_ -> Bool gcrPositive PNeg = False gcrPositive PZero = False gcrPositive _ = True gcrMult :: GCR_ -> GCR_ -> CycloMap -> (GCR_, CycloMap) gcrMult x POne m = (x, m) gcrMult POne y m = (y, m) gcrMult _ Pg0 m = (gcr0, m) gcrMult Pg0 _ m = (gcr0, m) gcrMult x@(xh@(xp,xe):xs) y@(yh@(yp,ye):ys) m = case compare xp yp of LT -> (xh:g, m') where (g, m') = gcrMult xs y m EQ -> if gcrNegative x || expSum == c0 then gcrMult xs ys m -- cancel double negs/exponents adding to zero else ((xp, expSum):gf, mf) where (expSum, m') = cAdd xe ye m (gf, mf) = gcrMult xs ys m' GT -> (yh:g, m') where (g, m') = gcrMult x ys m gcrMult x y _ = error e where e = "Non-exhaustive pattern error in gcrMult. Params: " ++ show x ++ "*" ++ show y gcr1, gcr0 :: GCR_ gcr1 = [] gcr0 = [(0, c1)] gcreN1 :: GCRE_ gcreN1 = (-1, c1) gcrToI :: GCR_ -> Integer -- ToDo: Add upperbound into conversaion step gcrToI g = product $ map f g where f (p, e) | ce > 0 = p ^ ce | otherwise = error negExpErr where ce = cToI e negExpErr = "gcrToI: Negative exponent found trying to convert " ++ show g gcrToD :: GCR_ -> Double gcrToD g = product $ map (\(p,e) -> (fromIntegral p) ** cToD e) g gcrCmp :: GCR_ -> GCR_ -> Ordering gcrCmp POne y = gcrCmpTo1 y True gcrCmp x POne = gcrCmpTo1 x False gcrCmp x y | x == y = EQ | xN && yN = compare (gcrToC $ tail y) (gcrToC $ tail x) | xN = LT | yN = GT | gcrIsZero x = LT | gcrIsZero y = GT | otherwise = case cSuperLogCmp (superLogCan x) (superLogCan y) of -- If equal: we have to break out the big guns and try to convert it to an integer EQ -> compare (gcrToI x) (gcrToI y) cmp -> cmp where (xN, yN) = (gcrNegative x, gcrNegative y) gcrIsZero Pg0 = True; gcrIsZero _ = False gcrCmpTo1 :: GCR_ -> Bool -> Ordering gcrCmpTo1 PNeg b = if b then GT else LT gcrCmpTo1 Pg0 b = if b then GT else LT gcrCmpTo1 _ b = if b then LT else GT -- | These internal functions should not be called directly. The definition of GCD and LCM -- are extended to handle non-Integers. gcrGCD, gcrLCM :: GCR_ -> GCR_ -> GCR_ gcrGCD Pg0 y = y gcrGCD x Pg0 = x gcrGCD x y | crFactCutoff > 0 && -- partial factorization workarounds can be disabled if <= 0 ((gcrIncompFact x && gcrUnsmooth y) || -- either has an imcomplete factorization and the other (gcrIncompFact y && gcrUnsmooth x)) = f spx spy -- in case of unfactored bits | otherwise = f x y where f Pg1 _ = gcr1 f _ Pg1 = gcr1 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 gcrLCM Pg0 _ = gcr0 gcrLCM _ Pg0 = gcr0 gcrLCM x y | gcrNegative x || gcrNegative y = f (gcrAbs x) (gcrAbs y) | crFactCutoff > 0 && -- partial factorization workarounds can be disabled if <= 0 ((gcrIncompFact x && gcrUnsmooth y) || -- either has an imcomplete factorization and the other (gcrIncompFact y && gcrUnsmooth x)) = f spx spy -- in case of unfactored bits | otherwise = f x y where f Pg1 y' = y' f x' Pg1 = 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) -- ToDo: This code as well as for gcd and lcm closely matches that in Internals.hs. Remove the duplicate code -- Note: This and related functions are only called when crFactCutoff > 0 (indicating a partial factorization is possible) spFactor :: GCR_ -> GCR_ -> (GCR_, GCR_) spFactor x y = (sx ++ (grp $ sortOn fst $ px ++ spF cx (py ++ cy)), sy ++ (grp $ sortOn fst $ py ++ spF cy (px ++ cx)) ) where spl n = (s, p, c) where (s, r) = partition gcrSmooth n (p, c) = partition (\ce -> not $ gcrBigComposite 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 gcrUnsmooth, gcrIncompFact :: GCR_ -> Bool gcrUnsmooth = any (\ce -> not $ gcrSmooth ce) gcrIncompFact = any gcrBigComposite gcrBigComposite, gcrSmooth :: GCRE_ -> Bool gcrSmooth (b,_) = b <= crSmallFactCutoff gcrBigComposite (b,_) = bigComposite b bigComposite :: Integer -> Bool bigComposite b = crFactCutoff > 0 && b > crFactCutoff && (not $ isPrime b) {- Use cSuperLog instead gcrLogDouble :: GCR_ -> Double gcrLogDouble g = sum $ map (\(p,e) -> (log $ fromIntegral p) * (cToD e)) g -} divisionError :: String divisionError = "gcrDiv: As requested per param, the dividend must be a multiple of the divisor." divByZeroError :: String divByZeroError = "gcrDiv: Division by zero error!" zeroDivZeroError :: String zeroDivZeroError = "gcrDiv: Zero divided by zero is undefined!" gcrDivStrict :: GCR_ -> GCR_ -> GCR_ gcrDivStrict x y = case gcrDiv x y of Left errorMsg -> error errorMsg Right results -> results gcrDiv :: GCR_ -> GCR_ -> Either String GCR_ gcrDiv Pg0 Pg0 = Left zeroDivZeroError gcrDiv Pg0 _ = Right gcr0 gcrDiv _ Pg0 = Left divByZeroError gcrDiv n d = g' n d where g' x POne = Right x g' POne _ = Left divisionError g' x y | gcrNegative y = g' (gcrNegate x) (gcrAbs y) | otherwise = case compare xp yp of LT -> case (g' xs y) of Left _ -> Left divisionError Right res -> Right ((xp, xe) : res) EQ | xe > ye -> case (g' xs ys) of Left _ -> Left divisionError Right res -> Right ((xp, xe - ye) : res) EQ | xe == ye -> gcrDiv xs ys _ -> Left divisionError where ((xp,xe):xs) = x ((yp,ye):ys) = y -- GCR functions (GCR is an acronym for generalized canonical representation) gcrAbs :: GCR_ -> GCR_ gcrAbs x | gcrNegative x = tail x | otherwise = x gcrToCR :: GCR_ -> CR_ gcrToCR c = map (\(p,e) -> (p, cToI e)) c gcrBare :: GCR_ -> Bool gcrBare PBare = True gcrBare POne = True gcrBare _ = False gcrPrime :: GCR_ -> Bool gcrPrime PgPrime = True gcrPrime _ = False gcrPrimePower :: GCR_ -> Bool gcrPrimePower PgPPower = True gcrPrimePower _ = False gcrOdd, gcrEven :: GCR_ -> Bool gcrOdd Pg0 = False gcrOdd POne = True gcrOdd c | gcrNegative c = gcrOdd (gcrAbs c) | otherwise = cp /= 2 where (cp,_):_ = c gcrEven g = not (gcrOdd g) gcrEqCheck :: GCR_ -> GCR_ -> Bool gcrEqCheck POne POne = True gcrEqCheck POne _ = False gcrEqCheck _ POne = False gcrEqCheck ((xp,xe):xs) ((yp,ye):ys) | xp /= yp || xe /= ye = False | otherwise = gcrEqCheck xs ys gcrEqCheck x y = error e where e = "Non-exhaustive patterns in gcrEqCheck comparing " ++ show x ++ " to " ++ show y gcrSimplified :: GCR_ -> Bool gcrSimplified = all (\(b,e) -> cSimplified e && check b) where check n = crFactCutoff <= 0 || n < crFactCutoff || (n > crFactCutoff && isPrime n) -- | Return the base b from a Canon Element (equivalent to b^e) cGetBase :: CanonElement -> Canon cGetBase (b, _) = b -- | Return the exponent e from a Canon Element (equivalent to b^e) cGetExponent :: CanonElement -> Canon cGetExponent (_, e) = e -- | Return the list of bases from a Canon (conceptually of the form [b^e])> cGetBases :: Canon -> [Canon] cGetBases = cGetBases' False False False -- don't check if in range cGetBases' :: Bool -> Bool -> Bool -> Canon -> [Canon] cGetBases' f d h c -- if f flag True, only keep the "True" matches based on allowed hyper op level range -- if d flag True, only return the bases that occur more than once (for grouping later) -- if h flag True, only include the "hyper bases" = if d then rMultiple else nub r' where rMultiple = map fst $ filter (\(_,ct) -> ct /= 1) $ map (\l -> (head l, length l)) $ group $ sort r' r' = map fst $ if f then (filter (\(_, f') -> f') r) else r r = g' (abs c) h g' b@(Bare _ _) h' = if h' || b == c1 then [] else [(b,True)] g' (Can g _) h' = map (\c' -> (cGetBase $ convGCREToCE c', True)) g2 where g2 = if not h' then g else (filter (\(_,e) -> cHyperExprAny e) g) g' c'@(HX PoA _ _) _ = [(c', not $ (f && (cHyOpLvlOutOfRange $ cGetHyperOp c')))] -- return sum itself, only option. g' (HX PoM cL _) h' = concat $ (g' (head cL) h'):(map (\e -> g' e False) (tail cL)) -- ToDo: weave in mult check even though in range g' (HX y cL _) _ = map (\(c',_) -> (c',pF)) $ g' (head cL) False where pF = not $ f && cHyOpLvlOutOfRange y -- If true, we drill down -- e.g 3^4 or 3<^>4. First member of list is 3. Could also be a composite -- | Similar to cGetBases except that it will do trial factoring of any hyper sums. So, for obvious reasons, this is not a complete factorization. cGetBasesDeep :: Canon -> [Canon] cGetBasesDeep c@(HX PoA l _) = sort $ nub ((c':i) ++ sF) where i = foldl1 intersect $ map cGetBasesDeep l iP = product i c' = simpleHX cAddOpLevel (map (\a -> div a iP) l) sF = filter (\p -> mod c p == 0) $ smallPrimeCanons cGetBasesDeep c = cGetBases c -- | Return the list of exponents from a Canon (conceptually of the form [b^e]). cGetExponents :: Canon -> [Canon] cGetExponents (Bare _ _) = [c1] -- always just one cGetExponents (Can g _) = map (cGetExponent . convGCREToCE) g cGetExponents (HX PoA _ _) = [c1] -- this is a sum so the exponent is just one cGetExponents hx@(HX _ _ _) = map (cGetExponent . convGCREToCE) $ cToGCR hx -- | Return the list of CanonElements from a Canon (conceptually of the form [b^e]). cGetElements:: Canon -> [CanonElement] cGetElements b@(Bare _ _) = [(b, c1)] cGetElements (Can g _) = map convGCREToCE g cGetElements hx@(HX _ _ _) = map convGCREToCE $ cToGCR hx -- | Convert a generalized canon rep element to a CanonElement convGCREToCE :: GCRE_ -> CanonElement convGCREToCE (b, e) = (makeCanon b, e) -- ToDo: Optimize as b is already known to be a prime here -- | Divisor functions should be called with integral Canons. Restricted to positive divisors. Returns Either String Canon cNumDivisors, cTau :: Canon -> Canon cNumDivisors c = case (cNumDivisors' c) of Left s -> error s Right v -> v cTau = cNumDivisors -- | Underlying divisor function that can return value or (error) message. cNumDivisors' :: Canon -> Either String Canon cNumDivisors' c | not (canComputeDivs c') = Left $ "Canon was zero, not integral or not completely factored. Can't compute." | otherwise = case cToGCR' c' of Just g -> Right $ product $ map (\(_,e) -> 1 + e) g Nothing -> Left "Unknown issue converting to underlying GCR structure." where c' = abs c -- | Compute the nth divisor of a Canon. It operates on the absolute value of the Canon and is zero based. -- Note: This is deterministic but it's not ordered by the value of the divisor. cNthDivisor :: Canon -> Canon -> Either String Canon cNthDivisor _ (Bare _ NSim) = Left "cNthDivisor: Bare integer has not been simplified." cNthDivisor n c | cNegative n || not (cIntegral n) = Left "cNthDivisor: Both n must be integral and >= 0" | not (canComputeDivs aC) = Left "cNthDivisor: Canon was zero, not integral or not completely factored. Can't compute" | cHyperExpr c = nthHyper $ cToGCR aC | otherwise = nth aC where aC = abs c nth Pc0 = Right n -- Zero has an infinite set of divisors. The nth divisor is just n as a Canon nth cn = case f (cAbs n) (cToGCR cn) of Right r -> Right $ gcrToC r Left e -> Left e where f Pc0 _ = Right gcr1 f _ Pg1 = Left "cNthDivisor: Bad dividend number requested." f n' c' = case f (div n' (e + 1)) cs of -- First param is the next n Right r -> Right $ if m == c0 then r else ((b,m):r) e' -> e' -- Return error message where (b,e):cs = c' m = mod n' (e + 1) nthHyper :: GCR_ -> Either String Canon nthHyper x | n >= nd = Left "cNthDivisor: Bad dividend number requested of hyper expression." | otherwise = Right $ compute n divL c1 where nd = product $ map snd divL -- Number of dividends divL = map (\(p,e) -> (Bare p Simp, e + 1)) x compute n' ((p,d):ds) wP | nq == 0 = newWp -- at 0, no need to continue | otherwise = compute nq ds newWp where (nq, m) | n' < d = (c0, n') | cHyperExpr d && n' < 100 * d = proc n' (0 :: Int) | otherwise = (quotRem n' d) -- limited -- ToDo: make it smarter newWp = if m == c0 then wP else (wP * p <^ m) proc n2 rc | n2 < d = (makeCanon $ toInteger rc, n2) | otherwise = proc (n2 - d) (rc + 1) compute _ _ wP = wP -- | Consider this to be the inverse of the cNthDivisor function. This function ignores signs -- but both parameters must be integral. cWhichDivisor :: Canon -> Canon -> Either String Canon cWhichDivisor d c | not (cIntegral d) = Left "cWhichDivisor: d must be integral" | d == c0 = Left "cWhichDivisor: Zero is not a proper divisor" | d == c1 && aC /= c0 = Right c0 | not (canComputeDivs aC) = Left "cWhichDivisor: Canon was either zero or not completely factored. Can't compute" | d == c = Right $ cNumDivisors c - 1 | otherwise = case gcrDiv gAc gAd of Right _ -> Right $ wD gAd divLProg c0 -- valid divisor so it's safe to compute Left _ -> Left "cWhichDivisor: d is not a divisor of c" where (aD, aC, gAd, gAc) = (abs d, abs c, cToGCR aD, cToGCR aC) divLProg = zip (map fst gAc) (c1:(init $ scanl1 (*) $ map (\(_,e) -> e + 1) gAc)) wD a@((aP,aE):as) ((bP,bC):bs) wS | aP < bP = error "Logic error: Unexpected factors found in divisor" | aP > bP = wD a bs wS | otherwise = wD as bs (wS + aE * bC) wD a _ wS | null a = wS | otherwise = error "Logic error: Unexpected factors found in divisor" -- | Efficiently compute all of the divisors based on the canonical representation. -- | Returns Either an error message or a list of Canons. cDivisors :: Canon -> Either String [Canon] cDivisors c | not (canComputeDivs c') = Left "cWhichDivisor: Canon was either zero or not completely factored. Can't compute" | otherwise = divs c' where c' = abs c divs (Bare x _) | x == 1 = Right [c1] | otherwise = Right [c1, makeCanon x] divs (Can g _) = Right $ map gcrToC $ foldr1 cartProd $ map pwrDivList g where cartProd xs ys = [x ++ y | y <- ys, x <- xs] pwrDivList (n,e) = [if y == 0 then gcr1 else [(n, makeCanon y)] | y <- [0 .. cToI e]] divs _ = error "cDivisors can't return all of the divisors for hyper expressions!" -- | Return the first N divisors of a hyper expression (if possible) cGetFirstNDivisors :: Int -> Canon -> Either String [Canon] cGetFirstNDivisors n c@(HX _ _ IntC) | any errPred divList = Left "cGetNDivisors: Canon was either zero or not completely factored. Can't compute" | otherwise = Right $ map clean divList where divList = map (\i -> cNthDivisor (makeCanon $ toInteger i) c) [0..n-1] errPred (Left _) = True errPred (Right _) = False clean (Left _) = error "'Dirty list' passed to clean function" clean (Right v) = v cGetFirstNDivisors n c = case cDivisors c of Right ds -> Right $ take n ds msg -> msg -- Assumes unsigned input canComputeDivs :: Canon -> Bool canComputeDivs c | cBare c && (cToI c == 0) = False | not (cSimplified c) || not (cIntegral c) = False | not (cHyperExpr c) = True | cHyperSum c = False | cGetHyperOp c > maxHyperOpDelveLevel = False | cHyperProd c && not (all canComputeDivs cL) = False | otherwise = canComputeDivs b where cL@(b:_) = cGetHyperList c smallPrimeCanons :: [Canon] smallPrimeCanons = map (\p -> Bare p Simp) $ take 1000 primes -- | This will determine if two arbitary expressions are relatively prime or not (if possible). Goes deep. cRelativelyPrime :: Canon -> Canon -> Maybe Bool cRelativelyPrime x y | x == c1 || y == c1 = Just True | cEven x && cEven y = Just False | not (null iBs) = Just False -- Intersection between the bases | xHs == c1 && yHs == c1 = Just True | xHs /= c1 && f xHs yB = Just False | yHs /= c1 && f yHs xB = Just False | xHs /= c1 && yHs /= c1 = if hpi then (Just False) else Nothing | otherwise = Nothing where (xB, yB) = (cGetBases x, cGetBases y) iBs = intersect xB yB hs b = product $ filter cHyperSum $ b (xHs, yHs) = (hs xB, hs yB) f x' l = any (\p -> mod x' p == 0) l -- does anything in l divide the hyperexpression x hpi = any (\p -> mod x p == 0 && mod y p == 0) $ smallPrimeCanons -- hyper prod intersection -- Super Log Logic: -- This section allows for comparison of hyperoperations by converting to a super log -- For instance: 9 ^ (8 ^ 7) = 10 ** 10 ** 6.301288668477042. We represent the super log as (2, 6.301288668477042) -- | Used when computing "Super Log" type SuperPairC = (Canon, Double) type SuperPairI = (Integer, Double) -- Trick to handle floating point issues (https://stackoverflow.com/questions/2354707/in-haskell-is-there-infinity-num-a-a) infinity :: Double infinity = read "Infinity" integerDblCutoff, manualLB :: Integer integerDblCutoff = 2 ^ ((snd $ floatRange (1.0 :: Double)) - 1) -- integral upper bound for converting to doubles. -- Note: I reduced this by one to be safe. When testing in 8.4.3, I noticed some inconsistent behavior close to 2^1024 manualLB = lBI ^ (16 :: Int) -- manually walk down number to this lower bound to get significant 'digits' dblLogMantissaUB :: Double dblLogMantissaUB = logD $ fromInteger integerDblCutoff -- we can safely compute lB ** lB ** doubleLogMantissaCutoff -- This is the cutoff for converting numbers to integers superLogICutoff :: SuperPairI superLogICutoff = (2, 5.0) -- 10^10^5 cSuperLogCutoff :: SuperPairC cSuperLogCutoff = (makeCanon $ fst superLogICutoff, snd superLogICutoff) -- 10^10^5 -- | This is the super or iterated log function. A level and mantissa is returned along with the number's sign. cSuperLog :: Canon -> (SuperPairC, Int) -- cSuperLog c | trace ("cSuperLog: Processing: (c=" ++ show c ++ ")") False = undefined cSuperLog (Bare n _) = (makeSuperLogC $ superLogI n, if n > 0 then 1 else (-1)) cSuperLog c@(Can g _) = if d == c1 then (superLogCan g, s) -- ToDo: if contains a hyperexpr, convert all? else (slDiv (superLogCan $ cToGCR n) (superLogCan $ cToGCR d), s) where (n, d) = cSplit c s = if cPositive c then 1 else (-1) cSuperLog (HX PoA s _) = cSuperLogSum s cSuperLog (HX PoM p _) = cSuperLogProd p cSuperLog (HX PoE e _) = (cSuperLogExp e, 1) -- ToDo: always positive? -- beyond exponentiation, get the tower height from the tail and adjust by offset cSuperLog c@(HX h (b:cs) _) | h > maxHyperOpDelveLevel = error $ "Unable to take cSuperLog of massive hyper expression: " ++ show c | h == cTetrOpLevel = ((sv1 + offset, m), 1) -- in case the cNestExpTail is not a hyper expr. | otherwise = ((c1 + offset + (head $ tail $ cGetHyperList $ cNestExpTail c False), m), 1) where (offset, m) = getTowerMantissa b sv1 sv1 = cApplyHy h cs True cSuperLog _ = error "Logic error in Super Log: Default Canon configuration unexpectedly reached" -- This implemented in the comparison logic. This reveals that instead of comparing (3 <<^>> 7) to (5 <<^>> 6) with super log -- we can compare superlog of 3 <<^>> 3 to 5 <<^>> 2 -- How the function would look: cSuperLogIter x n = if (n < 2) then (cSuperLog x) else (cSuperLog $ fst $ fst $ cSuperLogIter x (n-1)) superLogCan :: GCR_ -> SuperPairC -- don't really need to keep track of the sign superLogCan g = fst $ slProd $ map f g -- convert gcr to powerTower of sls for each base and then their products where f (b, e) = if e == 1 then cSuperLog (bFmt b) else (cSuperLogExp ((bFmt b):[e]), if b > 0 then 1 else (-1)) bFmt b = (Bare b Simp) -- for a power tower, the mantissa doesn't change (or more than negiligibly) after this height: 2 <^> x is the most "volatile". towerHeight :: Integer towerHeight = 8 towerHeightC :: Canon towerHeightC = makeCanon towerHeight -- get the tower mantissa (allowing from something shorter than the "towerHeight" above getTowerMantissa :: Canon -> Canon -> (Canon, Double) getTowerMantissa b h = (lvl - makeCanon htu, m) where (lvl, m) = slExp $ replicate (fromInteger htu) slb -- replicate the superlog "hgt" times slb = fst $ cSuperLog b htu = if h < towerHeightC then (cToI h) else towerHeight integerShowCutoff :: Integer -- numbers larger than this won't be fully printed by default integerShowCutoff = (10 :: Integer) ^ (1000 :: Integer) superLogI :: Integer -> SuperPairI superLogI i | i > 0 = suL i | i < 0 = suL $ abs i | otherwise = (0, 0.0) -- integral superlog rep. for zero where suL n | n > integerDblCutoff = refine $ spLog' (1, mC 0 n) -- manually compute above: save the last few significant digits | otherwise = refine $ spLog' (0, iToD) where spLog' (lvl', n') = if ln > lB then (spLog' (lvl' + 1, ln)) else (lvl', n') where ln = logB n' iToD | d == infinity = error "Raise bug: Number below cutoff still returning infinity." | otherwise = d d = fromIntegral n mC :: Integer -> Integer -> Double mC c n' | n' <= manualLB = fromIntegral c + (logB $ fromIntegral n') | otherwise = mC (c+1) (div n' lBI) refine (l,v) = if v > lB -- similar to slRefine then refine (l+1, logB v) else (if (v < 1.0 || (v == 1 && l > 0)) then (l-1, powB v) else (l, v)) makeSuperLogC :: SuperPairI -> SuperPairC -- Promote to a Canon makeSuperLogC (spi, d) = (makeCanon spi, d) cSuperLogProd, cSuperLogSum :: [Canon] -> (SuperPairC, Int) cSuperLogProd cL = slProd $ map cSuperLog cL cSuperLogSum cL = slSum $ map cSuperLog cL slEpsilon :: Double slEpsilon = 2.0 * 10.0 ** (-15.0) -- necessary for testing equality. Some rounding error can occur slProd, slSum :: [(SuperPairC, Int)] -> (SuperPairC, Int) slProd sL | null sL = error "Null list passed to slProd" -- ((c0, 1.0), 1) -- rep for nullary product aka 1 | otherwise = foldl1 (\(s, b) (ws, wb) -> (slMult s ws, b * wb)) sL slSum l | length posCL == 0 = (negSL, -1) -- all negative | length negCL == 0 = (posSL, 1) -- all positive | pl == nl && pm == nm = ((c0, 0.0), 0) -- superlog representation for zero | (pl == nl && pm > nm) || pl > nl = (cmpSl posSL negSL, 1) | otherwise = (cmpSl negSL posSL, -1) where (posCL, negCL) = partition (\(_, s) -> s == 1) l -- partition by positive then negative posSL = foldl1 slAdd $ map fst posCL negSL = foldl1 slAdd $ map fst negCL (pl,pm) = posSL (nl,nm) = negSL cmpSl (l1, m1) (l2, m2) = case (l1 - l2) of 0 -> case l1 of 0 -> (l1, m1-m2) -- literally m1 - m2 1 -> slRefine (l1, logB $ powB m1 - powB m2) -- lB^3 - lB^2 _ | l1 == c2 && m1 < dblLogMantissaUB -> -- lB^lB^2 - lB^lB^1.5 slRefine (l1, logD $ powD m1 - powD m2) _ | otherwise -> (l1, m1) -- ToDo: This is not accurate if l1 = 2 and m1 - m2 is v. small -- This needs a more robust solution 1 -> case l1 of 1 -> slRefine (l1, logB $ powB m1 - m2) -- lB^2 - 3 _ | l1 == c2 && m1 < dblLogMantissaUB -> -- lB^lB^2 - lB^3 slRefine (l1, logD $ powD m1 - powB m2) _ | otherwise -> (l1, m1) _ | l1 == c2 && m1 < dblLogMantissaUB -> -- lB^lB^2 - 3 slRefine (l1, logD $ powD m1 - m2) _ | otherwise -> (l1, m1) slAdd, slMult, slDiv :: SuperPairC -> SuperPairC -> SuperPairC slAdd e1 e2 = if cSuperLogGT e1 e2 then (add2 e1 e2) else (add2 e2 e1) -- larger first where add2 (l1, m1) (l2, m2) = case l1 - l2 of 0 -> case l1 of 0 -> slRefine (l1, m1 + m2) -- just m1 + m2 1 -> slRefine (l1, logB $ powB m1 + powB m2) -- lB^3 + lB^2 -- lB^lB^3 + lB^lB^2 and take double logB if in range _ | l1 == c2 && m1 < dblLogMantissaUB && tot < integerDblCutoff -> slRefine (l1, logD $ m1pd + m2pd) where (m1pd, m2pd) = (powD m1, powD m2) tot = round m1pd + round m2pd + 1 -- insurance _ | otherwise -> (l1, m1) 1 -> case l1 of 1 -> slRefine (l1, logB $ powB m1 + m2) -- lB^2 + 3 and take log -- lB^lB^2 + lB^3 and take double log if in range _ | l1 == c2 && m1 < dblLogMantissaUB -> slRefine (l1, logD $ powD m1 + powB m2) _ | otherwise -> (l1, m1) _ | l1 == c2 && m1 < dblLogMantissaUB -> slRefine (l1, logD $ powD m1 + powB m2) -- lB^lB^2 + 3 and dbl logB _ | otherwise -> (l1, m1) slMult e1 e2 = if cSuperLogGT e1 e2 then (mul2 e1 e2) else (mul2 e2 e1) -- larger first where mul2 (l1, m1) (l2, m2) = case l1 - l2 of 0 -> case l1 of 0 -> slRefine (l1, m1 * m2) -- 2 * 3 1 -> slRefine (l1, m1 + m2) -- lB^3 * lB^2 = lB ^ (3+2) -- lB^lB^3 * lB^lB^2 = lB ^(lB^3 + lB^2) and take logB 2 -> slRefine (l1, logB $ powB m1 + powB m2) -- lB^lB^lB^3 * lB^lB^lB^2 = lB ^(lB^lB^3 + lB^lB^2) and take logD _ | l1 == c3 && m1 < dblLogMantissaUB && tot < integerDblCutoff -> slRefine (l1, logD $ m1pd + m2pd) where (m1pd, m2pd) = (powD m1, powD m2) tot = round m1pd + round m2pd + 1 -- insurance _ | otherwise -> (l1, m1) 1 -> case l1 of 1 -> slRefine (l1, m1 + logB m2) -- lB^2 * 3 = lB^(2 + logB 3) -- lB^lB^2 * lB^3 = lB^(lB^2 + 3) and take logB 2 -> slRefine (l1, logB $ powB m1 + m2) -- lB^lB^lB^2 * lB^lB^3 = lB^(lB^lB^2 + lB^3) and take dbl logB if ... _ | l1 == c3 && m1 < dblLogMantissaUB -> slRefine (l1, logD $ powD m1 + powB m2) _ | otherwise -> (l1, m1) 2 -> case l1 of -- lB^lB^2 * 3 = lB^(lB^2 + logB 3) and take logB 2 -> slRefine (l1, logB $ powB m1 + logB m2) -- lB^lB^lB^2 * lB^3 = lB^(lB^lB^2 + 3) and take dbl log if ... _ | l1 == c3 && m1 < dblLogMantissaUB -> slRefine (l1, logD $ powD m1 + m2) _ | otherwise -> (l1, m1) _ | l1 == c3 && m1 < dblLogMantissaUB -> slRefine (l1, logD $ powD m1 + logB m2) _ | otherwise -> (l1, m1) slDiv a@(l1', m1') b@(l2', m2') | l1' < 0 && l2' < 0 = slDiv (slInvert b) (slInvert a) | l2' < 0 = slMult a (slInvert b) | l1' > l2' || (l1' == l2' && m1' > m2') = slRefine $ slDiv' a b | otherwise = slInvert $ slRefine $ slDiv b a where slInvert (l, m) = (negate l, m) slDiv' (l1, m1) (l2, m2) -- The 1st is > than the 2nd at this point and both levels are non-negative. | l1 == 0 = (c0, m1 / m2) -- simple division | l1 == 1 = case l2 of 0 -> (l1, m1 - logB m2) -- lB ^ (m1 - logB m2) _ -> (l1, m1 - m2) -- lB ^ (m1 - m2) | l1 == 2 = case l2 of 0 -> (l1, logB (powB m1 - logB m2)) -- lB ^ (lB^m1 - logB m2) 1 -> (l1, logB (powB m1 - m2)) -- lB ^ (lB^m1 - m2) _ -> (l1, logB (powB m1 - powB m2)) -- lB ^ (lB^m1 - lB^m2) | l1 == 3 && m1 < dblLogMantissaUB = case l2 of 0 -> (l1, logD (powD m1 - logB m2)) -- lB ^ (lB^lB^m1 - logB m2) 1 -> (l1, logD (powD m1 - m2)) -- lB ^ (lB^lB^m1 - m2) 2 -> (l1, logD (powD m1 - powB m2)) -- lB ^ (lB^lB^m1 - lB^m2) _ -> (l1, logD (powD m1 - powD m2)) -- lB ^ (lB^lB^m1 - lB^lB^m2) | otherwise = (l1, m1) -- Exp short for exponent. Interpret as a power tower (e.g. 2 ^ (3 ^ (5 ^ 7))) cSuperLogExp :: [Canon] -> SuperPairC cSuperLogExp cL = slExp $ map (fst . cSuperLog) cL -- ignore the sign returned. slExp :: [SuperPairC] -> SuperPairC slExp (x:r@(_:_)) = combineS x (slExp r) -- ToDo: change this to a fold where combineS (nl, nm) (rl, rm) -- Note in the comments below, ^ is a bit more legible than ** and internally, doubles are used = case rl of 0 -> case nl of 0 -> slRefine (nl, nm ** rm) -- 2 ^ 3 1 -> slRefine (nl, nm * rm) -- (lB^2) ^ 3 = lB^(2 * 3) -- (lB^lB^2) ^ 3 = lB^(3 * lB^2) = lB^lB^(2 + logB 3) 2 -> slRefine (nl, nm + logB rm) -- (lB^lB^lB^2) ^ 3 = lB^(3 * lB^lB^2) = -- lB ^ (lB^logB 3 * lB^lB^2) = lB^lB^(lB^2 + logB 3) then take the log 3 -> slRefine (nl, logB $ powB nm + logB rm) -- at this point the remainder has minimal to no effect _ | nl == 4 && nm < dblLogMantissaUB -> slRefine (nl, logD $ powD nm + logB rm) | otherwise -> (nl, nm) -- rm has no effect 1 -> case nl of -- 2 ^ (lB^3) = lB ^ (logB 2 * (lB ^ 3)) -- then progression from 0 to 2 0 -> slRefine (rl, logB nm * powB rm) 1 -> slRefine (rl, nm * powB rm) 2 -> slRefine (rl, powB nm * powB rm) -- (lB^lB^lB^2) ^ (lB^3) = lB ^ (lB^lB^2 * lB^3) = -- lB^lB^(lB^2 + 3) = then take double log 3 -> slRefine (nl, logB $ powB nm + rm) -- at this point the remainder has minimal to no effect. Note: We go down a level _ | nl == 4 && nm < dblLogMantissaUB -> slRefine (nl,logD $ powD nm + rm) | otherwise -> (nl, nm) -- rm has no effect 2 -> case nl of -- 2 ^ (lB^lB^3) = lB ^ (log2 * lB^lB^3) = lB ^ (lB^logLog2 * lB^lB^3) = -- lB ^ lB ^ (logLog2 + lB^3). Then, there's a progression 0 -> slRefine (rl, logD nm + powB rm) 1 -> slRefine (rl, logB nm + powB rm) 2 -> slRefine (rl, nm + powB rm) 3 -> slRefine (rl, powB nm + powB rm) -- at this point the remainder has minimal to no effect. Note: We go down 2 levels -- (lB^lB^lB^lB^2) ^ (lB^lB^3) = lB ^ (lB^lB^lB^2 * lB^lB^3) = -- lB^lB^(lB^lB^2 + lB^3) then take double log and reduce a level _ | nl == 4 && nm < dblLogMantissaUB -> slRefine (rl, powD nm + powB rm) | otherwise -> (nl, nm) -- rm has no effect _ | rl == 3 && rm < dblLogMantissaUB -> -- 2 ^ (lB^lB^lB^3) = lB^(log2 * lB^lB^lB^3) = -- lB^(lB^logD 2 * lB^lB^lB^3) = lB^lB^(logD 2 + lB^lB^3) -- Then we use similar logic for the progression of nl from 0 to 3 case nl of -- note the progression 0 -> slRefine (rl, logB $ logD nm + powD rm) 1 -> slRefine (rl, logB $ logB nm + powD rm) 2 -> slRefine (rl, logB $ nm + powD rm) 3 -> slRefine (rl, logB $ powB nm + powD rm) 4 -> if nm < dblLogMantissaUB && tot < integerDblCutoff then slRefine (rl, logB $ nmpd + rmpd) else (nl, max nm rm) where (nmpd, rmpd) = (powD nm, powD rm) tot = round nmpd + round rmpd + 1 -- insurance _ -> (nl, nm) -- just use nl | otherwise -> -- beyond we just compare the current and remainder case compare nl (rl + c1) of EQ -> (nl, max nm rm) GT -> (nl, nm) LT -> (rl + c1, rm) slExp (x:_) = slRefine x slExp _ = (0, 1.0) -- nullary product / poewr tower slRefine :: SuperPairC -> SuperPairC slRefine (lvl, v) = if v > lB then slRefine (lvl + c1, logB v) else (if v <= 1.0 then (lvl - c1, powB v) else (lvl, v)) -- | Compare the level and the "mantissa" cSuperLogCmp :: SuperPairC -> SuperPairC -> Ordering cSuperLogCmp (l1, m1) (l2, m2) | l1 > l2 = GT | l1 < l2 = LT | abs (m1 - m2) < slEpsilon = EQ | otherwise = if l1 >= 0 then (compare m1 m2) else (compare m2 m1) cSuperLogGT :: SuperPairC -> SuperPairC -> Bool cSuperLogGT x y = case cSuperLogCmp x y of GT -> True _ -> False lBI :: Integer lBI = 10 lB :: Double lB = fromIntegral lBI logB, logD, powB, powD :: Double -> Double logB a = log a / log lB logD a = logB $ logB a powB a = lB ** a powD a = lB ** (lB ** a) -- | Instance of CanonConv class instance CanonConv Canon where toSC c = toSC $ cToCR c toRC c = toRC $ cToCR c -- Canon of form x^1. (Does not match on 1 itself) pattern PBare :: forall t. [(t, Canon)] pattern PBare <- [(_, Bare 1 _)] -- Canon of form p^e where e >= 1. p has already been verified to be prime. pattern PgPPower :: forall t a. (Num a, Ord a) => [(a, t)] pattern PgPPower <- [(compare 1 -> LT, _ )] -- Canon of form p^1 where p is prime pattern PgPrime :: forall a. (Num a, Ord a) => [(a, Canon)] pattern PgPrime <- [(compare 1 -> LT, Bare 1 _)] -- Pattern looks for Canons beginning with negative 1 pattern PgNeg :: forall a. (Num a, Eq a) => [(a, Canon)] pattern PgNeg <- ((-1, Bare 1 _):_) -- Pattern for "generalized" zero pattern Pg0 :: forall a. (Num a, Eq a) => [(a, Canon)] pattern Pg0 <- [(0, Bare 1 _)] -- internal pattern for zero -- Pattern for "generalized" 1 pattern Pg1 :: forall t. [t] pattern Pg1 = [] -- Patterns for 0-2 pattern Pc0 :: Canon pattern Pc0 <- Bare 0 _ pattern Pc1 :: Canon pattern Pc1 <- Bare 1 _ pattern PoA :: Canon pattern PoA <- Pc1 -- addition operator pattern Pc2 :: Canon pattern Pc2 <- Bare 2 _ pattern PoM :: Canon pattern PoM <- Pc2 -- multiplication operator pattern Pc3 :: Canon pattern Pc3 <- Bare 3_ pattern PoE :: Canon pattern PoE <- Pc3 -- exponentiation operator pattern PcN1 :: Canon -- this pattern is only used in the "bad" function pattern PcN1 <- Can [(-1, Bare 1 _)] _ -- | Maximum exponent (of a polynomial) to distribute into a sum of terms. cMaxExpoToExpand :: Canon cMaxExpoToExpand = c4 -- need to finesse this to get the right operation returned. If sortByHypo has more than one entry, then it's a sum. -- Do we just create a hyper expr distHyperExpr, distSum, distProduct, distExpo :: Canon -> Canon -> ([Canon], Bool) distHyperExpr c m | h == cAddOpLevel = (dS, dS /= hL) | h == cMultOpLevel && fP = (dP, fP) -- distributed product (and maybe poly) | h == cExpOpLevel && fE = (dE, fE) -- distributed exponential | otherwise = ([c], False) where (h, hL) = (cGetHyperOp c, cGetHyperList c) ((dS,_), (dP,fP), (dE,fE)) = (distSum c m, distProduct c m, distExpo c m) distSum (HX PoA l _) m = (concat dM, any id fM) where (dM, fM) = unzip $ map (\x -> distHyperExpr x m) l -- check if any of the flags are set distSum c _ = ([c], False) distProduct c@(HX PoM l _) m | not (null sums) = (dist (head sums') (tail sums') [computeExpr cMultOpLevel nonSums], True) | otherwise = ([c], False) -- nothing to distribute where sums' = map exP sums exP s = if polyPred s m then simpleHX c1 (fst $ distExpo s m) else s -- clean up? (sums, nonSums) = partition (\e -> sumPred e || polyPred e m) l distProduct c _ = ([c], False) -- distribute sum against list of canons dist :: Canon -> [Canon] -> [Canon] -> [Canon] dist x y wL | length y > 0 = dist (head y) (tail y) cartProd | otherwise = cartProd where cartProd = if null wL then hLx else (concat $ map (\a -> map (* a) wL) hLx) hLx = cGetHyperList x distExpo c m | polyPred c m = cExpandPoly (cGetHyperList $ head cL) (fromInteger $ cToI (cL !! 1)) m | otherwise = ([c], False) where cL = cGetHyperList c sumPred :: Canon -> Bool sumPred c = cGetHyperOp c == cAddOpLevel polyPred :: Canon -> Canon -> Bool -- The 2nd param is a cutoff for the exponent of the polynomial polyPred (HX h (b:e:r) _) m = h == cExpOpLevel && sumPred b && e <= m && null r polyPred _ _ = False cExpandPoly :: [Canon] -> Int -> Canon -> ([Canon], Bool) -- e.g. (1 + x + y) ^ 3 cExpandPoly (a:r@(_:_:_)) e m = (eP, True) where eP = concat $ map (\x -> fst $ distHyperExpr x m) $ concat $ map (\x -> fst $ distHyperExpr x m) $ -- two dists needed to flatten it out? To Do: Investigate fst $ cExpandPoly [a, simpleHX cAddOpLevel r] e m cExpandPoly (a:b:_) e _ = (map (\i -> product [makeCanon $ binomial e i, raise a i, raise b (e-i)]) [0..e], True) where raise b' e' | e' == 0 = c1 | e' == 1 = b' | otherwise = expH b' (makeCanon $ toInteger e') cExpandPoly c _ _ = (c, False) -- No-op now. ToDo: Should this be an error condition factorial :: [Integer] factorial = (1 :: Integer) : 1 : zipWith (*) [2..] (tail factorial) binomial :: Int -> Int -> Integer binomial n m | n < 0 = error "Binomial: n must be >= 0" | m < 0 || m > n = error "Binomial: m must be >= 0 and <= n" | otherwise = div (factorial !! n) (factorial !! m * factorial !! (n-m)) -- This is essentially a wrapper to create a hyper expression from a hyper op and hyper list. -- This is lower level than calling cApplyHy / cHyperOp. Must be used with care computeExpr :: Canon -> [Canon] -> Canon -- computeExpr hy l | trace ("computeExpr: Processing: (hy=" ++ show hy ++ ", l=" ++ show l ++ ")") False = undefined computeExpr hy l | null nL = dV | length nL == 1 = head nL | otherwise = simpleHX hy nL where nL | hy == cAddOpLevel = filter (/= dV) (hE ++ [f nHe]) | hy == cMultOpLevel = filter (/= dV) ((f nHe):hE) | otherwise = l (hE, nHe) = partition cHyperExprAny l (dV, f) | hy == cAddOpLevel = (c0, sum) | otherwise = (c1, product) simplifySum, simplifyProd :: [Canon] -> Canon -- simplifySum l | trace ("simplifySum: Processing: (" ++ show l ++ ")") False = undefined simplifySum l = checkToFlipSum $ combineSum $ computeExpr cAddOpLevel l -- simplifyProd l | trace ("simplifyProd: Processing: (" ++ show l ++ ")") False = undefined simplifyProd l = if (any (== 0) l) then c0 else (combineProd $ computeExpr cMultOpLevel l) checkToFlipSum :: Canon -> Canon checkToFlipSum r = if cGetHyperOp r == cAddOpLevel && cNegative r then simpleHX cMultOpLevel [cN1, negate r] -- flip the signs and mult by negative one else r addH, multH, expH, addH', multH', expH' :: Canon -> Canon -> Canon (addH, multH, expH) = (prep addH', prep multH', prep expH') prep :: (Canon -> Canon -> Canon) -> Canon -> Canon -> Canon -- prep _ x y | trace ("prep: Processing: (" ++ show x ++ ") and (" ++ show y ++ ")") False = undefined prep f a b | (cHyperExprAny a || cHyperExprAny b) && (not (cIntegral a) || not (cIntegral b)) = error "Can't have sums or products with non-integers and hyper expressions" | otherwise = f a b -- addH' x y | trace ("addH': Processing: (" ++ show x ++ ") and (" ++ show y ++ ")") False = undefined addH' (HX PoA l1 _) (HX PoA l2 _) = simplifySum $ l1 ++ l2 addH' (HX PoA lA _) b | cHyperSum b = simplifySum $ lA ++ (negSumList b) | otherwise = simplifySum (b:lA) addH' b a@(HX PoA _ _) = addH' a b -- flip the terms addH' a@(HX PoM _ _) b@(HX PoM _ _) | aHs && bHs = simplifySum $ aNs ++ bNs | aHs = simplifySum (b:aNs) | bHs = simplifySum (a:bNs) | otherwise = simplifySum [a,b] where (aHs, bHs) = (cHyperSum a, cHyperSum b) (aNs, bNs) = (negSumList a, negSumList b) addH' a@(HX PoM _ _) b | cHyperSum a = simplifySum (b:negSumList a) | otherwise = simplifySum [a,b] addH' b a@(HX PoM _ _) = addH' a b addH' a b | cHyperExprAny a || cHyperExprAny b = simplifySum [a,b] | otherwise = a + b -- call the underlying function. Shouldn't happen in practice. -- multH' a b | trace ("multH': Processing: (" ++ show a ++ ") and (" ++ show b ++ ")") False = undefined multH' a b | nhs a && nhs b = simplifyProd [negate a, negate b] | cHyperExpr a || cHyperExpr b = if snp == cN1 then negate hp else hp -- flip sign if needed (abs value only internally) | otherwise = fst $ cMult a b crCycloInitMap where nhs x = cHyperSum x && cGetHyperOp x == cMultOpLevel -- negative hyper sum snp = signum a * signum b hp = simplifyProd $ cGetFactors (abs a) ++ cGetFactors (abs b) expH' a b | cHyperExprAny a || cHyperExprAny b = fst $ cHyperOp cExpOpLevel [a,b] crCycloInitMap | otherwise = fst $ cExp a b False crCycloInitMap negSumList :: Canon -> [Canon] negSumList c = map negate $ cGetHyperList ((cGetHyperList c) !! 1) -- e.g. -1 * (3 + 5) -> (-3 -5) -- | Convert a hyperexpression to a sum if possible. Useful in comparison. Will expand polynomials to a limited degree. cConvertToSum :: Canon -> (Canon, Bool) cConvertToSum x = cConvertToSum' x cMaxExpoToExpand cConvertToSum' :: Canon -> Canon -> (Canon, Bool) -- cConvertToSum' x m | trace ("cConvertToSum': (x = " ++ show x ++ ", m = " ++ show m ++ ")") False = undefined cConvertToSum' x m | b = (checkToFlipSum $ computeExpr cAddOpLevel (sortByHpo $ map f d), b) | otherwise = (x, b) -- just return itself. "dist" didn't do change the expression where f c = if cNegative c then (negate f') else f' where f' = product $ cFlattenAndGroup' (abs c) m (d,b) = distHyperExpr x m -- qL = cGetHyperList q, rL = cGetHyperList r, g = intersect qL rL, q' = qL \\ g, r' = rL' -- if there's a flag, filter out hyperexpressions from the gcd' -- this is the first level the next level is -- Note: These "simple" functions promote Canons to an exponent level before manipulating them then it demotes them. -- Example Transforms: a * b <^3 => [(a,1),(b,3)], x <^ y => [(x,y)], x + y => [(x+y, 1)], x <^> y => [(x <^> y, 1)] -- Pass ax and ay and return a * (x + y) where x and y are "relatively prime" simpleFactorPair :: Canon -> Canon -> Bool -> (Canon, Canon) -- simpleFactorPair n d hF | trace ("simpleFactorPair: (n = " ++ show n ++ ", d = " ++ show d ++ ", hF = " ++ show hF ++ ")") False = undefined simpleFactorPair n d hF = (applyFcnForHy v cMultOpLevel cFlattenAndGroup, gcd') where (gcd', x', y') = simpleReduce n d hF v = computeExpr cMultOpLevel [gcd',x'+y'] -- Promote each param so it can be easily manipulated and then demote it and return the answers. simpleReduce :: Canon -> Canon -> Bool -> (Canon, Canon, Canon) -- simpleReduce n d hF | trace ("simpleReduce: (n = " ++ show n ++ ", d = " ++ show d ++ ", hF = " ++ show hF ++ ")") False = undefined simpleReduce n d hF = (gcd', x', y') where proc xL@(x@(xP, xE):xs) yL wL wX wY | null yL = (wL, (wX ++ xL), wY) | null z = proc xs yL wL (wX ++ [x]) wY | length z > 1 = error "Entry occured was not unique in y!" | xE == yE = proc xs (yL \\ z) (wL ++ [(xP,xE)]) wX wY | xE < yE = proc xs (yL \\ z) (wL ++ [(xP,xE)]) wX (wY ++ [(xP,yE -xE)]) | otherwise = proc xs (yL \\ z) (wL ++ [(xP,yE)]) (wX ++ [(xP,xE - yE)]) wY where z = filter (\(p,_) -> xP == p) yL yE = snd $ head z proc _ yL wL wX wY = (wL, wX, (wY ++ yL)) ((nN, nH'), (dN, dH')) = (simplePrep n, simplePrep d) nHq = if hF then 1 else gcd nN dN -- if hyperFlag set, don't do non-hyper gcd dem (a, b, c) = (map expDemote a, map expDemote b, map expDemote c) (wL', wX', wY') = dem $ proc nH' dH' [] [] [] gcd' = computeExpr cMultOpLevel (nHq :wL') x' = computeExpr cMultOpLevel ((nN/nHq):wX') y' = computeExpr cMultOpLevel ((dN/nHq):wY') -- Split into non-hyper and hyper list and "expPromote" the hyper list simplePrep :: Canon -> (Canon, [(Canon, Canon)]) simplePrep c = (product nN, map expPromote nH) where (nH, nN) = partition cHyperExpr $ cGetFactors c -- Promote to exponential form : e.g. Two examples(1 + 3 <^> 4) -> (1 + 3<^>4, 1) and 7 <^ (5 <^> 4) -> (7, 5 <^> 4) expPromote :: Canon -> (Canon, Canon) expPromote v | cGetHyperOp v == cExpOpLevel = (head h, computeExpr cExpOpLevel (tail h)) | otherwise = (v, c1) where h = cGetHyperList v expPromoteFull :: Canon -> Canon expPromoteFull c | cGetHyperOp c > maxHyperOpDelveLevel = error "expPromoteFull: Can't perform this action. Max hyper op at base level exceeded." | otherwise = simpleHX cMultOpLevel newFactors where (hE, nonHe) = partition cHyperExpr $ cGetFactors $ cQuasiCanonize c prmNonHe :: [(Canon, Canon)] prmNonHe = map (\(p,e) -> (makeCanon p, e)) $ concat $ map cToGCR nonHe newFactors :: [Canon] newFactors = map (\(p,e) -> simpleHX cExpOpLevel [p,e]) $ sortOn fst $ prmNonHe ++ (concat $ map (\c' -> map expPromote $ cGetFactors c') hE) -- Demote a pair to an the canon itself or exponential hyper expression expDemote :: (Canon, Canon) -> Canon expDemote (p, e) = if e == c1 then p else computeExpr cExpOpLevel [p,e] -- | cFactor : Factor simple terms from a sum. If the flag is true, only factor by the gcd if the gcd is a hyper expression cFactorSum :: Canon -> Bool -> Canon cFactorSum c@(HX PoA hL _) hF | gcdL == c1 || (hF && not (cHyperExpr gcdL)) = c -- return as is | otherwise = computeExpr cMultOpLevel ((cGetFactors gcdL) ++ [computeExpr cAddOpLevel (map (\a -> div a gcdL) hL)]) where gcdL = foldl1 cGCD hL cFactorSum c _ = c -- To do: any additional poly factorizations -- | cFactorHorizon: Good for polynomial-like expressions like: (1 + 3<^>4 + 3<^>5) <^ 3 - 1, where there's a mixture of "canonical" and hE exponents. cFactorHorizon :: Canon -> Canon cFactorHorizon c | gcdL == c1 || length hL' == 1 = c -- return as is | otherwise = computeExpr cMultOpLevel ((cGetFactors gcdL) ++ [simpleHX cAddOpLevel aL]) where gcdL = foldl1 cGCD hL' hL' = map expPromoteFull $ cGetAddends $ fst $ cConvertToSum c mL = map (\a -> cCleanup $ div a gcdL) hL' (b,nB) = partition cBare mL sB = sum b aL = if sB == 0 then nB else (nB ++ [sB]) {- Cleanup / Hyperize / QuasiCanonize examples: Run cCleanup which is cHyperize . cQuasiCanonize Identity found as a result: (a <^> x) <^ (a <^> y) = (a <^> (y+1)) <^ (a <^> (y-1)) testsGood = [ -- Worked despite P3 bug (7 <^> (2<^2)) <^ 7 <^ (7 <<^>> 5), (2 <<^>> (7<^2 * 25303)) * (2 <<<^>>> (17 * 23 * 317)), (2 <<^>> (7<^2 * 25303)) * ((2 * 5) <<<^>>> (17 * 23 * 317)), 77 <<^>> 1239847 * 15 <<<<^>>>> 123947, 77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7, (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (79 ~~^~~ 101), (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (79 ~~^~~ 101), (1 + 3 <<^>> 7) <^> 5, (7 <<^>> 5) <^> 11, (35 <^> 5) <^ (35 <^> 3), 30 ~~^~~ 5, 3 * 3 <<^>> 5 * 6 <^> 5, 13 * (((3 + 5 <^> 7) <^> 4) <^> 5) <^> 7, (7 <^ (3 <<^>> 5)) * 7 <^> 4, (3 * 7 <^> 4) <^ 2, (28 <^ (3 <<^>> 5)) * 7 <^> 4, (((7 * 11) <<^>> (7<^2 * 25303)) <^ (13 ~~^~~ 101)) * (((3 * 5) <<<<^>>>> (17 * 23 * 317)) <^ (7 ~~^~~ 103)), (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7<^2) <^ (7 ~~^~~ 101), (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 ~~^~~ 101), (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 <<^>> 101), 10 <<^>> 1239847 * 10 <<<^>>> 123947, 14 <<^>> 1239847 * 10 <<<^>>> 123947, 60 <^> 4, 60 ~~^~~ 5, 90 ~~^~~ 5, 7 <<^>> 7 * 2 <<<<^>>>> 7 * 3 <^> 1234234 * 5 ~~~~^~~~~ 19 ~~~~^~~~~ 23982315987 * ((11*13*19*23*37*41) <^> 19) <^ 3, 75 <<^>> 11 * 45 <^> 5, 3 <^ 139478 * 9 <^> 6 * 27 <^> 77, (17 <^> 102) <^ (17 <^> 100),(7 <^> 5) <^ (7 <^> 3)] -- -> (17 <^> 101) <^> 2,(7 <^> 4) <^ (7 <^> 4), testHangDisplay0526Solved = [ (1 + 3 <<^>> 7) *((7 <^> 4) <^> 5 <^> 7), ((11 * 7 <<^>> 107) <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 <^> 101), (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4 * 3 <<^>> 1897 * 5 ~^~ 237, (7 * (1 + 3 <<^>> 7) <^> 5) <<^>> 17 * (1 + 3 <<^>> 7) ~~~~^~~~~ (1 + 7 ~~^~~ 1239847230), (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4 * 3 <<^>> 1897, (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4 * 3 <<^>> 1897 * 3 <^> 237, (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4 * 3 <<^>> 1897 * 3 ~^~ 237 * (cApplyHy (makeCanon 1001) [3,5] True), (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4 * 3 <<^>> 1897 * 3 ~^~ 237 * (cApplyHy (makeCanon 1001) [7,5] True), ((11 * 7 <<^>> 107) <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 <^ 101), ((3 * ((3 * 5) <<^>> 7) * (7 ~^~ 7)) <^> (2<^2)) * (3 <<^>> (7 * 271)) * (3 ~^~ (3 * 79)) * (cApplyHy (7 * 11 * 13) [2<^3, 5] True), (7 <<<^>>> 4) <<^>> 3 ] -- P3 testsQCHang0526Solved = [ (7<^ 2 * (7 <^> 4) <<^>> 13) <<<^>>>19, 13 * ((7 <^> 4) <^> 5 <^> 7), (7 <<^>> 3) <<^>> 3, (7 <^> 4) <<^>> 3, ((7 <<^>> (2<^2)) * (7 ~^~ (2<^2))) <<^>> 5, -- doing several comparisons (7 <<^>> 3) <<^>> 3, (7 <<^>> 4) <<^>> 3, (3 * 3 <<^>> 7) <^> 4, (30 <^> 7) ~~^~~ 5, (7 <<<^>>> 4) <<^>> 3, -- P3 ((3<^>17*7 ~^~ 4)) <<^>> 5, ((3*7 ~^~ 4)) <<^>> 5, 13 <^ (11 <<^>> 4) *((7 <^> 4) <^> 5 <^> 7), 13 <^ 4 *((7 <^> 4) <^> 5 <^> 7), ((3*7 <^> 4)) <<^>> 5, ((3 <<^>> 5) <^> 7) ~~^~~ 5, ((3 <^> 5) <<^>> 7) ~~^~~ 5, (3 * 15 <<^>> 7) <^> 4, (3 * 15 <<^>> 7 * 7 ~^~ 7) <^> 4, (4 <^> 7) ~~^~~ 5, (60 <^> 7) ~~^~~ 5, (12 <^> 7) ~~^~~ 5 ] -- Works, too but needs special functions (1 + 17 ~~~~~^~~~~~ 1329847 + 37 ~~~|^|~~~ grahamsNumber) ~~~^~~~ 8951237 testsU_P3 = [ -- Does not map them back (77 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 <^> 101), (539 <<^>> 1239847 * 15 <<<<^>>>> 123947 * 7) <^ (7 <^> 101)] Correctly doesn't convert: 13 <^> 3 Hangs: 3 ^ 3 * 6 <^> 4 -- Canonical issue 3 * 3 <<^>> 5 * 6 <^> 4 -- Canonical meeets Hyper expression issue (12 * 28 <^> 5) <^> 7 Utility func for verifying: v c = map hypMap $ map (\l -> (l !! 0, l !! 1)) $ map cGetHyperList $ cGetFactors $ cQuasiCanonize c -}