{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LINE 1 "Quipper/Libraries/Arith.hs" #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeSynonymInstances #-}

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

-- | This library provides a type of quantum integers, as well as
-- basic arithmetic functions on them.
-- 
-- The type 'QDInt' holds a fixed-size, ℓ-qubit quantum integer,
-- considered modulo 2[sup ℓ]. The integers may be regarded as
-- signed or unsigned, depending on the operation. If signs are used,
-- they are assumed to be in two's complement.
-- 
-- Some of the arithmetic operations are adapted from the GFI for the
-- Triangle Finding algorithm. Most algorithms used are, for now, very
-- naïve (ripple adders, etc).  Gate count estimates are given in the
-- Toffoli gatebase.

-- Documentation note: the table of contents is organized slightly
-- differently than the source code. The guiding idea is that the
-- documentation should expose the interface that is accessible to
-- user-level code, and changeable implementation details are hidden.

module Quipper.Libraries.Arith (
  -- * Quantum integers
  -- ** Data type definitions
  -- $DATATYPES
  XInt,  -- constructors not exported
  QDInt,
  CInt,
  IntM,
  -- ** Operations on QDInt
  qulist_of_qdint_bh,
  qdint_of_qulist_bh,
  qulist_of_qdint_lh,
  qdint_of_qulist_lh,
  qdint_length,
  qdint_extend_unsigned,
  qdint_extend_signed,
  -- ** Operations on CInt
  bitlist_of_cint_bh,
  cint_of_bitlist_bh,
  bitlist_of_cint_lh,
  cint_of_bitlist_lh,
  cint_length,
  cint_extend_unsigned,
  cint_extend_signed,
  -- ** Operations on IntM
  -- $INTMCLASSES
  boollist_of_intm_bh,
  intm_of_boollist_bh,
  intm_length,
  integer_of_intm_unsigned,
  integer_of_intm_signed,
  intm_with_length,
  intm_of_integer,
  intm,
  intm_promote,
  intm_interval_signed,
  intm_interval_unsigned,
  intm_extend_unsigned,
  intm_extend_signed,
  -- ** Shape parameters
  qdint_shape,
  cint_shape,
  -- ** Operations on XInt
  xint_maybe_length,
  list_of_xint_bh,
  xint_of_list_bh,
  list_of_xint_lh,
  xint_of_list_lh,
  -- * Quantum arithmetic operations
  -- ** The QNum type class
  QNum(..),
  -- ** In-place increment and decrement
  q_increment,
  q_decrement,
  -- ** In-place addition and subtraction
  q_add_in_place,
  q_sub_in_place,
  q_negate_in_place,
  -- ** Arithmetic with classical parameter
  q_add_param,
  q_sub_param,
  q_add_param_in_place,
  q_sub_param_in_place,
  q_mult_param,
  -- ** Comparison
  q_le_unsigned,
  q_le_signed,
  q_lt_signed,
  q_negative,
  -- ** Division and remainder
  q_moddiv_unsigned_in_place,
  q_mod_unsigned,
  q_divrem_unsigned,
  q_div_unsigned,
  q_div,
  q_quot,
  q_div_exact_unsigned,
  q_div_exact,
  -- ** Specialized functions
  q_ext_euclid,
  -- * Lifting of arithmetic functions
  -- $LIFTING
  template_symb_plus_,
  ) where

import Quipper
import Quipper.Internal

import Quipper.Utils.Sampling
import Quipper.Utils.Auxiliary

import Control.Monad
import Data.Typeable

-- ======================================================================
-- * Quantum integers

-- ** Data type definitions

-- $DATATYPES 
-- We define three versions of the fixed-length integer type: quantum,
-- classical input, and classical parameter. The triple ('IntM',
-- 'QDInt', 'CInt') forms an instance of the 'QShape' class.  All three
-- types are special cases of the type 'XInt' /x/.

-- | 'XInt' /x/ is the type of fixed-length integers, but using
-- elements of type /x/ instead of bits. It is an abstract type, and
-- details of its implementation is not exposed to user-level code.

-- Implemenation notes: The integer types are currently implemented as
-- big-headian bit lists. However, the details of the implementation
-- are not exposed to user-level code, and are subject to change.
-- 
-- The form ('XInt_indet' /n/ /id/) is permitted as a special case, to
-- represent an integer of indeterminate length. Such a value can only
-- be used when the length is deducible from the context. This form is
-- only allowed in the special case /x/ = 'Bool', and we use an
-- identity type to enforce this.
-- 
-- Integers of indeterminate length may only be used in certain
-- operations where the shape information is available from other
-- data. It can be used, for example, for terminating or controlling a
-- 'QDInt', but not for initializing a 'QDInt'.
data XInt x = XInt [x] | XInt_indet Integer (Identity Bool x)
  deriving (Show, Typeable)

-- | The type of fixed-length /m/-qubit quantum integers. This is a
-- circuit execution time value.
type QDInt = XInt Qubit

instance Show QDInt where
  show (XInt l) = "#" ++ show l
  show (XInt_indet n id) = error "IntM: internal error"

-- A better implementation might be something like:
--   show (XInt l) = "qdint_of_qulist_bh " ++ show l
-- However, as 'CInt' and 'QDInt' are currently synonyms, it seems best
-- to keep the instance agnostic between them, while still distinguishing
-- them somehow from naked lists.

-- | The type of fixed-length /m/-bit classical integer inputs. This
-- is a circuit execution time value.
type CInt = XInt Bit

-- Currently, 'CInt' is literally a synonym of 'QDInt', since 'Qubit' is
-- a synonym of 'Bit'.  If this changes, the following instance will be
-- required:
--
-- instance Show CInt where
--   show (XInt l) = show l
--   show (XInt_indet n id) = error "IntM: internal error"

-- | The type of fixed-length /m/-bit integer parameters.  Values of
-- this type are parameters, i.e., they are classical and known at
-- circuit generation time.
-- 
-- Unlike values of type 'QDInt' and 'CInt', a value of type 'IntM'
-- may have an indeterminate length. This happens, for example, if the
-- value is specified by means of an integer literal (e.g., 17), which
-- does not carry length information. In such cases, the value can
-- only be used when it can be deduced from the context. For example,
-- such values may be used for terminating or controlling a 'QDInt',
-- but not for initializing a 'QDInt'.
type IntM = XInt Bool

instance Show IntM where
  show (XInt l) = "intm " ++ show (length l) ++ " " ++ show (int_of_boollist_signed_bh l)
  show (XInt_indet n id) = show n

-- Note: for the purpose of forming intervals, we regard 'IntM' as an
-- /unsigned/ type.  This disagrees with the 'Enum' instance, where
-- 'IntM' is regarded as a /signed/ type. So for example, when /x/ =
-- 'intm' 4 (-2) = 'intm' 4 14 and /y/ = 'intm' 4 10, then [/x/../y/]
-- and 'interval' /y/ /x/ are non-empty, but [/y/../x/] and 'interval'
-- /x/ /y/ are empty. 
-- 
-- This is confusing but for the time being there is no better
-- alternative, unless we create distinct signed and unsigned types.

instance Interval IntM where
  interval x y = intm_interval_unsigned x y

instance Zero IntM where
  zero x = intm_with_length (intm_length x) 0
  
-- ----------------------------------------------------------------------
-- ** Primitive combinators on XInt

-- $ 'XInt' is intended to be an abstract data type, and all access to
-- it should pass through the three access functions of this
-- section.

-- *** Constructors

-- | Create a 'XInt' /x/ from a list of /x/s. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
xint_of_list_bh :: [x] -> XInt x
xint_of_list_bh xs = XInt xs

-- | Create an 'IntM' of indeterminate length from an 'Integer'. 
intm_of_integer :: Integer -> IntM
intm_of_integer n = XInt_indet n reflexivity

-- *** Destructor

-- | If the 'XInt' is of determinate length, return its list of digits
-- as a big-headian list, i.e., the head of the list holds the most
-- significant digit. If the 'XInt' is of indeterminate length, return
-- (/n/, /id/), where /n/ is the underlying 'Integer' and /id/ is the
-- witness proving that /x/ = 'Bool'.
-- 
-- This is a lowest-level access function not intended to be used by
-- user-level code.
xint_case :: XInt x -> Either [x] (Integer, Identity Bool x)
xint_case (XInt xs) = Left xs
xint_case (XInt_indet n id) = Right (n, id)

-- ----------------------------------------------------------------------
-- ** Other low-level operations

-- $ The operations in this section are the only ones intended to use
-- 'xint_case' directly.

-- | Set the length of an 'XInt' to /m/ ≥ 0. This operation is only
-- legal if the input (a) has indeterminate length or (b) has
-- determinate length already equal to /m/. In particular, it cannot
-- be used to change the length from anything other than from
-- indeterminate to determinate.
-- 
-- If both arguments already have determinate lengths, and they do not
-- coincide, throw an error. The 'String' argument is used as an error
-- message in that case.
-- 
-- Note that if /x/ ≠ 'Bool', the input is guaranteed to have
-- determinate length. However, we cannot test for equality of types
-- in a polymorphic function. This is where the 'id' argument to
-- 'XInt_indet' is used.
xint_set_length :: Int -> XInt x -> String -> XInt x
xint_set_length m x errmsg | m < 0 = 
  error "xint_set_length: negative length not permitted"  
xint_set_length m x errmsg =
  case xint_case x of
    Left xs | m == length xs -> x
            | otherwise -> error errmsg
    Right (n, id) -> XInt xs where
      xs = [ identity id b | b <- boollist_of_int_bh m n ]

-- | Return 'True' if the 'XInt' is of determinate length, and 'False'
-- if it is of indeterminate length. 
xint_is_determinate :: XInt x -> Bool
xint_is_determinate x = 
  case xint_case x of
    Left _ -> True
    Right _ -> False

-- | From a 'XInt', which must be of determinate length, extract a
-- list of /x/s. The conversion is big-headian, i.e., the head of the
-- list holds the most significant digit. It is an error to call this
-- function with an 'XInt' of indeterminate length.
list_of_xint_bh :: XInt x -> [x]
list_of_xint_bh x = 
  case xint_case x of
    Left xs -> xs
    Right _ -> error "list_of_xint_bh: integer has indeterminate length"

-- | Return the size of a 'XInt', or 'Nothing' if indeterminate.
xint_maybe_length :: XInt x -> Maybe Int
xint_maybe_length x =
  case xint_case x of
    Left xs -> Just (length xs)
    Right _ -> Nothing

-- | Convert an 'IntM' of length /m/ to an 'Integer' in the range {0,
-- …, 2[sup /m/]-1}. If the 'IntM' has indeterminate length, return the
-- original 'Integer'.
integer_of_intm_unsigned :: IntM -> Integer
integer_of_intm_unsigned x = 
  case xint_case x of
    Left xs -> int_of_boollist_unsigned_bh xs
    Right (n, id) -> n

-- | Convert an 'IntM' of length /m/ to an 'Integer' in the range
-- {-2[sup /m/-1], …, 2[sup /m/-1]-1}. If the 'IntM' has indeterminate
-- length, return the original 'Integer'.
integer_of_intm_signed :: IntM -> Integer
integer_of_intm_signed x =
  case xint_case x of
    Left xs -> int_of_boollist_signed_bh xs
    Right (n, id) -> n

-- | Equality test. If at least one argument has determinate length,
-- test equality modulo 2[sup /m/]. If both have indeterminate length,
-- check equality of the underlying integers.
xint_equals :: (Eq x) => XInt x -> XInt x -> Bool    
xint_equals x y =
  case (xint_case x, xint_case y) of
    (Left xs, Left ys) 
      | length xs == length ys -> xs == ys
      | otherwise -> error "Equality test on XInt: operands must be of equal length"
    (_, Left ys) -> xint_equals (xint_set_length m x "xint_equals") y
        where m = length ys
    (Left xs, _) -> xint_equals x (xint_set_length m y "xint_equals")
        where m = length xs
    (Right (n, _), Right (n', _)) -> n == n'
    
-- ----------------------------------------------------------------------
-- ** Derived operations on XInt

-- | Get the nominal length of an integer (in bits). It is an error to
-- apply this function to an integer of indeterminate length.
xint_length :: XInt x -> Int
xint_length x = 
  case xint_maybe_length x of
    Just m -> m
    Nothing -> error "xint_length: integer has indeterminate length"

-- | Convert an integer to a bit list. The conversion is
-- little-headian, i.e., the head of the list holds the least
-- significant digit.
list_of_xint_lh :: XInt x -> [x]
list_of_xint_lh = reverse . list_of_xint_bh

-- | Convert a bit list to an integer. The conversion is
-- little-headian, i.e., the head of the list holds the least
-- significant digit.
xint_of_list_lh :: [x] -> XInt x
xint_of_list_lh = xint_of_list_bh . reverse

-- | Extend a 'XInt' to the given length without changing its
-- (unsigned) value. This is done by adding the required number of
-- high bits initialized to /zero/. It is an error to call this
-- function when the new length is shorter than the old one.
xint_extend_unsigned :: (Monad m) => Int -> m x -> XInt x -> m (XInt x)
xint_extend_unsigned len zero x 
  | len < m =
    error "pad_xint: requested length is shorter than current length"
  | otherwise = do
    pad <- sequence (replicate extra zero)
    return $ xint_of_list_bh (pad ++ digits)
  where
    digits = list_of_xint_bh x
    m = length digits
    extra = len - m

-- | Extend a 'XInt' to the given length without changing its (signed)
-- value. This is done by adding the required number of high bits
-- initialized to copies of the sign bit (or to /zero/ if the original
-- integer was of length 0). It is an error to call this function when
-- the new length is shorter than the old one.
xint_extend_signed :: (Monad m) => Int -> m x -> (x -> m x) -> XInt x -> m (XInt x)
xint_extend_signed len zero copy x
  | len < m 
    = error "pad_xint: requested length is shorter than current length"
  | m == 0
    = xint_extend_unsigned len zero x
  | otherwise = do
    pad <- sequence (replicate extra (copy sign))
    return $ xint_of_list_bh (pad ++ digits)
  where
    digits = list_of_xint_bh x
    m = length digits
    extra = len - m
    sign = head digits

-- ----------------------------------------------------------------------
-- ** Operations on IntM

-- | Return the size of an 'IntM', or 'Nothing' if indeterminate.
intm_length :: IntM -> Maybe Int
intm_length = xint_maybe_length
  
-- | Convert an 'IntM' to a list of booleans. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit. As usual, 'False' is 0 and 'True' is 1. It is an error to
-- apply this operation to an 'IntM' whose length is indeterminate.
boollist_of_intm_bh :: IntM -> [Bool]
boollist_of_intm_bh = list_of_xint_bh

-- | Convert a boolean list to an 'IntM'. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
intm_of_boollist_bh :: [Bool] -> IntM
intm_of_boollist_bh = xint_of_list_bh

-- | Create an 'IntM' of the specified length (first argument) and
-- value (second argument).
intm :: Int -> Integer -> IntM
intm m n = intm_set_length m (intm_of_integer n) "intm: internal error"

-- | Set the length of an 'IntM' to /m/ ≥ 0. This operation is only
-- legal if the input (a) has indeterminate length or (b) has
-- determinate length already equal to /m/. In particular, it cannot
-- be used to change the length from anything other than from
-- indeterminate to determinate. (Use 'intm_extend_unsigned' or
-- 'intm_extend_signed' to increase a determinate length).
-- 
-- If both arguments already have determinate lengths, and they do not
-- coincide, throw an error. The 'String' argument is used as an error
-- message in that case.
intm_set_length :: Int -> IntM -> String -> IntM
intm_set_length = xint_set_length

-- | Extend an 'IntM' to the given length without changing its
-- (unsigned) value. This is done by adding the required number of
-- high bits initialized to 0. It is an error to call this function
-- when the new length is shorter than the old one.
intm_extend_unsigned :: Int -> IntM -> IntM
intm_extend_unsigned len x =
  getId $ xint_extend_unsigned len (return False) x

-- | Extend an 'IntM' to the given length without changing its
-- (signed) value. This is done by adding the required number of
-- high bits initialized to copies of the sign bit. It is an error to
-- call this function when the new length is shorter than the old one.
intm_extend_signed :: Int -> IntM -> IntM
intm_extend_signed len x =
  getId $ xint_extend_signed len (return False) (return) x

-- ----------------------------------------------------------------------
-- ** Operations on CInt

-- | Convert a 'CInt' to a list of qubits. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
bitlist_of_cint_bh :: CInt -> [Bit]
bitlist_of_cint_bh = list_of_xint_bh

-- | Convert a list of qubits to a 'CInt'. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
cint_of_bitlist_bh :: [Bit] -> CInt
cint_of_bitlist_bh = xint_of_list_bh

-- | Convert a 'CInt' to a list of bits. The conversion is
-- little-headian, i.e., the head of the list holds the least
-- significant digit.
bitlist_of_cint_lh :: CInt -> [Bit]
bitlist_of_cint_lh = list_of_xint_lh

-- | Convert a list of bits to a 'CInt'. The conversion is
-- little-headian, i.e., the head of the list holds the least significant
-- digit.
cint_of_bitlist_lh :: [Bit] -> CInt
cint_of_bitlist_lh = xint_of_list_lh

-- | Return the length of a 'CInt', in bits.
cint_length :: CInt -> Int
cint_length = xint_length

-- | Extend a 'CInt' to the given length without changing its
-- (unsigned) value. This is done by adding the required number of
-- high bits initialized to 0. It is an error to call this function
-- when the new length is shorter than the old one.
cint_extend_unsigned :: Int -> CInt -> Circ CInt
cint_extend_unsigned len x =
  xint_extend_unsigned len (cinit False) x

-- | Extend a 'CInt' to the given length without changing its
-- (signed) value. This is done by adding the required number of
-- high bits initialized to copies of the sign bit. It is an error to
-- call this function when the new length is shorter than the old one.
cint_extend_signed :: Int -> CInt -> Circ CInt
cint_extend_signed len x =
  xint_extend_signed len (cinit False) (qc_copy) x

-- ----------------------------------------------------------------------
-- ** Operations on QDInt

-- | Convert a 'QDInt' to a list of qubits. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
qulist_of_qdint_bh :: QDInt -> [Qubit]
qulist_of_qdint_bh = list_of_xint_bh

-- | Convert a list of qubits to a 'QDInt'. The conversion is
-- big-headian, i.e., the head of the list holds the most significant
-- digit.
qdint_of_qulist_bh :: [Qubit] -> QDInt
qdint_of_qulist_bh = xint_of_list_bh

-- | Convert a 'QDInt' to a list of qubits. The conversion is
-- little-headian, i.e., the head of the list holds the least significant
-- digit.
qulist_of_qdint_lh :: QDInt -> [Qubit]
qulist_of_qdint_lh = list_of_xint_lh

-- | Convert a list of qubits to a 'QDInt'. The conversion is
-- little-headian, i.e., the head of the list holds the least significant
-- digit.
qdint_of_qulist_lh :: [Qubit] -> QDInt
qdint_of_qulist_lh = xint_of_list_lh

-- | Return the length of a 'QDInt', in bits.
qdint_length :: QDInt -> Int
qdint_length = xint_length

-- | Extend a 'QDInt' to the given length without changing its
-- (unsigned) value. This is done by adding the required number of
-- high bits initialized to 0. It is an error to call this function
-- when the new length is shorter than the old one.
qdint_extend_unsigned :: Int -> QDInt -> Circ QDInt
qdint_extend_unsigned len x =
  xint_extend_unsigned len (qinit False) x

-- | Extend a 'QDInt' to the given length without changing its
-- (signed) value. This is done by adding the required number of
-- high bits initialized to copies of the sign bit. It is an error to
-- call this function when the new length is shorter than the old one.
qdint_extend_signed :: Int -> QDInt -> Circ QDInt
qdint_extend_signed len x =
  xint_extend_signed len (qinit False) (qc_copy) x

-- ----------------------------------------------------------------------
-- ** Shape parameters

-- | Return a piece of shape data to represent an /m/-qubit quantum
-- integer. Please note that the data can only be used as shape; it
-- will be undefined at the leaves.
qdint_shape :: Int -> QDInt
qdint_shape m = xint_of_list_bh (replicate m qubit)

-- | Return a piece of shape data to represent an /m/-bit 'CInt'.
-- Please note that the data can only be used as shape; it will be
-- undefined at the leaves.
cint_shape :: Int -> CInt
cint_shape m = xint_of_list_bh (replicate m bit)

-- ======================================================================
-- ** Type class instances

-- Note: instance declarations do not show up in the documentation

-- | Input the lengths of two inputs to a binary operation, and return
-- the common length. If one length is indeterminate, return the other;
-- if both are indeterminate, return 'Nothing'; if both are determinate but
-- not equal, throw an error with the given error message.
combine_length :: String -> Maybe Int -> Maybe Int -> Maybe Int
combine_length s Nothing m = m
combine_length s m Nothing = m
combine_length s (Just m) (Just m') | m == m' = Just m
                                    | otherwise = error s

-- | Try to set the length of an 'IntM' to that of another 'XInt'
-- value (which could be a 'QDInt', a 'CInt', or another 'IntM'). This
-- will fail with an error if both numbers already have determinate
-- lengths that don't coincide. In this case, the string argument is
-- used as an error message.
intm_promote :: IntM -> XInt x -> String -> IntM
intm_promote bi xi errmsg = 
  case xint_maybe_length xi of
    Nothing -> bi
    Just m -> intm_set_length m bi errmsg
    
type instance QCType x y (XInt z) = XInt (QCType x y z)
type instance QTypeB IntM = QDInt

instance QCLeaf x => QCData (XInt x) where
  qcdata_mapM shape f g xs = 
    mmap xint_of_list_bh $ qcdata_mapM (list_of_xint_bh shape) f g (list_of_xint_bh xs)
  qcdata_zip shape q c q' c' xs ys e = 
    xint_of_list_bh $ qcdata_zip (list_of_xint_bh shape) q c q' c' (list_of_xint_bh xs) (list_of_xint_bh ys) (const $ e "XInt length mismatch")
  qcdata_promote b q e = intm_promote b q (e "IntM length mismatch")

-- Labeling of QDInt is s[m-1], ..., s[0], with the least significant
-- bit at index 0.
instance QCLeaf x => Labelable (XInt x) String where
  label_rec qa = label_rec (list_of_xint_lh qa)

instance CircLiftingUnpack (Circ QDInt) (Circ QDInt) where
  pack x = x
  unpack x = x

-- ======================================================================
-- * Classical arithmetic on 'IntM'

-- ----------------------------------------------------------------------
-- ** Auxiliary functions
  
-- | A useful auxiliary function that converts a list of 'IntM's to a
-- list of 'Integer's, while also returning their common length. All
-- arguments whose length is not indeterminate must have the same common
-- length. If the lengths don't match, throw an error with the error
-- message specified by the 'String' argument.
integers_of_intms_signed :: [IntM] -> String -> (Maybe Int, [Integer])
integers_of_intms_signed xs s = (m, is) where
  m = foldl (combine_length s) Nothing [ intm_length x | x <- xs ]
  is = [ integer_of_intm_signed x | x <- xs ]  

-- | Like 'integers_of_intms_signed', but regards the 'IntM's as
-- unsigned integers.
integers_of_intms_unsigned :: [IntM] -> String -> (Maybe Int, [Integer])
integers_of_intms_unsigned xs s = (m, is) where
  m = foldl (combine_length s) Nothing [ intm_length x | x <- xs ]
  is = [ integer_of_intm_unsigned x | x <- xs ]  

-- | Create an 'IntM' of the given length and value. Leave the length
-- indeterminate if it is given as 'Nothing'.
intm_with_length :: Maybe Int -> Integer -> IntM
intm_with_length (Just m) n = intm m n
intm_with_length Nothing n = intm_of_integer n

-- | Auxiliary function for lifting a binary operator from 'Integer'
-- to 'IntM'. The string argument is the name of the operator, for
-- error messages.
intm_binop :: (Integer -> Integer -> Integer) -> String -> IntM -> IntM -> IntM 
intm_binop op opname x y = intm_with_length m (op x' y') where
  (m, [x',y']) = integers_of_intms_signed [x, y] ("Binary operation " ++ opname ++ " on IntM: operands must be of equal length")

-- | Auxiliary function for lifting a unary operator from 'Integer' to
-- 'IntM'.
intm_unop :: (Integer -> Integer) -> IntM -> IntM
intm_unop op x = intm_with_length xm (op x') where
  xm = intm_length x
  x' = integer_of_intm_signed x
  
-- ----------------------------------------------------------------------
-- ** Type class instances

-- $INTMCLASSES 
-- 
-- 'IntM' is an instance of Haskell's 'Eq', 'Num', 'Ord', 'Real',
-- 'Enum', and 'Integral' type classes. This means that integer
-- literals (e.g., 17), and the usual arithmetic functions, such as
-- '+', '-', '*', 'abs', 'succ', 'pred', 'mod', 'div', and others, can
-- be used for values of type 'IntM'. In general, we treat 'IntM' as a
-- signed integer type. Use 'fromIntegral' to convert an integer to an
-- 'IntM' of indeterminate length.
-- 
-- The general convention for binary operations (such as
-- multiplication) is: both operands must have the same length,
-- except: if one operand has indeterminate length, it takes on the
-- length of the other; if both operands have indeterminate length, the
-- result will have indeterminate length.

instance Eq x => Eq (XInt x) where
  x == y = xint_equals x y
    
instance Num IntM where
  (+) = intm_binop (+) "+"
  (*) = intm_binop (*) "*"
  (-) = intm_binop (-) "-"
  abs = intm_unop abs
  signum = intm_unop signum
  fromInteger = intm_of_integer

instance Ord IntM where
  compare x y = compare (toInteger x) (toInteger y)

instance Real IntM where
  toRational = toRational . integer_of_intm_signed

instance Enum IntM where
  succ = intm_unop succ
  pred = intm_unop pred
  toEnum = intm_of_integer . fromIntegral
  fromEnum = fromIntegral . integer_of_intm_signed
  enumFrom x = map (intm_with_length m) [x'..] where
    (m, [x']) = integers_of_intms_signed [x] "enumeration: IntM"
  enumFromThen x y = map (intm_with_length m) [x',y'..] where
    (m, [x',y']) = integers_of_intms_signed [x,y] "enumeration: IntM operands must be of equal length"
  enumFromTo x y = map (intm_with_length m) [x'..y'] where
    (m, [x',y']) = integers_of_intms_signed [x,y] "enumeration: IntM operands must be of equal length"
  enumFromThenTo x y z = map (intm_with_length m) [x',y'..z'] where
    (m, [x',y',z']) = integers_of_intms_signed [x,y,z] "enumeration: IntM operands must be of equal length"

-- | Return the interval @[/x/../y/]@, with /x/ and /y/ regarded as
-- signed values of type 'IntM'.
intm_interval_signed :: IntM -> IntM -> [IntM]
intm_interval_signed x y = [x..y]  -- this comes from the 'Enum' instance defined above

-- | Return the interval @[/x/../y/]@, with /x/ and /y/ regarded as
-- unsigned values of type 'IntM'.
intm_interval_unsigned :: IntM -> IntM -> [IntM]
intm_interval_unsigned x y = map (intm_with_length m) [x'..y'] where
  (m, [x',y']) = integers_of_intms_unsigned [x,y] "intm_interval: operands must be of equal length"

instance Integral IntM where
  toInteger = integer_of_intm_signed
  quotRem x y = (intm_with_length m q', intm_with_length m r') where
    (m, [x',y']) = integers_of_intms_signed [x, y] "Division on IntM: operands must be of equal length"
    (q',r') = quotRem x' y'

-- ======================================================================
-- * Quantum arithmetic on 'QDInt'

-- Developer note: each quantum arithmetic function is implemented by
-- an underlying low-level function operating on qubit lists. However,
-- these lower-level functions are *not* exported. If any module needs
-- them, it should define them explicitly using 'qulist_of_qdint_bh'
-- and 'qdint_of_qulist_bh'. In particular, the low-level functions
-- may not perform error checking if it is already performed by their
-- high-level counterparts.

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

-- | @'common_value' error_str [n1,n2,n3]@: if /n1/, /n2/, /n3/ all equal, return this value; else throw an error.  Useful for checking that sizes of 'QDInt's match.
common_value :: (Eq a) => String -> [a] -> a
common_value _ [] = error "common_value: no inputs given"
common_value _ [n] = n
common_value error_str (n:ns) = if common_value error_str ns == n then n else error error_str

-- | @'common_length' error_str [x1,x2,x3]@: if /x1/, /x2/, /x3/ all have the same length, return this value; else throw an error. 
common_length :: String -> [QDInt] -> Int
common_length error_str xs = common_value error_str $ map qdint_length xs

-- ----------------------------------------------------------------------
-- ** Incrementing and decrementing

-- | Increment a 'QDInt' in place.  /O/(ℓ) gates. 
-- 
-- Implementation note: currently tries to minimize gate count, at the cost of a rather long Quipper description.  Can the latter be reduced without increasing the former?  
q_increment :: QDInt -> Circ QDInt
q_increment = mmap xint_of_list_bh . q_increment_qulist . list_of_xint_bh

-- | Low-level implementation of 'q_increment': represents integers as
-- big-headian qubit lists.
q_increment_qulist :: [Qubit] -> Circ [Qubit]
q_increment_qulist [] = return []
q_increment_qulist [x_0] = do x_0 <- qnot x_0; return [x_0]
q_increment_qulist [x_1,x_0] = do x_0 <- qnot x_0; x_1 <- qnot x_1 `controlled` (x_0 .==. 0); return [x_1,x_0]
q_increment_qulist [x_2,x_1,x_0] = do
  x_0 <- qnot x_0
  x_1 <- qnot x_1 `controlled` (x_0 .==. 0)
  x_2 <- qnot x_2 `controlled` (x_0 .==. 0) .&&. (x_1 .==. 0)
  return [x_2,x_1,x_0]
q_increment_qulist x_bits = do
  let x_0 = last x_bits
      x_1 = last $ init x_bits
      x_higher = init $ init $ x_bits
  x_0 <- qnot x_0
  x_1 <- qnot x_1 `controlled` (x_0 .==. 0)
  (x_higher,x_1,x_0) <- with_ancilla (\c -> do
    c <- qnot c `controlled` (x_0 .==. 0) .&&. (x_1 .==. 0)
    (c,rev_x_higher) <- q_increment_qulist_aux c (reverse x_higher)
    c <- qnot c `controlled` (x_0 .==. 0) .&&. (x_1 .==. 0)
    return (reverse rev_x_higher, x_1, x_0))
  return (x_higher ++ [x_1,x_0])
  where
  -- increment a LITTLE-endian bit-string, controlled by a single qubit.
    q_increment_qulist_aux :: Qubit -> [Qubit] -> Circ (Qubit,[Qubit])
    q_increment_qulist_aux b [] = return (b,[])
    q_increment_qulist_aux b [x_0] =
      do x_0 <- qnot x_0 `controlled` b; return (b,[x_0])
    q_increment_qulist_aux b [x_0,x_1] = do
      x_0 <- qnot x_0 `controlled` b
      x_1 <- qnot x_1 `controlled` b .&&. (x_0 .==. 0)
      return (b,[x_0,x_1])
    q_increment_qulist_aux b (x_0:x_higher) = do
      x_0 <- qnot x_0 `controlled` b 
      (b,x_0,x_higher) <- with_ancilla (\c -> do
        c <- qnot c `controlled` b .&&. (x_0 .==. 0)
        (c,x_higher) <- q_increment_qulist_aux c x_higher
        c <- qnot c `controlled` b .&&. (x_0 .==. 0)
        return (b,x_0,x_higher))
      return (b, (x_0:x_higher))
    
-- | Decrement a 'QDInt' in place.  Inverse of 'q_increment'. /O/(ℓ).
q_decrement :: QDInt -> Circ QDInt
q_decrement = reverse_generic_endo q_increment

-- ----------------------------------------------------------------------
-- ** Addition and subtraction

-- | Add two 'QDInt's into a fresh one.  Arguments are assumed to be of equal size.  /O/(ℓ) gates, both before and after transformation to Toffoli.
q_add_qdint :: QDInt -> QDInt -> Circ (QDInt, QDInt, QDInt)
q_add_qdint x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y', z') <- q_add_qulist x' y'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  let z = xint_of_list_bh z'
  return (x, y, z)

-- | Low-level implementation of 'q_add_qdint': represents integers as
-- big-headian qubit lists.
q_add_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit], [Qubit])
q_add_qulist x y = do
  let l = length x
  when (l /= length y) $ do
    error "q_add_qdint: cannot add QDInts of different lengths"

  ((x,y),s_out) <- with_computed_fun (x,y)
    (\(x,y) -> do
      s <- qinit (replicate l False) -- holds the eventual sum
      c <- qinit (replicate l False) -- holds the carries

      -- bring the lists’ indexing in line with usual numeric indexing of binary expansons
      x <- return $ reverse x
      y <- return $ reverse y
      s <- return $ reverse s
      c <- return $ reverse c

      (x,y,s,c) <- loop_with_indexM (l-1) (x,y,s,c) (\j (x,y,s,c) -> do
        let c_j1 = c !! (j+1)
        let s_j = s !! j
        -- If any two of the current bits x_j, y_j, and the carry c_j are true, then set the next carry:
        -- (Note: we use that (a & b) xor (a & c) xor (b & c) == (a & b) or (a & c) or (b & c).)
        c_j1 <- qnot c_j1 `controlled` (x!!j .&&. y!!j)
        c_j1 <- qnot c_j1 `controlled` (x!!j .&&. c!!j)
        c_j1 <- qnot c_j1 `controlled` (y!!j .&&. c!!j)
        -- Now the current sum digit is computed mod 2:
        s_j <- qnot s_j `controlled` (x!!j)
        s_j <- qnot s_j `controlled` (y!!j)
        s_j <- qnot s_j `controlled` (c!!j) 
        c <- return $ overwriteAt (j+1) c_j1 c
        s <- return $ overwriteAt j s_j s

        return (x,y,s,c))

      -- Final sum digit; no carry required:
      let s_l1 = s !! (l-1)
      s_l1 <- qnot s_l1 `controlled` (x!!(l-1))
      s_l1 <- qnot s_l1 `controlled` (y!!(l-1))
      s_l1 <- qnot s_l1 `controlled` (c!!(l-1)) 
      s <- return $ overwriteAt (l-1) s_l1 s
     
      x <- return $ reverse x
      y <- return $ reverse y
      s <- return $ reverse s
      c <- return $ reverse c

      return (x,y,s,c))
    
    (\(x,y,s,c) -> do
      (s,s_out) <- qc_copy_fun s
      return ((x,y,s,c), s_out))
  return (x, y, s_out)

-- | Subtract two 'QDInt's, into a fresh one.  Arguments are assumed to be of equal size.  /O/(ℓ).
q_sub_qdint :: QDInt -> QDInt -> Circ (QDInt,QDInt,QDInt)
q_sub_qdint x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y', z') <- q_sub_qulist x' y'
  let x = xint_of_list_bh x'  
  let y = xint_of_list_bh y'
  let z = xint_of_list_bh z'
  return (x, y, z)

-- | Low-level implementation of 'q_sub_qdint': represents integers as
-- big-headian qubit lists.
q_sub_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit], [Qubit])
q_sub_qulist x y = do
  let l = length x
  when (l /= length y) $ do
    error "q_sub_qdint: cannot subtract QDInts of different lengths"
  
  ((x,y),d_out) <- with_computed_fun (x,y)
    (\(x,y) -> do
      d <- qinit (replicate l False) -- holds the eventual difference
      b <- qinit (replicate l False) -- holds the borrows

      -- bring the lists’ indexing in line with usual numeric indexing of binary expansons
      x <- return $ reverse x
      y <- return $ reverse y
      d <- return $ reverse d
      b <- return $ reverse b

      (x,y,d,b) <- loop_with_indexM (l-1) (x,y,d,b) (\j (x,y,d,b) -> do
        let b_j1 = b !! (j+1)
        let d_j = d !! j
        -- If any two of (not x_j), y_j, and the borrow c_j hold, then set the next borrow.
        -- (Note: we use that (a & b) xor (a & c) xor (b & c) == (a & b) or (a & c) or (b & c).)
        b_j1 <- qnot b_j1 `controlled` (x!!j .==. 0) .&&. (y!!j .==. 1)
        b_j1 <- qnot b_j1 `controlled` (x!!j .==. 0) .&&. (b!!j .==. 1)
        b_j1 <- qnot b_j1 `controlled` (y!!j .==. 1) .&&. (b!!j .==. 1)
        -- Now the current difference digit is computed mod 2:
        d_j <- qnot d_j `controlled` (x!!j)
        d_j <- qnot d_j `controlled` (y!!j)
        d_j <- qnot d_j `controlled` (b!!j) 
        b <- return $ overwriteAt (j+1) b_j1 b
        d <- return $ overwriteAt j d_j d
        return (x,y,d,b))

      -- Final difference digit; no carry required:
      let d_l1 = d !! (l-1)
      d_l1 <- qnot d_l1 `controlled` (x!!(l-1))
      d_l1 <- qnot d_l1 `controlled` (y!!(l-1))
      d_l1 <- qnot d_l1 `controlled` (b!!(l-1)) 
      d <- return $ overwriteAt (l-1) d_l1 d

      -- bring the lists’ indexing in line with usual numeric indexing of binary expansons
      x <- return $ reverse x
      y <- return $ reverse y
      d <- return $ reverse d
      b <- return $ reverse b

      return (x,y,d,b))
    
    (\(x,y,d,b) -> do
      (d,d_out) <- qc_copy_fun d
      return ((x,y,d,b), d_out))
  return (x, y, d_out)

-- | Add one 'QDInt' onto a second, in place; i.e. (/x/,/y/) ↦ (/x/,/x/+/y/).  Arguments are assumed to be of equal size.  /O/(ℓ) gates.
q_add_in_place :: QDInt -> QDInt -> Circ (QDInt,QDInt)
q_add_in_place x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y') <- q_add_in_place_qulist x' y'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  return (x, y)

-- | Low-level implementation of 'q_add_in_place': represents integers
-- as big-headian qubit lists.
q_add_in_place_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit])
q_add_in_place_qulist [] [] = return ([], [])
q_add_in_place_qulist [x0] [y0] = do
  y0 <- qnot y0 `controlled` x0
  return ([x0], [y0])
q_add_in_place_qulist x y = do
  let (x0:x_higher) = reverse x  
      (y0:y_higher) = reverse y

  y0 <- qnot y0 `controlled` x0

  ((x0,y0),(x_higher,y_higher)) <- with_computed_fun (x0,y0)
    (\(x0,y0) -> do
      c <- qinit False
      c <- qnot c `controlled` (x0 .==. 1) .&&. (y0 .==. 0)
      return (x0,y0,c))
    (\(x0,y0,c) -> do
      (x_higher,y_higher,c) <- q_add_aux (x_higher) (y_higher) c
      return ((x0,y0,c),(x_higher,y_higher)))
  return (reverse (x0:x_higher), reverse (y0:y_higher))
  where
  -- Aux: add two LITTLE-endian bit strings, and an optional extra 1.
    q_add_aux :: [Qubit] -> [Qubit] -> Qubit -> Circ ([Qubit],[Qubit],Qubit)
    q_add_aux [] [] c = return ([],[],c)
    q_add_aux [x0] [y0] c = do
      y0 <- qnot y0 `controlled` x0
      y0 <- qnot y0 `controlled` c
      return ([x0],[y0],c)
    q_add_aux (x0:xs) (y0:ys) c = do
      y0 <- qnot y0 `controlled` x0
      y0 <- qnot y0 `controlled` c
      ((x0,y0,c),(xs,ys)) <- with_computed_fun (x0,y0,c)
        (\(x0,y0,c) -> do
          c' <- qinit False
          c' <- qnot c' `controlled` (x0 .==. 1) .&&. (y0 .==. 0)
          c' <- qnot c' `controlled` (x0 .==. 1) .&&. (c .==. 1)
          c' <- qnot c' `controlled` (y0 .==. 0) .&&. (c .==. 1)
          return (x0,y0,c,c'))
        
        (\(x0,y0,c,c') -> do
          (xs,ys,c') <- q_add_aux xs ys c'
          return ((x0,y0,c,c'),(xs,ys)))
      return (x0:xs,y0:ys,c)
    q_add_aux _ _ _ = error "q_add_in_place: cannot add integers of different sizes."

-- | Subtract one 'QDInt' from a second, in place; i.e. (/x/,/y/) ↦ (/x/,/y/–/x/).  Arguments are assumed to be of equal size.  /O/(ℓ) gates.
q_sub_in_place :: QDInt -> QDInt -> Circ (QDInt,QDInt)
q_sub_in_place x y = reverse_generic_endo (\(x,d) -> q_add_in_place x d) (x,y)

-- | Low-level version of 'q_sub_in_place': represents integers
-- as big-headian qubit lists.
q_sub_in_place_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit],[Qubit])
q_sub_in_place_qulist x y = reverse_generic_endo (\(x,d) -> q_add_in_place_qulist x d) (x,y)

-- ----------------------------------------------------------------------
-- ** Arithmetic with parameters

-- | Add a parameter 'IntM' and a 'QDInt', into a fresh 'QDInt': (/x/, /y/) ↦ (/y/, /x/+/y/).  The parameter /x/ must be of the same length as /y/, or /x/ can also be of undetermined length.  /O/(ℓ).
q_add_param :: IntM -> QDInt -> Circ (QDInt,QDInt)
q_add_param x1 y = do
  let x = intm_promote x1 y "q_add_param: inputs must have equal length"
  let x' = boollist_of_intm_bh x
  let y' = list_of_xint_bh y
  (y', z') <- q_add_param_qulist x' y'
  let y = xint_of_list_bh y'
  let z = xint_of_list_bh z'
  return (y, z)
  
-- | Low-level implementation of 'q_add_param': represents integers as
-- big-headian qubit lists. Precondition: /x/ and /y/ have the same length.
q_add_param_qulist :: [Bool] -> [Qubit] -> Circ ([Qubit], [Qubit])
q_add_param_qulist x y = do 
  -- Implementation note: compare with q_add_qdint.  The code is almost verbatim identical, but since x consists of parameter Bools, `controlled` will omit gates when appropriate.

  -- Further slight optimisation would be possible, since if x doesn’t have low bits then no low carry qubits are required. 

  -- Implementation note: once we have a general-purpose classical
  -- circuit optimizer (even a primitive one), then this kind of
  -- source-code level optimization will hopefully become moot.
  
  let l = length x

  (y,s_out) <- with_computed_fun y
    (\y -> do
      s <- qinit (replicate l False) -- holds the eventual sum
      c <- qinit (replicate l False) -- holds the carries

      -- bring the lists’ indexing in line with usual numeric indexing of binary expansons
      x <- return $ reverse x
      y <- return $ reverse y
      s <- return $ reverse s
      c <- return $ reverse c

      (y,s,c) <- loop_with_indexM (l-1) (y,s,c) (\j (y,s,c) -> do
        let c_j1 = c !! (j+1)
        let s_j = s !! j
        -- If any two of the current bits x_j, y_j, and the carry c_j are true, then set the next carry:
        -- (Note: we use that (a & b) xor (a & c) xor (b & c) == (a & b) or (a & c) or (b & c).)
        c_j1 <- qnot c_j1 `controlled` (x!!j .&&. y!!j)
        c_j1 <- qnot c_j1 `controlled` (x!!j .&&. c!!j)
        c_j1 <- qnot c_j1 `controlled` (y!!j .&&. c!!j)
        -- Now the current sum digit is computed mod 2:
        s_j <- qnot s_j `controlled` (x!!j)
        s_j <- qnot s_j `controlled` (y!!j)
        s_j <- qnot s_j `controlled` (c!!j) 
        c <- return $ overwriteAt (j+1) c_j1 c
        s <- return $ overwriteAt j s_j s

        return (y,s,c))

      -- Final sum digit; no carry required:
      let s_l1 = s !! (l-1)
      s_l1 <- qnot s_l1 `controlled` (x!!(l-1))
      s_l1 <- qnot s_l1 `controlled` (y!!(l-1))
      s_l1 <- qnot s_l1 `controlled` (c!!(l-1)) 
      s <- return $ overwriteAt (l-1) s_l1 s
     
      x <- return $ reverse x
      y <- return $ reverse y
      s <- return $ reverse s
      c <- return $ reverse c

      return (y,s,c))
    
    (\(y,s,c) -> do
      (s,s_out) <- qc_copy_fun s
      return ((y,s,c), s_out))
  return (y, s_out)

-- | Subtract a parameter 'IntM' from a 'QDInt', into a fresh 'QDInt'.  The 'IntM' cannot be shorter than the 'QDInt' (that would give ill-defined behavior), but can be of undetermined length.  /O/(ℓ).
q_sub_param :: IntM -> QDInt -> Circ (QDInt,QDInt)
q_sub_param x y = q_add_param (-x) y

-- | Add a parameter 'IntM' onto a 'QDInt', in place; i.e. (/x/,/y/) ↦ /x/+/y/. The parameter /x/ must be of the same length as /y/, or /x/ can also be of undetermined length.  /O/(ℓ).
q_add_param_in_place :: IntM -> QDInt -> Circ QDInt
q_add_param_in_place x1 y = do
  let x = intm_promote x1 y "q_add_param_in_place: inputs must have equal length"
  let x' = boollist_of_intm_bh x
  let y' = list_of_xint_bh y
  y' <- q_add_param_in_place_qulist x' y'
  let y = xint_of_list_bh y'
  return y

-- | Low-level implementation of 'q_add_param_in_place': represents
-- integers as big-headian qubit lists. Precondition: /xlist/ and /y/
-- have the same length. Precondition: /x/ and /y/ have the same length.
q_add_param_in_place_qulist :: [Bool] -> [Qubit] -> Circ [Qubit]
q_add_param_in_place_qulist [] [] = return []
q_add_param_in_place_qulist [False] [y0] = return [y0]
q_add_param_in_place_qulist [True] [y0] = do
  y0 <- qnot y0
  return [y0]
q_add_param_in_place_qulist x y = do
  let l = length x

  let (x0:x_higher) = reverse x  
      (y0:y_higher) = reverse y

  y0 <- qnot y0 `controlled` x0

  (y0,y_higher) <- with_computed_fun y0
    (\y0 -> do
      c <- qinit False
      c <- qnot c `controlled` (x0 == 1) .&&. (y0 .==. 0)
      return (y0,c))
    (\(y0,c) -> do
      (y_higher,c) <- q_add_aux x_higher y_higher c
      return ((y0,c),y_higher))
  return (reverse (y0:y_higher))
  where
  -- Aux: add two LITTLE-endian bit strings, and an optional extra 1.
    q_add_aux :: [Bool] -> [Qubit] -> Qubit -> Circ ([Qubit],Qubit)
    q_add_aux [] [] c = return ([],c)
    q_add_aux [x0] [y0] c = do
      y0 <- qnot y0 `controlled` x0
      y0 <- qnot y0 `controlled` c
      return ([y0],c)
    q_add_aux (x0:xs) (y0:ys) c = do
      y0 <- qnot y0 `controlled` x0
      y0 <- qnot y0 `controlled` c
      ((y0,c),ys) <- with_computed_fun (y0,c)
        (\(y0,c) -> do
          c' <- qinit False
          c' <- qnot c' `controlled` (x0 == 1) .&&. (y0 .==. 0)
          c' <- qnot c' `controlled` (x0 == 1) .&&. (c .==. 1)
          c' <- qnot c' `controlled` (y0 .==. 0) .&&. (c .==. 1)
          return (y0,c,c'))
        
        (\(y0,c,c') -> do
          (ys,c') <- q_add_aux xs ys c'
          return ((y0,c,c'),ys))
      return (y0:ys,c)
    q_add_aux _ _ _ = error "q_add_in_place: cannot add integers of different sizes."

-- | Subtract a parameter 'IntM' from a 'QDInt', in place; i.e. (/x/,/y/) ↦ (/x/,/x/-/y/).  /x/ cannot be shorter than /y/.  /O/(/l/) gates.
q_sub_param_in_place :: IntM -> QDInt -> Circ QDInt
q_sub_param_in_place x = q_add_param_in_place (-x)

-- | Multiply a parameter 'IntM' by a 'QDInt', into a fresh 'QDInt'.  The 'IntM' cannot be shorter than the 'QDInt' (that would give ill-defined behavior), but can be of undetermined length.  /O/(ℓ).
q_mult_param :: IntM -> QDInt -> Circ (QDInt,QDInt)
q_mult_param x1 y = do
  let x = intm_promote x1 y "q_add_param: inputs must have equal length"
  let x' = boollist_of_intm_bh x
  let y' = list_of_xint_bh y
  (y', z') <- q_mult_param_qulist x' y'
  let y = xint_of_list_bh y'
  let z = xint_of_list_bh z'
  return (y, z)

-- | Low-level implementation of 'q_mult_param': represents integers as
-- big-headian (qu)bit lists.
q_mult_param_qulist :: [Bool] -> [Qubit] -> Circ ([Qubit],[Qubit])
q_mult_param_qulist [] [] = return ([], [])
q_mult_param_qulist xs ys = do
  let x0 = last xs 
      x_higher = init xs
      y_high = head ys
      y_lower = tail ys
  (y_lower, p_higher)
    <- q_mult_param_qulist x_higher y_lower
  p0 <- qinit False
  let y = y_high:y_lower
      p = p_higher ++ [p0]
  -- Once we have some optimisation, the following line should be replaced by:
  -- (y,p) <- q_add_in_place_qulist y p `controlled` x0
  (y,p) <- if x0 then q_add_in_place_qulist y p else return (y,p)
  return (y, p)  

-- ----------------------------------------------------------------------
-- ** Sign and negation

-- | Negate a (signed) 'QDInt' in place. /O/(ℓ).
q_negate_in_place :: QDInt -> Circ QDInt
q_negate_in_place x = do
  x <- mapUnary qnot x
  x <- q_increment x
  return x

-- | Low-level version of 'q_negate_in_place': represents integers as
-- big-headian qubit lists.
q_negate_in_place_qulist :: [Qubit] -> Circ [Qubit]
q_negate_in_place_qulist = mmap list_of_xint_bh . q_negate_in_place . xint_of_list_bh

-- | Compute the negation of a (signed) 'QDInt'. /O/(ℓ).
-- 
-- (Fixes the minimum value, consistently with Haskell’s @'negate' 'minBound' :: 'Int'@.) 
q_negate_qdint :: QDInt -> Circ (QDInt,QDInt)
q_negate_qdint x = do
  (x,nx) <- qc_copy_fun x
  nx <- q_negate_in_place nx
  return (x,nx)

-- | Compute the absolute value of a (signed) 'QDInt'.  /O/(ℓ).
-- 
-- (Fixes the minimum value, consistently with Haskell’s @'abs' 'minBound' :: 'Int'@.) 
q_abs_qdint :: QDInt -> Circ (QDInt,QDInt)
q_abs_qdint x = do
  let x' = list_of_xint_bh x
  (x', a') <- q_abs_qulist x'
  let x = xint_of_list_bh x'
  let a = xint_of_list_bh a'
  return (x, a)

-- | Low-level implementation of 'q_abs_qdint': represents integers as
-- big-headian qubit lists.
q_abs_qulist :: [Qubit] -> Circ ([Qubit],[Qubit])
q_abs_qulist [] = return ([], [])
q_abs_qulist (x_high:x_lower) = do
  a_high <- qinit False
  (x_lower,a_lower) <- qc_copy_fun x_lower
  a_lower <- mapUnary qnot a_lower `controlled` x_high
  let a = a_high:a_lower
  a <- q_increment_qulist a `controlled` x_high
  return (x_high:x_lower, a)

-- | Compute a 'QDInt' of the same length as the input, with value 1, 0, or –1, depending on the sign of the input.  Analogous to Haskell’s 'signum'.  /O/(ℓ).
q_signum_qdint :: QDInt -> Circ (QDInt,QDInt)
q_signum_qdint x = do
  let x' = list_of_xint_bh x
  (x', a') <- q_signum_qulist x'
  let x = xint_of_list_bh x'
  let a = xint_of_list_bh a'
  return (x, a)

-- | Low-level implementation of 'q_abs_qdint': represents integers as
-- big-headian qubit lists.
q_signum_qulist :: [Qubit] -> Circ ([Qubit],[Qubit])
q_signum_qulist [] = return ([], [])
q_signum_qulist x = do
  let l = length x
  (s_higher, s_low) <- qinit (replicate (l-1) False, False)
  -- first set s = 1
  s_low <- qnot s_low
  -- case x < 0: set s to -1 
  s_higher <- mapUnary qnot s_higher `controlled` (head x)
  -- case x = 0: reset s to 0
  s_low <- qnot s_low `controlled` (x .==. replicate l 0)
  return (x,s_higher ++ [s_low])

-- ----------------------------------------------------------------------
-- ** Comparison

-- | Comparison of two 'QDInt's, considered as unsigned.  /O/(ℓ).
q_le_unsigned :: QDInt -> QDInt -> Circ (QDInt,QDInt,Qubit)
q_le_unsigned x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x',y',le_out) <- q_le_unsigned_qulist x' y'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  return (x,y,le_out)

-- | Low-level implementation of 'q_le_unsigned': represents integers
-- as big-headian qubit lists.
q_le_unsigned_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit],[Qubit],Qubit)
q_le_unsigned_qulist x y = do
  ((x,y),le_out) <- with_computed_fun (x,y) q_le_unsigned_aux
    (\(x,y,(le,garbage)) -> do
      (le,le_out) <- qc_copy_fun le
      return ((x,y,(le,garbage)),le_out))
  return (x,y,le_out)

-- | Same as 'q_le_unsigned', but uncurried, represents integers as
-- big-headian qubit lists, and doesn’t clean up its garbage.  (If
-- each recursive call cleans up its own garbage, gate use becomes
-- /O/(2[sup /n/]).)
q_le_unsigned_aux :: ([Qubit], [Qubit]) -> Circ ([Qubit],[Qubit],(Qubit,[Qubit]))
q_le_unsigned_aux ([], []) = do q <- qinit True; return ([], [], (q,[]))
q_le_unsigned_aux (x_high:x_lower, y_high:y_lower) = do
  -- Classically, one would make the recursive call only if needed; but quantumly, it’s cheaper to do it unconditionally. 
  (x_lower, y_lower, (le_lower, garbage)) <- q_le_unsigned_aux (x_lower, y_lower)
  (x_high,y_high,eq_high) <- q_is_equal x_high y_high
  le <- qinit False
  le <- qnot le `controlled` (x_high .==. 0) .&&. (y_high .==. 1)
  le <- qnot le `controlled` eq_high .&&. le_lower
  return (x_high:x_lower, y_high:y_lower, (le, eq_high:le_lower:garbage))
q_le_unsigned_aux _ = error "q_le // QDInt: cannot compare QDInt’s of different lengths."

-- | Comparison of two 'QDInt's, considered as signed.  Used in @instance 'QOrd' 'QDInt'@.  /O/(ℓ).
q_le_signed :: QDInt -> QDInt -> Circ (QDInt,QDInt,Qubit)
q_le_signed x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x',y',le_out) <- q_le_signed_qulist x' y'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  return (x,y,le_out)

-- | Low-level implementation of 'q_le_signed': represents integers
-- as big-headian qubit lists.

-- Implementation note: defined in terms of q_le_unsigned_aux
q_le_signed_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit],[Qubit],Qubit)
q_le_signed_qulist [] [] = do q <- qinit True; return ([], [], q)
q_le_signed_qulist (x_high:x_lower) (y_high:y_lower) = do
  ((x,y), x_le_y) <- with_computed_fun (x_high:x_lower, y_high:y_lower)
    (\(x_high:x_lower, y_high:y_lower) -> do 
  -- Classically, one would make the recursive call only if needed; but quantumly, it’s cheaper to do it unconditionally. 
      (x_lower, y_lower, (le_lower, garbage)) <- q_le_unsigned_aux (x_lower, y_lower)
      (x_high,y_high,eq_high) <- q_is_equal x_high y_high
      return (x_high, y_high, x_lower, y_lower, le_lower, eq_high, garbage))
    (\(x_high, y_high, x_lower, y_lower, le_lower, eq_high, garbage) -> do
      le <- qinit False
      le <- qnot le `controlled` (x_high .==. 1) .&&. (y_high .==. 0)
      le <- qnot le `controlled` eq_high .&&. le_lower
      return ((x_high, y_high, x_lower, y_lower, le_lower, eq_high, garbage), le))
  return (x,y,x_le_y)
q_le_signed_qulist _ _ = error "q_le // QDInt: cannot compare QDInt’s of different lengths."

-- | Comparison of two 'QDInt's, considered as signed.  Used in @instance 'QOrd' 'QDInt'@. /O/(ℓ).
q_lt_signed :: QDInt -> QDInt -> Circ (QDInt,QDInt,Qubit)
q_lt_signed x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (y',x',y_le_x) <- q_le_signed_qulist y' x'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  x_lt_y <- qnot y_le_x
  return (x,y,x_lt_y)

-- | Text whether a 'QDInt' is nonnegative.   /O/(1)
q_negative :: QDInt -> Circ (QDInt,Qubit)
q_negative x = do
  let (x_high:x_lower) = list_of_xint_bh x
  (x_high,x_neg) <- qc_copy_fun x_high
  return (xint_of_list_bh (x_high:x_lower),x_neg)

instance QOrd QDInt where
  q_less qx qy = do (qx,qy,q) <- q_lt_signed qx qy; return q
  q_leq qx qy = do (qx,qy,q) <- q_le_signed qx qy; return q

-- ----------------------------------------------------------------------
-- ** Multiplication

-- | Multiply two 'QDInt's into a fresh third.  Arguments are assumed to be of equal size.  /O/(ℓ[sup 2]).
q_mult_qdint :: QDInt -> QDInt -> Circ (QDInt,QDInt,QDInt)
q_mult_qdint x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y', z') <- q_mult_qulist x' y'
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  let z = xint_of_list_bh z'
  return (x, y, z)

-- | Low-level implementation of 'q_mult_qdint': represents integers as
-- big-headian qubit lists.
q_mult_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit],[Qubit],[Qubit])
q_mult_qulist [] [] = return ([], [], [])
q_mult_qulist xs ys = do
  let x0 = last xs 
      x_higher = init xs
      y_high = head ys
      y_lower = tail ys
  (x_higher, y_lower, p_higher)
    <- q_mult_qulist x_higher y_lower
  p0 <- qinit False
  let y = y_high:y_lower
      p = p_higher ++ [p0]
  (y,p) <- q_add_in_place_qulist y p `controlled` x0
  return (x_higher ++ [x0], y, p)

-- ----------------------------------------------------------------------
-- ** Division and remainder

-- | Reduce one 'QDInt' modulo another, in place, returning also the
-- integer quotient, i.e. (/x/, /y/) ↦ (/x/ mod /y/, /y/, /x/ div /y/).
-- All inputs and outputs are considered unsigned.  Inputs are assumed
-- to have the same length.  Division by zero returns (/x/,0,-1).
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_moddiv_unsigned_in_place :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt)
q_moddiv_unsigned_in_place x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (r', y', q') <- q_moddiv_unsigned_in_place_qulist x' y'
  let r = xint_of_list_bh r'
  let y = xint_of_list_bh y'
  let q = xint_of_list_bh q'
  return (r, y, q)

-- | Low-level implementation of 'q_moddiv_unsigned_in_place':
-- represents integers as big-headian qubit lists.
-- 
-- Implementation note: this algorithm can be optimized significantly:
-- [scratch] is always 0, so can be optimized away.  This optimisation
-- is much clearer at the circuit level than the code level, however,
-- so is deferred for now.
q_moddiv_unsigned_in_place_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit],[Qubit])
q_moddiv_unsigned_in_place_qulist x y = do
  let l = common_value "q_divrem_unsigned_in_place: arguments must be same length" $ map length [x,y]
  -- Set quot = 0, rem = x.  At each stage, we will have y * quot + rem = x.
  quot_bits <- qinit (replicate l False)
  let rem = x
  with_ancilla_init (replicate (l-1) False) (\scratch -> do
    -- For j from (l-1) to 0, compute the jth bit of the quotient, and update the remainder 
    (y,scratch,quot_bits,rem) <- loop_with_indexM l (y,scratch,quot_bits,rem)
      (\i (y,scratch,quot_bits,rem) -> do
        let j = l-1-i
        -- Getting y*(2^j) requires no quantum operations, just formal re-bracketing of the data:
        let y_init = take j y
            y_2j = (drop j y) ++ (take j scratch)
            scratch_tail = drop j scratch
        -- Compare (2^j) y with rem.  If (2^j)y <= rem, then subtract (2^j)y from rem and flip the jth bit of quot.
        ((y_init,y_2j,rem,quot_bits),()) <- with_computed_fun 
          (y_init,y_2j,rem,quot_bits)
          (\(y_init,y_2j,rem,quot_bits) -> do
            -- first, test that (2^j)y has not overflowed the register, i.e. that y_init is still all 0;
            le1 <- qinit False
            le1 <- qnot le1 `controlled` (y_init .==. replicate (length y_init) False)
            -- next, test that (2^j)y <= rem (this will be correct if (2^j)y has not overflowed yet);
            (y_2j,rem,le2) <- q_le_unsigned_qulist y_2j rem
            return (y_init,y_2j,rem,quot_bits,le1,le2))
          (\(y_init,y_2j,rem,quot_bits,le1,le2) -> do
            -- flip jth bit of quot, if appropriate
            -- (leave rem unchanged for now, since it’s needed to uncompute le1)
            q_j <- qnot (quot_bits !! (l-1-j)) `controlled` le1 .&&. le2
            return ((y_init,y_2j,rem,(overwriteAt (l-1-j) q_j quot_bits),le1,le2),()))
        -- now, subtract (2^j) from rem, if appropriate.
        (y_2j,rem) <- q_sub_in_place_qulist y_2j rem `controlled` quot_bits !! (l-1-j)
        -- Finally, undo the formal re-bracketing:
        let y = y_init ++ (take (l-j) y_2j)
            scratch = (drop (l-j) y_2j) ++ scratch_tail
        return (y,scratch,quot_bits,rem))
    return (rem, y, quot_bits))

-- | Reduce one 'QDInt' modulo another, i.e. (/x/, /y/) ↦ (/x/, /y/, /x/ mod /y/).
-- All inputs and outputs are considered unsigned.  Inputs are assumed
-- to have the same length.  If /y/ = 0, returns (/x/,0,/x/).
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_mod_unsigned :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt)
q_mod_unsigned x y = do
  ((x,y),x_mod_y) <- with_computed_fun (x,y)
    (\(x,y) -> q_moddiv_unsigned_in_place x y)
    (\(x_mod_y, y, x_div_y) -> do
      (x_mod_y, x_mod_y_out) <- qc_copy_fun x_mod_y
      return ((x_mod_y, y, x_div_y), x_mod_y_out))
  return (x,y,x_mod_y)

-- | Integer division of two 'QDInt's, returning the quotient and remainder,
-- i.e. (/x/,/y/) ↦ (/x/,/y/,/x/ div /y/,/x/ mod /y/). 
-- All inputs and outputs are considered unsigned.
-- Inputs are assumed to have the same length.
-- Division by zero returns (-1,/x/).
--
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_divrem_unsigned :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt,QDInt)
q_divrem_unsigned x y = do
  (x,rem) <- qc_copy_fun x
  (rem,y,quot) <- q_moddiv_unsigned_in_place rem y
  return (x,y,quot,rem)

-- | Integer division of two 'QDInt's, returning just the quotient.  
-- All inputs/outputs considered unsigned.
-- Inputs are assumed to have the same length.
-- Division by zero returns –1.
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.

-- Implementation note: We use a multiplication to do the
-- uncomputation of the remainder, because this uses fewer total gates
-- than uncomputing divrem would.
q_div_unsigned :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt)
q_div_unsigned x y = do
  (x,y,quot,rem) <- q_divrem_unsigned x y
  ((x,y,quot),()) <- with_computed_fun (x,y,quot) 
    (\(x,y,quot) -> do
      (y,quot,y_quot) <- q_mult y quot
      (x,y_quot,rem_copy) <- q_sub x y_quot
      return (x,y,quot,y_quot,rem_copy))
    (\(x,y,quot,y_quot,rem_copy) -> do
      rem_copy <- qc_uncopy_fun rem_copy rem
      return ((x,y,quot,y_quot,rem_copy),()))
  return (x,y,quot)

-- | Low-level version of 'q_div_unsigned': represents integers as
-- big-headian qubit lists.
q_div_unsigned_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit], [Qubit])
q_div_unsigned_qulist x' y' = do
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  (x, y, q) <- q_div_unsigned x y            
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  let q' = list_of_xint_bh q
  return (x', y', q')

-- | Integer division of two 'QDInt's into a fresh third, rounding towards –∞. 
-- The first argument is the numerator, and is assumed to be signed.
-- The second argument is the denominator, and is assumed to be unsigned.
-- The output is signed.
-- Inputs are assumed to have the same length.
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_div :: QDInt -> QDInt -> Circ (QDInt, QDInt, QDInt)
q_div x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y', q') <- q_div_qulist x' y'            
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  let q = xint_of_list_bh q'
  return (x, y, q)

-- | Low-level implementation of 'q_div': represents integers as
-- big-headian qubit lists.
q_div_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit], [Qubit])
q_div_qulist [] [] = return ([], [], [])
q_div_qulist (x_high:x_lower) y = do
  -- Writing x/y for unsigned division, we have x `div` y = 
  -- if x>0 then x/y else ~((~x)/y)
  -- where ~x := -1-x --- i.e. the bitflip of x.
  --
  -- We do this with an ancilla (which can easily be optimized away): 
  with_ancilla_init False (\fake_x_high -> do
    -- flip x if x<0:
    x_lower <- mapUnary qnot x_lower `controlled` x_high
    let x' = fake_x_high:x_lower
    -- x' now holds (if x > 0 then x else ~x) 
    (x',y,quot) <- q_div_unsigned_qulist x' y
    -- restore x and flip quotient, if necessary
    let fake_x_high:x_lower = x'
    x_lower <- mapUnary qnot x_lower `controlled` x_high
    quot <-  mapUnary qnot quot `controlled` x_high
    return (x_high:x_lower, y, quot))
q_div_qulist _ _ = error "q_div: arguments must have same length."

-- | Integer division of two 'QDInt's into a fresh third, rounding towards 0. 
-- The first argument is the numerator, and is assumed to be signed.
-- The second argument is the denominator, and is assumed to be unsigned.
-- The output is signed.
-- Inputs are assumed to have the same length.
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_quot :: QDInt -> QDInt -> Circ (QDInt, QDInt, QDInt)
q_quot x y = do
  let x' = list_of_xint_bh x
  let y' = list_of_xint_bh y
  (x', y', q') <- q_quot_qulist x' y'            
  let x = xint_of_list_bh x'
  let y = xint_of_list_bh y'
  let q = xint_of_list_bh q'
  return (x, y, q)

-- | Low-level implementation of 'q_quot': represents integers as
-- big-headian qubit lists.
q_quot_qulist :: [Qubit] -> [Qubit] -> Circ ([Qubit], [Qubit], [Qubit])
q_quot_qulist [] [] = return ([], [], [])
q_quot_qulist (x_high:x_lower) y = do
  -- Writing x/y for unsigned division, we have x `div` y = 
  -- if x>0 then x/y else -((-x)/y)
  --
  -- We do this with an ancilla (which can easily be optimized away): 
  with_ancilla_init False (\fake_x_high -> do
    -- flip x if x<0:
    x_lower <- q_negate_in_place_qulist x_lower `controlled` x_high
    let x' = fake_x_high:x_lower
    -- x' now holds (if x > 0 then x else ~x) 
    (x',y,quot) <- q_div_unsigned_qulist x' y
    -- restore x and flip quotient, if necessary
    let fake_x_high:x_lower = x'
    x_lower <- q_negate_in_place_qulist x_lower `controlled` x_high
    quot <- q_negate_in_place_qulist quot `controlled` x_high
    return (x_high:x_lower, y, quot))
q_quot_qulist _ _ = error "q_quot: arguments must have same length."

-- | Integer division of two 'QDInt's, returning the quotient,
-- assuming that the second exactly divides the first and throwing an error otherwise.  
-- All inputs/outputs considered unsigned.
-- Inputs are assumed to have the same length.
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_div_exact_unsigned :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt)
q_div_exact_unsigned x y = do
  (x,y,quot,rem) <- q_divrem_unsigned x y
  qterm 0 rem
  return (x,y,quot)

-- | Integer division of two 'QDInt's, returning the quotient,
-- assuming that the second exactly divides the first.  
-- The first argument is the numerator, considered signed.
-- The second argument is the denominator, considered unsigned.
-- The output is signed.
-- 
-- /O/(ℓ[sup 2]) gates, /O/(ℓ) qubits.
q_div_exact :: QDInt -> QDInt -> Circ (QDInt, QDInt,QDInt)
q_div_exact x y = do
  (x,(y,quot)) <- with_computed_fun x
    (\x -> do
      (x,x_neg) <- q_negative x
      x' <- q_negate_in_place x `controlled` x_neg
      return (x',x_neg))
    (\(x',x_neg) -> do
      (x',y,quot) <- q_div_exact_unsigned x' y
      quot <- q_negate_in_place quot `controlled` x_neg
      return ((x',x_neg),(y,quot)))
  return (x,y,quot)

-- ----------------------------------------------------------------------
-- * The QNum type class

-- $ This section defines a quantum analogue of Haskell’s 'Num' type
-- class.  See "Quipper.Internal.QClasses" for the general philosophy behind
-- quantum type classes.
-- 
-- The type class 'QNum' corresponds to Haskell’s 'Num', with methods
-- including 
-- 
-- > add_q :: (QNum qa) => qa -> qa -> Circ (qa,qa,qa), 
-- 
-- and so on.
-- 
-- Note: type class instances don't show up in the documentation. View
-- the source code to see them.

-- | Quantum analogue of Haskell’s 'Num' type class. This provides
-- basic addition, subtraction, multiplication, sign operations, and
-- conversion from integers.

class (QData qa) => QNum qa where
  -- | Add two quantum numbers into a fresh one. The arguments are
  -- assumed to be of equal size. The 'QDInt' instance uses /O/(ℓ)
  -- gates, both before and after transformation to Toffoli.
  q_add :: qa -> qa -> Circ (qa,qa,qa)
  
  -- | Multiply two quantum numbers into a fresh third. The arguments
  -- are assumed to be of equal size.  The 'QDInt' instance is
  -- /O/(ℓ[sup 2]).
  q_mult :: qa -> qa -> Circ (qa,qa,qa)
  
  -- | Subtract two quantum numbers into a fresh one. The arguments
  -- are assumed to be of equal size. The 'QDInt' instance uses /O/(ℓ)
  -- gates, both before and after transformation to Toffoli.
  q_sub :: qa -> qa -> Circ (qa,qa,qa)
  
  -- | Compute the absolute value of a (signed) quantum number. The
  -- 'QDInt' instance is /O/(ℓ).
  q_abs :: qa -> Circ (qa,qa)
  
  -- | Compute the negation of a (signed) quantum number. The 'QDInt'
  -- instance is /O/(ℓ).
  q_negate :: qa -> Circ (qa,qa)
  
  -- | Compute a quantum number of the same precision as the input,
  -- with value 1, 0, or –1, depending on the sign of the
  -- input. Analogous to Haskell’s 'signum'. The 'QDInt' instance is
  -- /O/(ℓ).
  q_signum :: qa -> Circ (qa,qa)
  
  -- | Convert a 'QDInt' to a quantum number. For the 'QDInt'
  -- instance, this is just a copy operation.
  q_fromQDInt :: QDInt -> Circ (QDInt,qa)
   
instance QNum QDInt where
  q_add = q_add_qdint
  q_mult = q_mult_qdint
  q_sub = q_sub_qdint
  q_abs = q_abs_qdint
  q_negate = q_negate_qdint
  q_signum = q_signum_qdint
  q_fromQDInt = qc_copy_fun

-- ----------------------------------------------------------------------
-- ** Specialized functions

-- | Euclid's extended algorithm. 'q_ext_euclid' /a/ /b/: returns
-- (/a/,/b/,/x/,/y/,/d/) such that /ax/ + /by/ = /d/ = gcd(/a/,/b/).
q_ext_euclid :: QDInt -> QDInt -> Circ (QDInt,QDInt,QDInt,QDInt,QDInt)
q_ext_euclid a b = do
  let l = common_length "q_ext_euclid: inputs must have equal length" [a,b]

  -- Extended Euclidean algorithm:
  --
  -- compute three series of numbers (a_i), (x_i), (y_i), such that
  --
  -- a_0 = a, a_1 = b,
  -- a_{i+2} = a_i `mod` a_{i+1},
  -- and always a_i = (x_i * a) + (y_i * b),
  --
  -- When we first get a_{i+2} = 0, we know that a_{i+1} is gcd(a_0,a_1),
  -- so we save a_{i+1}, x_{i+1} and y_{i+1} for output.
  --
  -- The bound on the length of this loop is standard.

  ((a,b),(x,y)) <- with_computed_fun (a,b) (\(a,b) -> do
    x0 <- qinit (intm l 1)
    y0 <- qinit (intm l 0)
    x1 <- qinit (intm l 0)
    y1 <- qinit (intm l 1)
    done_yet <- qinit False
    x_final <- qinit (intm l 0)
    y_final <- qinit (intm l 0)

    (stuff1, stuff2, stuff3, (done_yet, x_final, y_final))
            <- loopM (euclid_bound l) ((a,b,x0,x1,y0,y1),[],[],(done_yet,x_final,y_final))
      (\((a_i, a_i1, x_i, x_i1, y_i, y_i1), quots_scratch, tests_scratch, (done_yet, x_final, y_final)) -> do
        (a_i2, a_i1, q_i) <- q_moddiv_unsigned_in_place a_i a_i1
        ((x_i1, y_i1, q_i), (x_i2, y_i2)) <- with_computed_fun
          (x_i1, y_i1, q_i)
          (\(x_i1, y_i1, q_i) -> do
            (q_i, x_i1, qx) <- q_mult q_i x_i1
            (q_i, y_i1, qy) <- q_mult q_i y_i1
            return (x_i1, y_i1, q_i, qx, qy))
          (\(x_i1, y_i1, q_i, qx, qy) -> do
            (qx,x_i2) <- q_sub_in_place qx x_i
            (qy,y_i2) <- q_sub_in_place qy y_i
            return ((x_i1, y_i1, q_i, qx, qy),(x_i2,y_i2)))
        done_this_time <- qinit False
        done_this_time <- qnot done_this_time `controlled` (a_i2 .==. 0) .&&. (done_yet .==. False)
        (x_i1,x_final) <- controlled_not x_final x_i1 `controlled` done_this_time
        (y_i1,y_final) <- controlled_not y_final y_i1 `controlled` done_this_time
        done_yet <- qnot done_yet `controlled` done_this_time
        return ((a_i1, a_i2, x_i1, x_i2, y_i1, y_i2),
                (q_i:quots_scratch), (done_this_time:tests_scratch), 
                (done_yet, x_final, y_final)))
    qterm True done_yet -- Assert that we really did reach the end of the algorithm. 
    return (x_final,y_final,(stuff1, stuff2, stuff3)))

    (\(x_final, y_final, stuff) -> do
      (x_final,x) <- qc_copy_fun x_final
      (y_final,y) <- qc_copy_fun y_final
      return ((x_final, y_final, stuff),(x,y)))

  ((a,b,x,y),gcd) <- with_computed_fun
    (a,b,x,y)
    (\(a,b,x,y) -> do
      (a,x,ax) <- q_mult a x
      (b,y,by) <- q_mult b y
      return (a,b,x,y,ax,by))
    (\(a,b,x,y,ax,by) -> do
      (ax,by,gcd) <- q_add ax by
      return ((a,b,x,y,ax,by),gcd))

  return (a,b,x,y,gcd)

  where
  -- Classical bound on running time (Lamé’s theorem), assuming b < a:
  -- t <= 1 + log_{phi} b <= 1 + l / (log_2 phi)
  --  
  -- We allow one extra step since we may have a < b, in which case first step swaps them.
    euclid_bound :: Int -> Int
    euclid_bound l = 2 + ceiling ( (fromIntegral l) / (logBase 2 phi) )
    phi = (1+sqrt(5))/2

-- ----------------------------------------------------------------------
-- * Lifting of arithmetic functions

-- $LIFTING 
-- This sections provides templates for lifting various arithmetic
-- functions in connection with the @build_circuit@ keyword. This
-- extends the liftings given in "Quipper.Internal.CircLifting" to operations
-- of the 'Num' type class.

-- | Quantum lifting of the '+' operator.
template_symb_plus_ :: (QNum qa) => Circ (qa -> Circ (qa -> Circ qa))
template_symb_plus_ = return $ \qx -> return $ \qy -> do (qx,qy,qz) <- q_add qx qy; return qz

-- TODO: complete templates of methods.