{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LINE 1 "Quipper/Algorithms/CL/Test.hs" #-}
-- | Test the Class Number algorithm, and its components, using classical computation

module Quipper.Algorithms.CL.Test where

import Quipper
import Quipper.Libraries.Arith
import Quipper.Libraries.FPReal
import Quipper.Algorithms.CL.Auxiliary
import Quipper.Algorithms.CL.Types
import Quipper.Algorithms.CL.RegulatorClassical
import Quipper.Algorithms.CL.CL
import Quipper.Algorithms.CL.SmithReduction
import Data.Ratio
import Data.List

-- * Sample data

-- $ Some fairly arbitrarily chosen sample elements of various types, for convenience in testing functions.

-- ** Matrices

-- | A sample square matrix
sample_matrix :: CLMatrix Integer
sample_matrix = matrix_from_list [
    [  8,  16, 16 ],
    [ 32,  6,  12 ],
    [  8, -4, -16 ]
 ]

-- | A sample non-square matrix
sample_matrix_2 :: CLMatrix Integer
sample_matrix_2 = matrix_from_list [
    [ 5, 1, 5, 253, 15, -725, 1 ],
    [ 253,2,1001,11,23,273,14079 ],
    [ 1,-185861,-28,11,91,29,-2717 ],
    [ -319,1,-19,11,3146,1,-1 ],
    [ 19285,-493,145,25,-1482,1,6647]
 ]

-- | Another sample non-square matrix
sample_matrix_3 :: CLMatrix Integer 
sample_matrix_3 = matrix_from_list [
    [ 4, 8, 4 ],
    [ 8, 4, 8 ]
 ]

-- ** Ideals and related types

-- | A sample 'CLReal'.
sample_CLReal :: Int -> FPReal
sample_CLReal l = (fprealx 0 (intm l 0))

-- | A sample 'Ideal'.
sample_Ideal :: CLIntP -> Ideal
sample_Ideal bigD =
  let l = max (length_for_ab bigD) (length_for_ml bigD)
      x = (intm l 0)
  in (Ideal bigD x x x x)

-- | A sample 'IdealQ'.
sample_IdealQ :: CLIntP -> IdealQ
sample_IdealQ = qshape . sample_Ideal

-- | A sample 'IdealRed'.
sample_IdealRed :: CLIntP -> IdealRed
sample_IdealRed bigD =
  let l = max (length_for_ab bigD) (length_for_ml bigD)
      x = (intm l 0)
  in (IdealRed bigD x x)

-- | A sample 'IdealRedQ'.
sample_IdealRedQ :: CLIntP -> IdealRedQ
sample_IdealRedQ = qshape . sample_IdealRed

-- | A sample 'IdDist'.
sample_IdDist :: CLIntP -> IdDist
sample_IdDist bigD = (sample_Ideal bigD, sample_CLReal (length_for_ab bigD))

-- | A sample 'IdDistQ'.
sample_IdDistQ :: CLIntP -> IdDistQ
sample_IdDistQ = qshape . sample_IdDist

-- | A sample 'IdRedDist'.
sample_IdRedDist :: CLIntP -> IdRedDist
sample_IdRedDist bigD = (sample_IdealRed bigD, sample_CLReal (length_for_ab bigD))

-- | A sample 'IdRedDistQ'.
sample_IdRedDistQ :: CLIntP -> IdRedDistQ
sample_IdRedDistQ = qshape . sample_IdRedDist

-- * Testing routines

-- ** Smith reduction

-- | Test the Smith Normal Form code.
test_SNF :: IO ()
test_SNF = do
    flip mapM_ [sample_matrix,sample_matrix_2,sample_matrix_3] $ \m -> do
      putStrLn $ show $ m
      putStrLn $ show $ structure_constants_from_matrix m
      putStrLn $ show $ group_order_from_matrix m
      putStrLn ""

-- ** Class group functions

-- | Classical period finding (just compare the \"next\" ideal to /O/ and see if
--   it is the same). Takes in the /O/ ideal with appropriate Δ, and returns
--   the circle length (sum δ(I)) and the list of ideals in the first iteration.
period_of_ideals :: (IdDist->IdDist) -> IdDist -> (CLReal, [IdDist])
period_of_ideals func o = (delta $ last list, list)
    where
        list = takePeriod False (iterate (\i -> func i) o)
        takePeriod :: Bool -> [IdDist] -> [IdDist]
        takePeriod got_first_o [] = undefined -- not reached
        takePeriod got_first_o (x:xs) =
            if (fst x == fst o) then
                if (got_first_o) then
                    [x]                     -- Have two O's, stop iterating here
                else
                    x : takePeriod True xs  -- This was first O, mark as such
            else
                x : takePeriod got_first_o xs

-- | Show period string for a given Δ.
show_period_for_bigD :: CLIntP -> String
show_period_for_bigD bigD =
    let (delta, ideals) = period_of_ideals rho_d $ (unit_ideal bigD, 0)
     in "For bigD=" ++ (show bigD) ++ " the period has "
                 ++ (show $ (length ideals) - 1) ++ " ideals and sum delta is "
                 ++ (show delta)

-- | Show the period for the first /n/ valid Δ's.
show_period_for_many_bigDs :: Int -> IO()
show_period_for_many_bigDs n = do
   putStrLn $ unlines $ map (\bigD -> show_period_for_bigD bigD) $ sort $ take n all_bigDs

-- | Show period string and the list of ideals for a given Δ.
show_period_for_some_bigD :: CLIntP -> IO()
show_period_for_some_bigD bigD = do
    putStrLn $ show_period_for_bigD bigD
    putStrLn "Fwd rho_d:"
    putStrLn $ unlines $ map printIdeal ideals
    putStrLn "Inv rho_d:"
    putStrLn $ unlines $ map printIdeal invideals
    where
        (delta,    ideals)    = period_of_ideals rho_d    $ (unit_ideal bigD, 0)
        (invdelta, invideals) = period_of_ideals rho_inv_d $ (unit_ideal bigD, 0)
        printIdeal ideal =
            (show ideal)
            ++ " Reduced: "
            ++ if (is_reduced $ fst ideal) then "true" else "false"

-- | Show a list of valid Δ's.
show_bigDs :: Int -> IO()
show_bigDs n = do
    putStrLn $ show $ take n all_bigDs

-- | Explicitly compute first few ideals for some Δ.
first_few :: IO()
first_few = do
    putStrLn $ "O   :" ++ show j_0
    putStrLn $ "j1/2:" ++ show j_05
    putStrLn $ "j1  :" ++ show j_1
    where
        bigD = 17
        j_0  = (unit_ideal bigD, 0)
        j_05 = rho_d j_0
        j_1  = rho_d j_05

-- | Perform an operation on all ideal pairs that are generated by Δ.
op_all_ideals :: (IdDist -> IdDist -> IdDist) -> String -> CLIntP -> IO()
op_all_ideals op opString bigD = do
    putStrLn $ unlines $ [ doOp i j | i <- ideals, j <- ideals ]
    where
        (delta, ideals_with_o) = period_of_ideals rho_d $ (unit_ideal bigD, 0)
        ideals = init ideals_with_o
        doOp i j = "(" ++ (show i) ++ ")" ++ opString ++ "(" ++ (show j) ++ ") = "
                    ++ (show (i_op_j))
                    ++ " Reduced:"
                    ++ (if (is_reduced $ fst i_op_j) then "true" else "false")
--                    ++ " rho_d of:" ++ (show $ rho_d i_op_j)
                    where i_op_j = i `op` j

-- | The the product of all pairs of ideals for a given Δ.
dot_all_ideals :: CLIntP -> IO()
dot_all_ideals bigD = op_all_ideals dot "." bigD

-- | Take the star product of all pairs of ideals for a given Δ.
star_all_ideals :: CLIntP -> IO()
star_all_ideals bigD = op_all_ideals star "*" bigD

-- | Test the 'bounded_while' functionality.
test_bounded_while :: (Show int, Integral int) => int -> int -> IO()
test_bounded_while bound start = do
    putStrLn $ show $
        bounded_while (\k -> k > 0) bound
            (\k -> k-1) start

-- | Run classical tests for Class Number algorithm.
main :: IO()
main = do
--    test_bounded_while 10 5
--    putStrLn $ "a=23, b=-41, bigD=28, tau =" ++ show (tau (-41) 23 28)
--    putStrLn $ "a=23, b=-41, bigD=28, itau=" ++ show (itau (-41) 23 28)
--    putStrLn $ unlines $ testTauForDelta 28 tau
--    first_few
--    showDs 50
--    showPeriodForManyBigDs 400
--    show_period_for_some_bigD 28
--    dot_all_ideals 28
    star_all_ideals 28
--    putStrLn $ show $ rho_d $ (Ideal 28 1 1 9 8, 0)

--    putStrLn $ show $ take 100 all_small_ds
--    putStrLn $ show $ sort $ take 100 all_bigDs

-- For bigD=2524 the period has 48 ideals and sum delta is 41.3199021281136

--    putStrLn $ show $ continued_list 649 200
--    putStrLn $ show $ convergents $ continued_list 649 200

-- | Test the primes code.
test_primes :: IO ()
test_primes = do
--    putStrLn $ show $ jacobi_symbol 1001 9907
--    putStrLn $ show $ jacobi_symbol 14 7
    putStrLn $ show $ primes_to 8000