{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ViewPatterns #-}
module Numeric.QuadraticIrrational
(
QI, qi, qi', runQI, runQI', unQI, unQI'
,
_qi, _qi', _qiABD, _qiA, _qiB, _qiC, _qiD
,
qiIsZero
, qiToFloat
, qiAddI, qiSubI, qiMulI, qiDivI
, qiAddR, qiSubR, qiMulR, qiDivR
, qiNegate, qiRecip, qiAdd, qiSub, qiMul, qiDiv, qiPow
, qiFloor
,
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
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)
qi :: Integer
-> Integer
-> Integer
-> Integer
-> 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 #-}
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 #-}
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 #-}
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)
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 #-}
qi' :: Rational
-> Rational
-> Integer
-> QI
qi' :: Rational -> Rational -> Integer -> QI
qi' Rational
a Rational
b Integer
c = QI
n
where
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' #-}
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 #-}
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' #-}
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 #-}
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' #-}
_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 #-}
_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' #-}
_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 #-}
_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
_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
_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
_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
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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 #-}
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))
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))
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))
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))
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))
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)))
qiFloor :: QI -> Integer
qiFloor :: QI -> Integer
qiFloor (QI -> (Integer, Integer, Integer, Integer)
unQI -> ~(Integer
a,Integer
b,Integer
c,Integer
d)) =
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)
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)
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)
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)
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"
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
~(((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 = []
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 #-}