Haskell code

{-# 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.