{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LINE 1 "Quipper/Algorithms/CL/Types.hs" #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}

{-# LANGUAGE IncoherentInstances #-}
{-# LANGUAGE DeriveDataTypeable  #-}

-- | This module defines the specialized datatypes of the Class Number algorithm, and basic utility functions on these types.
module Quipper.Algorithms.CL.Types where

import Quipper
import Quipper.Internal
import Data.Typeable
import Data.Ratio

import Quipper.Libraries.Arith hiding (q_mult_param)
import Quipper.Libraries.FPReal
import Quipper.Algorithms.CL.Auxiliary

-- ===========================================
-- * Type synonyms

-- $ First, we define some type synonyms for arithmetic types, selecting which will be used in the functions for the Class Number algorithm.
--
-- We use three different integer types.  For interfacing with quantum computation, we use 'CLInt' := 'IntM'.  For efficient classical (i.e. circuit-generation time) computation on potentially large integers, we use 'CLIntP' := 'Integer', Haskell’s arbitrary-precision integers.  (Δ, for instance, is taken to be a 'CLIntP').  For small classical integers (typically for register sizes), we use 'Int', Haskell’s bounded-precision integers.
--
-- For the first two of these, we define type synonyms, so that they can be swapped out to other types if desired (they are to a large extent modular).  For 'Int' we do not, since we make un-coerced use of built-in Haskell functions like 'length' which give it specifically.
--
-- Where not dictated by these conventions, integer types are generalized, i.e., @(Integral a) =>@ …
--
-- Rational and real numbers have not yet been similarly stratified.

-- | Integers that may be passed into or received out of quantum computations.
type CLInt      = IntM

-- | Integers that will be used for parameter computation only, potentially large.
type CLIntP     = Integer

-- | Rational numbers for the Class Number code.
type CLRational = Rational

-- | Real numbers for the Class Number code.
type CLReal     = FPReal

-- ===========================================
-- * Algebraic number fields

-- ===========================================
-- ** Discriminants

-- $ The functions of this subsection are needed only for circuit-generation-time classical computation, not for quantum circuit computation.

-- | Compute Δ, given /d/.
-- (See [Jozsa 2003], Prop. 6 et seq.  We use Δ, or in code /bigD/,
-- where Jozsa uses /D/.)
bigD_of_d :: Integral a => a -> a
bigD_of_d d = case (d `mod` 4) of
  1 -> d
  _ -> 4*d

-- | Compute /d/, given Δ.
-- (Again, see [Jozsa 2003], Prop. 6 et seq.)
d_of_bigD :: Integral a => a -> a
d_of_bigD bigD = case (bigD `mod` 4) of
  0 -> bigD `div` 4
  _ -> bigD

-- | Check if /d/ is a valid input to Hallgren’s algorithm,
-- i.e. correctly defines a real quadratic number field.
is_valid_d :: (Integral a) => a -> Bool
is_valid_d d = d > 1 && is_square_free d

-- | Check if /Δ/ is a valid input to Hallgren’s algorithm,
-- i.e. is the discriminant of a real quadratic number field.
-- (Cf. <http://en.wikipedia.org/wiki/Fundamental_discriminant>)
is_valid_bigD :: (Integral a) => a -> Bool
is_valid_bigD bigD = bigD > 1 && case (bigD `mod` 4) of
  1 -> is_square_free bigD
  0 -> (d `mod` 4 == 2 || d `mod` 4 == 3) && is_square_free d
    where d = bigD `div` 4
  _ -> False

-- | The (infinite, lazy) list of all valid inputs /d/, 
-- i.e. of all square-free integers above 2.
all_small_ds :: (Integral int) => [int]
all_small_ds = filter (\n -> is_valid_d n)  [2..]

-- | The (infinite, lazy) list of all valid inputs Δ, 
-- i.e. of all discriminants of real quadratic number fields.
all_bigDs :: (Integral int) => [int]
all_bigDs = map bigD_of_d all_small_ds

-- ===========================================
-- ** Field elements

-- | A data type describing a number in the algebraic number field K =
-- ℚ[√Δ]: @'AlgNum' /a/ /b/ Δ@ represents /a/ + /b/√Δ.
--
-- In general, the type of coefficients may be any type of (classical
-- or quantum) numbers, i.e. an instance of the 'Num' or 'QNum' class.
-- Given this, the algebraic numbers with a fixed Δ will in turn be an
-- instance of 'Num' or 'QNum'.
-- 
-- A value @/a/ :: /x/@ may also be used as an @'AlgNumGen' /x/@, with
-- no Δ specified, to represent simply /a/ + 0√Δ; this can be
-- considered polymorphic over all possible values of Δ.
-- 
-- This is similar to the use of 'IntM's or 'FPReal's of indeterminate
-- size, although unlike for them, we do not restrict this to the
-- classical case.  However, the question of whether an 'AlgNum' has
-- specified √Δ is (like e.g. the length of a list) is a parameter
-- property, known at circuit generation time, not a purely quantum
-- property.
data AlgNumGen a = AlgNum a a CLIntP | AlgNum_indet a deriving (Show)

-- | The specific instance of 'AlgNumGen' used for classical
-- (parameter) computation.
type AlgNum = AlgNumGen CLRational

-- | Extract the first co-ordinate of an 'AlgNumGen'
fst_AlgNum :: AlgNumGen a -> a
fst_AlgNum (AlgNum u _ _) = u
fst_AlgNum (AlgNum_indet u) = u

-- | Extract the second co-ordinate of an 'AlgNumGen'
snd_AlgNum :: (Num a) => AlgNumGen a -> a
snd_AlgNum (AlgNum _ v _) = v
snd_AlgNum (AlgNum_indet _) = 0

instance (Eq a, Num a) => Eq (AlgNumGen a) where
  (AlgNum a b bigD) == (AlgNum a' b' bigD') =
    if bigD == bigD' then a == a' && b == b'
    else error "Operation = on AlgNum: operands must have same Δ."
  (AlgNum a b bigD) == (AlgNum_indet a') = (AlgNum a b bigD) == (AlgNum a' 0 bigD)
  (AlgNum_indet a) == (AlgNum a' b' bigD') = (AlgNum a 0 bigD') == (AlgNum a' b' bigD')
  (AlgNum_indet a) == (AlgNum_indet a') = a == a'

-- | Print an algebraic number in human-readable (though not
-- Haskell-readable) format, as e.g. /a/ + /b/√Δ.
pretty_show_AlgNum :: Show a => AlgNumGen a -> String
pretty_show_AlgNum (AlgNum a b bigD) = (show a) ++ " + " ++ (show b) ++ " √" ++ show bigD
pretty_show_AlgNum (AlgNum_indet a) = show a

-- | Realize an algebraic number as a real number (of any 'Floating' type).
floating_of_AlgNum :: (Real a, Floating b) => AlgNumGen a -> b
floating_of_AlgNum (AlgNum a b bigD) = (realToFrac a) + (realToFrac b) * (sqrt $ fromIntegral bigD)
floating_of_AlgNum (AlgNum_indet a) = (realToFrac a)

-- | Coerce one algebraic number into the field of a second, if possible.  If not possible (i.e. if their Δ’s mismatch), throw an error.
number_promote :: Num a => AlgNumGen a -> AlgNumGen b -> ErrMsg -> AlgNumGen a
number_promote (AlgNum a b bigD) (AlgNum _ _ bigD') e =
  if bigD == bigD' then AlgNum a b bigD
  else error $ e "mismatched Δ."
number_promote (AlgNum_indet a) (AlgNum _ _ bigD') _ = AlgNum a 0 bigD'
number_promote n (AlgNum_indet _) _ = n

instance (Ord a, Num a) => Ord (AlgNumGen a) where
  compare (AlgNum a b bigD) (AlgNum a' b' bigD') =
    if bigD == bigD' then
      case (compare a a', compare b b') of
        (EQ,y) -> y
        (x,EQ) -> x
        (GT,GT) -> GT
        (LT,LT) -> LT
        (GT,LT) -> compare ((a-a')^2) ((b-b')^2 * fromInteger bigD)
        (LT,GT) -> compare ((b-b')^2 * fromInteger bigD) ((a-a')^2)
    else
      error "compare // AlgNumGen: mismatched Δ."
  compare (AlgNum a b bigD) (AlgNum_indet a') = compare (AlgNum a b bigD) (AlgNum a' 0 bigD)
  compare (AlgNum_indet a) (AlgNum a' b' bigD') = compare (AlgNum a 0 bigD') (AlgNum a' b' bigD')
  compare (AlgNum_indet a) (AlgNum_indet a') = compare a a'

instance (Ord a, Num a) => Num (AlgNumGen a) where
  (AlgNum a b bigD) + (AlgNum a' b' bigD') =
    if bigD == bigD' then AlgNum (a+a') (b+b') bigD
    else error "Operation + on AlgNum: operands must have same Δ."
  (AlgNum a b bigD) + (AlgNum_indet a') = (AlgNum a b bigD) + (AlgNum a' 0 bigD)
  (AlgNum_indet a) + (AlgNum a' b' bigD') = (AlgNum a 0 bigD') + (AlgNum a' b' bigD')
  (AlgNum_indet a) + (AlgNum_indet a') = (AlgNum_indet (a + a'))

  (AlgNum a b bigD) * (AlgNum a' b' bigD') =
    if bigD == bigD' then AlgNum (a*a' + b*b'*(fromIntegral bigD)) (a*b' + a'*b) bigD
    else error "Operation * on AlgNum: operands must have same Δ."
  (AlgNum a b bigD) * (AlgNum_indet a') = (AlgNum a b bigD) * (AlgNum a' 0 bigD)
  (AlgNum_indet a) * (AlgNum a' b' bigD') = (AlgNum a 0 bigD') * (AlgNum a' b' bigD')
  (AlgNum_indet a) * (AlgNum_indet a') = (AlgNum_indet (a * a'))

  (AlgNum a b bigD) - (AlgNum a' b' bigD') =
    if bigD == bigD' then AlgNum (a-a') (b-b') bigD
    else error "Operation - on AlgNum: operands must have same Δ."
  (AlgNum a b bigD) - (AlgNum_indet a') = (AlgNum a b bigD) - (AlgNum a' 0 bigD)
  (AlgNum_indet a) - (AlgNum a' b' bigD') = (AlgNum a 0 bigD') - (AlgNum a' b' bigD')
  (AlgNum_indet a) - (AlgNum_indet a') = (AlgNum_indet (a - a'))

  abs n = if (n >= 0) then n else -n
  signum n = number_promote (if n > 0 then 1 else if n == 0 then 0 else (-1)) n
               (const "CL.Types: internal error (signum // AlgNum)")
  fromInteger = AlgNum_indet . fromInteger

instance (Real a) => Real (AlgNumGen a) where
  toRational = toRational . floating_of_AlgNum

instance (Ord a, Fractional a) => Fractional (AlgNumGen a) where
  fromRational = AlgNum_indet . fromRational

  recip (AlgNum a b bigD) =
    let c = (a^2) - (b^2 * (fromIntegral bigD))
    in assert (c /= 0) (if (a == 0 && b == 0) then "CL.Types: divide-by-zero error"
                        else if is_valid_bigD bigD then "CL.Types: internal error (AlgNum // recip)"
                        else error "CL.Types: " ++ show bigD ++ " not a valid discriminant")
       (AlgNum (a/c) (-b/c) bigD)
  recip (AlgNum_indet a) = AlgNum_indet $ recip a

instance (RealFrac a) => RealFrac (AlgNumGen a) where
  properFraction x = (x',x - fromIntegral x')
    where x' = truncate $ floating_of_AlgNum x

-- | The algebraic conjugate: sends /a/ + /b/ √Δ to /a/ - /b/ √Δ.
conjugate :: (Num a) => AlgNumGen a -> AlgNumGen a
conjugate (AlgNum a b bigD) = (AlgNum a (-b) bigD)
conjugate (AlgNum_indet a) = (AlgNum_indet a)

-- | Test whether an algebraic number is an algebraic integer.
--
-- (A number is an algebraic integer iff it can be written in the form /m/ + /n/(Δ + √Δ)\/2, where /m/, /n/ are integers.
-- See [Jozsa 2003], proof of Prop. 14.)
is_alg_int :: (Ord a, RealFrac a) => AlgNumGen a -> Bool
is_alg_int (AlgNum a b bigD) = is_int n && is_int m
  where
-- solve for m, n in the equation [a + b √D = m + n(Δ + √Δ)/2]
    n = 2 * b
    m = a - b * fromIntegral bigD
is_alg_int (AlgNum_indet a) = is_int a

-- | Test whether an algebraic number is a unit of the ring of algebraic integers.
is_unit :: (Ord a, RealFrac a) => AlgNumGen a -> Bool
is_unit n = if n == 0 then False else (is_alg_int n) && (is_alg_int (recip n))

-- | The number ω associated to the field /K/.
omega_of_bigD :: CLIntP -> AlgNum
omega_of_bigD bigD =
  if (bigD `mod` 4 == 1)
  then (AlgNum (1/2) (1/2) bigD)
  else (AlgNum 0 1 bigD)

-- ===========================================
-- * Ideals

-- | Data specifying an ideal in an algebraic number field.  An ideal is described by a tuple
-- (Δ,/m/,/l/,/a/,/b/), representing the ideal
-- 
-- /m/\//l/ (/aZ/ + (/b/+√Δ)\/2 /Z/),
-- 
-- where moreover we assume and ensure always that the ideal is in /standard form/ ([Jozsa 2003], p.11, Prop. 16).  Specifically,
-- 
-- * /a/,/k/,/l/ > 0;
--
-- * 4/a/ | /b/[sup 2] – Δ; 
--
-- * /b/ = τ(/a/,/b/);
--
-- * gcd(/k/,/l/) = 1
-- 
-- In particular, this gives us bounds on the size of /a/ and /b/, 
-- and hence tells us the sizes needed for these registers (see 'length_for_ab' below).
data IdealX x = Ideal CLIntP (XInt x) (XInt x) (XInt x) (XInt x)
  deriving (Show, Eq, Typeable)

-- | Classical parameter specifying an ideal.
type Ideal = IdealX Bool

-- | Quantum circuit-type counterpart of 'Ideal'.
type IdealQ = IdealX Qubit

-- | Classical circuit-type counterpart of 'Ideal'.
type IdealC = IdealX Bit

type instance QCType x y (IdealX z) = IdealX (QCType x y z)
type instance QTypeB Ideal = IdealQ

instance Show Ideal where
  show (Ideal bigD m l a b) =
    "Ideal "
    ++ show bigD ++ " "
    ++ show m ++ " "
    ++ show l ++ " "
    ++ show a ++ " "
    ++ show b

instance QCLeaf x => QCData (IdealX x) where

  qcdata_mapM ~(Ideal _ msh lsh ash bsh) f g (Ideal bigD m l a b) = do
    m' <- qcdata_mapM msh f g m
    l' <- qcdata_mapM lsh f g l
    a' <- qcdata_mapM ash f g a
    b' <- qcdata_mapM bsh f g b
    return (Ideal bigD m' l' a' b')

  qcdata_zip ~(Ideal _ msh lsh ash bsh) q c q' c' (Ideal bigD m l a b) (Ideal bigD' m' l' a' b') e
    | bigD /= bigD'
      = error (e "Ideal exponent mismatch")
    | otherwise
      = (Ideal bigD m'' l'' a'' b'')
      where
        m'' = qcdata_zip msh q c q' c' m m' errmsg
        l'' = qcdata_zip lsh q c q' c' l l' errmsg
        a'' = qcdata_zip ash q c q' c' a a' errmsg
        b'' = qcdata_zip bsh q c q' c' b b' errmsg
        errmsg x = e ("in Ideal: " ++ x)

  qcdata_promote (Ideal bigD m l a b) (Ideal bigD' m' l' a' b') e
    | bigD /= bigD'
      = error (e "Ideal exponent mismatch")
    | otherwise
      = (Ideal bigD m'' l'' a'' b'')
    where
      m'' = qcdata_promote m m' errmsg
      l'' = qcdata_promote l l' errmsg
      a'' = qcdata_promote a a' errmsg
      b'' = qcdata_promote b b' errmsg
      errmsg x = e ("in Ideal: " ++ x)

-- Labeling of IdealQ is (m,l,a,b).
instance QCLeaf x => Labelable (IdealX x) String where
  label_rec (Ideal _ qm ql qa qb) s = do
    label_rec qm s `dotted_indexed` "m"
    label_rec ql s `dotted_indexed` "l"
    label_rec qa s `dotted_indexed` "a"
    label_rec qb s `dotted_indexed` "b"

-- We also provide an alternate labeling by a 4-tuple of strings, in
-- case this is ever useful (maybe for an ideal where the components
-- are called something other than /m/, /l/, /a/, and /b/).
instance Labelable IdealQ (String, String, String, String) where
  label_rec (Ideal _ qm ql qa qb) (sm, sl, sa, sb) = do
    label_rec qm sm
    label_rec ql sl
    label_rec qa sa
    label_rec qb sb

instance Eq Ideal where
  i1@(Ideal bigD m l a b) == i2@(Ideal bigD' m' l' a' b')
    = if (bigD /= bigD')
      then error error_string
      else (m == m' && l' == l' && a == a' && b == b')
    where error_string = "Comparing two ideals of different Δ: " ++ (show i1) ++ "," ++ (show i2)

-- | Data specifying a reduced ideal, by a tuple (Δ,/a/,/b/); this 
-- corresponds to the ideal specified by (Δ,1,/a/,/a/,/b/), i.e.,
-- /Z/ + (/b/+√Δ)\/2/a/ /Z/.
data IdealRedX x = IdealRed CLIntP (XInt x) (XInt x)
  deriving (Show, Typeable)

-- | Classical parameter specifying a reduced ideal.
type IdealRed = IdealRedX Bool

-- | Quantum circuit-type counterpart of 'IdealRed'.
type IdealRedQ = IdealRedX Qubit

-- | Classical circuit-type counterpart of 'IdealRed'.
type IdealRedC = IdealRedX Bit

instance Show IdealRed where
  show (IdealRed bigD a b) =
    "IdealRed "
    ++ show bigD ++ " "
    ++ show a ++ " "
    ++ show b

instance Eq IdealRed where
  i1@(IdealRed bigD a b) == i2@(IdealRed bigD' a' b')
    = if (bigD /= bigD')
      then error error_string
      else (a == a' && b == b')
    where error_string = "Comparing two reduced ideals of different Δ: "
                         ++ (show i1) ++ "," ++ (show i2)

type instance QCType x y (IdealRedX z) = IdealRedX (QCType x y z)
type instance QTypeB IdealRed = IdealRedQ

instance QCLeaf x => QCData (IdealRedX x) where

  qcdata_mapM ~(IdealRed _ ash bsh) f g (IdealRed bigD a b) = do
    a' <- qcdata_mapM ash f g a
    b' <- qcdata_mapM bsh f g b
    return (IdealRed bigD a' b')

  qcdata_zip ~(IdealRed _ ash bsh) q c q' c' (IdealRed bigD a b) (IdealRed bigD' a' b') e
    | bigD /= bigD'
      = error (e "IdealRed exponent mismatch")
    | otherwise
      = (IdealRed bigD a'' b'')
      where
        a'' = qcdata_zip ash q c q' c' a a' errmsg
        b'' = qcdata_zip bsh q c q' c' b b' errmsg
        errmsg x = e ("in IdealRed: " ++ x)

  qcdata_promote (IdealRed bigD a b) (IdealRed bigD' a' b') e
    | bigD /= bigD'
      = error (e "IdealRed exponent mismatch")
    | otherwise
      = (IdealRed bigD a'' b'')
      where
        a'' = qcdata_promote a a' errmsg
        b'' = qcdata_promote b b' errmsg
        errmsg x = e ("in IdealRed: " ++ x)

-- Labeling of IdealRedQ is (a,b).
instance QCLeaf x => Labelable (IdealRedX x) String where
  label_rec (IdealRed _ qa qb) s = do
    label_rec qa s `dotted_indexed` "a"
    label_rec qb s `dotted_indexed` "b"

-- We also provide an alternate labeling by a pair of strings, in case
-- this is ever useful (maybe for an ideal where the two components
-- are called something other than /a/ and /b/).
instance QCLeaf x => Labelable (IdealRedX x) (String, String) where
  label_rec (IdealRed _ qa qb) (sa, sb) = do
    label_rec qa sa
    label_rec qb sb

-- | An ideal /I/, together with a distance δ for it — that is, /some/ representative, mod /R/, for δ(/I/) as defined on /G/ p.4.  
-- Most functions described as acting on ideals need in fact to be seen as a pair of an ideal and a distance for it. 
type IdDist = (Ideal,FPReal)

-- | Quantum analogue of 'IdDist'. 
type IdDistQ = (IdealQ,FPRealQ)

-- | A reduced ideal /I/, together with a distance δ for it.
type IdRedDist = (IdealRed,FPReal)

-- | Quantum analogue of 'IdRedDist'. 
type IdRedDistQ = (IdealRedQ,FPRealQ)

-- ===========================================
-- ** Trivial access functions

-- | Extract the /d/ component from an 'IdealQ'.
d_of_Ideal :: IdealX a -> CLIntP
d_of_Ideal (Ideal bigD _ _ _ _) = d_of_bigD bigD

-- | Extract the /d/ component from an 'IdealRedQ'.
d_of_IdealRed :: IdealRedX a -> CLIntP
d_of_IdealRed (IdealRed bigD _ _) = d_of_bigD bigD

-- | Extract Δ from an 'IdealQ'.
bigD_of_Ideal :: IdealX a -> CLIntP
bigD_of_Ideal (Ideal bigD _ _ _ _) = bigD

-- | Extract Δ from an 'IdealRedQ'.
bigD_of_IdealRed :: IdealRedX a -> CLIntP
bigD_of_IdealRed (IdealRed bigD _ _) = bigD

-- | Extract the delta part from an ideal/distance pair.
delta :: IdDist -> CLReal
delta (Ideal _ _ _ _ _, del) = del

-- ===========================================
-- ** Assertions, coercions

-- $ Elements of the types 'Ideal', 'IdealRed', etc are assumed to
-- satisfy certain extra conditions.  This section includes functions
-- for checking that these conditions are satisfied, and for safely
-- coercing between these types.

-- | @'tau' Δ /b/ /a/@: the function τ(/b/,/a/).  Gives the representative for /b/ mod /2a/, in a range dependent on /a/ and √Δ.  
--
-- (This doesn't quite belong here, but is included as a prerequisite of the assertions).
tau :: (Integral int, Integral int') => int' -> int -> int -> int
tau bigD b a = mod_with_max b (2*a) max
  where
    max = if a > root_bigD then a else root_bigD
    root_bigD = floor $ sqrt $ fromIntegral bigD

-- | Return 'True' if the given ideal is in standard form.  (Functions should /always/ keep ideals in standard form).
is_standard :: Ideal -> Bool
is_standard (Ideal bigD m l a b) =
  (a > 0) && (l > 0) && (m > 0)
  && ((bigD - (fromIntegral b)^2) `mod` (4 * (fromIntegral a)) == 0)
  && b == tau bigD b a

-- | Test whether an 'Ideal' is reduced.  (An ideal \</m/,/l/,/a/,/b/> is reduced iff /m/ = 1, /l/ = /a/, /b/ ≥ 0 and /b/ + √Δ > 2/a/ ([Jozsa 2003], Prop. 20)). 
is_reduced :: Ideal -> Bool
is_reduced (Ideal bigD m l a b) = (m == 1) && (l == a) && (b >= 0) && (b + root_bigD > 2 * a)
  where root_bigD = ceiling $ sqrt $ fromIntegral bigD

-- | Test whether an 'IdealRed' is really reduced.  (An ideal \<1,/a/,/a/,/b/> is reduced iff /b/ ≥ 0 and /b/ + √Δ > 2/a/ ([Jozsa 2003], Prop. 20)). 
is_really_reduced :: IdealRed -> Bool
is_really_reduced (IdealRed bigD a b) = (b >= 0) && (b + root_bigD > 2 * a)
  where root_bigD = ceiling $ sqrt $ fromIntegral bigD

-- | Coerce an 'IdealRed' to an 'Ideal'.
forget_reduced :: IdealRed -> Ideal
forget_reduced (IdealRed bigD a b) = (Ideal bigD 1 a a b)

-- | Coerce an 'Ideal' to an 'IdealRed', if it is reduced, or throw an error otherwise.  Cf. [Jozsa 2003], Prop. 20.
to_reduced :: Ideal -> IdealRed
to_reduced ii@(Ideal bigD m l a b) =
  if is_reduced ii then (IdealRed bigD a b)
                   else error $ "to_reduced: (" ++ (show ii) ++ ") is not reduced."

-- | Throw an error if an 'Ideal' is not reduced; otherwise, the identity function.
assert_reduced :: Ideal -> a -> a
assert_reduced ii =
  assert (is_reduced ii) ("assert_reduced: (" ++ (show ii) ++ ") is not reduced.")

-- | Throw an error if an 'IdealRed' is not really reduced; otherwise, the identity function.
assert_really_reduced :: IdealRed -> a -> a
assert_really_reduced ii =
  assert (is_really_reduced ii) ("assert_really_reduced: (" ++ (show ii) ++ ") is not reduced.")

-- | Quantum analogue of 'tau'.  @'q_tau' Δ /qb/ /qa/@: compute the representative for /qb/ mod 2/qa/, in a range dependent on /qa/ and √Δ.
q_tau :: CLIntP -> QDInt -> QDInt -> Circ (QDInt, QDInt, QDInt)
q_tau bigD = box ("tau, Δ = " ++ show bigD) $ \a b -> do
  let root_bigD = floor $ sqrt $ fromIntegral bigD
  t <- with_computed
    (do
      (_, a_gt_rtD) <- q_gt_param a root_bigD
      max <- qinit $ qc_false a
      (max, _) <- controlled_not max a `controlled` a_gt_rtD
      max <- bool_controlled_not max (intm_promote root_bigD max "q_tau: internal error") `controlled` (a_gt_rtD .==. False)
      (_, twice_a) <- q_mult_param 2 a
      (_, _, _, t) <- q_mod_with_max b twice_a max
      return t)
    qc_copy
  return (a,b,t)

-- | Test whether a given 'IdealQ' is reduced.  \</m/,/l/,/a/,/b/> is reduced iff /m/ = 1, /l/ = /a/, /b/ ≥ 0 and /b/ + √Δ > 2/a/ ([Jozsa 2003], Prop. 20).  
q_is_reduced :: IdealQ -> Circ (IdealQ,Qubit)
q_is_reduced = box "is_reduced" $ \qii ->
  let bigD = bigD_of_Ideal qii in
  with_computed_fun qii
    (\(Ideal bigD qm ql qa qb) -> do
      test1 <- qinit False
      test1 <- qnot test1 `controlled` qm .==. 1
      (ql,qa,test2) <- q_is_equal ql qa
      (qb,test3) <- q_ge_param qb 0
      (qa, q2a) <- q_mult_param 2 qa
      qx <- q_sub_param_in_place (ceiling $ sqrt $ fromIntegral bigD) q2a
      (qb, qx, test4) <- q_gt qb qx
      return ([test1,test2,test3,test4], (qm,ql,qa,qb,qx)))
    (\(tests, rest) -> do
      test_out <- qinit False
      test_out <- qnot test_out `controlled` tests
      return ((tests, rest), test_out))

-- | Test whether a given 'IdealQ' is really reduced (as it should always be, if code is written correctly).  An ideal \<1,/a/,/a/,/b/> is reduced iff /b/ ≥ 0 and /b/ + √Δ > 2/a/ ([Jozsa 2003], Prop. 20). 
q_is_really_reduced :: IdealRedQ -> Circ (IdealRedQ,Qubit)
q_is_really_reduced = box "is_really_reduced" $ \qii ->
  let bigD = bigD_of_IdealRed qii in
  with_computed_fun qii
    (\(IdealRed bigD qa qb) -> do
      (qb,test1) <- q_ge_param qb 0
      (qa, q2a) <- q_mult_param 2 qa
      qx <- q_sub_param_in_place (ceiling $ sqrt $ fromIntegral bigD) q2a
      (qb, qx, test2) <- q_gt qb qx
      return ([test1,test2], (qa,qb,qx)))
    (\(tests, rest) -> do
      test_out <- qinit False
      test_out <- qnot test_out `controlled` tests
      return ((tests, rest), test_out))

-- | Coerce an 'IdealRedQ' to an 'IdealQ', initializing the extra components appropriately.
q_forget_reduced :: IdealRedQ -> Circ IdealQ
q_forget_reduced = box "forget_reduced" $ \(IdealRed bigD a b) -> do
   let a_bits = qulist_of_qdint_bh a
       n = length a_bits
   m <- qinit (intm n 1)
   (a,l) <- qc_copy_fun a
   return (Ideal bigD m l a b)

-- | Coerce an 'IdealQ' to an 'IdealRedQ', assertively terminating the extra components 
-- (and hence throwing an error at quantum runtime if the input is not reduced).
q_assert_reduced :: IdealQ -> Circ IdealRedQ
q_assert_reduced = box "assert_reduced" $ \x@(Ideal bigD m l a b) -> do
  x_red <- reverse_generic q_forget_reduced (IdealRed bigD a b) x
  q_assert_really_reduced x_red

-- | Throw a (quantum-runtime) error if an 'IdealRedQ' is not really reduced; otherwise, do nothing.
--
-- Compare 'assert_reduced', 'q_is_really_reduced' in "Quipper.Algorithms.CL.RegulatorQuantum", and [Jozsa 2003] Prop. 20.
q_assert_really_reduced :: IdealRedQ -> Circ IdealRedQ
q_assert_really_reduced = box "assert_really_reduced" $ \ii -> do
  (ii,test) <- q_is_really_reduced ii
  qterm True test
  return ii

-- ======================================================================
-- ** Bounds on coefficient sizes

-- $ Given Δ, how much space should be allocated for the coefficients of ideals?  Most of these bounds are currently missing or uncertain, as documented below.  Note these bounds are intended to be sufficient for the calculations occurring in this algorithm, /not/ for representing arbitrary ideals.

-- | Given Δ, return the size of integers to be used for the coefficients /a/, /b/ of reduced ideals.
--
-- Note: can we bound this more carefully?  In reduced ideals, we always have 0 ≤ /a/,/b/ ≤ √Δ (see notes on 'is_standard', 'is_reduced'), and the outputs of ρ, ρ[sup –1] and dot-products of reduced ideals always keep |/a/| ≤ Δ.  However, intermediate calculations may involve larger values, so we allocate a little more space.  For now, this padding is a seat-of-the-pants estimate.  
length_for_ab :: CLIntP -> Int
length_for_ab bigD = 3 + (ceiling $ logBase 2 $ fromIntegral bigD)

-- | Given Δ, return the size of integers to be used for the coefficients /m/, /l/ of general ideals.
--
-- TODO: bound this!  Neither Hallgren nor [Jozsa 2003] discusses bounds on the values of /m/ and /l/ that will appear, and we do not yet have a bound.  For now we use the same length as for /a/ and /b/, for convenience; this should be considered a dummy bound, quite possibly not sufficient in general.
length_for_ml :: CLIntP -> Int
length_for_ml = length_for_ab

-- | Given Δ, return the precision /n/ = log[sub 2]/N/ to be used for
-- discretizing the quasi-periodic function /f/ to /f/[sub /N/].
--
-- (“Precision” here means the number of binary digits after the point).
-- 
-- Taken to ensure 1\//N/ < 3/(32 Δ log Δ).  (Cf. [Jozsa 2003], Prop. 36 (iii).)
n_of_bigD :: (Integral int) => CLIntP -> int
n_of_bigD bigD =
  ceiling $ logBase 2 $ 32 * (fromIntegral bigD) * (log $ fromIntegral bigD) / 3

-- | Given Δ, /n/, /l/ (as for
-- 'Quipper.Algorithms.CL.RegulatorClassical.fN',
-- 'Quipper.Algorithms.CL.RegulatorQuantum.q_fN'), return the
-- precision required for intermediate distance calculations during
-- the computation of /f/[sub /N/].
--
-- TODO: bound this more carefully.  [Jozsa 2003] asks for the final output to be precision /n/, but does not discuss intermediate precision, and we have not yet got a confident answer.  For now, just a back-of-the-envelope estimate, which should be sufficient and /O/(correct), but is almost certainly rather larger than necessary.
precision_for_fN :: CLIntP -> Int -> Int -> Int
precision_for_fN _ n l = 2 * (n + l)

-- | Set the 'IntM' coefficients of an 'Ideal' to the standard lengths, if they are not already fixed incompatibly.  The standard lengths are determined by 'length_for_ml', 'length_for_ab'.   (Compare 'intm_promote', etc.)
fix_sizes_Ideal :: Ideal -> Ideal
fix_sizes_Ideal (Ideal bigD m l a b) = (Ideal bigD (f m) (f l) (f a) (f b))
  where
    f x = intm_promote x (intm n 0) "set_sizes_Ideal: lengths already fixed incompatibly"
    n = max (length_for_ml bigD) (length_for_ab bigD)

-- | Set the 'IntM' coefficients of an 'IdealRed' to the standard lengths, if they are not already fixed incompatibly.  The standard lengths are determined by 'length_for_ml', 'length_for_ab'.   (Compare 'intm_promote', etc.)
fix_sizes_IdealRed :: IdealRed -> IdealRed
fix_sizes_IdealRed (IdealRed bigD a b) = (IdealRed bigD (f a) (f b))
  where
    f x = intm_promote x (intm n 0) "set_sizes_Ideal: lengths already fixed incompatibly"
    n = max (length_for_ml bigD) (length_for_ab bigD)