-- |
-- Module:      Math.NumberTheory.Canon
-- Copyright:   (c) 2015-2018 Frederick Schneider
-- Licence:     MIT
-- Maintainer:  Frederick Schneider <frederick.schneider2011@gmail.com> 
-- Stability:   Provisional
--
-- A Canon is an exponentation-based representation for arbitrarily massive numbers, including prime towers.

{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, PatternSynonyms, ViewPatterns, RankNTypes #-}

module Math.NumberTheory.Canon ( 
  makeCanon, makeC,
  canonToGCR, cToGCR,

  cMult, cDiv, cAdd, cSubtract, cExp,
  cReciprocal,
  cGCD, cLCM, cMod, cOdd, cTotient, cPhi,
  cLog, cLogDouble,
  cNegative, cPositive,
  cIntegral, cRational, cIrrational,
  cSimplify, cSimplified, 
  cDepth, 
  cSplit, cNumerator, cDenominator,
  cCanonical, cBare, cBareStatus, cValueType,
  cIsPrimeTower, cPrimeTowerLevel,
                                               
  cTetration, cPentation, cHexation, cHyperOp,
  (>^), (<^), (%), (<^>), (<<^>>), (<<<^>>>)         
)
where

import Math.NumberTheory.Primes.Testing (isPrime)
import Data.List (intersperse)
import GHC.Real (Ratio(..))
import Math.NumberTheory.Canon.Internals
import Math.NumberTheory.Canon.Additive
import Math.NumberTheory.Canon.AurifCyclo (CycloMap, crCycloInitMap)
import Math.NumberTheory.Canon.Simple (CanonConv(..))

-- | CanonValueType: 3 possibilities for this GADT.  Imaginary/complex numbers are not supported
data CanonValueType = IntegralC | NonIntRationalC | IrrationalC deriving (Eq, Ord, Show)

-- | GCR_ stands for Generalized Canonical Representation
type GCR_  = [GCRE_]
type GCRE_ = (Integer, Canon)

-- | Canon: GADT for either Bare or some variation of a canonical form.
data Canon = Bare Integer BareStatus | Canonical GCR_ CanonValueType 

-- | BareStatus: A "Bare Simplified" number means a prime number, +/-1 or 0.  The code must set the flag properly
--               A "Bare NotSimplified" number is an integer that has not been checked (to see if it can be factored).
data BareStatus = Simplified | NotSimplified deriving (Eq, Ord, Show)

makeCanon, makeC, makeCanonFull, makeDefCanonForExpnt :: Integer -> Canon

-- | Create a Canon from an integer.  This may involve expensive factorization.
makeCanon n = makeCanonI n False

-- | Shorthand for makeCanon
makeC       = makeCanon

-- | Make a Canon and attempt a full factorization
makeCanonFull n = makeCanonI n True

makeCanonI :: Integer -> Bool -> Canon
makeCanonI n b = crToC (crFromI n) b 
-- TODO: next step: enhance this once we can partially factor numbers

cCutoff :: Integer
cCutoff = 1000000

-- | Create type of Canon based on whether it exceeds a cutoff
makeDefCanonForExpnt e | e > cCutoff = Bare e (getBareStatus e) 
                       | otherwise   = makeCanonFull e

-- | 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              Simplified
crToC c    b | crSimplified c = Bare (fst $ head c) Simplified -- a little ugly
             | otherwise      = Canonical g (gcrCVT g)
                                where g          = map (\(p,e) -> (p, convPred e)) c
                                      convPred e | b =         makeCanonFull e        -- do complete factorization
                                                 | otherwise = makeDefCanonForExpnt e
                                                  -- Leave exponents "Bare" with flag based on if whether it's "simplified"
                                                 -- (can't be reduced any further)

-- | Instances for Canon
instance Eq Canon where
  x == y = cEq x y

instance Show Canon where 
  show (Bare n NotSimplified) = "(" ++ (show n) ++ ")" -- Note the extra characters.  This does not mean the figure is negative.
  show (Bare n Simplified)    = show n
  show c                      | denom == c1 = s numer False 
                              | otherwise   = s numer True ++ " / " ++ s denom True
                              where (numer, denom)      = cSplit c  
                                    s (Bare n f) _ = show $ Bare n f
                                    s v          w | w         = "(" ++ catList ++ ")" 
                                                   | otherwise = catList               -- w = with(out) parens
                                                   where catList   = concat $ intersperse " * " $ map sE $ cToGCR v -- sE means showElem
                                                         sE (p, e) | ptLevel > 2 = sp ++ " <^> " ++ s (makeCanonFull $ ptLevel) True
                                                                   | otherwise   = case e of
                                                                                     Bare 1 _ -> sp 
                                                                                     Bare _ _ -> sp ++ "^" ++ se
                                                                                     _        -> sp ++ " ^ (" ++ se ++ ")"
                                                                   where ptLevel = cPrimeTowerLevelI e p 1 
                                                                         sp      = show p
                                                                         se      = show e
                                                                          -- TODO: try for add'l collapse into <<^>>
                                                                                     
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 
  -- ToDo: mod n m     = fst $ cModBAD n m crCycloInitMap -- fix this "bad" logic and use this instead of the original function
  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

-- | Is the Canon a more complex expression? 
cCanonical :: Canon -> Bool
cCanonical (Canonical _ _ ) = 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      _ _ ) = IntegralC
cValueType (Canonical _ 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 (Bare x _ )            (Bare y _ )            = x == y
cEq (Bare _ Simplified)    (Canonical _ _ )       = False
cEq (Canonical _ _ )       (Bare _ Simplified)    = False

cEq (Bare x NotSimplified) y                      | cValueType y /= IntegralC = False
                                                  | otherwise                 = cEq (makeCanon x) y

cEq x                      (Bare y NotSimplified) | cValueType x /= IntegralC = False
                                                  | otherwise                 = cEq x (makeCanon y)

cEq (Canonical x a )       (Canonical y b)        = if a /= b then False else gcrEqCheck x y

-- | Check if a Canon is an odd integer.  Note: If the Canon is not integral, return False 
cOdd :: Canon -> Bool
cOdd (Bare x _)               = mod x 2 == 1
cOdd (Canonical c IntegralC ) = gcrOdd c
cOdd (Canonical _ _ )         = False

-- | GCD and LCM functions for Canon
cGCD, cLCM :: Canon -> Canon -> Canon
cGCD x y = cLGApply gcrGCD x y
cLCM x y = cLGApply gcrLCM x y

-- | Compute log as a Rational number.
cLog :: Canon -> Rational
cLog x = gcrLog $ cToGCR x 

-- | Compute log as a Double.
cLogDouble :: Canon -> Double
cLogDouble x = gcrLogDouble $ cToGCR x 

-- | Compare Function
cCmp :: Canon -> Canon -> Ordering
cCmp (Bare x _) (Bare y _) = compare x y
cCmp x          y          = gcrCmp (cToGCR x) (cToGCR y)

-- | QuotRem Function
cQuotRem :: Canon -> Canon -> CycloMap -> ((Canon, Canon), CycloMap)
cQuotRem x y m | 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
                                        -- ToDo: Left  _        -> (q,        md, mm) -- fix "cModBAD" and stop pointing to orig fcn
                                        Left  _        -> (q,        md, m')
                                        Right quotient -> (quotient, c0, m)
                                      where gy       = cToGCR y
                                            -- ToDo: fix  (md, mm) = cModBAD x y m'  -- Better to compute quotient this way .. to take adv. of alg. forms
                                            md       = cMod x y
                                            q        = gcrDivStrict (cToGCR d) gy  -- equivalent to: (x - x%y) / y.
                                            (d, m')  = cSubtract x md m

-- | Mod function
cMod :: Canon -> Canon -> Canon
cMod c m = if (cIntegral c) && (cIntegral m) then (makeCanon $ cModI c (cToI m))
                                             else error "cMod: Must both parameters must be integral"

cModI :: Canon -> Integer -> Integer
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 * cModI (cAbs c) (abs m)
                  | (cn && not mn) ||
                    (mn && not cn) = (signum m) * ( (abs m) - (cModI' (cAbs c) (abs m)) )
                  | otherwise        = cModI' c m
                    where cn         = cNegative c
                          mn         = m < 0

cModI' :: Canon -> Integer -> Integer
cModI' (Bare      n _         ) m = mod n m
cModI' (Canonical c IntegralC ) m = mod (product $ map (\x -> pmI (fst x) (mmt $ snd x) m) c) m
                                    where tm    = totient m
                                          mmt e = cModI e tm -- optimization
cModI' (Canonical _ _         ) _ = error "cModI': Logic error:  Canonical var has to be integral at this point" 

-- | 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"
             | 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           = makeC p -- "Canon-ize" this.  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

-- | Hyperoperations (including tetration and beyond): https://en.wikipedia.org/wiki/Hyperoperation
-- | The thinking around the operators is that they should look progressively scarier :)
infixr 9 <^>, <<^>>, <<<^>>>
(<^>), (<<^>>), (<<<^>>>) :: Canon -> Integer -> Canon
a <^>     b = fst $ cTetration a b crCycloInitMap
a <<^>>   b = fst $ cPentation a b crCycloInitMap
a <<<^>>> b = fst $ cHexation  a b crCycloInitMap

cTetration, cPentation, cHexation :: Canon -> Integer -> CycloMap -> (Canon, CycloMap)

-- | Tetration function
cTetration = cHyperOp 4 

-- | Pentation Function
cPentation = cHyperOp 5

-- | Hexation Function
cHexation  = cHyperOp 6

-- | Generalized Hyperoperation Function
cHyperOp :: Integer -> Canon -> Integer -> CycloMap -> (Canon, CycloMap)
cHyperOp n a b m | b < -1                       = error "Hyperoperations not defined when b < -1"
                 | n < 0                        = error "Hyperoperations require the level n >= 0"
                 | a /= c0 && a /= c1 && 
                   b > 1 && (a /= c2 && b == 2) = c n cb m
                 | otherwise                    = (sp n a b, m)
                 where cb = makeCanon b
                       -- Function for regular cases
                       c 1  b' m'  = cAdd a b' m'  -- addition
                       c 2  b' m'  = cMult a b' m' -- multiplication
                       c 3  b' m'  = (a <^ b', m')  -- exponentiation: ToDo: Plug in the CycloMap logic for expo.
                       -- Tetration and beyond
                       c _  Pc1 m' = (a, m')
                       c n' b'  m' = c (n'-1) r m'' -- TODO: Find a way to optimize this
                                     where (r, m'') = c n' (b'-1) m'

                       -- Function for special cases:                 
                       -- Note: When n (first param) is zero, that is known as "succession"
                       --   Cases when a is zero ...
                       sp 0 Pc0 b'   = makeCanon (b' + 1)
                       sp 1 Pc0 b'   = makeCanon b'
                       sp 2 Pc0 _    = c0
                       sp 3 Pc0 b'   = if b' == 0 then c1 else c0
                       sp _ Pc0 b'   = if (mod b' 2) == 1 then c0 else c1
                       --   Cases when b is zero ...
                       sp 0 _   0    = c1 
                       sp 1 a'  0    = a'
                       sp 2 _   0    = c0 
                       sp _ _   0    = c1 
                       --   Cases when b is -1 ...
                       sp 0 _   (-1) = c0
                       sp 1 a'  (-1) = a' - 1
                       sp 2 a'  (-1) = cNegate a'
                       sp 3 a'  (-1) = cReciprocal a'
                       sp _ _   (-1) = c0
                       --   Other Cases ...
                       sp h Pc2 2    | h == 0    = makeCanon 3
                                     | otherwise = makeCanon 4 -- recursive identity
                       sp _ Pc1 _    = c1
                       sp _ a'  1    = a'
                       sp _ _   _    = error "Can't compute this hyperoperation.  b must be >= -1"                

infixl 7 %
-- | Mod operator
(%) :: (Integral a) => a -> a -> a
n % m = mod n m 

-- | Exponentation operator declaration
infixr 9 <^   
-- Note: Even with Flexible Contexts switched on, it doesn't infer a bare number to be an Integer

-- | Dedicated multi-param typeclass for exponentiation operator.
class CanonExpnt a b c | a b -> c where 
  -- | Exponentiation operator
  (<^) :: a -> b -> c

instance CanonExpnt Canon Canon Canon where
  p <^ e = fst $ cExp p e True crCycloInitMap 
  
instance CanonExpnt Integer Integer Canon where
  p <^ e = fst $ cExp (makeCanon p) (makeDefCanonForExpnt e) True crCycloInitMap

instance CanonExpnt Canon Integer Canon where
  p <^ e = fst $ cExp p (makeDefCanonForExpnt e) True crCycloInitMap

instance CanonExpnt Integer Canon Canon where
  p <^ e = fst $ cExp (makeCanon p) e True crCycloInitMap

-- | Operator declaration: r >^ n means: attempt to take the rth root of n 
infixr 9 >^ 

-- | Dedicated multi-param typeclass for radical or root operator.
class CanonRoot a b c | a b -> c where 
  -- | Root operator
  (>^) :: a -> b -> c

instance CanonRoot Canon Canon Canon where
  r >^ n = cRoot n r
  
instance CanonRoot Integer Integer Canon where
  r >^ n = cRoot (makeCanon n) (makeCanon r)
  
instance CanonRoot Integer Canon Canon where
  r >^ n = cRoot n (makeCanon r) 

instance CanonRoot Canon Integer Canon where
  r >^ n = cRoot (makeCanon n) r 

-- | Check if underlying rep is simplified
crSimplified :: CR_ -> Bool
crSimplified POne  = True
crSimplified PZero = True                                                        
crSimplified PN1   = True  
crSimplified c     = crPrime c

-- | Convert a Canon back to its underlying rep (if possible).
cToCR :: Canon -> CR_
cToCR (Canonical c v)         | v /= IrrationalC = gcrToCR c 
                              | otherwise        = error "cToCR: Cannot convert irrational canons to underlying data structure"
cToCR (Bare 1 _ )             = cr1
cToCR (Bare n NotSimplified)  = crFromI n
cToCR (Bare n Simplified)     = [(n,1)] -- not ideal

-- | Convert generalized canon rep to Canon.
gcrToC :: GCR_ -> Canon
gcrToC g | gcrBare g = Bare (gcrToI g) Simplified
         | otherwise = Canonical g (gcrCVT g)

-- | For generalized canon rep, determine the CanonValueType.   
gcrCVT :: GCR_ -> CanonValueType         
gcrCVT POne = IntegralC
gcrCVT g    = g' g IntegralC -- start Integral, can only get "worse"
              where g' _           IrrationalC = IrrationalC -- short-circuits once IrrationalC 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 IrrationalC _                             = IrrationalC
                    dcv _           (Canonical _ IrrationalC)     = IrrationalC
                    dcv _           (Canonical _ NonIntRationalC) = IrrationalC
                    dcv IntegralC   (Bare      n _ )              = if n < 0 then NonIntRationalC else IntegralC
                    dcv v           (Bare      _ _ )              = v
                    dcv v           c                             = if cNegative c then NonIntRationalC else v

c1, c0, cN1, c2 :: Canon
c1  = makeCanon 1
c0  = makeCanon 0
cN1 = makeCanon (-1)
c2  = makeCanon 2

-- | Convert Canon to Integer if possible.
cToI :: Canon -> Integer
cToI (Bare i _ )     = i
cToI (Canonical c v) | v == IntegralC = gcrToI c 
                     | otherwise      = error "Can't convert non-integral canon to an integer"

-- | Convert Canon To Double.
cToD :: Canon -> Double
cToD (Bare i _ )      = fromIntegral i
cToD (Canonical c _ ) = gcrToD c 

-- | Multiply Function: Generally speaking this will be much cheaper.
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 = (gcrToC g, 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.
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   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 = (gcd' * r, m')
                               where gcd'    = cGCD x y 
                                     x'      = x / gcd'
                                     y'      = y / gcd'
                                     r       = crToC c False
                                     (c, m') = crApplyAdtvOptConv b (cToCR x') (cToCR y') m -- costly bit

-- | Exponentiation: This does allow for negative exponentiation if the Bool flag is True.
cExp :: Canon -> Canon -> Bool -> CycloMap -> (Canon, CycloMap)
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 m
                         where cExp' Pc0 e'  m' | cPositive e' = (c0, m')
                                                | otherwise    = error "0^e where e <= 0 is either undefined or illegal"
                               cExp' Pc1 _   m' = (c1, m')
                               cExp' _   Pc0 m' = (c1, m')
                               cExp' c'   e' m' = (gcrToC g, mg)
                                                  where (g, mg) = gE (cToGCR c') e' m' 
                               gE :: GCR_ -> Canon -> CycloMap -> (GCR_, CycloMap)
                               gE g' e' m' | gcrNegative g' 
                                             = case cValueType e' of  -- gcr exponentiation
                                                 IntegralC       -> if cOdd e' then (gcreN1:absTail, m'')
                                                                               else (absTail, m'')
                                                 NonIntRationalC -> if cOdd d then (gcreN1:absTail, m'')
                                                                              else error "gE: Imaginary numbers not supported"
                                                 IrrationalC     -> 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 denominator means you will have 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 (Bare n      _ ) = n < 0
cNegative (Canonical c _ ) = gcrNegative c

cPositive (Bare n _      ) = n > 0
cPositive (Canonical c _ ) = gcrPositive c

-- | Functions for negation, absolute value and signum
cNegate, cAbs, cSignum :: Canon -> Canon 

cNegate (Bare 0 _)             = c0
cNegate (Bare 1 _)             = cN1
cNegate (Bare x Simplified)    = Canonical (gcreN1 : [(x, c1)]) IntegralC -- prepend a "-1", not ideal
cNegate (Bare x NotSimplified) = Bare (-1 * x) NotSimplified 
cNegate (Canonical x v)        = gcrNegateCanonical x v
      
cAbs x | cNegative x = cNegate x
       | otherwise   = x

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 _ Pc0 y   = y
cLGApply _ x   Pc0 = x
cLGApply f x   y   | cNegative x || 
                     cNegative y = gcrToC $ f (cToGCR $ cAbs x) (cToGCR $ cAbs y)
                   | otherwise   = gcrToC $ f (cToGCR x)        (cToGCR y)

-- | Div function : Multiply by the reciprocal.
cDiv :: Canon -> Canon -> CycloMap -> (Canon, CycloMap)
cDiv _ Pc0 _ = error "cDiv: Division by zero error"
cDiv x y   m = cMult (cReciprocal y) x m -- multiply by the reciprocal

-- | Compute reciprocal (by negating exponents).
cReciprocal :: Canon -> Canon
cReciprocal x = fst $ cExp x cN1 True crCycloInitMap  -- raise number to (-1)st power

-- | Functions to check if a Canon is integral, (ir)rational, "simplified" or a prime tower
cIntegral, cIrrational, cRational, cSimplified, cIsPrimeTower :: Canon -> Bool

cIntegral (Bare      _ _ ) = True
cIntegral (Canonical _ v ) = v == IntegralC

cIrrational (Canonical _ IrrationalC ) = True
cIrrational _                          = False

cRational c = not $ cIrrational c

cSimplified (Bare      _ Simplified)    = True
cSimplified (Bare      _ NotSimplified) = True
cSimplified (Canonical c _ )            = gcrSimplified c

cIsPrimeTower c = cPrimeTowerLevel c > 0 -- x^x would not be, but x^x^x would be

-- | cNumerator and cDenominator are for processing "rational" canon reps.
cNumerator, cDenominator :: Canon -> Canon

cNumerator (Canonical c _ ) = gcrToC $ filter (\x -> cPositive $ snd x) c -- filter positive exponents
cNumerator b                = b 

cDenominator (Canonical c _ ) = gcrToC $ map (\(p,e) -> (p, cN1 * e)) $ filter (\(_,e) -> cNegative e) c -- negate negative expnts
cDenominator _                = c1 

-- | Determines the depth/height of maximum prime tower in the Canon.
cDepth :: Canon-> Integer
cDepth (Bare      _ _ ) = 1
cDepth (Canonical c _ ) = 1 + gcrDepth c

-- | Force the expression to be simplified.  This could potentially be very expensive.
cSimplify :: Canon -> Canon
cSimplify (Bare      n NotSimplified) = makeCanonFull n
cSimplify (Canonical c _ )            = gcrToC $ gcrSimplify c
cSimplify g                           = g  -- Bare number already simplified : Fix when expr come into play

-- | Compute the rth-root of a Canon.
cRoot :: Canon -> Canon -> Canon 
cRoot c r | cNegative r = error "r-th roots are not allowed when r <= 0" 
          | otherwise   = gcrToC $ gcrRootI (cToGCR c) (cToGCR r)

-- | This is used for tetration, etc.  It defaults to zero for non-integral reps.
cPrimeTowerLevel :: Canon -> Integer                  
cPrimeTowerLevel (Bare      _ Simplified) = 1
cPrimeTowerLevel (Canonical g IntegralC)  = case gcrPrimePower g of
                                              False -> 0
                                              True  -> cPrimeTowerLevelI (snd $ head g) (fst $ head g) (1 :: Integer)
cPrimeTowerLevel _                        = 0

-- | Internal workhorse function
cPrimeTowerLevelI :: Canon -> Integer -> Integer -> Integer
cPrimeTowerLevelI (Bare b _ )             n l | b == n    = l + 1 
                                              | otherwise = 0
cPrimeTowerLevelI (Canonical g IntegralC) n l | gcrPrimePower g == False = 0 
                                              | n /= (fst $ head g)      = 0
                                              | otherwise                = cPrimeTowerLevelI (snd $ head g) n (l+1)
cPrimeTowerLevelI _                       _ _ = 0

-- | Functions to Convert Canon to Generalized Canon Rep
canonToGCR, cToGCR :: Canon -> GCR_
canonToGCR (Canonical x _) = x
canonToGCR (Bare x NotSimplified) = canonToGCR $ makeCanon x -- ToDo: Thread in CycloMap?
canonToGCR (Bare x Simplified)    | x == 1    = gcr1 
                                  | otherwise = [(x, c1)]
cToGCR = canonToGCR

-- 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) Simplified
                                            False -> Canonical (tail g) v             
                        | otherwise     = Canonical (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
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 compare (gcrLogDouble x) (gcrLogDouble y) of
                              -- If equal: we have to break out the big guns, both evaluate to infinity
                             EQ  -> compare (gcrLog'' x) (gcrLog'' y) 
                             cmp -> cmp

           where xN          = gcrNegative x
                 yN          = gcrNegative y  

                 -- This is much more expensive but accurate. You have an "infinity" result issue potentially with gcrLogDouble
                 gcrLog'' g = sum $ map f g
                 f (p,e)    = (toRational $ logDouble $ fromIntegral p) * (toRational e)
                 logDouble :: Double -> Double
                 logDouble n = log n
                
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 

gcrLog :: GCR_ -> Rational
gcrLog g = crLog $ gcrToCR g   

-- | 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 POne _    = gcr1
gcrGCD _    POne = gcr1
gcrGCD x    y    = case compare xp yp of
                      LT -> gcrGCD xs y
                      EQ -> (xp, min xe ye) : gcrGCD xs ys                    
                      GT -> gcrGCD x ys
                    where ((xp,xe):xs) = x
                          ((yp,ye):ys) = y    
gcrLCM POne y    = y
gcrLCM x    POne = x        
gcrLCM x    y    = case compare xp yp of
                      LT -> xh : gcrLCM xs y
                      EQ -> (xp, max xe ye) : gcrLCM xs ys
                      GT -> yh : gcrLCM x ys
                    where (xh@(xp,xe) : xs)  = x
                          (yh@(yp,ye) : ys)  = y  

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 

gcrIsZero :: GCR_ -> Bool
gcrIsZero Pg0 = True;
gcrIsZero _   = 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)

gcrDepth :: GCR_ -> Integer
gcrDepth g = maximum $ map (\(_,e) -> cDepth e) g

gcrSimplified :: GCR_ -> Bool
gcrSimplified g = all (\(_,e) -> cSimplified e) g               

gcrSimplify :: GCR_ -> GCR_
gcrSimplify g = map (\(p,e) -> (p, cSimplify e)) g

gcrRootI :: GCR_ -> GCR_ -> GCR_
gcrRootI POne _ = gcr1   
gcrRootI c    r | not $ gcrNegative c = case gcrDiv (cToGCR ce) r of
                                          Left  _        -> error e 
                                          Right quotient -> (cp, gcrToC quotient) : gcrRootI cs r
                | gcrEven r           = error "Imaginary numbers not allowed: Even root of negative number requested."
                | otherwise           = gcreN1 : gcrRootI (gcrAbs c) r
                where ((cp,ce):cs) = c  
                      e            = "gcrRootI: All expnts must be multiples of " ++ (show r) ++ ".  Not so with " ++ (show c)

-- | Check if the number is simplified rather than factoring it.  Simplified is equivalent to having one term in the list.
getBareStatus :: Integer -> BareStatus
getBareStatus n | n < -1              = NotSimplified 
                | n <= 1 || isPrime n = Simplified 
                | otherwise           = NotSimplified

-- | 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

-- | Patterns for 0 and 1
pattern Pc0 :: Canon
pattern Pc0 <- Bare 0 _

pattern Pc1 :: Canon
pattern Pc1 <- Bare 1 _ 

pattern Pc2 :: Canon
pattern Pc2 <- Bare 2 _

-- ToDo: Fix this Mod function.  "Proper" rewrite has terrible performance
{-
pattern PcN1 :: Canon  -- this pattern is only used in the "bad" function
pattern PcN1 <- Canonical [(-1, Bare 1 _)] _

cModBAD :: Canon -> Canon -> CycloMap -> (Canon, CycloMap)
cModBAD c m cm | cIntegral c && cIntegral m = f c m cm
            | otherwise                  = error "cModBAD: Must both parameters must be integral"
     where f _   Pc0  _   = error "cModBAD: Divide by zero error when computing n mod 0"
           f _   Pc1  cm' = (0, cm')
           f _   PcN1 cm' = (0, cm')
           f Pc0 _    cm' = (0, cm')
           f c'  m'   cm' | m' == c0         = error "cModBAD: Divide by zero error when computing n mod 0"
                          | ma == c1         = (c0, cm')
                          | ca == ma         = (c0, cm')
                          | cn && mn         = (cNegate mrn, cmn) -- both (n)egative
                          | (not cn) && (not mn) &&
                            ca < ma          = (ca, cm')
                          | (cn && not mn) ||
                            (mn && not cn) = ((cSignum m') * (makeC $ maI - mrm), cmm) -- (m)ixed sign: TODO: CycloMap threading
                          | otherwise        = (makeC io, mo)
                          where (cn, mn)     = (cNegative c', cNegative m')
                                (ca, ma)     = (cAbs c', cAbs m')
                                (mrn, cmn)   = f ca ma cm'
                                (mrm, cmm)   = f' ca maI cm'
                                maI          = cToI ma
                                (io, mo)     = f' c' (cToI m') cm'
           f' (Bare      n  _           ) mI cm' = (mod n mI, cm')
           f' ic@(Canonical c' IntegralC) mI cm' | cNegative ic = error "The canon must be positive here"
                                                 | otherwise    = (mod ip mI, cmf)
                                                 where (ip, cmf) = i c' cm'' (1 :: Integer) -- performs fold-like product
                                                       i []     cmi pri = (pri, cmi)          -- with CycloMap threading
                                                       i (l:ls) cmi pri | pri == 0   = (pri, cmi)
                                                                        | otherwise = i ls cmv (pri * v)
                                                                        where (v, cmv) = pf l cmi
                                                       pf (p,e) mp      = (pmI p (cToI v) mI, cmv)
                                                                          where (v, cmv) = f e tm mp
                                                       (tm, cm'')       = cTotient ic cm'
           f' (Canonical _  _           ) _  _   = error "cModBAD: Logic error:  Canonical var has to be integral at this point"
-}