-- | This module implements the classical operations on ideals used in Hallgren’s
-- algorithm (including also classical versions of the quantum operations required).
 
module Quipper.Algorithms.CL.RegulatorClassical where

import Data.Maybe

import Quipper.Libraries.Arith
import Quipper.Libraries.FPReal

import Quipper.Algorithms.CL.Auxiliary
import Quipper.Algorithms.CL.Types

-- ======================================================================
-- * Basic operations on ideals

-- | @'unit_ideal' /bigD/ /l/@: the unit ideal /O/, for Δ = /bigD/, with the ideal’s coefficients given as /l/-bit integers.  
-- ([Jozsa 2003], Prop. 14.)
unit_ideal :: CLIntP -> Ideal
unit_ideal = forget_reduced . unit_idealRed 

-- | Like 'unit_ideal', but considered as a reduced ideal.
unit_idealRed :: CLIntP -> IdealRed
unit_idealRed bigD = assert (is_valid_bigD bigD) ("unit_idealRed: " ++ show bigD ++ " not valid discriminant") 
                  $ IdealRed bigD 1 (tau bigD (fromIntegral bigD) 1)

-- | The integer constant /c/ of an ideal.
--   ([Jozsa 2003], page 14 bottom: \"Since 4/a/ divides /b/[sup 2]-/D/
--   (cf. proposition 16) we introduce the integer /c/ = |/D/ − /b/[sup 2]|\/(4/a/)\")
c_of_ideal :: Ideal -> CLInt
c_of_ideal i@(Ideal bigD m l a b) =
    if (num `mod` denom == 0) then num `div` denom
                              else error error_string
    where num   = abs ((fromIntegral bigD) - b*b)
          denom = 4*a
          error_string = "rho of ideal [" ++ show i ++ "] produces non-integer c"

-- | γ(/I/) = (/b/+√Δ)\/(2/a/) for a given ideal /I/.
--   ([Jozsa 2003], Sec. 6.2.)
gamma_of_ideal :: Ideal -> AlgNum
gamma_of_ideal (Ideal bigD m l a b) = AlgNum a' b' bigD
    where
        a' = (fromIntegral b / fromIntegral (2*a)) :: CLRational
        b' = (fromIntegral 1 / fromIntegral (2*a)) :: CLRational
  -- Recall: @'AlgNum' u v bigD@ represents (u + v * √bigD).

-- | The reduction function ρ on ideals.
-- ([Jozsa 2003], Sec. 6.2.)
rho :: Ideal -> Ideal
rho ii@(Ideal bigD m l a b) = (Ideal bigD m'' l'' a' b')
    where
        -- With little algebra, one can derive these:
        m'   = m * a
        l'   = l * a'
        m''  = m' `div` (gcd m' l')
        l''  = l' `div` (gcd m' l')
        -- From  [Jozsa 2003], p.15, equation (10)
        a'   = c_of_ideal ii
        b'   = tau bigD (-b) a'

-- | The ρ[sup -1] function on ideals.  Inverse to 'rho'.
-- ([Jozsa 2003], Sec. 6.4.)
rho_inv :: Ideal -> Ideal
rho_inv (Ideal bigD m l a b) = (Ideal bigD m'' l'' a'' b'')
    where
        -- Create m and l
        m'    = m * a
        l'    = l * a''
        -- Reduce them to have no common denominator
        m''   = m' `div` (gcd m' l')
        l''   = l' `div` (gcd m' l')
        -- Calculate b* (b' below)
        b'    = tau bigD (-b) a
        -- Calculate new a and b
        a''   = ((fromIntegral bigD) - b'*b') `divchk` (4*a)
        b''   = tau bigD b' a''

-- | The ρ operation on reduced ideals.
rho_red :: IdealRed -> IdealRed
rho_red = to_reduced . rho . forget_reduced

-- | The ρ[sup –1] operation on reduced ideals.
rho_inv_red :: IdealRed -> IdealRed
rho_inv_red = to_reduced . rho_inv . forget_reduced

-- | The ρ operation on ideals-with-distance.
rho_d :: IdDist -> IdDist
rho_d (ii, del) = (rho ii, del + del_change)
  where
    gamma = gamma_of_ideal ii
    gamma_bar_by_gamma = (floating_of_AlgNum $ conjugate gamma) / (floating_of_AlgNum gamma)
    del_change = (log $ abs gamma_bar_by_gamma) / 2

-- | The ρ[sup –1] operation on ideals-with-distance.
rho_inv_d :: IdDist -> IdDist
rho_inv_d (ii, del) = (ii', del - del_change)
  where
    ii' = rho_inv ii
    gamma = gamma_of_ideal ii'
    gamma_bar_by_gamma = (floating_of_AlgNum $ conjugate gamma) / (floating_of_AlgNum gamma)
    del_change = (log $ abs gamma_bar_by_gamma) / 2

-- | The ρ operation on ideals-with-generator (i.e. pairs of an ideal /I/ and an 'AlgNum' /x/ such that /I/ is the principal ideal (/x/)).
rho_num :: (Ideal,AlgNum) -> (Ideal,AlgNum)
rho_num (ii, gen) = (rho ii, gen / (gamma_of_ideal ii))

-- | Apply ρ to an reduced-ideals-with-generator
rho_red_num :: (IdealRed,AlgNum) -> (IdealRed,AlgNum)
rho_red_num (ii, gen) = (rho_red ii, gen / (gamma_of_ideal $ forget_reduced ii))

-- ======================================================================
-- * Ideal reductions (bounded)

-- | Reduce an ideal, by repeatedly applying ρ.
reduce :: Ideal -> IdealRed
reduce ii = to_reduced $ while (not . is_reduced) rho ii

-- | Reduce an ideal within a bounded loop. Applies the ρ function
--   until the ideal is reduced. Used in 'star' and 'fJN' algorithms.
bounded_reduce :: IdDist -> IdDist
bounded_reduce k@(ideal@(Ideal bigD m l a b),dist) = 
    fst $ bounded_while condition bound func (k,0)
    where
        -- NOTE: some uncertainty regarding loop bound.
        bound             = ceiling $ bound_log + 1
        bound_log         = (logBase 2 ((fromIntegral bound_a) 
                            / (sqrt $ fromIntegral $ bigD_of_Ideal $ ideal)))
        bound_a           = 2 ^ (fromJust (intm_length $ a))

        condition (k,itr) = not $ is_reduced $ fst k
        func      (k,itr) = ((rho_d k),(itr+1))

-- | Apply a function (like ρ,ρ[sup -1],ρ[sup 2]) to an ideal, bounded
--   by 3*ln(Δ)\/2*ln(2). Execute while satisfies condition function.
bounded_step :: (IdDist -> Bool) -> (IdDist -> IdDist) -> IdDist -> IdDist
bounded_step condition step_function ideal =
    -- Execute a bounded while loop
    bounded_while condition bound step_function ideal
    where
        bound = ceiling $ (3 * (log $ fromIntegral $ bigD_of_Ideal $ fst ideal)) / (2 * log 2)

-- | Like 'bounded_step', but the condition is checked against delta of the
--   current ideal.
bounded_step_delta :: (CLReal -> Bool) -> (IdDist -> IdDist) -> IdDist -> IdDist
bounded_step_delta condition step_function ideal =
    bounded_step new_condition step_function ideal
    where new_condition = (\k -> condition (delta k))

-- ======================================================================
-- * Products of ideals

-- | The ordinary (not necessarily reduced) product of two reduced fractional ideals.
-- 
-- /I/⋅/J/ of [Jozsa 2003], Sec 7.1, following the description
-- given in Prop. 34.

-- NOTE: assumes I, J reduced.  Type should reflect this!
dot :: IdDist -> IdDist -> IdDist
dot i1@(Ideal bigD1 m1 l1 a1 b1, delta1) i2@(Ideal bigD2 m2 l2 a2 b2, delta2) =
    assert_reduced (fst i1) $
        assert_reduced (fst i2) $
            (Ideal bigD1 m l a3 b3, del)
    where
        sqrtd :: CLReal
        sqrtd = sqrt (fromIntegral bigD1)
        (k', u', v') = extended_euclid a1 a2
        (k,  x,  w ) = extended_euclid k' ((b1 + b2) `divchk` 2)
        a3 = (a1 * a2) `divchk` (k * k)
        t1 = x * u' * a1 * b2
        t2 = x * v' * a2 * b1
        t3 = w*(b1 * b2 + (fromIntegral bigD1)) `divchk` 2
        t  = (t1 + t2 + t3) `divchk` k
        b3 = tau bigD1 t a3
        m = k
        l = a3
        del = ((delta i1) + (delta i2))

-- | The dot-square /I/⋅/I/ of an ideal-with-distance /I/.
dot' :: IdDist -> IdDist
dot' i1@(Ideal bigD1 m1 l1 a1 b1, delta1) =
    assert_reduced (fst i1) $
            (Ideal bigD1 m l a3 b3, del)
    where
        sqrtd :: CLReal
        sqrtd = sqrt (fromIntegral bigD1)
        (k, u, w) = extended_euclid (abs a1) (abs b1)
        a3 = (a1 * a1) `divchk` (k * k)
        t1 = u * a1 * b1
        t3 = w*(b1 * b1 + (fromIntegral bigD1)) `divchk` 2
        t  = (t1 + t3) `divchk` k
        b3 = tau bigD1 t a3
        m = k
        -- l = a3
        l = 1
        del = ((delta i1) + (delta i1))

-- | The star-product of two ideals-with-distance.
--
-- This is /I/*/J/ of [Jozsa 2003], Sec. 7.1, defined as the first reduced
-- ideal-with-distance following /I/⋅/J/.

-- NOTE: assumes I, J reduced.  Type should reflect this!
star :: IdDist -> IdDist -> IdDist
star i j =
    if (delta k1 <= delta i_dot_j)
        then       bounded_step_delta (<= delta i_dot_j) rho_d k1
        else rho_d $ bounded_step_delta (>= delta i_dot_j) rho_inv_d k1
    where
        i_dot_j = i `dot` j
        -- FIX: Over bound (infinite loop) in reducing the following:
        --      <m:1 l:3 a:3 b:2 bigD:28 del:1.1>*<m:1 l:3 a:3 b:4 bigD:28 del:1.6>
        k1      = bounded_reduce i_dot_j

-- ======================================================================
-- * The function f[sub /N/], and variants

-- | Compute the expression i\/N + j\/L.
compute_injl :: (Integral int) => CLInt -> int -> CLInt -> int -> CLReal
compute_injl i nn j ll =
    (fromIntegral i)/(fromIntegral nn) + (fromIntegral j)/(fromIntegral ll)

-- |  @'fN' /i/ /j/ /n/ /l/ Δ@: find the minimal ideal-with-distance (/J/,δ[sub /J/]) such that δ[sub /J/] > /x/, where /x/ = /i/\//N/ + /j/\//L/, where /N/ = 2[sup /n/], /L/ = 2[sup /l/]. Return (/i/,/J/,δ[sub /J/]–/x/).  Work under the assumption that /R/ < 2[sup /s/].
--
-- This is the function /h/ of [Jozsa 2003, Section 9], discretized with precision 1\//N/ = 2[sup −/n/],
-- and perturbed by the jitter parameter /j/\//L/.
fN :: (Integral int) => CLInt -> CLInt -> int -> int -> CLIntP -> (Ideal, CLInt)
fN i j nn ll bigD =
  let ((ideal_J, _), diff) = fN_d i j nn ll bigD
  in (ideal_J, diff)

-- | Like 'fN', but returning an ideal-with-distance not just an ideal.
fN_d :: (Integral int) => CLInt -> CLInt -> int -> int -> CLIntP -> (IdDist, CLInt)
fN_d i j nn ll bigD =
    (j_star_19, floor $ (fromIntegral nn) * (toRational $ injl - (snd j_star_19)))
    where
        -- Expression "i/N + j/L" used repeatedly
        injl     = compute_injl i nn j ll

        -- Generate J1
        j1       = rho_d $ rho_d $ (unit_ideal bigD, 0)

        -- Generate Jk's (make an infinite list and take only what is needed)
        jks      = takeWhile (\jk -> (delta jk) <= injl) $
                        bounded_iterate bound_jks
                            (\jk -> jk `star` jk) j1
                   where
                       -- Bound for jk generation
                       max_i     = 2^(fromJust (intm_length $ i))
                       bound_jks = ceiling $
                                    (log $ fromIntegral max_i) /
                                        ((fromIntegral nn) * (delta j1))

        -- Apply all Jk's to J* using '*' in reverse (remember that the last
        -- element is J* itself)
        j_star_14   = foldr1 applyJkIfConditionHolds jks

        -- Apply Jk to J* if a condition holds.
        applyJkIfConditionHolds :: IdDist -> IdDist -> IdDist
        applyJkIfConditionHolds jk jstar =
            if (delta (jstar `star` jk) <= injl) then jstar `star` jk else jstar

        -- Go forward one step at a time as much as needed
        j_star_17   = bounded_step_delta (< injl) (rho_d.rho_d) j_star_14

        -- Go back one step if needed
        j_star_19   = if (delta (rho_inv_d j_star_17) >= injl)
                      then rho_inv_d j_star_17
                      else j_star_17

-- | Analogue of 'fN', working within the cycle determined by a given ideal /J/.
--   ([Hallgren 2006], Section 5.)
fJN :: IdDist -> CLInt -> CLInt -> CLInt -> CLInt -> CLIntP -> (Ideal, CLInt)
fJN ideal_J i j nn ll bigD =
  let ((ideal_J', _), diff) = fJN_d ideal_J i j nn ll bigD
  in (ideal_J', diff)

-- | Like 'fJN', but returning an ideal-with-distance not just an ideal.
fJN_d :: IdDist -> CLInt -> CLInt -> CLInt -> CLInt -> CLIntP -> (IdDist, CLInt)
fJN_d ideal_J i j nn ll bigD =
    (ideal_KFinal, floor $ (fromIntegral nn) * (injl - (delta ideal_KFinal)))
    where
        -- Expression "i/N + j/L"
        injl      = compute_injl i nn j ll

        -- Generate I and K
        ideal_I = fst (fN_d i j nn ll bigD)
        ideal_K = bounded_reduce (ideal_I `dot` ideal_J)

        -- Step forward/backward as much as needed
        ideal_KFinal =
            if (delta ideal_K <= injl)
            then        bounded_step_delta (<  injl) (rho_d)    ideal_K
            else rho_d $ bounded_step_delta (>= injl) (rho_inv_d) ideal_K

-- ======================================================================
-- * Classical period-finding

-- $ Functions for classically finding the regulator and fundamental unit of a field using the period of /f_N/.

-- ======================================================================
-- ** Auxiliary functions

-- | Find the order of an endofunction on an argument.  That is, @'order' /f/ /x/@ returns the first /n/ > 0 such that /f/[sup /n/]\(/x/) = /x/. 
--
--  Method: simple brute-force search/comparison.
order :: (Eq a) => (a -> a) -> a -> Int
order = order_with_projection id

-- | Given a function /p/, an endofunction /f/, and an argument /x/, returns the first /n/ > 0 such that /p/(/f/[sup /n/]\(/x/)) = /p/(/x/).  
--
-- Method: simple brute-force search/comparison.
order_with_projection :: (Eq b) => (a -> b) -> (a -> a) -> a -> Int
order_with_projection p f x = 1 + (length $ takeWhile (\y -> p y /= p x) (tail $ iterate f x))

-- | Given a function /p/, an endofunction /f/, and an argument /x/, return /f/[sup /n/]\(/x/), for the first /n/ > 0 such that /p/(/f/[sup /n/]\(/x/)) = /p/(/x/).  
--
-- Method: simple brute-force search/comparison.
first_return_with_projection :: (Eq b) => (a -> b) -> (a -> a) -> a -> a
first_return_with_projection p f x = head $ dropWhile (\y -> p y /= p x) $ tail $ iterate f x

-- | Given a bound /b/, a function /p/, an endofunction /f/, and an argument /x/, return /f/[sup /n/]\(/x/), for the first /n/ > 0 such that /p/(/f/[sup /n/]\(/x/)) = /p/(/x/), if there exists such an /n/ ≤ /b/.
--
-- Method: simple brute-force search/comparison.
first_return_with_proj_bdd :: (Eq b) => Int -> (a -> b) -> (a -> a) -> a -> Maybe a
first_return_with_proj_bdd b p f x = listToMaybe $ dropWhile (\y -> p y /= p x) $ take b $ tail $ iterate f x

-- | Find the period of a function on integers, assuming that it is periodic and injective on its period.  That is, @'period' /f/@ returns the first /n/ > 0 such that /f/(/n/) = /f/(0).  Method: simple brute-force search/comparison.
period :: (Eq a, Integral int) => (int -> a) -> int
period f = minimum [ n | n <- [1..], f n == f 0 ]
 
-- ======================================================================
-- ** Haskell native arithmetic

-- $ The functions of this section use Haskell’s native integer and floating computation.

-- | Find the regulator /R/ = log ε[sub 0] of a field, given the discriminant Δ, 
-- by finding (classically) the order of ρ.  
--
-- Uses 'IdDist' and 'rho_d'.
regulator :: CLIntP -> FPReal
regulator bigD = snd $ head $ dropWhile (\(ii,_) -> ii /= calO) $ tail $ iterate rho_d $ (calO,0)
  where calO = unit_ideal bigD

-- | Find the fundamental unit ε[sub 0] of a field, given the discriminant Δ, 
-- by finding (classically) the order of ρ.
--
-- Uses '(Ideal,Number)' and 'rho_num'.
fundamental_unit :: CLIntP -> AlgNum
fundamental_unit bigD = maximum [eps, -eps, 1/eps, -1/eps]
  where eps = snd $ first_return_with_projection fst rho_num (calO,1)
        calO = unit_ideal bigD

-- | Find the fundamental solution of Pell’s equation, given /d/.
--
-- Solutions of Pell’s equations are integer pairs (/x/,/y/) such that
-- /x/,/y/ > 0, and (/x/ + /y/√d)(/x/ – /y/√d) = 1.
--
-- In this situation, (/x/ + /y/√d) is a unit of the algebraic integers 
-- of /K/, and is >1, so we simply search the powers of ε[sub 0] for a
-- unit of the desired form.
fundamental_solution :: CLIntP -> (Integer,Integer)
fundamental_solution d = head pell_solutions
  where bigD = bigD_of_d d
        eps0 = fundamental_unit bigD
        pell_solutions = 
          [ (round x, round y) | n <- [1..],
                                 let eps@(AlgNum a b _) = eps0^n,
                                 a >= 0, b >= 0,
                                 let x = a,
                                 let y = if bigD == d then b else 2*b,
                                 is_int x, is_int y,
                                 eps * (conjugate eps) == 1]

-- ======================================================================
-- ** Fixed-precision arithmetic

-- $ The functions of this section perform period-finding using fixed-precision arithmetic.
-- This should parallel closely (though at present not exactly, due to the implementations 
-- of floating-point operations) the quantum circuit implementations, and hence allows
-- testing of whether the chosen precisions are accurate.

-- | Find the regulator /R/ = log ε[sub 0] of a field, given the discriminant Δ, 
-- by finding (classically) the order of ρ
-- using fixed-precision arithmetic: 'fix_sizes_Ideal' for 'Ideal's,
-- and given an assumed bound /b/ on log[sub 2] /R/.
--
-- Uses 'IdDist' and 'rho_d'.
regulator_fixed_prec :: Int -> CLIntP -> Maybe FPReal
regulator_fixed_prec b bigD = fmap snd (first_return_with_proj_bdd b' fst rho_d (calO,zero))
  where calO = fix_sizes_Ideal $ unit_ideal bigD
        n = n_of_bigD bigD
  -- S = 2RN, so this /i/ gives an assumed bound on log_2 S:  
        i = 1 + b + n
  -- Following precisions are as used in 'approximate_regulator_circuit', 'q_fN':
        q = 2 + 2 * i 
        l = 4 + q + n
        p = precision_for_fN bigD n l
        zero = fprealx (-p) (intm (q - n + p) 0)
  -- δ(I, ρ^2 (I)) ≥ √2 for any I, so the order of ρ is at most 2R / (√2 / 2):
        b' = ceiling $ 2 * sqrt(2) * (fromIntegral 2^b)


{-
Waiting for 'AlgNum' to be rebased using 'IntM' instead of 'Integer'.

-- | Find the fundamental unit ε[sub 0] of a field, given the discriminant Δ, 
-- by finding (classically) the order of ρ,
-- using fixed-precision arithmetic: 'fix_sizes_Ideal' for 'Ideal's,
-- and a given /l/ for the 'AlgNum's generating them.
--
-- Uses '(Ideal,Number)' and 'rho_num'.
fundamental_unit_with_fixed :: Int -> CLIntP -> AlgNum
fundamental_unit_with_fixed l bigD = maximum [eps, -eps, 1/eps, -1/eps]
  where eps = snd $ head $ dropWhile (\(ii,_) -> ii /= calO) $ tail $ iterate rho_num $ (calO,one)
        calO = fix_sizes_Ideal $ unit_ideal bigD
        one = intm_with_length (Just l) 1
-}