{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE ViewPatterns        #-}

-- |
-- Module      : Numeric.QuadraticIrrational
-- Description : An implementation of quadratic irrationals
-- Copyright   : © 2014 Johan Kiviniemi
-- License     : MIT
-- Maintainer  : Johan Kiviniemi <devel@johan.kiviniemi.name>
-- Stability   : provisional
-- Portability : FlexibleContexts, ViewPatterns
--
-- A library for exact computation with
-- <http://en.wikipedia.org/wiki/Quadratic_irrational quadratic irrationals>
-- with support for exact conversion from and to
-- <http://en.wikipedia.org/wiki/Periodic_continued_fraction (potentially periodic) simple continued fractions>.
--
-- A quadratic irrational is a number that can be expressed in the form
--
-- > (a + b √c) / d
--
-- where @a@, @b@ and @d@ are integers and @c@ is a square-free natural number.
--
-- Some examples of such numbers are
--
-- * @7/2@,
--
-- * @√2@,
--
-- * @(1 + √5)\/2@
--   (<http://en.wikipedia.org/wiki/Golden_ratio the golden ratio>),
--
-- * solutions to quadratic equations with rational constants – the
--   <http://en.wikipedia.org/wiki/Quadratic_formula quadratic formula> has a
--   familiar shape.
--
-- A simple continued fraction is a number expressed in the form
--
-- > a + 1/(b + 1/(c + 1/(d + 1/(e + …))))
--
-- or alternatively written as
--
-- > [a; b, c, d, e, …]
--
-- where @a@ is an integer and @b@, @c@, @d@, @e@, … are positive integers.
--
-- Every finite SCF represents a rational number and every infinite, periodic
-- SCF represents a quadratic irrational.
--
-- > 3.5      = [3; 2]
-- > (1+√5)/2 = [1; 1, 1, 1, …]
-- > √2       = [1; 2, 2, 2, …]

module Numeric.QuadraticIrrational
  ( -- * Constructors and deconstructors
    QI, qi, qi', runQI, runQI', unQI, unQI'
  , -- * Lenses
    _qi, _qi', _qiABD, _qiA, _qiB, _qiC, _qiD
  , -- * Numerical operations
    qiIsZero
  , qiToFloat
  , qiAddI, qiSubI, qiMulI, qiDivI
  , qiAddR, qiSubR, qiMulR, qiDivR
  , qiNegate, qiRecip, qiAdd, qiSub, qiMul, qiDiv, qiPow
  , qiFloor
  , -- * Continued fractions
    continuedFractionToQI, qiToContinuedFraction
  , continuedFractionApproximate
  , module Numeric.QuadraticIrrational.CyclicList
  ) where

import Control.Monad.Trans.State (evalState, gets, modify)
import Data.Foldable (toList)
import Data.List (foldl')
import Data.Maybe (fromMaybe)
import Data.Proxy
import Data.Ratio ((%), denominator, numerator)
import Data.Semigroup
import qualified Data.Set as Set (empty, insert, member)
import Data.Type.Equality
import Math.NumberTheory.Primes (Prime, unPrime, factorise)
import Math.NumberTheory.Roots (integerSquareRoot)
import Text.Read (Lexeme (Ident), Read (readListPrec, readPrec),
  lexP, parens, prec, readListPrecDefault, step)

import Numeric.QuadraticIrrational.CyclicList
import Numeric.QuadraticIrrational.Internal.Lens
import Numeric.QuadraticIrrational.Ring
import Numeric.QuadraticIrrational.SquareFree

-- $setup
-- >>> import Data.Number.CReal

-- | @(a + b √c)@
data QI where
  QI :: KnownSquareFree c => QuadExt c Rational -> QI

instance Eq QI where
  (QI (QuadExt c Rational
q1 :: QuadExt c1 Rational)) == :: QI -> QI -> Bool
== (QI (QuadExt c Rational
q2 :: QuadExt c2 Rational)) = case forall (a :: SquareFree) (b :: SquareFree).
(KnownSquareFree a, KnownSquareFree b) =>
Proxy a -> Proxy b -> Maybe (a :~: b)
sameSquareFree (forall {k} (t :: k). Proxy t
Proxy :: Proxy c1) (forall {k} (t :: k). Proxy t
Proxy :: Proxy c2) of
      Maybe (c :~: c)
Nothing   -> Bool
False
      Just c :~: c
Refl -> QuadExt c Rational
q1 forall a. Eq a => a -> a -> Bool
== QuadExt c Rational
q2

instance Show QI where
  showsPrec :: Int -> QI -> ShowS
showsPrec Int
p QI
q = forall a.
QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
runQI QI
q forall a b. (a -> b) -> a -> b
$ \Integer
a Integer
b Integer
c Integer
d -> Bool -> ShowS -> ShowS
showParen (Int
p forall a. Ord a => a -> a -> Bool
> Int
10)
                           forall a b. (a -> b) -> a -> b
$ String -> ShowS
showString String
"qi " forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 Integer
a
                           forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar   Char
' '   forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 Integer
b
                           forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar   Char
' '   forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 Integer
c
                           forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> ShowS
showChar   Char
' '   forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a. Show a => Int -> a -> ShowS
showsPrec Int
11 Integer
d

instance Read QI where
  readPrec :: ReadPrec QI
readPrec = forall a. ReadPrec a -> ReadPrec a
parens ReadPrec QI
rQI
    where
      rQI :: ReadPrec QI
rQI = forall a. Int -> ReadPrec a -> ReadPrec a
prec Int
10 forall a b. (a -> b) -> a -> b
$ do
        Ident String
"qi" <- ReadPrec Lexeme
lexP
        Integer -> Integer -> Integer -> Integer -> QI
qi forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. ReadPrec a -> ReadPrec a
step forall a. Read a => ReadPrec a
readPrec forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall a. ReadPrec a -> ReadPrec a
step forall a. Read a => ReadPrec a
readPrec forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall a. ReadPrec a -> ReadPrec a
step forall a. Read a => ReadPrec a
readPrec
           forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> forall a. ReadPrec a -> ReadPrec a
step forall a. Read a => ReadPrec a
readPrec

  readListPrec :: ReadPrec [QI]
readListPrec = forall a. Read a => ReadPrec [a]
readListPrecDefault

type QITuple = (Integer, Integer, Integer, Integer)

-- | Given @a@, @b@, @c@ and @d@ such that @n = (a + b √c)\/d@, constuct a 'QI'
-- corresponding to @n@.
--
-- >>> qi 3 4 5 6
-- qi 3 4 5 6
--
-- >>> qi 3 0 42 1
-- qi 3 0 42 1
--
-- >>> qi 1 1 (-1) 1
-- qi 1 1 (-1) 1
--
-- The fractions are reduced:
--
-- >>> qi 30 40 5 60
-- qi 3 4 5 6
-- >>> qi (-30) (-40) (-5) (-60)
-- qi 3 4 (-5) 6
--
-- If @c = 0@ then @b@ is zeroed and @c@ is put equal to 2:
--
-- >>> qi 3 42 0 1
-- qi 3 0 2 1
--
-- The @b √c@ term is simplified:
--
-- >>> qi 0 1 (5*5*6) 1
-- qi 0 5 6 1
--
-- If @c = 1@ (after simplification) then @b@ is moved to @a@:
--
-- >>> qi 1 5 (2*2) 1
-- qi 11 0 2 1
qi :: Integer  -- ^ a
   -> Integer  -- ^ b
   -> Integer  -- ^ c
   -> Integer  -- ^ d
   -> QI
qi :: Integer -> Integer -> Integer -> Integer -> QI
qi Integer
a Integer
_ Integer
0 (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qi" -> Integer
d) = Integer -> Integer -> Integer -> Integer -> QI
simplifyReduceCons Integer
a Integer
0 Integer
2 Integer
d
qi Integer
a Integer
b Integer
c (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qi" -> Integer
d) = Integer -> Integer -> Integer -> Integer -> QI
simplifyReduceCons Integer
a Integer
b Integer
c Integer
d
{-# INLINE qi #-}

-- Construct a 'QI' without simplifying @b √c@. Make sure it has already been
-- simplified.
qiNoSimpl :: Integer -> Integer -> Integer -> Integer -> QI
qiNoSimpl :: Integer -> Integer -> Integer -> Integer -> QI
qiNoSimpl Integer
a Integer
b (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qiNoSimpl" -> Integer
c) (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qiNoSimpl" -> Integer
d)
  | Integer
b forall a. Eq a => a -> a -> Bool
== Integer
0    = Integer -> Integer -> Integer -> Integer -> QI
reduceCons Integer
a Integer
0 Integer
2 Integer
d
  | Integer
c forall a. Eq a => a -> a -> Bool
== Integer
0    = Integer -> Integer -> Integer -> Integer -> QI
reduceCons Integer
a Integer
0 Integer
2 Integer
d
  | Integer
c forall a. Eq a => a -> a -> Bool
== Integer
1    = Integer -> Integer -> Integer -> Integer -> QI
reduceCons (Integer
a forall a. Num a => a -> a -> a
+ Integer
b) Integer
0 Integer
2 Integer
d
  | Bool
otherwise = Integer -> Integer -> Integer -> Integer -> QI
reduceCons Integer
a Integer
b Integer
c Integer
d
{-# INLINE qiNoSimpl #-}

-- Simplify @b √c@ before constructing a 'QI'.
simplifyReduceCons :: Integer -> Integer -> Integer -> Integer -> QI
simplifyReduceCons :: Integer -> Integer -> Integer -> Integer -> QI
simplifyReduceCons Integer
a Integer
b (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"simplifyReduceCons" -> Integer
c) Integer
d
  | Integer
c' forall a. Eq a => a -> a -> Bool
== Integer
1   = Integer -> Integer -> Integer -> Integer -> QI
reduceCons (Integer
a forall a. Num a => a -> a -> a
+ Integer
b') Integer
0 Integer
2 Integer
d
  | Bool
otherwise = Integer -> Integer -> Integer -> Integer -> QI
reduceCons Integer
a Integer
b' Integer
c' Integer
d
  where ~(Integer
b', Integer
c') = Integer -> Integer -> (Integer, Integer)
separateSquareFactors Integer
b Integer
c
{-# INLINE simplifyReduceCons #-}

-- | Given @b@ and @c@ such that @n = b √c@, return a potentially simplified
-- @(b, c)@.
separateSquareFactors :: Integer -> Integer -> (Integer, Integer)
separateSquareFactors :: Integer -> Integer -> (Integer, Integer)
separateSquareFactors Integer
b (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"separateSquareFactors" -> Integer
c) =
  case forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (Integer, Integer) -> (Prime Integer, Word) -> (Integer, Integer)
go (Integer
1,Integer
1) (forall a. UniqueFactorisation a => a -> [(Prime a, Word)]
factorise (forall a. Num a => a -> a
abs Integer
c)) of
    ~(Integer
bMul, Integer
c') -> (Integer
bforall a. Num a => a -> a -> a
*Integer
bMul, Integer
c' forall a. Num a => a -> a -> a
* forall a. Num a => a -> a
signum Integer
c)
  where
    go :: (Integer, Integer) -> (Prime Integer, Word) -> (Integer, Integer)
    go :: (Integer, Integer) -> (Prime Integer, Word) -> (Integer, Integer)
go ~(Integer
i, Integer
j) ~(Prime Integer
p, Word
pow) =
      Integer
i seq :: forall a b. a -> b -> b
`seq` Integer
j seq :: forall a b. a -> b -> b
`seq` Prime Integer
p seq :: forall a b. a -> b -> b
`seq` Word
pow seq :: forall a b. a -> b -> b
`seq`
        (Integer
i forall a. Num a => a -> a -> a
* forall a. Prime a -> a
unPrime Prime Integer
p forall a b. (Num a, Integral b) => a -> b -> a
^ (Word
pow forall a. Integral a => a -> a -> a
`div` Word
2), if forall a. Integral a => a -> Bool
even Word
pow then Integer
j else Integer
j forall a. Num a => a -> a -> a
* forall a. Prime a -> a
unPrime Prime Integer
p)

-- Reduce the @a@, @b@, @d@ factors before constructing a 'QI'.
reduceCons :: Integer -> Integer -> Integer -> Integer -> QI
reduceCons :: Integer -> Integer -> Integer -> Integer -> QI
reduceCons Integer
a Integer
b Integer
c (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"reduceCons" -> Integer
d) = case Integer -> Maybe SomeSquareFree
someSquareFreeVal Integer
c of
  Just (SomeSquareFree (Proxy c
Proxy :: Proxy c)) -> forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (forall (c :: SquareFree) t. t -> t -> QuadExt c t
QuadExt (Integer
a forall a. Integral a => a -> a -> Ratio a
% Integer
d) (Integer
b forall a. Integral a => a -> a -> Ratio a
% Integer
d) :: QuadExt c Rational)
  Maybe SomeSquareFree
Nothing -> forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"reduceCons: square root term is not square free and equals to " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Integer
c
{-# INLINE reduceCons #-}

-- | Given @a@, @b@ and @c@ such that @n = a + b √c@, constuct a 'QI'
-- corresponding to @n@.
--
-- >>> qi' 0.5 0.7 2
-- qi 5 7 2 10
-- >>> qi' 0.3 0.4 (-3)
-- qi 3 4 (-3) 10
qi' :: Rational  -- ^ a
    -> Rational  -- ^ b
    -> Integer   -- ^ c
    -> QI
qi' :: Rational -> Rational -> Integer -> QI
qi' Rational
a Rational
b Integer
c = QI
n
  where
    -- (aN/aD) + (bN/bD) √c = ((aN bD) + (bN aD) √c) / (aD bD)
    n :: QI
n = Integer -> Integer -> Integer -> Integer -> QI
qi (Integer
aNforall a. Num a => a -> a -> a
*Integer
bD) (Integer
bNforall a. Num a => a -> a -> a
*Integer
aD) Integer
c (Integer
aDforall a. Num a => a -> a -> a
*Integer
bD)
    (Integer
aN, Integer
aD) = (forall a. Ratio a -> a
numerator Rational
a, forall a. Ratio a -> a
denominator Rational
a)
    (Integer
bN, Integer
bD) = (forall a. Ratio a -> a
numerator Rational
b, forall a. Ratio a -> a
denominator Rational
b)
{-# INLINE qi' #-}

-- | Given @n@ and @f@ such that @n = (a + b √c)\/d@, run @f a b c d@.
--
-- >>> runQI (qi 3 4 5 6) (\a b c d -> (a,b,c,d))
-- (3,4,5,6)
-- >>> runQI (qi 1 1 (-1) 1) (\a b c d -> (a,b,c,d))
-- (1,1,-1,1)
runQI :: QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
runQI :: forall a.
QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
runQI (QI (QuadExt Rational
a Rational
b :: QuadExt c Rational)) Integer -> Integer -> Integer -> Integer -> a
f = Integer -> Integer -> Integer -> Integer -> a
f Integer
a' Integer
b' Integer
c Integer
d
  where
    c :: Integer
c = forall (c :: SquareFree). KnownSquareFree c => Proxy c -> Integer
squareFreeVal (forall {k} (t :: k). Proxy t
Proxy :: Proxy c)
    d :: Integer
d = forall a. Integral a => a -> a -> a
lcm (forall a. Ratio a -> a
denominator Rational
a) (forall a. Ratio a -> a
denominator Rational
b)
    a' :: Integer
a' = forall a. Ratio a -> a
numerator Rational
a forall a. Num a => a -> a -> a
* (Integer
d forall a. Integral a => a -> a -> a
`quot` forall a. Ratio a -> a
denominator Rational
a)
    b' :: Integer
b' = forall a. Ratio a -> a
numerator Rational
b forall a. Num a => a -> a -> a
* (Integer
d forall a. Integral a => a -> a -> a
`quot` forall a. Ratio a -> a
denominator Rational
b)
{-# INLINE runQI #-}

-- | Given @n@ and @f@ such that @n = a + b √c@, run @f a b c@.
--
-- >>> runQI' (qi' 0.5 0.7 2) (\a b c -> (a, b, c))
-- (1 % 2,7 % 10,2)
-- >>> runQI' (qi' 0.3 0.4 (-3)) (\a b c -> (a, b, c))
-- (3 % 10,2 % 5,-3)
runQI' :: QI -> (Rational -> Rational -> Integer -> a) -> a
runQI' :: forall a. QI -> (Rational -> Rational -> Integer -> a) -> a
runQI' (QI (QuadExt Rational
a Rational
b :: QuadExt c Rational)) Rational -> Rational -> Integer -> a
f = Rational -> Rational -> Integer -> a
f Rational
a Rational
b Integer
c
  where
    c :: Integer
c = forall (c :: SquareFree). KnownSquareFree c => Proxy c -> Integer
squareFreeVal (forall {k} (t :: k). Proxy t
Proxy :: Proxy c)
{-# INLINE runQI' #-}

-- | Given @n@ such that @n = (a + b √c)\/d@, return @(a, b, c, d)@.
--
-- >>> unQI (qi 3 4 5 6)
-- (3,4,5,6)
-- >>> unQI (qi 1 1 (-1) 1)
-- (1,1,-1,1)
unQI :: QI -> (Integer, Integer, Integer, Integer)
unQI :: QI -> (Integer, Integer, Integer, Integer)
unQI QI
n = forall a.
QI -> (Integer -> Integer -> Integer -> Integer -> a) -> a
runQI QI
n (,,,)
{-# INLINE unQI #-}

-- | Given @n@ such that @n = a + b √c@, return @(a, b, c)@.
--
-- >>> unQI' (qi' 0.5 0.7 2)
-- (1 % 2,7 % 10,2)
-- >>> unQI' (qi' 0.3 0.4 (-3))
-- (3 % 10,2 % 5,-3)
unQI' :: QI -> (Rational, Rational, Integer)
unQI' :: QI -> (Rational, Rational, Integer)
unQI' QI
n = forall a. QI -> (Rational -> Rational -> Integer -> a) -> a
runQI' QI
n (,,)
{-# INLINE unQI' #-}

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @(a, b, c, d)@.
--
-- >>> view _qi (qi 3 4 5 6)
-- (3,4,5,6)
-- >>> view _qi (qi 1 1 (-1) 1)
-- (1,1,-1,1)
--
-- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 3 4 5 6)
-- qi 13 14 15 16
-- >>> over _qi (\(a,b,c,d) -> (a+10, b+10, c+10, d+10)) (qi 1 1 (-1) 1)
-- qi 4 0 2 1
_qi :: Lens' QI (Integer, Integer, Integer, Integer)
_qi :: Lens' QI (Integer, Integer, Integer, Integer)
_qi (Integer, Integer, Integer, Integer)
-> f (Integer, Integer, Integer, Integer)
f QI
n = (\ ~(Integer
a',Integer
b',Integer
c',Integer
d') -> Integer -> Integer -> Integer -> Integer -> QI
qi Integer
a' Integer
b' Integer
c' Integer
d') forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Integer, Integer, Integer, Integer)
-> f (Integer, Integer, Integer, Integer)
f (QI -> (Integer, Integer, Integer, Integer)
unQI QI
n)
{-# INLINE _qi #-}

-- | Given a 'QI' corresponding to @n = a + b √c@, access @(a, b, c)@.
--
-- >>> view _qi' (qi' 0.5 0.7 2)
-- (1 % 2,7 % 10,2)
-- >>> view _qi' (qi' 0.3 0.4 (-3))
-- (3 % 10,2 % 5,-3)
--
-- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 3 4 5 6)
-- qi 9 10 15 90
-- >>> over _qi' (\(a,b,c) -> (a/5, b/6, c*3)) (qi 1 1 (-1) 1)
-- qi 6 5 (-3) 30
_qi' :: Lens' QI (Rational, Rational, Integer)
_qi' :: Lens' QI (Rational, Rational, Integer)
_qi' (Rational, Rational, Integer) -> f (Rational, Rational, Integer)
f QI
n = (\ ~(Rational
a',Rational
b',Integer
c') -> Rational -> Rational -> Integer -> QI
qi' Rational
a' Rational
b' Integer
c') forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Rational, Rational, Integer) -> f (Rational, Rational, Integer)
f (QI -> (Rational, Rational, Integer)
unQI' QI
n)
{-# INLINE _qi' #-}

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @(a, b, d)@.
-- Avoids having to simplify @b √c@ upon reconstruction.
--
-- >>> view _qiABD (qi 3 4 5 6)
-- (3,4,6)
-- >>> view _qiABD (qi 1 1 (-1) 1)
-- (1,1,1)
--
-- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 3 4 5 6)
-- qi 13 14 5 16
-- >>> over _qiABD (\(a,b,d) -> (a+10, b+10, d+10)) (qi 1 1 (-1) 1)
-- qi 1 1 (-1) 1
_qiABD :: Lens' QI (Integer, Integer, Integer)
_qiABD :: Lens' QI (Integer, Integer, Integer)
_qiABD (Integer, Integer, Integer) -> f (Integer, Integer, Integer)
f (QI -> (Integer, Integer, Integer, Integer)
unQI -> ~(Integer
a,Integer
b,Integer
c,Integer
d)) =
  (\ ~(Integer
a',Integer
b',Integer
d') -> Integer -> Integer -> Integer -> Integer -> QI
qiNoSimpl Integer
a' Integer
b' Integer
c Integer
d') forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Integer, Integer, Integer) -> f (Integer, Integer, Integer)
f (Integer
a,Integer
b,Integer
d)
{-# INLINE _qiABD #-}

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @a@. It is more
-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
--
-- >>> view _qiA (qi 3 4 5 6)
-- 3
-- >>> view _qiA (qi 1 1 (-1) 1)
-- 1
--
-- >>> over _qiA (+ 10) (qi 3 4 5 6)
-- qi 13 4 5 6
-- >>> over _qiA (+ 10) (qi 1 1 (-1) 1)
-- qi 11 1 (-1) 1
_qiA :: Lens' QI Integer
_qiA :: Lens' QI Integer
_qiA = Lens' QI (Integer, Integer, Integer)
_qiABD forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {f :: * -> *} {t} {a} {b} {c}.
Functor f =>
(t -> f a) -> (t, b, c) -> f (a, b, c)
go
  where go :: (t -> f a) -> (t, b, c) -> f (a, b, c)
go t -> f a
f ~(t
a,b
b,c
d) = (\a
a' -> (a
a',b
b,c
d)) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> f a
f t
a

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @b@. It is more
-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
--
-- >>> view _qiB (qi 3 4 5 6)
-- 4
-- >>> view _qiB (qi 1 1 (-1) 1)
-- 1
--
-- >>> over _qiB (+ 10) (qi 3 4 5 6)
-- qi 3 14 5 6
-- >>> over _qiB (+ 10) (qi 1 1 (-1) 1)
-- qi 1 11 (-1) 1
_qiB :: Lens' QI Integer
_qiB :: Lens' QI Integer
_qiB = Lens' QI (Integer, Integer, Integer)
_qiABD forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {f :: * -> *} {t} {b} {a} {c}.
Functor f =>
(t -> f b) -> (a, t, c) -> f (a, b, c)
go
  where go :: (t -> f b) -> (a, t, c) -> f (a, b, c)
go t -> f b
f ~(a
a,t
b,c
d) = (\b
b' -> (a
a,b
b',c
d)) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> f b
f t
b

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @c@. It is more
-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
--
-- >>> view _qiC (qi 3 4 5 6)
-- 5
-- >>> view _qiC (qi 1 1 (-1) 1)
-- -1
--
-- >>> over _qiC (+ 10) (qi 3 4 5 6)
-- qi 3 4 15 6
-- >>> over _qiC (+ 10) (qi 1 1 (-1) 1)
-- qi 4 0 2 1
_qiC :: Lens' QI Integer
_qiC :: Lens' QI Integer
_qiC = Lens' QI (Integer, Integer, Integer, Integer)
_qi forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {f :: * -> *} {t} {c} {a} {b} {d}.
Functor f =>
(t -> f c) -> (a, b, t, d) -> f (a, b, c, d)
go
  where go :: (t -> f c) -> (a, b, t, d) -> f (a, b, c, d)
go t -> f c
f ~(a
a,b
b,t
c,d
d) = (\c
c' -> (a
a,b
b,c
c',d
d)) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> f c
f t
c

-- | Given a 'QI' corresponding to @n = (a + b √c)\/d@, access @d@. It is more
-- efficient to use '_qi' or '_qiABD' when modifying multiple terms at once.
--
-- >>> view _qiD (qi 3 4 5 6)
-- 6
-- >>> view _qiD (qi 1 1 (-1) 1)
-- 1
--
-- >>> over _qiD (+ 10) (qi 3 4 5 6)
-- qi 3 4 5 16
-- >>> over _qiD (+ 10) (qi 1 1 (-1) 1)
-- qi 1 1 (-1) 11
_qiD :: Lens' QI Integer
_qiD :: Lens' QI Integer
_qiD = Lens' QI (Integer, Integer, Integer)
_qiABD forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {f :: * -> *} {t} {c} {a} {b}.
Functor f =>
(t -> f c) -> (a, b, t) -> f (a, b, c)
go
  where go :: (t -> f c) -> (a, b, t) -> f (a, b, c)
go t -> f c
f ~(a
a,b
b,t
d) = (\c
d' -> (a
a,b
b,c
d')) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> t -> f c
f t
d

-- | Check if the value is zero.
--
-- >>> map qiIsZero [qi 0 0 0 1, qi 1 0 0 1, qiSubR (qi 7 0 0 2) 3.5, qi 0 9 0 1, qi 0 0 (-1) 5]
-- [True,False,True,True,True]
qiIsZero :: QI -> Bool
qiIsZero :: QI -> Bool
qiIsZero (QI QuadExt c Rational
q) = QuadExt c Rational
q forall a. Eq a => a -> a -> Bool
== QuadExt c Rational
0
{-# INLINE qiIsZero #-}

-- | Convert a 'QI' number into a 'Floating' one.
-- Returns NaN for imaginary values.
--
-- >>> qiToFloat (qi 3 4 5 6) == ((3 + 4 * sqrt 5)/6 :: Double)
-- True
-- >>> qiToFloat (qi 2 3 (-5) 7)
-- NaN
qiToFloat :: Floating a => QI -> a
qiToFloat :: forall a. Floating a => QI -> a
qiToFloat (QI -> (Integer, Integer, Integer, Integer)
unQI -> ~(Integer
a,Integer
b,Integer
c,Integer
d)) =
  (forall a. Num a => Integer -> a
fromInteger Integer
a forall a. Num a => a -> a -> a
+ forall a. Num a => Integer -> a
fromInteger Integer
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt (forall a. Num a => Integer -> a
fromInteger Integer
c)) forall a. Fractional a => a -> a -> a
/ forall a. Num a => Integer -> a
fromInteger Integer
d
{-# INLINE qiToFloat #-}

-- | Add an 'Integer' to a 'QI'.
--
-- >>> qi 3 4 5 6 `qiAddI` 1
-- qi 9 4 5 6
-- >>> qi 1 1 (-1) 1 `qiAddI` 1
-- qi 2 1 (-1) 1
qiAddI :: QI -> Integer -> QI
qiAddI :: QI -> Integer -> QI
qiAddI (QI QuadExt c Rational
q) Integer
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
+ forall a. Num a => Integer -> a
fromInteger Integer
x)
{-# INLINE qiAddI #-}

-- | Add a 'Rational' to a 'QI'.
--
-- >>> qi 3 4 5 6 `qiAddR` 1.2
-- qi 51 20 5 30
-- >>> qi 1 1 (-1) 1 `qiAddR` 1.2
-- qi 11 5 (-1) 5
qiAddR :: QI -> Rational -> QI
qiAddR :: QI -> Rational -> QI
qiAddR (QI QuadExt c Rational
q) Rational
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
+ forall a. Fractional a => Rational -> a
fromRational Rational
x)
{-# INLINE qiAddR #-}

-- | Subtract an 'Integer' from a 'QI'.
--
-- >>> qi 3 4 5 6 `qiSubI` 1
-- qi (-3) 4 5 6
-- >>> qi 1 1 (-1) 1 `qiSubI` 1
-- qi 0 1 (-1) 1
qiSubI :: QI -> Integer -> QI
qiSubI :: QI -> Integer -> QI
qiSubI (QI QuadExt c Rational
q) Integer
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
- forall a. Num a => Integer -> a
fromInteger Integer
x)
{-# INLINE qiSubI #-}

-- | Subtract a 'Rational' from a 'QI'.
--
-- >>> qi 3 4 5 6 `qiSubR` 1.2
-- qi (-21) 20 5 30
-- >>> qi 1 1 (-1) 1 `qiSubR` 1.2
-- qi (-1) 5 (-1) 5
qiSubR :: QI -> Rational -> QI
qiSubR :: QI -> Rational -> QI
qiSubR (QI QuadExt c Rational
q) Rational
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
- forall a. Fractional a => Rational -> a
fromRational Rational
x)
{-# INLINE qiSubR #-}

-- | Multiply a 'QI' by an 'Integer'.
--
-- >>> qi 3 4 5 6 `qiMulI` 2
-- qi 3 4 5 3
-- >>> qi 1 1 (-1) 1 `qiMulI` 2
-- qi 2 2 (-1) 1
qiMulI :: QI -> Integer -> QI
qiMulI :: QI -> Integer -> QI
qiMulI (QI QuadExt c Rational
q) Integer
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
* forall a. Num a => Integer -> a
fromInteger Integer
x)
{-# INLINE qiMulI #-}

-- | Multiply a 'QI' by a 'Rational'.
--
-- >>> qi 3 4 5 6 `qiMulR` 0.5
-- qi 3 4 5 12
-- >>> qi 1 1 (-1) 1 `qiMulR` 0.5
-- qi 1 1 (-1) 2
qiMulR :: QI -> Rational -> QI
qiMulR :: QI -> Rational -> QI
qiMulR (QI QuadExt c Rational
q) Rational
x = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
* forall a. Fractional a => Rational -> a
fromRational Rational
x)
{-# INLINE qiMulR #-}

-- | Divide a 'QI' by an 'Integer'.
--
-- >>> qi 3 4 5 6 `qiDivI` 2
-- qi 3 4 5 12
-- >>> qi 1 1 (-1) 1 `qiDivI` 2
-- qi 1 1 (-1) 2
qiDivI :: QI -> Integer -> QI
qiDivI :: QI -> Integer -> QI
qiDivI (QI QuadExt c Rational
q) (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qiDivI" -> Integer
x) = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
* forall a. Fractional a => Rational -> a
fromRational (Integer
1 forall a. Integral a => a -> a -> Ratio a
% Integer
x))
{-# INLINE qiDivI #-}

-- | Divide a 'QI' by a 'Rational'.
--
-- >>> qi 3 4 5 6 `qiDivR` 0.5
-- qi 3 4 5 3
-- >>> qi 1 1 (-1) 1 `qiDivR` 0.5
-- qi 2 2 (-1) 1
qiDivR :: QI -> Rational -> QI
qiDivR :: QI -> Rational -> QI
qiDivR (QI QuadExt c Rational
q) (forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
"qiDivR" -> Rational
x) = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q forall a. Num a => a -> a -> a
* forall a. Fractional a => Rational -> a
fromRational (forall a. Fractional a => a -> a
recip Rational
x))
{-# INLINE qiDivR #-}

-- | Negate a 'QI'.
--
-- >>> qiNegate (qi 3 4 5 6)
-- qi (-3) (-4) 5 6
-- >>> qiNegate (qi 1 1 (-1) 1)
-- qi (-1) (-1) (-1) 1
qiNegate :: QI -> QI
qiNegate :: QI -> QI
qiNegate (QI QuadExt c Rational
q) = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (forall a. Num a => a -> a
negate QuadExt c Rational
q)
{-# INLINE qiNegate #-}

-- | Compute the reciprocal of a 'QI'.
--
-- >>> qiRecip (qi 5 0 0 2)
-- Just (qi 2 0 2 5)
--
-- >>> qiRecip (qi 0 1 5 2)
-- Just (qi 0 2 5 5)
--
-- >>> qiRecip (qi 1 1 (-1) 1)
-- Just (qi 1 (-1) (-1) 2)
--
-- >>> qiRecip (qi 0 0 0 1)
-- Nothing
qiRecip :: QI -> Maybe QI
qiRecip :: QI -> Maybe QI
qiRecip (QI QuadExt c Rational
0) = forall a. Maybe a
Nothing
qiRecip (QI QuadExt c Rational
q) = forall a. a -> Maybe a
Just (forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (forall a. Fractional a => a -> a
recip QuadExt c Rational
q))

-- | Add two 'QI's if the square root terms are the same or zeros.
--
-- >>> qi 3 4 5 6 `qiAdd` qi 1 0 5 1
-- Just (qi 9 4 5 6)
--
-- >>> qi 3 4 5 6 `qiAdd` qi 3 4 5 6
-- Just (qi 3 4 5 3)
--
-- >>> qi 3 4 (-5) 6 `qiAdd` qi 3 4 (-5) 6
-- Just (qi 3 4 (-5) 3)
--
-- >>> qi 0 1 5 1 `qiAdd` qi 0 1 6 1
-- Nothing
qiAdd :: QI -> QI -> Maybe QI
qiAdd :: QI -> QI -> Maybe QI
qiAdd (QI (QuadExt c Rational
q1 :: QuadExt c1 Rational)) (QI (QuadExt c Rational
q2 :: QuadExt c2 Rational)) =
  case forall (a :: SquareFree) (b :: SquareFree).
(KnownSquareFree a, KnownSquareFree b) =>
Proxy a -> Proxy b -> Maybe (a :~: b)
sameSquareFree (forall {k} (t :: k). Proxy t
Proxy :: Proxy c1) (forall {k} (t :: k). Proxy t
Proxy :: Proxy c2) of
    Maybe (c :~: c)
Nothing   -> forall a. Maybe a
Nothing
    Just c :~: c
Refl -> forall a. a -> Maybe a
Just (forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q1 forall a. Num a => a -> a -> a
+ QuadExt c Rational
q2))

-- | Subtract two 'QI's if the square root terms are the same or zeros.
--
-- >>> qi 3 4 5 6 `qiSub` qi 1 0 5 1
-- Just (qi (-3) 4 5 6)
--
-- >>> qi 3 4 5 6 `qiSub` qi 3 4 5 6
-- Just (qi 0 0 5 1)
--
-- >>> qi 3 4 (-5) 6 `qiSub` qi 3 4 (-5) 6
-- Just (qi 0 0 (-5) 1)
--
-- >>> qi 0 1 5 1 `qiSub` qi 0 1 6 1
-- Nothing
qiSub :: QI -> QI -> Maybe QI
qiSub :: QI -> QI -> Maybe QI
qiSub (QI (QuadExt c Rational
q1 :: QuadExt c1 Rational)) (QI (QuadExt c Rational
q2 :: QuadExt c2 Rational)) =
  case forall (a :: SquareFree) (b :: SquareFree).
(KnownSquareFree a, KnownSquareFree b) =>
Proxy a -> Proxy b -> Maybe (a :~: b)
sameSquareFree (forall {k} (t :: k). Proxy t
Proxy :: Proxy c1) (forall {k} (t :: k). Proxy t
Proxy :: Proxy c2) of
    Maybe (c :~: c)
Nothing   -> forall a. Maybe a
Nothing
    Just c :~: c
Refl -> forall a. a -> Maybe a
Just (forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q1 forall a. Num a => a -> a -> a
- QuadExt c Rational
q2))

-- | Multiply two 'QI's if the square root terms are the same or zeros.
--
-- >>> qi 3 4 5 6 `qiMul` qi 0 0 5 1
-- Just (qi 0 0 5 1)
--
-- >>> qi 3 4 5 6 `qiMul` qi 1 0 5 1
-- Just (qi 3 4 5 6)
--
-- >>> qi 3 4 5 6 `qiMul` qi 3 4 5 6
-- Just (qi 89 24 5 36)
--
-- >>> qi 3 4 (-5) 6 `qiMul` qi 3 4 (-5) 6
-- Just (qi (-71) 24 (-5) 36)
--
-- >>> qi 0 1 5 1 `qiMul` qi 0 1 6 1
-- Nothing
qiMul :: QI -> QI -> Maybe QI
qiMul :: QI -> QI -> Maybe QI
qiMul (QI (QuadExt c Rational
q1 :: QuadExt c1 Rational)) (QI (QuadExt c Rational
q2 :: QuadExt c2 Rational)) =
  case forall (a :: SquareFree) (b :: SquareFree).
(KnownSquareFree a, KnownSquareFree b) =>
Proxy a -> Proxy b -> Maybe (a :~: b)
sameSquareFree (forall {k} (t :: k). Proxy t
Proxy :: Proxy c1) (forall {k} (t :: k). Proxy t
Proxy :: Proxy c2) of
    Maybe (c :~: c)
Nothing   -> forall a. Maybe a
Nothing
    Just c :~: c
Refl -> forall a. a -> Maybe a
Just (forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q1 forall a. Num a => a -> a -> a
* QuadExt c Rational
q2))

-- | Divide two 'QI's if the square root terms are the same or zeros.
--
-- >>> qi 3 4 5 6 `qiDiv` qi 0 0 5 1
-- Nothing
--
-- >>> qi 3 4 5 6 `qiDiv` qi 1 0 5 1
-- Just (qi 3 4 5 6)
--
-- >>> qi 3 4 5 6 `qiDiv` qi 3 4 5 6
-- Just (qi 1 0 5 1)
--
-- >>> qi 3 4 5 6 `qiDiv` qi 0 1 5 1
-- Just (qi 20 3 5 30)
--
-- >>> qi 3 4 (-5) 6 `qiDiv` qi 0 1 (-5) 1
-- Just (qi 20 (-3) (-5) 30)
--
-- >>> qi 0 1 5 1 `qiDiv` qi 0 1 6 1
-- Nothing
qiDiv :: QI -> QI -> Maybe QI
qiDiv :: QI -> QI -> Maybe QI
qiDiv QI
_ (QI QuadExt c Rational
0) = forall a. Maybe a
Nothing
qiDiv (QI (QuadExt c Rational
q1 :: QuadExt c1 Rational)) (QI (QuadExt c Rational
q2 :: QuadExt c2 Rational)) =
  case forall (a :: SquareFree) (b :: SquareFree).
(KnownSquareFree a, KnownSquareFree b) =>
Proxy a -> Proxy b -> Maybe (a :~: b)
sameSquareFree (forall {k} (t :: k). Proxy t
Proxy :: Proxy c1) (forall {k} (t :: k). Proxy t
Proxy :: Proxy c2) of
    Maybe (c :~: c)
Nothing   -> forall a. Maybe a
Nothing
    Just c :~: c
Refl -> forall a. a -> Maybe a
Just (forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (QuadExt c Rational
q1 forall a. Fractional a => a -> a -> a
/ QuadExt c Rational
q2))

-- | Exponentiate a 'QI' to an 'Integer' power.
--
-- >>> qi 3 4 5 6 `qiPow` 0
-- qi 1 0 5 1
--
-- >>> qi 3 4 5 6 `qiPow` 1
-- qi 3 4 5 6
--
-- >>> qi 3 4 5 6 `qiPow` 2
-- qi 89 24 5 36
--
-- >>> qi 1 1 (-1) 1 `qiPow` 4
-- qi (-4) 0 (-1) 1
qiPow :: QI -> Integer -> QI
qiPow :: QI -> Integer -> QI
qiPow (QI QuadExt c Rational
q) Integer
pow = forall (c :: SquareFree).
KnownSquareFree c =>
QuadExt c Rational -> QI
QI (forall a. Product a -> a
getProduct (forall a b. (Semigroup a, Integral b) => b -> a -> a
stimes Integer
pow (forall a. a -> Product a
Product QuadExt c Rational
q)))

-- | Compute the floor of a 'QI'.
-- Throws an error for imaginary values.
--
-- >>> qiFloor (qi 10 0 0 2)
-- 5
--
-- >>> qiFloor (qi 10 2 2 2)
-- 6
--
-- >>> qiFloor (qi 10 2 5 2)
-- 7
--
-- >>> qiFloor (qi 1 1 (-1) 1)
-- *** Exception: ...
-- ...
qiFloor :: QI -> Integer
qiFloor :: QI -> Integer
qiFloor (QI -> (Integer, Integer, Integer, Integer)
unQI -> ~(Integer
a,Integer
b,Integer
c,Integer
d)) =
  -- n = (a + b √c)/d
  -- n d = a + b √c
  -- n d = a + signum b · √(b² c)
  Integer
n_d forall a. Integral a => a -> a -> a
`div` Integer
d
  where
    n_d :: Integer
n_d = Integer
a forall a. Num a => a -> a -> a
+ forall a. Ord a => a -> a -> a
min (forall a. Num a => a -> a
signum Integer
b forall a. Num a => a -> a -> a
* Integer
b2cLow) (forall a. Num a => a -> a
signum Integer
b forall a. Num a => a -> a -> a
* Integer
b2cHigh)

    ~(Integer
b2cLow, Integer
b2cHigh) = Integer -> (Integer, Integer)
iSqrtBounds (Integer
bforall a. Num a => a -> a -> a
*Integer
b forall a. Num a => a -> a -> a
* Integer
c)

-- | Convert a (possibly periodic) simple continued fraction to a 'QI'.
--
-- @[2; 2] = 2 + 1\/2 = 5\/2@.
--
-- >>> continuedFractionToQI (2,CycList [2] [])
-- qi 5 0 2 2
--
-- The golden ratio is @[1; 1, 1, …]@.
--
-- >>> showCReal 1000 (qiToFloat (continuedFractionToQI (1,CycList [] [1])))
-- "1.6180339887498948482045868343656381177203091798057628621354486227052604628189024497072072041893911374847540880753868917521266338622235369317931800607667263544333890865959395829056383226613199282902678806752087668925017116962070322210432162695486262963136144381497587012203408058879544547492461856953648644492410443207713449470495658467885098743394422125448770664780915884607499887124007652170575179788341662562494075890697040002812104276217711177780531531714101170466659914669798731761356006708748071013179523689427521948435305678300228785699782977834784587822891109762500302696156170025046433824377648610283831268330372429267526311653392473167111211588186385133162038400522216579128667529465490681131715993432359734949850904094762132229810172610705961164562990981629055520852479035240602017279974717534277759277862561943208275051312181562855122248093947123414517022373580577278616008688382952304592647878017889921990270776903895321968198615143780314997411069260886742962267575605231727775203536139362"
--
-- >>> continuedFractionToQI (0,CycList [83,78,65,75,69] [32,66,65,68,71,69,82])
-- qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886
continuedFractionToQI :: (Integer, CycList Integer) -> QI
continuedFractionToQI :: (Integer, CycList Integer) -> QI
continuedFractionToQI (Integer
i0_, CycList Integer
is_) = QI -> Integer -> QI
qiAddI (CycList Integer -> QI
go CycList Integer
is_) Integer
i0_
  where
    go :: CycList Integer -> QI
go (CycList [Integer]
as [Integer]
bs) = [Integer] -> QI -> QI
goNonCyc [Integer]
as (if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [Integer]
bs then Integer -> Integer -> Integer -> Integer -> QI
qi Integer
0 Integer
0 Integer
0 Integer
1 else [Integer] -> QI
goCyc [Integer]
bs)

    goNonCyc :: [Integer] -> QI -> QI
goNonCyc ((Integer -> Integer
pos -> Integer
i):[Integer]
is) QI
final = QI -> QI
sudoQIRecip (QI -> Integer -> QI
qiAddI ([Integer] -> QI -> QI
goNonCyc [Integer]
is QI
final) Integer
i)
    goNonCyc []              QI
final = QI
final

    goCyc :: [Integer] -> QI
goCyc [Integer]
is = QI -> QI
sudoQIRecip ([Integer] -> QI
solvePeriodic [Integer]
is)

    -- x = (a x + b) / (c x + d)
    -- x (c x + d) = a x + b
    -- c x² + d x = a x + b
    -- c x² + (d − a) x − b = 0
    -- Apply quadratic formula, positive solution only.
    solvePeriodic :: [Integer] -> QI
solvePeriodic [Integer]
is =
      case [Integer] -> (Integer, Integer, Integer, Integer)
solvePeriodic' [Integer]
is of
        ~(Integer
a,Integer
b,Integer
c,Integer
d) ->
          Integer
a seq :: forall a b. a -> b -> b
`seq` Integer
b seq :: forall a b. a -> b -> b
`seq` Integer
c seq :: forall a b. a -> b -> b
`seq` Integer
d seq :: forall a b. a -> b -> b
`seq`
            Integer -> Integer -> Integer -> QI
qfPos Integer
c (Integer
d forall a. Num a => a -> a -> a
- Integer
a) (forall a. Num a => a -> a
negate Integer
b)
      where
        qfPos :: Integer -> Integer -> Integer -> QI
qfPos Integer
i Integer
j Integer
k = Integer -> Integer -> Integer -> Integer -> QI
qi (forall a. Num a => a -> a
negate Integer
j) Integer
1 (Integer
jforall a. Num a => a -> a -> a
*Integer
j forall a. Num a => a -> a -> a
- Integer
4forall a. Num a => a -> a -> a
*Integer
iforall a. Num a => a -> a -> a
*Integer
k) (Integer
2forall a. Num a => a -> a -> a
*Integer
i)

    -- i + 1/((a x + b) / (c x + d))      =
    -- i + (c x + d)/(a x + b)            =
    -- ((a i x + b i + c x + d)/(a x + b) =
    -- ((a i + c) x + (b i + d))/(a x + b)
    solvePeriodic' :: [Integer] -> (Integer, Integer, Integer, Integer)
solvePeriodic' ((Integer -> Integer
pos -> Integer
i):[Integer]
is) =
      case [Integer] -> (Integer, Integer, Integer, Integer)
solvePeriodic' [Integer]
is of
        ~(Integer
a,Integer
b,Integer
c,Integer
d) ->
          Integer
a seq :: forall a b. a -> b -> b
`seq` Integer
b seq :: forall a b. a -> b -> b
`seq` Integer
c seq :: forall a b. a -> b -> b
`seq` Integer
d seq :: forall a b. a -> b -> b
`seq` Integer
i seq :: forall a b. a -> b -> b
`seq`
            (Integer
aforall a. Num a => a -> a -> a
*Integer
iforall a. Num a => a -> a -> a
+Integer
c, Integer
bforall a. Num a => a -> a -> a
*Integer
iforall a. Num a => a -> a -> a
+Integer
d, Integer
a, Integer
b)

    -- x = (1 x + 0) / (0 x + 1)
    solvePeriodic' [] = (Integer
1,Integer
0,Integer
0,Integer
1)

    sudoQIRecip :: QI -> QI
sudoQIRecip QI
n =
      forall a. a -> Maybe a -> a
fromMaybe (forall a. HasCallStack => String -> a
error String
"continuedFractionToQI: Divide by zero") (QI -> Maybe QI
qiRecip QI
n)

    pos :: Integer -> Integer
pos = forall a. (Num a, Ord a, Show a) => String -> a -> a
positive String
"continuedFractionToQI"

-- | Convert a 'QI' into a (possibly periodic) simple continued fraction.
-- Throws an error for imaginary values.
--
-- @5\/2 = 2 + 1\/2 = [2; 2]@.
--
-- >>> qiToContinuedFraction (qi 5 0 0 2)
-- (2,CycList [2] [])
--
-- The golden ratio is @(1 + √5)\/2@. We can compute the corresponding PCF.
--
-- >>> qiToContinuedFraction (qi 1 1 5 2)
-- (1,CycList [] [1])
--
-- >>> qiToContinuedFraction (qi 987601513930253257378987883 1 14116473325908285531353005 81983584717737887813195873886)
-- (0,CycList [83,78,65,75,69] [32,66,65,68,71,69,82])
--
-- >>> qiToContinuedFraction (qi 1 1 (-1) 1)
-- *** Exception: ...
-- ...
qiToContinuedFraction :: QI
                      -> (Integer, CycList Integer)
qiToContinuedFraction :: QI -> (Integer, CycList Integer)
qiToContinuedFraction QI
num
  | Just ((Integer, Integer, Integer, Integer), Integer) -> Bool
isLoopQI <- forall a. Maybe (((Integer, Integer, Integer, Integer), a) -> Bool)
loopQI =
      case forall a. (a -> Bool) -> [a] -> ([a], [a])
break ((Integer, Integer, Integer, Integer), Integer) -> Bool
isLoopQI [((Integer, Integer, Integer, Integer), Integer)]
cfs of
        ([((Integer, Integer, Integer, Integer), Integer)]
preLoop, ~(((Integer, Integer, Integer, Integer), Integer)
i:[((Integer, Integer, Integer, Integer), Integer)]
postLoop)) ->
          let is :: [((Integer, Integer, Integer, Integer), Integer)]
is = forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Bool -> Bool
not forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Integer, Integer, Integer, Integer), Integer) -> Bool
isLoopQI) [((Integer, Integer, Integer, Integer), Integer)]
postLoop
          in  (Integer
i0, forall a. [a] -> [a] -> CycList a
CycList (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd [((Integer, Integer, Integer, Integer), Integer)]
preLoop) (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd forall a b. (a -> b) -> a -> b
$ ((Integer, Integer, Integer, Integer), Integer)
i forall a. a -> [a] -> [a]
: [((Integer, Integer, Integer, Integer), Integer)]
is))
  | Bool
otherwise =
      (Integer
i0, forall a. [a] -> [a] -> CycList a
CycList (forall a b. (a -> b) -> [a] -> [b]
map forall a b. (a, b) -> b
snd [((Integer, Integer, Integer, Integer), Integer)]
cfs) [])
  where
    (Integer
i0, [((Integer, Integer, Integer, Integer), Integer)]
cfs) = QI -> (Integer, [((Integer, Integer, Integer, Integer), Integer)])
qiToContinuedFractionList QI
num

    loopQI :: Maybe ((QITuple,a) -> Bool)
    loopQI :: forall a. Maybe (((Integer, Integer, Integer, Integer), a) -> Bool)
loopQI = forall s a. State s a -> s -> a
evalState (forall {m :: * -> *} {b} {b} {b}.
(Monad m, Ord b) =>
[(b, b)] -> StateT (Set b) m (Maybe ((b, b) -> Bool))
go [((Integer, Integer, Integer, Integer), Integer)]
cfs) forall a. Set a
Set.empty
      where
        go :: [(b, b)] -> StateT (Set b) m (Maybe ((b, b) -> Bool))
go ((b
n,b
_) : [(b, b)]
xs) = do
          Bool
haveSeen <- forall (m :: * -> *) s a. Monad m => (s -> a) -> StateT s m a
gets (forall a. Ord a => a -> Set a -> Bool
Set.member b
n)
          forall (m :: * -> *) s. Monad m => (s -> s) -> StateT s m ()
modify (forall a. Ord a => a -> Set a -> Set a
Set.insert b
n)
          if Bool
haveSeen
            then forall (m :: * -> *) a. Monad m => a -> m a
return (forall a. a -> Maybe a
Just ((forall a. Eq a => a -> a -> Bool
== b
n) forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fst))
            else [(b, b)] -> StateT (Set b) m (Maybe ((b, b) -> Bool))
go [(b, b)]
xs
        go [] = forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Maybe a
Nothing

qiToContinuedFractionList :: QI -> (Integer, [(QITuple, Integer)])
qiToContinuedFractionList :: QI -> (Integer, [((Integer, Integer, Integer, Integer), Integer)])
qiToContinuedFractionList QI
num =
  case Maybe QI -> [((Integer, Integer, Integer, Integer), Integer)]
go (forall a. a -> Maybe a
Just QI
num) of
    -- There is always a first number.
    ~(((Integer, Integer, Integer, Integer)
_,Integer
i) : [((Integer, Integer, Integer, Integer), Integer)]
xs) -> (Integer
i, [((Integer, Integer, Integer, Integer), Integer)]
xs)
  where
    go :: Maybe QI -> [((Integer, Integer, Integer, Integer), Integer)]
go (Just QI
n) = (QI -> (Integer, Integer, Integer, Integer)
unQI QI
n, Integer
i) forall a. a -> [a] -> [a]
: Maybe QI -> [((Integer, Integer, Integer, Integer), Integer)]
go (QI -> Maybe QI
qiRecip (QI -> Integer -> QI
qiSubI QI
n Integer
i))
      where i :: Integer
i = QI -> Integer
qiFloor QI
n
    go Maybe QI
Nothing  = []

-- | Compute a rational partial evaluation of a simple continued fraction.
--
-- Rational approximations that converge toward φ:
--
-- >>> [ continuedFractionApproximate n (1, repeat 1) | n <- [0,3..18] ]
-- [1 % 1,5 % 3,21 % 13,89 % 55,377 % 233,1597 % 987,6765 % 4181]
continuedFractionApproximate :: Foldable f
                             => Int -> (Integer, f Integer) -> Rational
continuedFractionApproximate :: forall (f :: * -> *).
Foldable f =>
Int -> (Integer, f Integer) -> Rational
continuedFractionApproximate Int
n (Integer
i0, forall (t :: * -> *) a. Foldable t => t a -> [a]
toList -> [Integer]
is) =
  forall a. Num a => Integer -> a
fromInteger Integer
i0 forall a. Num a => a -> a -> a
+
    forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (\(Integer -> Integer
pos -> Integer
i) Rational
r -> forall a. Fractional a => a -> a
recip (forall a. Num a => Integer -> a
fromInteger Integer
i forall a. Num a => a -> a -> a
+ Rational
r)) Rational
0 (forall a. Int -> [a] -> [a]
take Int
n [Integer]
is)
  where
    pos :: Integer -> Integer
pos = forall a. (Num a, Ord a, Show a) => String -> a -> a
positive String
"continuedFractionApproximate"

iSqrtBounds :: Integer -> (Integer, Integer)
iSqrtBounds :: Integer -> (Integer, Integer)
iSqrtBounds Integer
n = (Integer
low, Integer
high)
  where
    low :: Integer
low = forall a. Integral a => a -> a
integerSquareRoot Integer
n
    high :: Integer
high | Integer
lowforall a. Num a => a -> a -> a
*Integer
low forall a. Eq a => a -> a -> Bool
== Integer
n = Integer
low
         | Bool
otherwise    = Integer
low forall a. Num a => a -> a -> a
+ Integer
1

positive :: (Num a, Ord a, Show a) => String -> a -> a
positive :: forall a. (Num a, Ord a, Show a) => String -> a -> a
positive String
name = forall a. Show a => String -> String -> (a -> Bool) -> a -> a
validate String
name String
"positive" (forall a. Ord a => a -> a -> Bool
> a
0)
{-# INLINE positive #-}

nonZero :: (Num a, Eq a, Show a) => String -> a -> a
nonZero :: forall a. (Num a, Eq a, Show a) => String -> a -> a
nonZero String
name = forall a. Show a => String -> String -> (a -> Bool) -> a -> a
validate String
name String
"non-zero" (forall a. Eq a => a -> a -> Bool
/= a
0)
{-# INLINE nonZero #-}

validate :: Show a => String -> String -> (a -> Bool) -> a -> a
validate :: forall a. Show a => String -> String -> (a -> Bool) -> a -> a
validate String
name String
expected a -> Bool
f a
a
  | a -> Bool
f a
a = a
a
  | Bool
otherwise =
      forall a. HasCallStack => String -> a
error (String
"Numeric.QuadraticIrrational." forall a. [a] -> [a] -> [a]
++ String
name forall a. [a] -> [a] -> [a]
++ String
": Got " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show a
a
              forall a. [a] -> [a] -> [a]
++ String
", expected " forall a. [a] -> [a] -> [a]
++ String
expected)
{-# INLINE validate #-}