{-# LANGUAGE ScopedTypeVariables #-}

-- | This module implements the approximate single-qubit synthesis
-- algorithm of
-- 
-- * N. J. Ross and P. Selinger, \"Optimal ancilla-free Clifford+/T/
-- approximation of /z/-rotations\". <http://arxiv.org/abs/1403.2975>.
-- 
-- The algorithm is near-optimal in the following sense: it produces
-- an operator whose expected /T/-count exceeds the /T/-count of the
-- second-to-optimal solution to the approximate synthesis problem by
-- at most /O/(log(log(1/ε))).

module Quantum.Synthesis.GridSynth where

import Quantum.Synthesis.Ring
import Quantum.Synthesis.Ring.FixedPrec
import Quantum.Synthesis.Matrix
import Quantum.Synthesis.CliffordT
import Quantum.Synthesis.SymReal
import Quantum.Synthesis.GridProblems
import Quantum.Synthesis.Diophantine
import Quantum.Synthesis.StepComp
import Quantum.Synthesis.QuadraticEquation

import System.Random
import Data.Number.FixedPrec

-- ----------------------------------------------------------------------
-- * Approximate synthesis

-- ----------------------------------------------------------------------
-- ** User-friendly functions

-- | Output a unitary operator in the Clifford+/T/ group that
-- approximates /R/[sub /z/](θ) = [exp −/i/θ/Z/\/2] to within ε in the
-- operator norm. This operator can then be converted to a list of
-- gates with 'to_gates'.
-- 
-- The parameters are:
-- 
-- * a source of randomness /g/;
-- 
-- * the angle θ;
--   
-- * the precision /b/ ≥ 0 in bits, such that ε = 2[sup -/b/];
-- 
-- * an integer that determines the amount of \"effort\" to put into
-- factoring. A larger number means more time spent on factoring. 
-- A good default for this is 25.
-- 
-- Note: the argument /theta/ is given as a symbolic real number. It
-- will automatically be expanded to as many digits as are necessary
-- for the internal calculation. In this way, the caller can specify,
-- e.g., an angle of pi\/128 @::@ 'SymReal', without having to worry
-- about how many digits of π to specify.
gridsynth :: (RandomGen g) => g -> Double -> SymReal -> Int -> U2 DOmega
gridsynth g prec theta effort = m where
  (m, _, _) = gridsynth_stats g prec theta effort

-- | A version of 'gridsynth' that returns a list of gates instead of a
-- matrix.
-- 
-- Note: the list of gates will be returned in right-to-left order,
-- i.e., as in the mathematical notation for matrix multiplication.
-- This is the opposite of the quantum circuit notation.
gridsynth_gates :: (RandomGen g) => g -> Double -> SymReal -> Int -> [Gate]
gridsynth_gates g prec theta effort = synthesis_u2 (gridsynth g prec theta effort)
    
-- | A version of 'gridsynth' that also returns some statistics:
-- log[sub 0.5] of the actual approximation error (or 'Nothing' if the
-- error is 0), and a data structure with information on the
-- candidates tried.
gridsynth_stats :: (RandomGen g) => g -> Double -> SymReal -> Int -> (U2 DOmega, Maybe Double, [(DOmega, DStatus)])
gridsynth_stats g prec theta effort = dynamic_fixedprec2 digits f prec theta where
  digits = ceiling (15 + 2 * prec * logBase 10 2)
  f prec theta = gridsynth_internal g prec theta effort
        
-- | Information about the status of an attempt to solve a Diophantine
-- equation. 'Success' means the Diophantine equation was solved;
-- 'Fail' means that it was proved that there was no solution;
-- 'Timeout' means that the question was not decided within the
-- allotted time.
data DStatus = Success | Fail | Timeout
             deriving (Eq, Show)
                                
-- ----------------------------------------------------------------------
-- * Implementation details

-- ----------------------------------------------------------------------
-- ** The ε-region

-- | The ε-/region/ for given ε and θ is a convex subset of the closed
-- unit disk, given by [nobr /u/ ⋅ /z/ ≥ 1 - ε²\/2], where [nobr /z/ =
-- [exp −/i/θ\/2]], and “⋅” denotes the dot product of ℝ² (identified
-- with ℂ).
-- 
-- \[center [image Re.png]]

epsilon_region :: (Floating r, Ord r, RootHalfRing r, Quadratic r) => r -> r -> ConvexSet r
epsilon_region epsilon theta = ConvexSet ell tst int where
  
  -- A bounding ellipse for the ε-region.
  ell = Ellipse mat ctr
  ctr = (d*zx, d*zy)
  mat = bmat * mmat * special_inverse bmat
  mmat = toOperator ((ev1, 0), (0, ev2))
  bmat = toOperator ((zx, -zy), (zy, zx))
  ev1 = 4 * (1 / epsilon)^4
  ev2 = (1 / epsilon)^2
  
  -- A line intersector for the ε-region.
  int p v
    | q == Nothing         = (1, 0)
    | vz == 0 && rhs <= 0  = (t0, t1)
    | vz == 0 && otherwise = (1, 0)
    | vz > 0               = (max t0 t2, t1)
    | otherwise            = (t0, min t1 t2)
    where
      a = iprod v v
      b = 2 * iprod v p
      c = iprod p p - 1
      q = quadratic (fromDRootTwo a) (fromDRootTwo b) (fromDRootTwo c)
      Just (t0, t1) = q
    
      -- solve (p + tv) * z >= d
      -- equivalently, t * vz >= d - pz
      vz = iprod (point_fromDRootTwo v) z
      rhs = d - iprod (point_fromDRootTwo p) z
      t2 = rhs / vz

  -- The characteristic function of the ε-region.
  tst (x,y) = x^2 + y^2 <= 1 && zx * fromDRootTwo x + zy * fromDRootTwo y >= d
  
  zx = cos (-theta/2)
  zy = sin (-theta/2)
  d = 1 - epsilon^2/2
  z = (zx, zy)
  
-- ----------------------------------------------------------------------
-- ** Main algorithm implementation
    
-- | The internal implementation of the ellipse-based approximate
-- synthesis algorithm. The parameters are a source of randomness /g/,
-- the angle θ, the precision /b/ ≥ 0 in bits, and an amount of
-- \"effort\" to put into factoring.
-- 
-- The outputs are a unitary operator in the Clifford+/T/ group that
-- approximates /R/[sub /z/](θ) to within ε in the operator norm;
-- log[sub 0.5] of the actual error, or 'Nothing' if the error is 0;
-- and the number of candidates tried.
-- 
-- Note: the parameter θ must be of a real number type that has enough
-- precision to perform intermediate calculations; this typically
-- requires precision O(ε[sup 2]).  A more user-friendly function that
-- selects the required precision automatically is 'gridsynth'.
gridsynth_internal :: forall r g.(RootHalfRing r, Ord r, Floating r, Adjoint r, Floor r, RealFrac r, Quadratic r, RandomGen g) => g -> r -> r -> Int -> (U2 DOmega, Maybe Double, [(DOmega, DStatus)])
gridsynth_internal g prec theta effort = (uU, log_err, candidate_info) where
  epsilon = 2 ** (-prec)
  region = epsilon_region epsilon theta
  candidates = gridpoints2_increasing region unitdisk
  (uU, log_err, candidate_info) = first_solvable [] g candidates
  
  first_solvable candidate_info g [] = error "gridsynth: internal error: finite list of candidates?"
  first_solvable candidate_info g (u : us) = case answer_t of
    Just (Just t) -> let (uU, log_err) = with_successful_candidate u t in (uU, log_err, ((u, Success) : candidate_info))
    Just Nothing -> first_solvable ((u, Fail) : candidate_info) g2 us
    Nothing -> first_solvable ((u, Timeout) : candidate_info) g2 us
    where
      (g1, g2) = split g
      xi = real (1 - adj u * u)
      answer_t = run_bounded effort $ diophantine_dyadic g1 xi
  
  with_successful_candidate u t = (uU, log_err) where
    uU | denomexp (u + t) < denomexp (u + omega * t)
               = matrix2x2 (u, -(adj t)) (t, adj u)
       | otherwise
               = matrix2x2 (u, -(adj (omega*t))) (omega*t, adj u)
    log_err 
      | err <= 0  = Nothing
      | otherwise = Just (logBase_double 0.5 err)
    err = sqrt (real (hs_sqnorm (uU_fixed - zrot_fixed)) / 2)
    uU_fixed = matrix_map fromDOmega uU
    zrot_fixed = zrot (theta :: r)