isQ Compiler Docs

Description	Document for isQ Programming language and isQ Compiler.
Author(s)	isQ Team
Repository	https://github.com/isQ-Team/isQ-Compiler
Version	0.2.8+f48ba067f

This box will disappear when printing mkdocs-print-site-plugin

This page has combined all site pages into one. You can export to PDF using File > Print > Save as PDF.

See also export to PDF and export to standalone HTML.

Welcome to isQ docs

This is a brief introduction to our quantum programming software, isQ, a part of quantum software toolkit of Institute of Software, Chinese Academy of Science. The compiler of isQ accepts isQ programs as input, and produces low-level codes written in several kinds of intermediate representations or instruction sets (depending on the specific compiler options). The compilation results could be run on one of several backends: either on real superconducting hardware or on isQ simulators.

As the programming software in a relatively early version, if you find any bugs, please feel free to inform us and we will be very grateful. We also appreciate any kind suggestions. :)

Please head to Installation to get started.

User Manual

Install

Binary Tarball

We provide a traditional binary tarball based on nix-user-chroot. Note that the tarball requires Linux kernel user namespaces to work.

A list of all released versions of isQ compiler binaries can be seen at https://www.arclightquantum.com/isq-releases/isqc-standalone/.

VERSION=0.2.8
ARCH=x86_64-unknown-linux-gnu
# Create empty directory for isQ installation.
mkdir isqc && cd isqc
# Check if user namespace is supported for your Linux kernel.
# If not, see FAQ below.
unshare --user --pid echo YES
# Download and unpack tarball.
TARBALL=isqc-${VERSION}-${ARCH}.tar.gz
wget https://www.arclightquantum.com/isq-releases/isqc-standalone/${VERSION}/${TARBALL}
wget https://www.arclightquantum.com/isq-releases/isqc-standalone/${VERSION}/${TARBALL}.sha256
sha256sum -c ${TARBALL}.sha256
tar -xvf ${TARBALL}
# Now isQ is here.
./isqc --version

Nix Flake (Recommended)

isQ is built with Nix Flakes, making it super easy to obtain when you have Nix installed:

# Add isQ binary cache to Cachix to prevent building from source.
nix-shell -p cachix --run "cachix use arclight-quantum"
# Enter the environment with isQ installed.
nix shell github:isQ-Team/isQ-Compiler
# Now isQ is placed in $PATH.
isqc --version

Or you may create a project folder pinned to a compiler version.

nix flake new --template github:isQ-Team/isQ-Compiler hello-quantum
cd hello-quantum && nix develop

Docker Container

We provide two Docker images with isQ compiler builtin: one for normal users providing a full Ubuntu environment, and the other for professional Docker users with only binary files necessary for isQ.

# Ubuntu-based Docker image.
docker run -it arclightquantum/isqc:ubuntu-0.0.1 bash
isqc --version # Run in container.
# Binary only Docker image.
docker run --rm -v $(pwd):/workdir arclightquantum/isqc:0.0.1 isqc --version

Frequently Asked Questions

unshare failed and nix-user-chroot cannot be used.

Q:

Error occurs while running unshare:

user@server:~$ unshare --user --pid echo YES
unshare: unshare failed: Operation not permitted

or error occurs while running isqc:

user@server:~/isqc$ ./isqc
thread 'main' panicked at src/main.rs:124:70:
unshare failed: Sys(EPERM)
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

A:

nix-user-chroot requires unpriviledege user namespaces to work.

• Kernel version must be >=3.8.
• Follow the guides here to enable user namespaces, roughly:
  • Make sure CONFIG_USER_NS=y is set in kernel compile options.
  • (Suggested for RedHat/CentOS users): make sure user.max_user_namespaces is not zero by running:

    cat /proc/sys/user/max_user_namespaces
    

If you still cannot get user namespaces to work (e.g. you're in a container environment), an alternative is to unpack the tarball at the root directory.

binary_path=/usr/bin/isqc # Set your installation path.
# First unpack the tarball.
tar -xvf isqc.tar.gz
# This will move the `nix` folder to the root directory.
# Note: this may conflict with your Nix installation if you already have Nix installed!
cp -r ./nix /
# isqc should be located at path like:
# /nix/store/hps1c4vap5zc8nkdq1yshpqg9mm3aqd2-isqc/bin/isqc
# This path resides in our `isqc` entrypoint.
# The line below extracts the path from our entry-point script.
isqc_path=$(perl -ne 'print "$1\n" if /\s(\/nix\/store\/.*\/bin\/isqc)/' ./isqc) 
# Do a traditional installation
echo '#!/usr/bin/env bash' > $binary_path
echo "$isqc_path \"\$@\"" >> $binary_path
chmod +x $binary_path
# Remove the unpacked files.
rm -rf ./isqc ./nix 
# Test
isqc --version

AppArmor constraining unprivileged user namespace.

Q:

Unshare succeeded, but error occurs while running isqc:

user@server:~/isqc$ ./isqc --version
thread 'main' panicked at src/main.rs:138:43:
failed to list /nix directory: Os { code: 13, kind: PermissionDenied, message: "Permission denied" }
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Kernel dmesg outputs be like:

[ 1804.496798] audit: type=1400 audit(1718885622.496:241): apparmor="AUDIT" operation="userns_create" class="namespace" info="Userns create - transitioning profile" profile="unconfined" pid=5921 comm="nix-user-chroot" requested="userns_create" target="unprivileged_userns"
[ 1804.497103] audit: type=1400 audit(1718885622.496:242): apparmor="DENIED" operation="open" class="file" profile="unprivileged_userns" name="/" pid=5921 comm="nix-user-chroot" requested_mask="r" denied_mask="r" fsuid=1000 ouid=0

This problem is first reported on Ubuntu 24.04.

A:

Default security profiles of some distros using AppArmor has put stricter restrictions on unprivileged user namespaces, including but not limited to:

• Denying capabilities required to mount /nix.
• Denying accessing / from user namespaces.
• Denying mounting.

Details can be seen here: https://gitlab.com/apparmor/apparmor/-/wikis/unprivileged_userns_restriction

The simplest way to fix this is to disable AppArmor restriction on unprivileged user namespaces:

echo 0 | sudo tee /proc/sys/kernel/apparmor_restrict_unprivileged_userns

Warning

While the line above provides a hands-on workaround by disabling some functionalities of AppArmor, you may want to rollout finer-grain AppArmor policies if you care about safety.

Usage

• Usage
  • Compile
    • compile to qir
    • compile to qcis
  • Simulate
    • simulate qir
    • simulate qcis
    • result
    • show possibilites
  • Run

Compile

The isQ compiler receives the isQ source file as input, and the instruction format is as follows:

isqc compile [options] <Input>

options:

• --emit <FORMAT> : output content format, format value can be mlir, mlirqir, llvm, mlir-optimized, binary, out [default: out]
• -o, --output <OUTFILE> : output file
• -O, --opt-level <OPT_LEVEL> : llvm opt-level such as -O1, -O2, -O3 etc.
• --target <TARGET_IR> : target ir, now support qir, open-qasm3, qcis [default: qir]
• --qcis-config <MAPPING_CONFIG_FILE> : qcis mapping config file
• -I, --inc-path <INC_PATH>: library path of isQ used in the source code, a default path is $ISQ_ROOT/lib
• -i <INT>: int type paramter, which is used when compiled to qcis
• -d <DOUBLE> double type paramter, which is used when compiled to qcis

compile to qir

QIR is a new Microsoft-developed intermediate representation for quantum programs. Users can compile source file to qir using the following command

# compile to qir, default output file is source.so
isqc compile source.isq

or may compile to intermediate results, such as mlir

# compile to mlir, default output file is source.mlir
isqc compile source.isq --emit mlir

compile to qcis

QCIS is specially tailored for the superconducting quantum hardware at the University of Science and Technology of China. Users can compile source file to qcis using the following command

# compile to qcis, default output file is source.qcis
isqc compile --target qcis source.isq

Qcis instructions can run on real superconducting hardware, so isQ compiler provides qubit-mapping function. If want to use this feature, users should feed a configuration file mapping_config_file. The file is usually in JSON format, and needs to include the following fields:

{
// required
    "qbit_num": 12, // qubit number on hardware
    "topo": [[1,2],[2,3],[3,4],[4,5],[5,6],[6,7],[7,8],[8,9],[9,10],[10,11],[11,12]], // topology of hardware
// option 
    "init_map": "simulated_annealing" // the way to get init mapping
}

and the configuration file needs to be specified through --qcis config

isqc compile --target qcis --qcis config mapping_config_file source.isq

QCIS hardware does not support feedback control, so there are some restrictions on the usage of isQ

• can not use reset
• can not print
• can not use measurement results as right value, like assignment, condition

procedure main(){
    // These operations will report errors
    qbit q;
    int a = M(q);
    print a;
}

isQ supports compiling with parameters. When compiling to QCIS, users need pass parameter values in order to generate a definite circuit.

qbit q;

procedure main(int x[], double d[]){
    if (x[0] == 1){
        X(q);
    }else{
        Rx(d[0], q);
    }
    M(q);
}

compile using the following command

# pass values when compiling
isqc compile --target qcis -i 0 -d 1.23 source.isq

Simulate

The simulator of isQ provides the simulation of qir or qcis. The instruction format is as follows

isqc simulate [options] <Input>

options:

• --shots <NUM>: simulate times, default is 1
• --debug: use debug mod and get the print result
• --cuda <NUM>: use gpu, NUM is the number of qubits to be simulated
• -i <INT>: int type paramter
• -d <DOUBLE>: double type paramter
• --qcis: simulate qcis
• --qn <QN>: maximum qubit number, default is 25
• --probs: show the probabilites of measurement results

simulate qir

The input must be a .so file generated by compiler, and can use as follows:

# simulate in cpu
isqc simulate ./source.so
# simulate in gpu
isqc simulate --cuda 10 ./source.so
# run 10 times and use debug mod
isqc simulate --shots 10 --debug ./source.so

isQ supports compiling with parameters and passing values at runtime. The following program is compiled as qir with parameters, and users can pass values during simulation through -i, -d. isQ will gather values and generate int and double arrays (can be empty).

qbit q;

procedure main(int x[], double d[]){
    if (x[0] == 1){
        X(q);
    }else{
        Rx(d[0], q);
    }
    M(q);
}

simulate using the following command

# pass values when simulating
isqc simulate -i 0 -d 1.23 source.so

simulate qcis

isQ simulator can also simulate qcis through --qcis. The input must be a .qcis file and isQ will check the format of the QCIS instructions.

# simulate qcis, 1000 times
isqc simulate --qcis --shots 1000 ./source.qcis

result

The output is the statistical result of measurement in source.isq, like {"00": 4, "11": 6}. The results of all measurement operations will be added to the result string, even if it is an assignment operation like int a = M(q);

import std;

procedure main(){
    qbit q[2];
    H(q[0]);
    ctrl X(q[0],q[1]);
    int a = M(q[0]);
    int b = M(q[1]);
    print a - b;
}

The result of above isQ program can be {"11": 47, "00": 53}. In this program, there is a print operation, the result of this operation will not output until using debug mode through --debug. In debug mode, print result of each simulation will output to stderr, the format like {ith simulate print: 0}. When the shot_num is large, there will be many print result. It is suitable that redirect the stderr to a file. Taking the above program as an example:

# redirect print result to a file
isqc simulate --shots 100 --debug ./bell.so 2>res.txt

output

{"11": 47, "00": 53}

res.txt

0th simulation print:
0
1th simulation print:
0
...
99th simulation print:
0

show possibilites

Instead of repeating the simulation with many shots, isQ also supports generating the theory possibilities of measurement results. This can be done with the --probs argument. Then, isQ prints the possibilities of measurement results in an array. The index of the array is the measurement result, whereas the earlier measured qubit is the higher bit. Note that the measured qubits must be global. The reason is that local qubits will be deallocated at the end of execution, so the simulator cannot find out the possibilities.

The following program, prob.isq, is a simple example.

import std;

qbit q[3];

procedure main(){      
    H(q[0]);
    X(q[1]);
    M(q[1]);
    M(q[0]);
}

Two qubits are measured, i.e., q[1] and q[0]. The measurement results have four combinations:

Measurement result	q[1]	q[0]	Possibility
00	0	0	0
01	0	1	0
10	1	0	0.5
11	1	1	0.5

Simulate the previous program with command

isqc simulate --probs ./prob.so

The result is

[0.0, 0.0, 0.5, 0.5]

Run

Users can use isQ’s run command to compile and simulate quantum programs (compile to qir)

    isqc run [options] <Input>

options

• --shots <NUM>: simulate times, default is 1
• --debug: use debug mod and get the print result
• --qn <QN>: maximum qubit number, default is 25
• --probs: show the probabilites of measurement results

isQ will compile source file to qir first, and then simulate. The options will be forwarded to isqc simulate.

# simulate in cpu
isqc run --shots 100 ./source.isq

isQ Programming Language

This section mainly introduces isQ's grammar which is divided into the following parts:

• Package
• Import
• Types
• Classical Operation
• Quantum Operation
• Procedure
• Control Flow
• Modifier
• Oracle
• Parameter

Package

Each isQ file belongs to a package. The root of the package is specified with keyword package at the begining of the file. For example, the first lines of /some/path/aaa/bbb/aaa/ccc.isq are

// /some/path/aaa/bbb/aaa/ccc.isq
package aaa;

It means that the root of the package is /some/path/aaa/bbb/aaa. The root is matched with file or folder name in the abosolute path of the isQ file, from lower to higher hierachies. In this example, there are two folders of the name aaa. The root is the one closer to the file. If the specified package name cannot match any file or folder name, the compiler will report an error. If no package is specified, the default package name is the file name.

The global variables and procedures can be referred to with qualified names from another isQ file, after being imported. Continue with the previous example,

// /some/path/aaa/bbb/aaa/ccc.isq
package aaa;
int g;

The global variable g can be referred to as aaa.ccc.g. This qulified name is constructed with the relative path from the package root (replacing / with .) plus the variable name.

Import

An isQ file can use the variables and procedures defined in another isQ file by importing it. The syntax is like this example:

// /some/path/aaa/bbb/ddd/eee.isq
package ddd;
import aaa.ccc; // /some/path/aaa/bbb/aaa/ccc.isq
import std;     // ${ISQ_ROOT}/lib/std.isq

The root of this file is /some/path/aaa/bbb/ddd which is placed in the same folder (i.e., /some/path/aaa/bbb/) as the example in the Package part. It specifies the imported file with the relative path from the parent folder of the root, replacing / with . and omiting the .isq suffix. Therefore, the above line imports /some/path/aaa/bbb/aaa/ccc.isq. The user can also uses the relative path from the isQ library, which is specified by environmental variable ${ISQ_ROOT}/lib. A standard library file, std.isq, is placed in this folder. It contains the defination of basic gates such as H and X. Therefore, it should be imported by nearly all the isQ files.

Once a file is imported, its global variables and procedures (except main) can be used directly. For example, the eee.isq mentioned above can use g direclty in its program. However, if multiple variables that are defined in different imported files share the same name, they have to be referred to using qualified names, e.g., aaa.ccc.g and some.other.g.

Types

primitive types

isQ mainly supports four primitive types:

• int: a 64-bit signed integer, e.g. -1, 1;
• bool: a boolean value that can be either true or false;
• double: a double precision floating point number, e.g. -0.1, 3.1415926;
• unit: the return type of a procedure that returns nothing;
• qbit: an opaque type that represents a quantum bit;

Users can use these type to define variables anywhere. All qubits are set to default value |0>. For example, we could define variables like:

// global variable
int a, b = 4; // int variables
double d = pi; // double variable with initial value 3.14159...
qbit q; // qbit

procedure main() {    
    // local variable with initial value 3 
    int c = 3;
    ...
}

Note that we have defined a keyword pi to represent the value of $\pi$.

array

All the four primitive type supports one-dimensional array. When the array is defined with intialization values, no length should be provided. Otherwise, a length must be provided. The length of a global array must be a positive integer, while for a local variable, any int expression can be used as a length. The length of an array can be obtained using the operator .length.

qbit q[3]; // the length of a global array must be a positive integer

unit main() {
    // when an array is initilized, the length should not be specified
    int a[] = [1, 2]; 

    // a local variable may use expressions as length
    int b[a.length];
    ...
}

A slice of an array can be obtained by appending a [start:end:step] structure after the array. It is a subview of the original array, starting from the start position, increasing with step, and till the end (excluding). For example, a[2:4:1] refers to a[2] and a[3]. The three parts can be omitted, and the default values will be used. The default values of start, end, and step are 0, the length of the array, and 1, respectively. When the step is omitted, the second : can be omitted as well. For example, a[:2] represents the first two elements of a. The step can be negtive integers. In that case, any field cannot be omitted. For example, a[2:0:-1] refers to a[2] and a[1].

Classical Operation

isQ supports the following operators for classical types:

• Arithmetic operators: +, -, *, /, ** (power), %
• Comparison operators: ==, !=, >, <, >=, <=
• Logical operators: &&, and, ||, or, !, not
• Bitwise operators: &, |||, ^(XOR)
• Shift operators: >>, <<

Their semantics, operand types and precedence are consistent with common programming languages such as C and Python.

isQ has built-in automatic type conversion, which can convert bool to int and int to double. Specifically, true is converted to 1 and false to 0. Moreover, isQ provides a print command to print int and double values. Users can do arithmetic and print like this:

unit main() {
    int a = 2 * (3 + 4) % 3;
    double d = 3.14 * a;

    print a;
    print d;
}

In addition, isQ provides an assert command, which operates with a bool value. If the value evaluates true during runtime, this command is omitted. Otherwise, the program aborts and reports this event. For example,

unit main() {
    assert true;    // OK
    assert(3 == 4); // Causing a program abort
}

Note that putting paratheses around a bool value results in the same bool value.

Pure classical expressions

Warning

This feature is highly experimental and subject to changes. In particular, evaluation of compile-time pure classical expressions may not terminate.

Since 0.2.7, isQ introduces a special type of pure (i.e. side-effect-free) classical computation, enabling users to write classical computation in a functional-programming style, in particular, enabling users to express more powerful compile-time computation.

The grammar of pure classical expressions is inspired by languages from ML family. Basic arithmetic operations are supported. Moreover, let-binding is allowed but the types must be annotated explicitly.

// isqc/tests/input/lambda/let_add.isq
unit main(){
    int a = let b: int = 1 in b+b;
    print a; // print 2.
}

Pure classical expressions allows users to define inline, recursive lambdas (a.k.a. functions) by using letrec. Again, all types must be annotated explicitly. By now lambdas are second-class and cannot capture any closure.

// isqc/tests/input/lambda/factorial.isq
unit main()
{
    int a = letrec factorial(n : int) -> int = 
            if n <= 0 then 1
            else n * factorial(n - 1)
        in factorial(4);
    print a;
}

What makes pure classical expressions powerful is that they can be treated as compilation-time constants, therefore be used in combination with templates.

unit fun<int N>(){
    print(N);
}
unit fun::<6>(){
    print(6666);
}
unit main() {
    fun::<(letrec collatz(x: int, n: int)-> int = 
        if x==1 
            then n 
            else if (x%2==0)
                then collatz(x/2, n+1)
                else collatz(x*3+1, n+1)
        in collatz(10, 0))>(); // template argument evaluates to 6.
    // 6666 will be printed.
}

Quantum Operation

The quantum operation in isQ is simple, users can apply a gate or do measurement on qubits in a way similar to function call. All quantum operations are inside procedures and all qubits used are defined beforehand. When doing measurement, the result must be stored in an int variable.

basic operation

isQ supports some basic gate: X, Y, Z, H, S, T, Rx, Ry, Rz, CNOT, Toffoli, U3 (see this link for details), and two non-unitary operation: M(measure), |0>(set qubit to |0>). Users can directly use these gates like this:

qbit q[2];
unit main() {
    H(q[0]);
    CNOT(q[0], q[1]);
    bool x = M(q[0]);
    bool y = M(q[1]);
}

qbit array and slice can be used as the parameters of quantum gates and measurement. In this case, it represents applying the gates or measurement to all the qubits in the array or slice. For multiple-qubit gates, it means applying gates to each group of qubits. The number of gates depends on the shortest qbit array or slice. For example:

    qbit p[3], q[3];
    H(p);           // = H(p[0]); H(p[1]); H(p[2]);
    CNOT(p, q[:2]); // = CNOT(p[0], q[0]); CNOT(p[1], q[1]);
    int x = M(q);

The measurement result of a qbit array or slice is an int value where the measurement result of the first qbit is used as the lowest-order bit. In the previous example, if the measurement result of q[0], q[1], and q[2] are true, true, and false, respectively, x would be 3 (i.e., 011).

defgate

isQ allows users to define gate by using keyword defgate and a matrix. An element in the matrix can be an int/double/complex value or an arithmetic expression. Complex value can be written in a way similar to Python. There are three points to note:

• The user-defined gate should be unitary and its size should be a power of 2.
• Users should define gate outside procedures and then use it inside.
• A user-defined gate's name can not conflict with the names of built-in gates

For example, we can define a 2-qubit gate and use it like this:

// define gate Rs
defgate Rs = [
    0.5+0.8660254j,0,0,0;
    0,1,0,0;
    0,0,1,0;
    0,0,0,1
];

qbit q[2];

unit main() {
    // apply Rs on q;
    Rs(q[0], q[1]);
}

If all the non-zero elements in the unitary is 1, i.e., the gate maps compuatational bases to a permutation of compuatational basis, you can use a simplified version with keyword perm. For example, the Toffoli gate can be defined as

defgate Tof(3) = perm [0, 1, 2, 3, 4, 5, 7, 6];

The syntax means that the 3-qubit gate Tof maps $|i\rangle$ to $|i\rangle$ for $0\leq i\leq5$, $|6\rangle$ to $|7\rangle$, and $|7\rangle$ to $|6\rangle$.

quantum state preparation

You can prepare a quantum state using an assignment. The left side of the assignment should be a qbit array with a fixed length. The right side can be two forms:

• A complex number array that represents the amplitude of the desired state. This array will be padded with 0 if its length is smaller than the Hilbert space dimension. After that, the amplitudes will be normalized.

• An expression consisting of kets such as (0.5+0.5j)|2>. Only computational basis states are supported currently.

import std;

procedure main()
{
    qbit q[2];

    // Initialize q in a Bell state using an array
    q = [1, 0, 0, -1];

    // or using ket expressions.
    q = |0> - |3>;
}

Note that the computational basis is coded in little endian. For example, qbit q[2]; q=|2>; results that q[1] is in $|1\rangle$ and q[0] in $|0\rangle$.

Procedure

The main body of isQ program is composed of procedures, which is similar to functions in C. The entry procedure is named main, which has no paramters and output.

Users can define their own procedures, and procedures may have no output or return with classical type. When a procedure has no output, it uses keyword procedure or unit at the beginning of definition, otherwise, it uses classical type. User-procedure can has parameters with any type mentioned above, and can be called by others or by itself.

// proceduer with two qbit as paramters and has no return
// paramtesr can also be written like (a: qbit, b: qbit)
unit swap(qbit a, qbit b) {
    ...
}

// proceduer with two classical paramters and return with a double
double compute(int a, double b[]) {
    double s = 0.0;
    ...
    return s;
}

// entry procedure
unit main() {
    qbit q[2];
    // call swap
    swap(q[0], q[1]);

    double b[2];
    // call compute
    double c = compute(1, b);
}

In addition, procedure itself can also serve as a parameter to another procedure. Such a parameter is specified as returnType identifier(TypeList), in which returnType and TypeList are the return type and the argument types of the procedure. For example,

// proceduer with two qubits and a two-qubit-gate procedure as paramters and has no return.
// parameters can also be written like (a: qbit, b: qbit, two_bit_gate: (qbit, qbit)->unit)
unit circuit(qbit a, qbit b, unit two_bit_gate(qbit, qbit)) {
    H(a);
    H(b);
    two_bit_gate(a, b);
}

// two qubit gate procedure1
unit two_bit_1(qbit a, qbit b){
    ...
}

// two qubit gate procedure2
unit two_bit_2(qbit a, qbit b){
    ...
}

// entry procedure
unit main() {
    qbit q[2];
    // call circuit with two_bit_1
    circuit(q[0], q[1], two_bit_1);
    // call circuit with two_bit_2
    circuit(q[0], q[1], two_bit_2);
}

deriving

User-procedure can be purely quantum. Such kind of procedures work as gates, and we can do some operations on these procedures like on ordinary quantum gates, like modifiers.

isQ provides a keyword deriving gate, converting a procedure to quantum gate, so that user can add modifiers when calling procedure. For example, if we need a controlled-swap, we can write codes like:

// pure quantum procedure that swap two qbit
unit swap(qbit a, qbit b) {
    CNOT(b, a);
    CNOT(a, b);
    CNOT(b, a);
} deriving gate

unit main() {
    qbit q[3];
    // set |q> -> |110>
    X(q[0]);
    X(q[1]);
    // controlled-swap, set |q> -> |101> 
    ctrl swap(q[0], q[1], q[2])
}

When a qbit procedure argument is of an arry type, its length must be fixed. For example, qbit a[3] is legal but qbit a[] is not.

template

A template represents a group of procedures. The syntax for defining and using a template is shown in the following example:

unit fun<int N>(){
    print N;
}

unit main(){
    fun::<2>();
    fun::<(3 + 1)>();
}

Note that expressions with operators must be wrapped with parenthesis ().

A normal procedure is generated based on the instantiation value before other compilation processes. Therefore, it may bypass some limitations of deriving gate, e.g.,

unit g<int N>(qbit a[N]){
    ...
} deriving gate

In addition, we can define some procedures for instantiating a template with particular values. This feature is helpful in defining the base case of recursive templates. For example, the factorial function can be implemented as follows.

int factorial<int N>() {
    return N * factorial::<(N - 1)>();
}

int factorial::<0>() {
    return 1;
}

When calling factorial::<0>(), the second procedure will be used instead of the first.

Control Flow

isQ provides four kinds of classical control flow:

• if-statement
• for-loop
• while-loop
• switch-case

if-else

If-else statement in isQ is similar to that in C. However, the braces {...} cannot be omitted, even if there is single statement following the condition.

unit main() {
    int a = 1;
    int b;
    if (a > 0) {
        b = 1;
    } else {
        b = -1;
    }

    bool c = true;
    if (c) {
        b = 2;
    }
}

for-loop

For-loop in isQ is similar to that in python. User can uses a non-defined variable as the iterating variable, and use two integer value a, b as loop range [a, b), where a must be less than b. The step size is optional and defaulted to 1. Like for-loop in other languages, keywords break and continue are supported.

unit main() {
    int a = 0;

    for i in 0:10 {
        a = a + 1;
        if (a > 5) {
            break;   
        }
    }

    for j in 0:10:3 {
        a = a - 1;
    }
}

isQ also supports array iteration. For example,

unit main() {
    int a[] = [2, 3, 4];
    for v in a {
        print v;
    }
}

Printed results should be "2", "3", and "4".

while-loop

While-loop in isQ is similar to that in C. However, the braces {...} cannot be omitted, even if there is single statement following the condition.

unit main() {
    bool a = true;
    qbit q;
    while (a) {
        H(q);
        a = M(q);
    }
}

switch-case

Based on the value of an int-type variable, a switch-case selects a correspoding branch. The execution terminates before another case or default keyword. For example,

unit main() {
    int a = 1, b;
    switch a {
    case 1:
        b = 3;
    case 2:
        b = 4;
    default:
        b = 5;
    }
}

That is, the switch-case statement does not have the fall through semantics in some other languages such as C. Therefore, the final value of b would be 3, not 5. The switch-case statement does not support break as it is not needed.

Modifiers

Controlled gates are common in quantum circuit. As we introduced before, users can define gate, so theoretically, any controlled gate can be defined. But it requires users to calculate the whole matrix. To simplify this work, isQ provide gate modifiers like some other quantum programming languages.

There are three gate modifiers in isQ:

• ctrl: modify gate to controlled-gate, can has a paramter representing the number of control qubits(defualt is 1).
• nctrl: modify gate to neg-controlled-gate, can has a paramter representing the number of control qubits(defualt is 1).
• inv: modify gate to conjugate transpose.

Controlled gates can be easily defined by ctrl and nctrl. And of course, you can compound these modifiers. For example:

unit main() {
    qbit p, q, r;

    // p ---●---
    //      |  
    // q ---S---
    ctrl S(p, q);

    // p ---●---
    //      |  
    // q ---●---
    //      |  
    // r ---S---
    ctrl<2> S(p, q, r);

    // q ---S^{-1}---
    inv S(q);

    // p --X--●--X---
    //        |  
    // q ---S^{-1}---
    nctrl inv S(p, q);
}

When a group of gates share the same control qubits, sometimes it is better to use a switch-case statement. For example

    qbit q[2], p, r;
    switch q {
    case |3>:
        ctrl Rx(pi, r, p);
    default:
        H(p);
    }

It is equivalent to the following statements

    nctrl nctrl H(q[1], q[0], p);
    nctrl ctrl H(q[1], q[0], p);
    ctrl nctrl H(q[1], q[0], p);
    ctrl ctrl Rx(pi, q[1], q[0], p);

Oracle

value table

isQ supports simple quantum oracle definitions. For those simple oracles whose function like $f: \{0,1\}^n \rightarrow \{0,1\}^m$, keyword oracle can be used like this:

oracle g(2, 1) = [0, 1, 0, 0];

qbit q[3];

unit main() {
    ...
    g(q[2], q[1], q[0]);
    ...
}

In the above example, g is the oracle name, and there are two integer parameters n, m, where n represents the number of work qubits and m represents the number of ancilla qubits. The value list is ${f(i)} {\ } for {\ } i {\ } in [0, 2^n) $, and every $f(i)$ must be in $[0, 2^m)$. Usually, $m$ is 1, and the list is a truth table like above.

oracle function

Alternatively, you can define an oracle using a Boolean function. This function accepts one or more parameters with bool array type, and returns a bool array. For example, we can rewrite the oracle g as follows:

oracle bool[1] g(bool x[2]) {
    bool res[] = [x[0] && !x[1]];
    return res; 
}

The derived quantum gate applies to qbit arrays. For example:

unit main() {
    qbit p[2], q[1];
    ...
    g(p, q);
    ...
}

Parameter

isQ supports compiling with parameters and passing values at runtime. If you want to do this, you can define two arrays in the parameter list of the procedure main, and you can use them in the function body, like this:

unit main(int i_par[], double d_par[]) {
    ...
    Rx(d_par[0], q);
    if (i_par[1] == 2){...}
    ...
}

When simulating, you can pass values by -i or -d. isQ will gather values and generate int and double arrays (can be empty). For the above example, you can pass values like this:

isqc simulate -i 1 -i 2 -d 1.3 xxx.isq

Basic Gates

isQ supports the following quantum gates.

X

X(qbit q)

Pauli X gate. $$X=\begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix}$$

Y

Y(qbit q)

Pauli Y gate. $$Y=\begin{bmatrix} 0 & -i\\ i & 0 \end{bmatrix}$$

Z

Z(qbit q)

Pauli Z gate. $$Z=\begin{bmatrix} 1 & 0\\ 0 & -1 \end{bmatrix}$$

H

H(qbit q)

Hardmard gate. $$H=\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & 1\\ 1 & -1 \end{bmatrix}$$

S

S(qbit q)

Phase gate. $$S=\begin{bmatrix} 1 & 0\\ 0 & i \end{bmatrix}$$ Note that $S^2=Z$.

T

T(qbit q)

(\pi/8) gate. $$T=\begin{bmatrix} 1 & 0\\ 0 & e^{i\pi/4} \end{bmatrix}$$ Note that $T^2=S$.

Rx

Rx(double theta, qbit q)

Rotation about the X axis. $$Rx(\theta)=e^{-i\theta X/2}=\begin{bmatrix} \textrm{cos}\frac{\theta}{2} & -i\textrm{sin}\frac{\theta}{2}\\ -i\textrm{sin}\frac{\theta}{2} & \textrm{cos}\frac{\theta}{2} \end{bmatrix}$$

Ry

Ry(double theta, qbit q)

Rotation about the Y axis. $$Ry(\theta)=e^{-i\theta Y/2}=\begin{bmatrix} \textrm{cos}\frac{\theta}{2} & -\textrm{sin}\frac{\theta}{2}\\ \textrm{sin}\frac{\theta}{2} & \textrm{cos}\frac{\theta}{2} \end{bmatrix}$$

Rz

Rz(double theta, qbit q)

Rotation about the Z axis. $$Rz(\theta)=e^{-i\theta Z/2}=\begin{bmatrix} e^{-i\theta/2} & 0\\ 0 & e^{i\theta/2} \end{bmatrix}$$

X2P

X2P(qbit q)

$\pi/2$ rotation about the X axis. $$X2P=Rx(\pi/2)=\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & -i\\ -i & 1 \end{bmatrix}$$

X2M

X2M(qbit q)

$-\pi/2$ rotation about the X axis. $$X2M=Rx(-\pi/2)=\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & i\\ i & 1 \end{bmatrix}$$

Y2P

Y2P(qbit q)

$\pi/2$ rotation about the Y axis. $$Y2P=Ry(\pi/2)=\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & -1\\ 1 & 1 \end{bmatrix}$$

Y2M

Y2M(qbit q)

$-\pi/2$ rotation about the Y axis. $$Y2M=Ry(-\pi/2)=\frac{1}{\sqrt{2}}\begin{bmatrix} 1 & 1\\ -1 & 1 \end{bmatrix}$$

U3

U3(double theta, double phi, double lambda, qbit q)

Generic single-qubit rotation gate with 3 Euler angles. $$U3(\theta,\phi,\lambda)=\begin{bmatrix} e^{-i(\phi+\lambda)/2}\textrm{cos}\frac{\theta}{2} & -e^{i(\lambda-\phi)/2}\textrm{sin}\frac{\theta}{2}\\ e^{i(\phi-\lambda)/2}\textrm{sin}\frac{\theta}{2} & e^{i(\phi+\lambda)/2}\textrm{cos}\frac{\theta}{2} \end{bmatrix}$$ Note that $$U3\left(\theta,-\frac{\pi}{2},\frac{\pi}{2}\right)=Rx(\theta)$$ $$U3(\theta,0,0)=Ry(\theta)$$

CNOT

CNOT(qbit c, qbit t)

Controlled-NOT gate. $$CNOT=\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 \end{bmatrix}$$

CZ

CZ(qbit c, qbit t)

Controlled-Z gate. $$CZ=\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & -1 \end{bmatrix}$$

Toffoli

Toffoli(qbit c1, qbit c2, qbit t)

Toffoli gate. $$Toffoli=\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \end{bmatrix}$$

GPhase

GPhase(double theta)

Global phase gate. $$GPhase(\theta)=e^{i\theta}$$ Note that this gate does not have any qbit-type parameter.

Ended: User Manual

Developer Manual

Development Setup

To build up development environment for isQ itself, you need to install Nix with Flakes support on a GNU/Linux distro first.

Install Nix Package Manager

We recommend using Nix installer by Determinate Systems to install Nix with Flakes support out of the box.

# Install Nix.
curl --proto '=https' --tlsv1.2 -sSf -L https://install.determinate.systems/nix | sh -s -- install
# Setup Cachix binary cache.
nix-shell -p cachix --run "cachix use arclight-quantum"

Enter Development Environment

We provide necessary tools to build isQ components locally.

git clone git@github.com:isQ-Team/isQ-Compiler.git && cd isQ-Compiler
# Virtual environment with necessary building tools. The bash prompt should become (nix-isqc:dev).
nix develop
# Virtual environment with a VSCodium. The bash prompt should become (nix-isqc:codium).
nix develop .#codium

Git Commit

We use pre-commit git hooks to: - prevent accidentally forgetting to update flake.lock; - format all Nix build scripts.

If a pre-commit git hook prevents you from committing, just re-add the files changed by git hooks and commit again.

Versioning

isQ follows Semantics Versioning 2.0.0 for versioning. Specifically, the version is given in format MAJOR.MINOR.PATCH or MAJOR.MINOR.PATCH+METADATA.

The main version of a build MAJOR.MINOR.PATCH is defined as-is:

• MAJOR version when you make incompatible API changes.
• MINOR version when you add functionality in a backward compatible manner.
• PATCH version when you make backward compatible bug fixes.

The metadata part might be empty, +<REV>, +<DIRTY>, +<REV>.<DIRTY> containing the git tree status of a given version:

• The <REV> part contains the short version of git commit hash, if the commit is not tagged as a released version.

• The <DIRTY> part is dirty if the current git tree is dirty, or empty if the current git tree is clean.

A few examples:

• Version 0.1.0 is a clean build of released version.
• Version 0.1.0+dirty is a dirty build of tagged version 0.1.0, by for example building after modifying some code.
• Version 0.1.0+8330d80 is a clean build of git commit 8330d80; the closest tag to this commit is 0.1.0.
• Version 0.1.0+8330d80.dirty is a dirty build of git commit 8330d80; the closest tag to this commit is 0.1.0.

Usage of Versioning

Single-sourced version

The single source of version is in version.json. The config now have two tags:

• version indicating the main version of the commit.
• frozen indicating whether this commit is a frozen released version. This flag should never be set by developer by hand.

Subproject metadata versioning

The main version of a build (e.g. 0.1.0) is propagated to the version of package metadata of different projects, e.g. Cargo.toml for Rust subprojects and package.yaml for Haskell subprojects.

This versioning scheme makes sure we only need to update the main version (for all projects) when we are ready for a release.

Displayed version

The full version should be used in the project, e.g. when invoking isqc --version.

Due to the limitation of Nix, currently the revision of a dirty build cannot be passed to the Nix-built binary (See Issue #4682 of Nix) For example, the version 0.1.0+8330d80.dirty can only show as 0.1.0+dirty.

Versioning Cycle

• Development of major and minor versions focuses on main branch. Patches are ported to release/MAJOR.MINOR.x branches.
• Any version bumping must go through a pull request. The commit can be as simple as updating version.json, since it will be rebased onto the source branch anyway.
• When the pull request is merged, the commits will be rebased onto main.
• When the major or minor version is bumped on main branch (patch version is 0), a new branch release/MAJOR.MINOR.x will be created.
• Whenever a new patch version is bumped on release/MAJOR.MINOR.x, a new commit setting the flag frozen to true will be created and pushed to corresponding git tag MAJOR.MINOR.PATCH.

An example of the main branch would be as follows: the X1 commit is added to main branch by a pull request bumping the version to 0.2.0, as well as branched as release/0.2.x. A new commit Y1 based on X1 setting frozen to true will be created and pushed as git tag 0.2.0. Suppose there is unfortunately a bug on 0.2.0, commits f through X3 fix the bug and bump the patch version to 0.2.1.

   PR: Bump 0.2.0  PR: Bump 0.3.0
 - a - b - X1 - c - d - X2 - e (main)
           \            \ 
           |-Y1(0.2.0)  |-Y2(0.3.0)
           |             \ i - j - k - l (release/0.3.x)
            \         PR: Bump 0.2.1
              - f - g - X3 - h (release/0.2.x)
                        \
                         Y3(0.2.1)

Roadmap

Currently isQ is still in early stage of development (0.x.y) and everything may evolve. Informally:

• MAJOR version is 0.
• MINOR version is bumped when some important feature is added to isQ.
• PATCH version is bumped when some patches or fixes are added to isQ.

Ended: Developer Manual

Examples

Basic examples

We will show you how to write isQ code with some basic examples. For readers interested in the algorithm details, please refer to the textbook Quantum Computation and quantum information written by M.A. Nilsen and I. Chuang. All the examples can be found in the examples folder in isQ open-source code repository.

• Bell state

• Quantum teleportation

• Deutch-Jozsa algorithm

• Bernstein-Vazirani algorithm

Bell state

The Bell's states or EPR pairs are specific quantum states of two qubits that represent the simplest (and maximal) examples of quantum entanglement. They can be created using the following circuit:

It can be easily expressed in isQ:

import std;

qbit a, b;

procedure main()
{
    H(a);
    CNOT(a, b);
    print M(a);
    print M(b);
}

The measured result should both be 0, or both be 1.

Quantum teleportation

Quantum teleportation is the transfer of a quantum state over a distance. Alice and Bob share an EPR pair and each took one qubit before they became separated. Alice must deliver a qubit of information to Bob, but she does not know the state of this qubit and can only send classical information to Bob.

It is performed step by step as the following:

1. Alice sends her qubits through a CNOT gate.
2. Alice then sends the first qubit through a Hadamard gate.
3. Alice measures her qubits, obtaining one of four results, and sends this information to Bob.
4. Given Alice's measurements, Bob performs one of four operations on his half of the EPR pair and recovers the original quantum state. The following quantum circuit describes teleportation:

The corresponding isQ program is as follows:

import std;
procedure transform(qbit a, qbit b, qbit c) {
    // Prepare the EPR pair
    H(b);
    CNOT(b, c);

    CNOT(a, b);
    H(a);

    // Fix up the state based on measurement result
    if (M(b)) { X(c); }
    if (M(a)) { Z(c); }
}
procedure main()
{
    qbit q[3];
    Rx(pi/3, q[0]); // Set initial state as 'sqrt(0.75)|0>-sqrt(0.25)i|1>'
    transform(q[0], q[1], q[2]);
    print M(q[2]);
}

The measurement result will be consistent with the probability distribution of measuring $\sqrt{0.75}|0\rangle-\sqrt{0.25}i|1\rangle$.

Deutsch-Jozsa algorithm

The Deutsch–Jozsa algorithm is a deterministic quantum algorithm proposed by David Deutsch and Richard Jozsa in 1992. Although of little current practical use, it is one of the first examples of a quantum algorithm that is exponentially faster than any possible deterministic classical algorithm.

We are given a black box quantum computer known as an oracle that implements some function: $$f:\{0,1\}^n\rightarrow\{0,1\}$$ The function takes n-bit binary values as input and produces either a 0 or a 1 as output for each such value. We are promised that the function is either constant (0 on all inputs or 1 on all inputs) or balanced (1 for exactly half of the input domain and 0 for the other half). The task then is to determine if $f$ is constant or balanced by using the oracle.

For a conventional deterministic algorithm, $2^{n-1}+1$ evaluations of $f$ will be required in the worst case. However, the quantum solution requires only one evaluation using the following circuit:

Its corresponding isQ program is as follows:

import std;

oracle bool[1] constant(bool a[4]) {
    bool res[] = [true];
    return res;
}

procedure main()
{
    qbit q[4], anc[1];
    X(anc[0]);
    H(q);
    H(anc[0]);
    constant(q, anc);
    H(q);
    print M(q);
}

Note that we defined a constant oracle that returns 1 for all inputs. Therefore, the printed result should be 0.

Bernstein-Vazirani algorithm

The Bernstein–Vazirani algorithm is a restricted version of the Deutsch–Jozsa algorithm where instead of distinguishing between two different classes of functions, it tries to learn a string encoded in a function. The oracle is defined as the function: $$f(x)=x\cdot s=x_1s_1\oplus x_2s_2\oplus\cdots\oplus x_ns_n$$ where $s\in\{0,1\}^n$ is a secret string and $\oplus$ represents addition modulo 2.

Following the same circuit as the Deutsch–Jozsa algorithm, the measurement result of the first register is a secret string.

We write an isQ program to demonstrate the Bernstein–Vazirani algorithm. We select $s=1011$ as the secret bitstring and use it to define the oracle g. Note that in its definition, $s$ is defined as a bool array, and the addition modulo 2 operation is implemented by the != operator.

import std;

oracle bool[1] g(bool a[4]) {
    bool s[] = [true, true, false, true];
    bool ba[] = s & a;  // bit-wise AND
    bool res[1];
    res[0] = ba[0] != ba[1] != ba[2] != ba[3];
    return res;
}

procedure main()
{
    qbit q[4], res[1];
    X(res);
    H(res);
    H(q);
    g(q, res);
    H(q); 
    print M(q); // should print 11 (=0x1011)
}

After execution, the program should print 11, i.e., the secret string $s$.

Quantum arithmetic operations

We demonstrate building some quantum arithmetic operations using basic gates. Some operations would be used in other quantum algorithms, such as Shor's algorithm.

• Adder

• Subtractor

• Carrier

• Multiplier

• Modular adder

• Modular multiplier

• Modular exponentiation

Adder

We will present two quantum adders based on the ripple-carry approach, i.e., starting from low-order bits of the input and working our way up to the high-order bits. We assume the inputs are $a=a_{n−1}\cdots a_0$ and $b=b_{n−1}\cdots b_0$, both with $n$ qubits and 0 is the lowest-order bit. The adder transfers $$|a\rangle|b\rangle\rightarrow|a\rangle|a+b\rangle$$

First, we show a basic design. It relies on two major components: a majority (MAJ) gate and an UnMajority and Add (UMA) gate. They make a 1-bit full-adder together.

These two components can be built with CNOT and Toffoli gates easily:

It follows that we can string together MAJ and UMA gates to build a ripple-carry adder. Below is a 3-bit adder

Note that the first qubit is ancillary. The isQ code implementing this adder is

/*
* Majority gate
*
* |a_i>|b_i>|c_i> -> |c_{i+1}>|b_i^a_i>|c_i^a_i>.
* when the inputs contain two or three 1s, c_{i+1} = 1; otherwise, 0.
*/
procedure maj(qbit a, qbit b, qbit c) {
    CNOT(c, a);
    CNOT(c, b);
    Toffoli(a, b, c);
} deriving gate

/*
* UMA (UnMaj-and-Add) gate
*
* |c_i^a_i>|b_i^a_i>|c_{i+1}> -> |c_i>|s_i>|a_i>
* s_i is the addition result of a_i, b_i, and c_i.
*/
procedure uma(qbit ca, qbit ba, qbit c)
{
    Toffoli(ca, ba, c);
    CNOT(c, ca);
    CNOT(ca, ba);
} deriving gate

/*
* 3-bit adder
* |a>|b> -> |a+b>|a>
*/
procedure adder3(qbit a0, qbit a1, qbit a2, qbit b0, qbit b1, qbit b2) {
    qbit anc;
    maj(anc, a0, b0);
    maj(b0, a1, b1);
    maj(b1, a2, b2);

    uma(b1, a2, b2);
    uma(b0, a1, b1);
    uma(anc, a0, b0);
} deriving gate

The keyword deriving gate transfers a procedure into a gate so that it can be decorated with ctrl and inv.

Above adder design can be optimized to improve scalability. One of the more efficient designs is shown below

Note that this design does not require any ancillary qubit. We write an isQ procedure for this design

procedure __add(qbit a[], qbit b[])
{
    int len = a.length;
    CNOT(a[1:len], b[1:]);
    CNOT(a[len-2:0:-1], a[len-1:0:-1]);
    Toffoli(a[0:], b, a[1:len]);
    for i in len-1:0:-1 {
        CNOT(a[i], b[i]);
        Toffoli(a[i - 1], b[i - 1], a[i]);
    }
    CNOT(a[1:], a[2:len]);
    CNOT(a, b);
}

We use two qbit arrays as parameters. This procedure works with an array of any length. Many gates used in this procedure are bundle operations, i.e., operations applied to array elements. For example, CNOT(a[1:len], b[1:]) means applying CNOT gates pair-wise to (a[1], b[1]), (a[2], b[2]), ..., (a[len-1], b[len-1]). This feature of isQ makes programs concise. However, this feature and classical control flow conflict with deriving gate because only pure quantum circuits can be decorated.

The second adder design has been adopted as the default adder in the current isQ compiler. It will be called with the form a += b; where a and b are both qubit arrays. For example, the following program will print 3, 6, and 1 in sequence.

import std;
procedure main()
{
    qbit p[3];
    qbit q[3];
    X(p[0:2]); // = X(p[0]); X(p[1]);
    q += p;
    print M(q); // should print 3
    q += p;
    print M(q); // should print 6
    q += p;
    print M(q); // should print 1
}

By setting an addend to a constant number, we can turn an adder into a constant adder. More formally, a constant adder maps $$|a\rangle\rightarrow|a+c\rangle$$ where $c$ is a constant number.

We can implement the constant adder using a regular adder and a simple component SET ($S_c$). The subscript $c$ represents the constant number. The set gate is simply some $X$ gates corresponding to the bit 1s in $c$. Therefore, its reverse is itself.

For example, an increment gate can be seen as a constant adder where $c=1$. We implement it in isQ:

procedure add1(qbit a0, qbit a1, qbit a2) {
    qbit b0, b1, b2;
    X(b0);
    adder3(a0, a1, a2, b0, b1, b2);
    X(b0);
} deriving gate

Note that $S_1$ is implemented with X(b0).

Subtractor

Subtraction is the dual operation of addition. Therefore, a subtractor is just the reverse of an adder. The following example demonstrates the way of subtracting using an adder.

import std;
import adder3;
procedure main()
{
    qbit a[3], b[3];

    // Set |a> = |7>
    X(a);

    // Set |b> = |1>
    X(b[0]);

    // |a>|b> -> |a-b>|b>
    inv adder3(a[0], a[1], a[2], b[0], b[1], b[2]);

    print M(a);     // should be 6
    print M(b);     // should be 1
}

Similarly to addition, isQ supports using -= as the short form of calling the subtraction procedure.

Carrier

A special gate named Carrier is helpful for gates such as the modular adder. It computes $$|a\rangle|b\rangle|0\rangle\rightarrow\left\{\begin{array}{l} |a\rangle|b\rangle|0\rangle\text{ if }a+b<2^n,\\ |a\rangle|b\rangle|1\rangle\text{ if }a+b\geq2^n, \end{array}\right.$$ i.e., it only computes whether $a+b$ leads to a carry without changing the values of $a$ and $b$. Therefore, we can use a circuit similar to the adder to implement the carrier gate.

Note that the gate on the right side of the figure is the inverse of $MAJ$. Thus, $a$ and $b$ change back to their initial values. The corresponding isQ code is

procedure carrier3(qbit a0, qbit a1, qbit a2, qbit b0, qbit b1, qbit b2, qbit c) {
    qbit anc;
    maj(anc, a0, b0);
    maj(b0, a1, b1);
    maj(b1, a2, b2);

    // Copy out the result
    CNOT(b2, c);

    // Uncompute
    inv maj(b1, a2, b2);
    inv maj(b0, a1, b1);
    inv maj(anc, a0, b0);
} deriving gate

Multiplier

We use isQ to build a naive multiplier in this section. It implements the following transformation:

$$|a\rangle|b\rangle|c\rangle\rightarrow|a\rangle|b\rangle|c+a\times b\rangle$$

where the lengths of $a$, $b$, and $c$ are $n$, $n$, and $2n$, respectively. The algorithm mimics the way of multiplying binary numbers on paper, e.g.,

We prepare an ancillary qubit array anc as the scratchpad, i.e., the memory storing the shifted numbers shown between the two lines in the above figure. The shift operation is implemented using Toffoli gates, which can also be seen as the multiplication of two bits. Next, anc is added to c using an adder. Finally, anc is uncomputed to $|0\rangle$, so the computation continues. The complete program is shown below.

// |a>|b>|c> -> |a>|b>|c+a*b>
procedure mul(qbit a[], qbit b[], qbit c[])
{
    // Check the length of qubit arrays
    int len = a.length;
    assert(len == b.length);
    assert(len * 2 == c.length);

    qbit anc[len * 2];
    for i in 0:len {
        // anc = b << i
        for j in 0:len {
            Toffoli(a[i], b[j], anc[i + j]);
        }

        c += anc;

        // anc = 0
        for j in 0:len {
            Toffoli(a[i], b[j], anc[i + j]);
        }
    }
}

Modular adder

Let us consider a type of modular adder that computes $$|a\rangle\rightarrow|a+c\bmod N\rangle$$ where $c$ and $N$ are constant numbers.

One way to implement a modular adder is first to test whether the result is larger than or equal to $N$. If $a+c\geq N$, then $a+(2^n+c-N)\geq2^n$. Therefore, this step can be done by setting ancillary qubits to $2^n+c-N$ and then applying the carrier gate. If $a+c\geq N$, we need to add $c$ to $a$ and subtract $N$. It equals adding $2^n+c-N$. If $a+c<N$, then we can add $c$ directly. Finally, we convert the carry bit back to $|0\rangle$ by calling the carrier gate again. The overall circuit is shown as follows.

As all the components require at most $O(n)$ quantum gates, the cost of a modular adder is $O(n)$.

We implement an example where $c=1$ and $N=15$ in isQ:

// |a> -> |a+1(mod 15)>
procedure add1mod15(qbit a0, qbit a1, qbit a2, qbit a3) {
    qbit b0, b1, b2, b3, c;

    // set b = |2>
    X(b1);

    // Test whether a + 1 >= 15
    carrier4(a0, a1, a2, a3, b0, b1, b2, b3, c);

    // If true, add 1 then subtract 15, i.e., subtract 14, i.e., add 2
    ctrl add2(c, a0, a1, a2, a3);

    // If false, add 1 directly
    X(c);
    ctrl add1(c, a0, a1, a2, a3);

    // Uncompute
    X(b0);
    X(b2);
    X(b3);
    carrier4(a0, a1, a2, a3, b0, b1, b2, b3, c);
    X(b3);
    X(b2);
    X(b1);
    X(b0);
} deriving gate

Modular multiplier

Based on the modular adder, we can build a modular multiplier. That is, given an $n$-qubit input $|a\rangle$, we calculate the result $|a\times c(\bmod N)\rangle$ where $c$ and $N$ are constant numbers.

If we express $a$ in binary form $a_{n-1}\cdots a_0$, i.e., $a=\sum_{j=0}^{n-1}a_j2^j$, we can see that $$a\times c=\sum_{j=0}^{n-1}a_j2^jc(\bmod N)$$ i.e., modular multiplication is a series of modular addition chained together.

This observation motivates us to implement the modular multiplier using the following circuit.

As a modular multiplier contains $n$ modular adders, the overall cost is $O(n^2)$.

The following is an isQ implementation where $c=7$ and $N=15$. Note that $7\times2^1=14(\bmod15)$, $7\times2^2=13(\bmod15)$, and $7\times2^3=11(\bmod15)$.

// |a>|0> -> |a>|a*7(mod 15)>
procedure mul7mod15(qbit a3, qbit a2, qbit a1, qbit a0, qbit b3, qbit b2, qbit b1, qbit b0) {
    ctrl add7mod15(a0, b0, b1, b2, b3);
    ctrl add14mod15(a1, b0, b1, b2, b3);
    ctrl add13mod15(a2, b0, b1, b2, b3);
    ctrl add11mod15(a3, b0, b1, b2, b3);
} deriving gate

The above implementation keeps input $a$ unchanged. That is, gate $MUL_c$ maps $$|a\rangle|0\rangle\rightarrow|a\rangle|a\times c(\bmod N)\rangle$$ However, a unitary gate $|a\rangle\rightarrow|a\times c(\bmod N)\rangle$ is preferred sometimes. It can be built with the following circuit

The box in the middle of the figure represents a $SWAP$ gate, i.e., $|a\rangle|b\rangle\rightarrow|b\rangle|a\rangle$ Note that $c\times c^{-1}=1(\bmod N)$, and $c^{-1}$ does not always exist for any $c$. Also, it does not always possible to build the unitary form since the modular multiplication may not be reversible. This circuit is implemented with the following isQ program

// |a>|b> -> |b>|a>
procedure swap4(qbit a0, qbit a1, qbit a2, qbit a3, qbit b0, qbit b1, qbit b2, qbit b3) {
    SWAP(a0, b0);
    SWAP(a1, b1);
    SWAP(a2, b2);
    SWAP(a3, b3);
} deriving gate

// |a> -> |a*7(mod 15)>
procedure mul7mod15_u(qbit a3, qbit a2, qbit a1, qbit a0) {
    qbit anc[4];
    mul7mod15(a3, a2, a1, a0, anc[3], anc[2], anc[1], anc[0]);
    swap4(a3, a2, a1, a0, anc[3], anc[2], anc[1], anc[0]);
    inv mul13mod15(a3, a2, a1, a0, anc[3], anc[2], anc[1], anc[0]);
} deriving gate

Note that $7^{-1}=13(\bmod15)$, i.e., $7\times13=1(\bmod15)$.

Modular exponentiation

In a similar way, we can also build modular exponentiation ($|p\rangle|a\rangle\rightarrow|p\rangle|x^pa(\bmod N)\rangle$) based on modular multiplication. Assume $p=\sum_{j=0}^{n-1}p_j2^j$. Then, $x^p=\prod_{j=0}^{n-1}x^{p_j2^j}$. The corresponding circuit is shown as follows.

As the above circuit contains $n$ modular multipliers, its cost is $O(n^3)$.

We implement such a gate where $x=7$ and $N=15$ with isQ:

// |p>|a> -> |a*7^p(mod 15)>
procedure pow7mod15(qbit p[], qbit a[4]) {
    // 7^(2^0) = 7 (mod 15)
    ctrl mul7mod15_u(p[0], a[3], a[2], a[1], a[0]);

    // 7^(2^1) = 4 (mod 15)
    ctrl mul4mod15_u(p[1], a[3], a[2], a[1], a[0]);
}

For bits other than the lowest two, the modular result is always 1, e.g., $7^{2^2}\bmod15=1$. Therefore, the corresponding modular multiplication gate would be a trivial identical gate.

Reference

1. Steven A. Cuccaro, Thomas G. Draper, et al. "A new quantum ripple-carry addition circuit." arXiv preprint quant-ph/0410184, 2004.

2. Y. Takahashi, S. Tani, and N. Kunihiro, "Quantum Addition Circuits and Unbounded Fan-Out." arXiv preprint 0910.2530, 2009.

3. D. Beckman, A. Chari, et al. "Efficient networks for quantum factoring." Physical Review A, 1996, 54: 1034-1063.

Hidden subgroup problems

Quantum computing can efficiently solve hidden subgroup problems (HSPs) for finite Abelian groups. Given a group $G$, a subgroup $H\leq G$, and a set $X$, we say a function $f:G\rightarrow X$ hides the subgroup $H$ if $f$ is constant on the cosets of $H$, while it is different between the different cosets of $H$. HSP is to determine a generating set of $H$ using $f$ as an oracle.

Simon's algorithm, quantum Fourier transformation, quantum phase estimation are all solutions to specific HSPs. We will demonstrate these algorithms using isQ examples.

• Simon's algorithm

• Quantum Fourier transformation

• Quantum phase estimation

• Shor's algorithm

Simon's algorithm

Given a function $f:\{0,1\}^n\rightarrow\{0,1\}^n$ with the promise that, for some unknown $s\in\{0,1\}^n$, for all $x,y\in\{0,1\}^n$, $$f(x)=f(y)\text{ if and only if }x\oplus y\in\{0^n,s\}$$ where $\oplus$ denotes bitwise XOR. The goal is to identify $s$ by making as few queries to $f(x)$ as possible. Note that $$a\oplus b=0^n\text{ if and only if }a=b$$ Furthermore, for some $x$ and $s$ in $x\oplus y=s$, $y$ is unique ($\neq x$) if and only if $s\neq0^n$. This means that $f$ is two-to-one when $s\neq0^n$, and one-to-one when $s=0^n$. It is also the case that $x\oplus y=s$ implies $y=s\oplus x$, meaning that $$f(x)=f(y)=f(x\oplus s)$$ which shows how $f$ is periodic.

The quantum circuit used in Simon's algorithm is:

The measurement result of the first register $y$ satisfies $$y\cdot s=0$$ Therefore, by executing the circuit $O(n)$ times, we can get a list of independent bitstrings $y_1,\cdots,y_{n-1}$. With the help of a classical computer, we can calculate $s$.

Next, we compose an example of Simon's algorithm in isQ. We choose $n=4$ and thus $0\leq x<16$. Let $f(x)=x$ (mod 8). Therefore, if $f(x) = f(y)$ and $x\neq y$, $x$ differs from $y$ only by the highest bit, i.e., $x=y\oplus8$ and $s=8$ (1000 as a bitstring). The complete code is as follows:

import std;

// g(x) = g(x + 8)
oracle bool[4] g(bool a[4]) {
    bool res[] = [a[0], a[1], a[2], false];
    return res;
}

procedure main(){
    qbit q[4], res[4];
    H(q);
    g(q, res);
    H(q); 
    print M(q);
}

We executed the program several times, and the first three outputs are 1, 3, and 6. It is easy to check that they are independent bitstrings. The only 4-bit bitstring $s$ that satisfies $s\oplus1=0$, $s\oplus3=0$, and $s\oplus6=0$ is 8, precisely the secret bitstring we defined.

Quantum Fourier transformation

Quantum Fourier transform (QFT) is the quantum analogue of the discrete Fourier transform. It is a key component of many HSP solutions, such as Shor's algorithm.

QFT is a linear operator with the following action on the basis states: $$|j\rangle\mapsto\frac{1}{N}\sum_{k=0}^{N-1}e^{2\pi ijk/N}|k\rangle$$ where $N=2^n$ and $n$ is the number of qubits. Equivalently, QFT can be given by the following product representation: $$|j_1,\cdots,j_n\rangle\mapsto\frac{\left(|0\rangle+e^{2\pi i0.j_n}\right)\left(|0\rangle+e^{2\pi i0.j_{n-1}j_n}\right)\cdots\left(|0\rangle+e^{2\pi i0.j_1j_2\cdots j_n}\right)}{2^{n/2}}$$ where $j_1,\cdots,j_n$ is the binary representation of $j$, and $0.j_lj_{l+1}\cdots j_m$ represents the binary fraction $j_l/2+j_{l+1}/4+\cdots+j_m/2^{m-l+1}$.

Inspired by the product representation, we use the following quantum circuit to implement QFT: where $$R_k\equiv\begin{bmatrix} 1 & 0\\ 0 & e^{2\pi i/2^k} \end{bmatrix}$$ Note that extra SWAP gates are needed at the end of the circuit to reverse the order of the qubits.

The isQ program describing the above circuit is:

import std;

procedure R(int k, qbit q) {
    double phase = pi / 2 ** (k - 1);
    ctrl GPhase(phase, q);
} deriving gate

procedure qft(qbit q[]) {
    int len = q.length;
    for i in len-1:-1:-1 {
        H(q[i]);
        for j in 0:i {
            ctrl R(i-j+1, q[j], q[i]);
        }
    }
    for i in 0:len/2 {
        SWAP(q[i], q[len-i-1]);
    }
}

// The inverse of qft
procedure qft_inv(qbit q[]) {
    int len = q.length;
    for i in 0:len/2 {
        SWAP(q[i], q[len-i-1]);
    }
    for i in 0:len {
        for j in 0:i {
            ctrl inv R(i-j+1, q[j], q[i]);
        }
        H(q[i]);
    }
}

Note that ctrl is a keyword that turns a gate into its controlled form, and its first qbit argument is the control qubit. deriving gate is a keyword that turns a procedure into a quantum gate. Measurement, branch, and loops cannot be used in such a procedure. .length is another keyword that acquires the length of an array. The inverse of QFT, named qft_inv, is also defined here. It will be used in the following examples.

Quantum phase estimation

The quantum phase estimation (QPE) algorithm is a quantum algorithm to estimate the phase (or eigenvalue) of an eigenvector of a unitary operator. More precisely, given a unitary matrix $U$ with an eigenvector $|u\rangle$, such that $$U|u\rangle=e^{2\pi i\theta}|u\rangle$$ The algorithm estimates the value of $\theta$ with high probability.

The overall structure of the QPE algorithm is as follows:

First, $H^{\otimes n}$ gates are applied to $|0\rangle^{\otimes n}$ to create the superposition of all numbers in $[0,2^n)$. Then, the superposition controls the gate $U^j$ so that $$\frac{1}{2^n}\sum_{j=0}^{2^n-1}|j\rangle|u\rangle\mapsto\frac{1}{2^n}\sum_{j=0}^{2^n-1}e^{2\pi ij\theta}|j\rangle|u\rangle$$

After that, applying the inverse of QFT ($FT^{-1}$) to the first register results in a state similar to $|2^n\theta\rangle$. For simplicity, we will not analyze the error.

The detailed circuit is shown as follows:

We choose $U=Rz(\theta)$ whose eigenvector is $|0\rangle$ and $\theta=2\pi\times23/2^6\approx2.258$. The corresponding isQ program is as follows:

import std;
import qft; // include qft_inv

int x = 23;
int n = 6;

double theta() {
    return 2 * pi * x / (2 ** n);
}

procedure U(double theta, qbit q) {
    X(q);
    ctrl GPhase(theta, q);
    X(q);
} deriving gate

procedure pow2_ctrlU(int n, qbit anc, qbit ev) {
    double t = theta() * (2 ** n);
    ctrl U(t, anc, ev);
}

int phase_estimation(int n, qbit ev) {
    qbit anc[n];
    for i in 0:n {
        H(anc[i]);
        pow2_ctrlU(i, anc[i], ev);
    }
    qft_inv(anc);
    return M(anc);
}

procedure main()
{
    qbit ev; /* ev is |0>, the eigen vector of U */
    print phase_estimation(n, ev);
}

Note that $FT^{-1}$ (qft_inv) is defined in qft.isq and not shown here. After execution, the printed result is 23, the phase we set.

Shor's algorithm

Shor’s algorithm is one of the most well-known quantum algorithms and finds the prime factors for input integer $N$ in polynomial time. This algorithm could break public-key schemes since they are based on the assumptions that factoring large integers is computationally intractable.

To factor a composite number $N$, Shor's algorithm is summarized below:

1. If $N$ is even, return the factor 2.

2. Determine whether $N=a^b$ for integers $a\geq1$ and $b\geq2$, and if so return the factor $a$.

3. Randomly choose $x$ in the range 1 to $N-1$. If gcd$(x,N)>1$ then return the factor gcd$(x,N)$.

4. Use the order-finding subroutine to find the order $r$ of $x$ modulo $N$. The \textit{order} $r$ is the minimal positive integer so that $x^r(\bmod N)=1$.

5. If $r$ is even and $x^{r/2}\neq-1(\bmod N)$ then compute gcd$(x^{r/2}-1, N)$ and gcd$(x^{r/2}+1, N)$, and test to see if one of these is a non-trivial factor, returning that factor if so. Otherwise, the algorithm fails.

Note that, in these steps, only Step 4 is conducted on a quantum machine. We will also focus on this part. The quantum circuit used for finding the order is

which is quite similar to the one used in QPE. The key reason that leads to this similarity is that the state defined by $$\left| {{u}_{s}} \right\rangle\equiv\frac{1}{\sqrt{r}}\sum\limits_{k=0}^{r-1}{\exp \left( \frac{-2\pi isk}{r} \right)}\left| {{x}^{k}}\bmod N \right\rangle$$ is an eigenvector of gate $x^j\bmod N$ with the eigenvalue $\exp(2\pi is/r)$. Therefore, we can infer $r$ by estimating the phase.

When implementing the quantum part of Shor's algorithm, the most challenging part is implementing the modular exponentiation gate, which is detailed in the previous tutorial. Its cost is $O(n^3)$. Therefore, the cost of order finding and Shor's algorithm is also $O(n^3)$.

Next, we use isQ to demonstrate finding the factors of 15 following Shor's algorithm. In the order-finding step, we select $x=7$ and try to find its order using a quantum circuit.

// |0> -> |r> where 7 ^ r = 1 (mod 15)
int find_order(qbit k[6]) {
    qbit a[4];
    H(k);
    X(a[0]);    // set |a> = |1>
    pow7mod15(k, a);
    qft_inv(k);
}

procedure main()
{
    qbit k[6];
    find_order(k);
    bp;     // check the amplitude of the state vector
    print M(k);
}

After executing the program, four results can be printed with equal possibilities: 0, 16, 32, and 48. When the returned value is 16 or 48, we can get the fraction $32/64=1/4$ or $48/64=3/4$. Thus, the order is 4, i.e., $7^4(\bmod15)=1$. Otherwise, we cannot infer the order, and the program needs to be executed again. Since gcd$(7^{4/2}-1, 15)=3$, we finally get 3 as a factor of 15.

Grover's algorithm finds with high probability the unique input to a black box function that produces a particular output value, using just $O(\sqrt {N})$ evaluations of the function, where $N$ is the size of the function's domain. As a comparison, the analogous problem in classical computation cannot be solved in fewer than $O(N)$ evaluations. Although Grover's algorithm provides only a quadratic speedup, it is considerable when $N$ is large, and Grover's algorithm can be applied to speed up broad classes of algorithms.

As input for Grover's algorithm, suppose we have a function $f:\{0,1,\cdots,N-1\}\rightarrow\{0,1\}$. In the "unstructured database" analogy, the domain represents indices to a database, and $f(x) = 1$ if and only if the data that $x$ points to satisfies the search criterion. We additionally assume that only one index satisfies $f(x) = 1$, and we call this index $\omega$. Our goal is to identify $\omega$.

We can access $f$ with an oracle in the form of a unitary operator $U_w$ that acts as $$U_\omega|x\rangle=(-1)^{f(x)}|x\rangle$$ i.e., $$\left\{\begin{array}{l} U_\omega|x\rangle=-|x\rangle\quad\text{for }x=\omega,\\ U_\omega|x\rangle=|x\rangle\quad\ \ \text{for }x\neq\omega \end{array}\right.$$ This uses the $N$-dimensional state space $\mathcal{H}$, which is supplied by a register with $n=\lceil\text{log}_2N\rceil$ qubits.

The steps of Grover's algorithm are

1. Initialize the system to uniform superposition over all states $$|s\rangle=\frac{1}{\sqrt{N}}\sum_{x=0}^{N-1}|x\rangle$$

2. Perform the following "Grover iteration" $O(\sqrt{N})$ times:

   1. Apply oracle $U_\omega$;

   2. Apply the Grover diffusion operator $U_s=2|s\rangle\langle s|-I$

3. Measure the resulting quantum state in the computational basis.

The circuit is shown as the following:

We write an example using isQ to demonstrate Grover's algorithm. In this example, $N=8$, corresponding to the Hilbert space of 3 qubits. With this size, we can calculate that 2 Grover iterations are needed. We select $\omega=1$, i.e., only $f(1)=1$. This oracle is defined using a truth table. Next, we define $U0=2|0\rangle\langle0|-I$ using a unitary matrix. The full program is shown below.

import std;

oracle U_omega(3,1) = [0,1,0,0,0,0,0,0];

defgate U0 = [
    1, 0, 0, 0, 0, 0, 0, 0;
    0, -1, 0, 0, 0, 0, 0, 0;
    0, 0, -1, 0, 0, 0, 0, 0;
    0, 0, 0, -1, 0, 0, 0, 0;
    0, 0, 0, 0, -1, 0, 0, 0;
    0, 0, 0, 0, 0, -1, 0, 0;
    0, 0, 0, 0, 0, 0, -1, 0;
    0, 0, 0, 0, 0, 0, 0, -1
];

qbit q[3];
qbit anc;

int grover_search() {
    // Initialize
    H(q);
    X(anc);
    H(anc);

    // Grover iterations
    for i in 0:2 {
        U_omega(q[2], q[1], q[0], anc);
        H(q);
        U0(q[2], q[1], q[0]);
        H(q);
    }

    // Finilize
    H(anc);
    X(anc);
    return M(q);
}

procedure main(){
    print grover_search();
}

We executed the program 100 times using isQ compiler's command bin/isqc run --shots 100 examples/grover.isq. The execution result is

{"010": 1, "110": 1, "011": 1, "001": 97}

Most of the results (97 out of 100) are 1, equal to $\omega$. It shows that Grover's algorithm finds the desired answer with a high probability.

Recursive Fourier sampling

To present the recursive Fourier sampling (RFS) problem, we begin by defining a type of tree. Let $n$ and $l$ be positive integers and consider a symmetric tree where each node, except the leaves, has $2^n$ children, and the depth is $l$. Let the root be labeled by $\varnothing$. The root’s children are labeled $(x_1)$ with $x_1\in\{0,1\}^n$. Each child of $(x_1)$ is, in turn, labeled $(x_1, x_2)$ with $x_2\in\{0,1\}^n$. We continue until we have reached the leaves, which are labeled by $(x_1,\cdots,x_l)$. Thus each node’s label can be thought of as a path describing how to find the node from the root.

Now we add the Fourier component to a tree. We begin by fixing an efficiently computable function $g:\{0,1\}^n\rightarrow\{0,1\}$. With each node of the tree $(x_1,\cdots,x_k)$, we associate a "secret" string $s_{(x_1,\cdots,x_k)}\in\{0,1\}^n$. These secrets are promised to obey $$g\left(s_{(x_1,\cdots,x_k)}\right)=s_{(x_1,\cdots,x_{k-1})}\cdot x_k$$ and the inner product has taken modulo 2. If $k=1$, we take $s_{(x_1,\cdots,x_{k-1})}$ to mean $s_{(\varnothing)}$. In this way, each node’s secret encodes one bit of information about its parent’s secret.

Now suppose that we are given an oracle $A:\left(\{0,1\}^n\right)^l\rightarrow\{0,1\}$ which behaves as $$A(x_1,\cdots,x_l)=g\left(s_{(x_1,\cdots,x_l)}\right)$$ Note that $A$ works for the leaves of the tree only. Our goal is to find $g\left(s_{(\varnothing)}\right)$. This is the recursive Fourier sampling problem.

The recursive nature of the RFS problem is obvious. First, note that the subtree rooted at any node obeys the same promises as the whole tree. Thus each subtree defines an RFS problem. (The subtree rooted at a leaf is a trivial problem, where the oracle just returns the solution.) Thus we have a methodology for solving RFS problems: solve subtrees in order to determine information about the root’s secret, then calculate the secret. Solving subtrees uses the same method, except when we reach a leaf, where the oracle is used instead. This is a type of top-down recursive structure, where the tree is built from a number of subtrees with the same structure.

Before presenting the quantum solution, we first define the behavior of oracles $A$ and $G$: $$A|x_1\rangle\cdots|x_l\rangle|y\rangle=|x_1\rangle\cdots|x_l\rangle|y\oplus g\left(s_{(x_1,\cdots,x_l)}\right)\rangle$$ $$G|s\rangle|y\rangle=|s\rangle|y\oplus g(s)\rangle$$ The main idea behind the algorithm is to use the fact that $H^{\otimes n}$ transforms $|y\rangle$ into $\sum_s(-1)^{x\cdot y}|x\rangle$ and vice versa. We have phase feedback and a call to $A$ to create the state $\sum_{x_l}^{s_{(x_1,\cdots,x_{l-1})}\cdot x_l}|x\rangle$, then apply $H^{\otimes n}$ to obtain $|s_{(x_1,\cdots,x_{l-1})}\rangle$. After calculating $g\left(s_{(x_1,\cdots,x_l)}\right)$, we uncompute $|s_{(x_1,\cdots,x_{l-1})}\rangle$ to disentangle this register from others.

The process of algorithm QRFS is as follows:

  Input: oracle $A$ and (G), $l, k$, quantum registers $x_1,...,x_k, y$

​    1. If $k = l$ then apply $A$ to $(x_1 ...x_l, y)$, then return.

​   2. Introduce ancilla $x_{k+1}$ in the state $\frac{1}{\sqrt{2^n}}\sum_{x \in\{0,1\}^n}|x>$

​   3. Introduce ancilla $y'$ in state $\frac{1}{\sqrt{2}}(|0>-|1>)$

​   4. Call $QRFS(A,G,l,k+1,(x_1...x_{k+1}),y′)$

​   5. Apply $H^{⊗n}$ on register $x_{k+1}$.

​   6. Apply $G$ to $(x_{k+1}, y)$.

​   7. Apply $H^{⊗n}$ on register $x_{k+1}$.

​   8. Call $QRFS(A,G,l,k+1,(x_1...x_{k+1}),y′)$

​   9. Discard $x_{k+1}$ and $y'$.

Now, we use isQ to compose an example of solving RFS. We set $n=l=2$. So there are two labels $x_1, x_2$, each one is 2 bits. $A$ is a unitary gate acting on 5 bits, and $G$ acts on 3 bits. Implement the function recursive_fourier_sampling according to the above algorithm process.

In this example, according to the previous formula, we can calculate $s_{00} = 11, s_{01} = 00, s_{10} = 10, s_{11} = 01 => s_{(\varnothing)} = 10 =>g\left(s_{(\varnothing)}\right) = 1$. So the measurement result of $p[0]$ should be 1

import std;

oracle A(4,1) = [0,1,1,0,0,0,0,0,0,0,1,1,0,1,0,1];
oracle g(2,1) = [0,1,1,0]; 

int a;
qbit q[4], p[3];

procedure hadamard_layer(int k){
    H(q[2*k]);
    H(q[2*k+1]);
}

procedure recursive_fourier_sampling(int k){
    if (k == 2){
        A(q[0], q[1], q[2], q[3], p[2]);
    }else{
        hadamard_layer(k);
        X(p[k+1]);
        H(p[k+1]);

        recursive_fourier_sampling(k+1);

        hadamard_layer(k);
        g(q[2*k],q[2*k+1], p[k]);
        hadamard_layer(k);

        recursive_fourier_sampling(k+1);

        hadamard_layer(k);
        H(p[k+1]);
        X(p[k+1]);
    }
}

procedure main(){
    a = 0;
    recursive_fourier_sampling(a);
    print M(p[0]);
}

Reference

1. McKague, Matthew. "Interactive proofs with efficient quantum prover for recursive Fourier sampling." arXiv preprint arXiv:1012.5699, 2010.

2. Ethan Bernstein and Umesh Vazirani. Quantum complexity theory. SIAM Journal on Computing, 26(5):1411–1473, 1997.

Harrow-Hassidim-Lloyd (HHL) quantum algorithm [1] can be used to solve linear system problems and can provide exponential speedup over the classical conjugate gradient method chosen for this purpose. HHL is the basis of many more advanced algorithms and is important in various applications such as machine learning and modeling of quantum systems.

A linear system problem (LSP) can be represented as the following $$A\vec{x}=\vec{b}$$ where $A$ is an $N_b\times N_b$ Hermitian matrix, and $\vec{x}$ and $\vec{b}$ are $N_b$-dimensional vectors. For simplicity, it is assumed $N_b = 2^{n_b}$. $A$ and $\vec{b}$ are known, and $\vec{x}$ is the unknown to be solved, i.e., $$\vec{x}=A^{−1}\vec{b}$$

Assume the spectral decomposition of $A$ is $$A=\sum_{i=0}^{N_b}\lambda_i|u_i\rangle\langle u_i|$$ where $\{\lambda_i\}$ and $\{|u_i\rangle\}$ are the eigenvalues and the eigenvectors of $A$. Since $\{|u_i\rangle\}$ are orthogonal, $|b\rangle$ can also be decomposed onto these basis: $$|b\rangle=\sum_{i=0}^{N_b}b_i|u_i\rangle$$

Therefore, $$|x\rangle=A^{-1}|b\rangle=\sum_{i=0}^{N_b}\lambda_i^{-1}b_i|u_i\rangle$$

The HHL algorithm consists of four major steps, as shown as follows

• Step 1: Phase Estimation. Let $U=e^{iAt}$. Then $$U=\sum_{j=0}^{N_b}e^{i\lambda_jt}|u_j\rangle\langle u_j|$$ i.e., the eigenvalue corresponding to eigenvector $|u_j\rangle$ is $e^{i\lambda_jt}$. Recall that the quantum phase estimation (QPE) algorithm obtains $N\theta$ for eigenvalue $e^{2\pi i\theta}$. Therefore, applying QPE with $U$ gets $|N_b\lambda_it/2\pi\rangle|u_j\rangle$. We should choose a suitable $t$ so that $\tilde{\lambda}_j=N_b\lambda_it/2\pi$ is (close to) an integer. The qubits that hold $|\tilde{\lambda}_j\rangle$ are usually named the clock register.

• Step 2: Controlled Rotation. This step consists of a gate $R$ that maps $$\left\{\begin{array}{l} |\lambda\rangle|0\rangle\mapsto|\lambda\rangle|0\rangle,\text{ if }\lambda=0\\ |\lambda\rangle|0\rangle\mapsto|\lambda\rangle\left(\sqrt{1-\frac{c^2}{ \lambda^2}}|0\rangle+\frac{c}{\lambda}|1\rangle\right),\text{ if }\lambda\neq0 \end{array}\right.$$

• Step 3: Inverse Phase Estimation. We need to uncompute Step 1 to change the clock register back to $|0\rangle$. Therefore, this step is precisely the inverse of Step 1.

• Step 4: Repeat the Circuit until the measurement result of the ancilla qubit is 1. In this case, the amplitude $c/\lambda$ merges into the output register, whose final state is $$\sum_{i=0}^{N_b}\frac{cb_i}{\lambda}|u_i\rangle$$ without normalization. It is the same with $|x\rangle$ except for a constant $c$. Thus, the algorithm achieves the goal.

To implement the algorithm, we need to further compose two key gates: Hamiltonian simulation gate $e^{iAt}$ and rotation gate $R$. First, we implement the former. Section 4.7.3 of [2] introduces a way of simulating Hamiltonian $H=Z_1\otimes Z_2\otimes Z_3$:

Based on this idea, we can simulate any Pauli matrices product by changing a Pauli gate to $Z$. For example, since $Z=HXH$, we can apply an $H$ gate before the circuit to change $X$ to $Z$. Since any matrix can be decomposed as the sum of Pauli matrices products, i.e., $A=\sum_jP_j$, then $$e^{iAt}=\prod_je^{iP_jt}$$ It means that we can simulate each $P_j$ in sequence, and the overall result is $e^{iAt}$.

To implement $R$, we separate it into $2^{n_c}-1$ control gates where $n_c$ is the number of qubits in the clock register. More specifically, when the clock register is $|\lambda$, we apply a $R_y(\theta)$ gate to the ancilla qubit, where $\theta=2\textrm{arcsin}(c/\lambda)$. Note that $$R_y(\theta)=\begin{bmatrix} \textrm{cos}\frac{\theta}{2} & -\textrm{sin}\frac{\theta}{2}\\ \textrm{sin}\frac{\theta}{2} & \textrm{cos}\frac{\theta}{2} \end{bmatrix}=\begin{bmatrix} \sqrt{1-\frac{c^2}{\lambda^2}} & -\frac{c}{\lambda}\\ \frac{c}{\lambda} & \sqrt{1-\frac{c^2}{\lambda^2}} \end{bmatrix}$$ Therefore, $R_y(\theta)|0\rangle=\sqrt{1-c^2/\lambda^2}|0\rangle+c/\lambda|1\rangle$, precisely as we intended.

Next, we write an isQ program to demonstrate the HHL algorithm. In this example, we will solve the problem $$\begin{bmatrix} -1 & 4\\ 4 & 8 \end{bmatrix}|x\rangle=\begin{bmatrix} 5\\ 16 \end{bmatrix}$$ It can be calculated that the eigenvalues of $A$ are $\lambda=(7\pm\sqrt{145})/2$. We set $t=0.078$ so that $\tilde{\lambda}_j=N_b\lambda_it/2\pi$ are roughly integers. In addition, we set $c=1$ in gate $R$ to maximize the possibility of measuring $|1\rangle$.

The complete source code is as follows.

import std;
import qft;

/*
* Set qInput as B, i.e., |qInput> = c(5|0> + 16|1>)
* where c = \sqrt{5^2 + 16^2} is the normalization factor
*/
procedure SetVectorB(qbit qInput[]) {
    double angle = 2.53582291684;   // 2 * arctan(16/5)
    H(qInput[0]);
    Rz(angle, qInput[0]);
    H(qInput[0]);
    S(qInput[0]);
}

/*
* We encode the tensor product of Paulis as an int.
*
* Two bits are used for each qubit:
*   00 - I, 01 - X, 10 - Y, 11 - Z
*/
int DecodePauliTensor(int N, int code, int i) {
    return code / 4 ** (N - i - 1) % 4;
}

/*
* Simulate controlled e^{iPt} where P is a Pauli product
*
* Based on Section 4.7.3 of [Nielson & Chuang, QCQI, 2010]
*/
procedure ctrlExp(qbit qctrl, int PTcode, double dt, qbit qInput[]) {
    int N = qInput.length;
    qbit qanc;
    for i in 0 : N {
        int p = DecodePauliTensor(N, PTcode, i);
        if (p == 1) { // X
            ctrl H(qctrl, qInput[i]);
            ctrl CNOT(qctrl, qInput[i], qanc);
        }
        if (p == 2) {  // Y
            ctrl S(qctrl, qInput[i]);
            ctrl H(qctrl, qInput[i]);
            ctrl CNOT(qctrl, qInput[i], qanc);
        }
        if (p == 3) {  // Z
            ctrl CNOT(qctrl, qInput[i], qanc);
        }
    }
    ctrl Rz(2 * dt, qctrl, qanc);
    for i in 0 : N {
        int p = DecodePauliTensor(N, PTcode, i);
        if (p == 1) {   // X
            ctrl CNOT(qctrl, qInput[i], qanc);
            ctrl H(qctrl, qInput[i]);
        }
        if (p == 2) {   // Y
            ctrl CNOT(qctrl, qInput[i], qanc);
            ctrl H(qctrl, qInput[i]);
            ctrl inv S(qctrl, qInput[i]);
        }
        if (p == 3) {   // Z
            ctrl CNOT(qctrl, qInput[i], qanc);
        }
    }
}

/*
* Simulate controlled e^{iAt} where t is a full time step
*/
procedure ctrlHSim(qbit qctrl, int Ps[], double Fs[], double t, qbit qInput[], int HSimPrecision) {
    int N = Ps.length;
    assert(N == Fs.length);

    // Trotter formula: lim(e^{iA/n})^n = e^{iA}
    double dt = t / HSimPrecision;
    for j in 0 : HSimPrecision {
        for i in 0 : N {
            double D = dt * Fs[i];
            ctrlExp(qctrl, Ps[i], D, qInput);
        }
    }
}

/*
* |b>|0> -> \sum{b_j|u_j>|v'_j>}
* where v_j and |u_j> are the eigenvalue and eigenvector of A,
*   |b> = \sum{b_j|u_j>},
*   v'_j = 2^n * v_j * t / (2 * \pi)
*/
procedure QPEforHSim(int Ps[], double Fs[], double t, qbit qsclock[], qbit qInput[]) {
    int Nclock = qsclock.length;
    H(qsclock);
    for i in 0 : Nclock {
        double actT = t * 2.0 ** i;
        ctrlHSim(qsclock[i], Ps, Fs, actT, qInput, Nclock);
    }
    qft_inv(qsclock);
}

// The uncompute process of QPEforHSim
procedure invQPEforHSim(int Ps[], double Fs[], double t, qbit qsclock[], qbit qInput[]) {
    int Nclock = qsclock.length;
    qft(qsclock);
    for i in 0 : Nclock {
        double actT = -t * 2.0 ** i;
        ctrlHSim(qsclock[i], Ps, Fs, actT, qInput, Nclock);
    }
    H(qsclock);
}

/*
* |v>|0> -> \sqrt{1-c^2/v^2}|0> + c/v|1>
* where we set c = 1
*/
procedure CRot(qbit qsctrl[], qbit qtar) {
    H(qtar);

    // theta = 2*arcsin(c/v)
    switch qsctrl {
    case |1>:
        Ry(3.141592654, qtar);
    case |2>:
        Ry(1.047197551, qtar);
    case |3>:
        Ry(0.679673819, qtar);
    case |4>:
        Ry(0.50536051, qtar);
    case |5>:
        Ry(0.402715842, qtar);
    case |6>:
        Ry(0.334896158, qtar);
    case |7>:
        Ry(0.286695138, qtar);
    case |8>:
        Ry(0.250655662, qtar);
    case |9>:
        Ry(0.222682029, qtar);
    case |10>:
        Ry(0.200334842, qtar);
    case |11>:
        Ry(0.182069556, qtar);
    case |12>:
        Ry(0.166860173, qtar);
    case |13>:
        Ry(0.153998281, qtar);
    case |14>:
        Ry(0.1429789, qtar);
    case |15>:
        Ry(0.133432297, qtar);
    case |16>:
        Ry(-0.125081524, qtar);
    case |17>:
        Ry(-0.133432297, qtar);
    case |18>:
        Ry(-0.1429789, qtar);
    case |19>:
        Ry(-0.153998281, qtar);
    case |20>:
        Ry(-0.166860173, qtar);
    case |21>:
        Ry(-0.182069556, qtar);
    case |22>:
        Ry(-0.200334842, qtar);
    case |23>:
        Ry(-0.222682029, qtar);
    case |24>:
        Ry(-0.250655662, qtar);
    case |25>:
        Ry(-0.286695138, qtar);
    case |26>:
        Ry(-0.334896158, qtar);
    case |27>:
        Ry(-0.402715842, qtar);
    case |28>:
        Ry(-0.50536051, qtar);
    case |29>:
        Ry(-0.679673819, qtar);
    case |30>:
        Ry(-1.047197551, qtar);
    case |31>:
        Ry(-3.141592654, qtar);
    }

    H(qtar);
    S(qtar);
}

/*
* HHL algorithm body
*
* Solve the problem A|x>=|b> where
*   |x> = x_0|0> + x_1|1>
* x_0 and x_1 are expressed with @NSimPrecision bits
*/
int HHL(int NSimPrecision) {
    // A = 3.5I + 4X - 4.5Z
    //               I  X  Z
    int Apaulis[] = [0, 1, 3];
    double Afactors[] = [3.5, 4, -4.5];

    // The simulation time step
    double facT = 0.078;

    while (true) {
        qbit qAncilla;
        qbit qInput[1];
        qbit qClock[NSimPrecision];

        SetVectorB(qInput);
        QPEforHSim(Apaulis, Afactors, facT, qClock, qInput);
        CRot(qClock, qAncilla);
        invQPEforHSim(Apaulis, Afactors, facT, qClock, qInput);
        if (M(qAncilla)) {
            return M(qInput);
        }
    }
}

procedure main()
{
    int NRepeat = 130;
    int NSimPrecision = 5;
    int Statistics[] = [0, 0];

    for i in 0 : NRepeat {
        int m = HHL(NSimPrecision);
        Statistics[m] += 1;
    }
    for s in Statistics {
        print s;
    }
}

The solution of $|x\rangle$ is $|0\rangle+3/2|1\rangle$ without normalization. Therefore, the printed statistics should be $4:9$ for 0s and 1s.

Reference

1. A. W. Harrow, A. Hassidim, and S. Lloyd, "Quantum Algorithm for Linear Systems of Equations," Phys. Rev. Lett. 103, 150502, 2009.
2. M. A. Nielsen, and I. L. Chuang. "Quantum Computation and Quantum Information: 10th Anniversary Edition." Cambridge University Press, 2010.

The eigenvalues of the Hamiltonian determine almost all properties in a molecule or material of interest. The ground state for molecule Hamiltonian is of particular interest since the energy gap between the ground state and the first excited state of electrons at room temperature is usually larger. Most molecules are in the ground state.

Here, molecular electronic Hamiltonian is represented as $\hat{H}$. A trial wave function $\left| \varphi \left( \theta \right) \right\rangle$ is parameterized with $\overrightarrow{\theta }$, which is called Ansatz. VQE is represented as follows:

$$\frac{\left\langle \varphi \left( {\vec{\theta }} \right) \right|\hat{H}\left| \varphi \left( {\vec{\theta }} \right) \right\rangle }{\left\langle \varphi \left( {\vec{\theta }} \right) \right|\left. \varphi \left( {\vec{\theta }} \right) \right\rangle }\ge {{E}_{0}}$$

${{E}_{0}}$ is the lowest energy of molecular electronic Hamiltonian $\hat{H}$. To estimate the ${{E}_{0}}$, the left-hand side of the equation above is minimized by updating the parameters of the Ansatz $\left| \varphi \left( \theta \right) \right\rangle$.

The molecular electronic Hamiltonian $\hat{H}$ has the second quantized form:

$$\hat{H}=\sum\limits_{pq}{{{h}_{pq}}}\hat{a}_{p}^{+}{{\hat{a}}_{q}}+\frac{1}{2}\sum\limits_{pqrs}{{{h}_{pqrs}}}\hat{a}_{p}^{+}\hat{a}_{q}^{+}{{\hat{a}}_{r}}{{\hat{a}}_{s}}$$

It can be mapped to the linear combination of Pauli operator by Jordan-Wigner transformation or Bravyi-Kitaev transformation as follows:

$$\hat{H}=\sum\limits_{\alpha }{{{\omega }_{\alpha }}{{{\hat{P}}}_{\alpha }}}$$

Then, we can calculate the expectation of Hamitonian $\hat{H}$ with respect to Ansatz $\left| \varphi \left( \theta \right) \right\rangle$ by designing quantum circuit. Finally, updating the parameters in the Ansatz $\left| \varphi \left( \theta \right) \right\rangle$ quantum circuit to get the ${{E}_{0}}$ the lowest energy of molecular electronic Hamiltonian.

The overview of VQE schematic diagram is as follows:

VQE belongs to hybrid quantum-classical algorithms, in which quantum computer is responsible for executing quantum circuit and classical computer is responsible for updating the parameters of quantum gates in the Ansatz $\left| \varphi \left( \theta \right) \right\rangle$. The lowest energy of molecular electronic Hamiltonian $\hat{H}$ can be obtained and the wave function corresponding to ${{E}_{0}}$.

To show the algorithm flow, we take the solution for molecule hydrogen's lowest energy as an example in the following code. After Jordan-Wigner transformation in a minimal basis molecule hydrogen Hamiltonian is

$$\begin{align}& {{{\hat{H}}}_{JW}}=-0.81261I+0.171201\sigma _{0}^{z}+0.171201\sigma _{1}^{z}-0.2227965\sigma _{2}^{z}-0.2227965\sigma _{3}^{z} \\ & +0.16862325\sigma _{1}^{z}\sigma _{0}^{z}+0.12054625\sigma _{2}^{z}\sigma _{0}^{z}+0.165868\sigma _{2}^{z}\sigma _{1}^{z}+0.165868\sigma _{3}^{z}\sigma _{0}^{z} \\ & +0.12054625\sigma _{3}^{z}\sigma _{1}^{z}+0.17434925\sigma _{3}^{z}\sigma _{2}^{z}-0.04532175\sigma _{3}^{x}\sigma _{2}^{x}\sigma _{1}^{y}\sigma _{0}^{y} \\& +0.04532175\sigma _{3}^{x}\sigma _{2}^{y}\sigma _{1}^{y}\sigma _{0}^{x}+0.04532175\sigma _{3}^{y}\sigma _{2}^{x}\sigma _{1}^{x}\sigma _{0}^{y}-0.04532175\sigma _{3}^{y}\sigma _{2}^{y}\sigma _{1}^{x}\sigma _{0}^{x}\\ \end{align}$$

We used two qubit in this example，the Ansatz $\left| \varphi \left( \theta \right) \right\rangle$ is:

Coefficients of the H2 Hamiltonian.
Hamiltonian gates: $Z_0$, $Z_1$, $Z_0Z_1$, $Y_0Y_1$, $X_0X_1$.
For more information, see 2.

coeffs = [-0.4804, +0.3435, -0.4347, +0.5716, +0.0910, +0.0910]
gates_group =  [
    ((0, "Z"),),
    ((1, "Z"),),
    ((0, "Z"), (1, "Z")),
    ((0, "Y"), (1, "Y")),
    ((0, "X"), (1, "X")),
]

h2.isq file.

import std;

int num_qubits = 2; 
int pauli_gates[] = [
    3, 0,
    0, 3,
    3, 3,
    2, 2,
    1, 1
];
// Hamiltonian gates: Z0, Z1, Z0Z1, Y0Y1, X0X2 

unit vqe_measure(qbit q[], int idx) {
    // using arrays for pauli measurement
    // I:0, X:1, Y:2, Z:3 
    int start_idx = num_qubits*idx;
    int end_idx = num_qubits*(idx+1);

    for i in start_idx:end_idx {
        if (pauli_gates[i] == 0) {
            // I:0
            continue;
        }
        if (pauli_gates[i] == 1) {
            // X:1
            H(q[i%num_qubits]);
            M(q[i%num_qubits]);
        }
        if (pauli_gates[i] == 2) {
            // Y:2
            X2P(q[i%num_qubits]);
            M(q[i%num_qubits]);
        }
        if (pauli_gates[i] == 3) {
            // Z:3
            M(q[i%num_qubits]);
        }
    }
}

unit main(int i_par[], double d_par[]) {

    qbit q[2];
    X(q[1]);

    Ry(1.57, q[0]);
    Rx(4.71, q[1]);
    CNOT(q[0],q[1]);
    Rz(d_par[0], q[1]);
    CNOT(q[0],q[1]);
    Ry(4.71, q[0]);
    Rx(1.57, q[1]);

    vqe_measure(q, i_par[0]);
}

Using isqtools to complie isq file.

from isqtools import compile, simulate
compile("h2.isq", target="qir")
# compile "h2.isq" and generate the qir file "h2.so"

Define hyperparameters.

import numpy as np 

shots = 100
learning_rate = 1.
energy_list = []
epochs = 20
theta = np.array([0.0])

Define functions.

def get_expectation(theta):
    measure_results = np.zeros(len(gates_group) + 1)
    measure_results[0] = 1.0
    # The first does not require quantum measurement, which is constant
    # As a result, the other 5 coefficients need to be measured
    for idx in range(len(gates_group)):
        result_dict = simulate(
            "h2.so",
            shots=shots,
            int_param=idx,
            double_param=theta,
        )
        result = 0
        for res_index, frequency in result_dict.items():
            parity = (-1) ** (res_index.count("1") % 2)
            # parity checking to get expectation value
            result += parity * frequency / shots
            # frequency instead of probability
        measure_results[idx + 1] = result
        # to accumulate every expectation result
        # The result is multiplied by coefficient
    return np.dot(measure_results, coeffs)

def parameter_shift(theta):
    # to get gradients
    theta = theta.copy()
    theta += np.pi / 2
    forwards = get_expectation(theta)
    theta -= np.pi
    backwards = get_expectation(theta)
    return 0.5 * (forwards - backwards)

Run VQE.

import time

energy = get_expectation(theta)
energy_list.append(energy)
print(f"Initial VQE Energy: {energy_list[0]} Hartree")

start_time = time.time()
for epoch in range(epochs):
    theta -= learning_rate * parameter_shift(theta)
    energy = get_expectation(theta)
    energy_list.append(energy)
    print(f"Epoch {epoch+1}: {energy} Hartree")

print(f"Final VQE Energy: {energy_list[-1]} Hartree")
print("Time:", time.time() - start_time)

Execution results:

Initial VQE Energy: -0.28835999999999995 Hartree
Epoch 1: -0.33809399999999995 Hartree
Epoch 2: -0.444432 Hartree
Epoch 3: -0.772314 Hartree
Epoch 4: -1.387832 Hartree
Epoch 5: -1.8272059999999999 Hartree
Epoch 6: -1.824372 Hartree
Epoch 7: -1.8609519999999997 Hartree
Epoch 8: -1.8728760000000002 Hartree
Epoch 9: -1.8528600000000002 Hartree
Epoch 10: -1.842748 Hartree
Epoch 11: -1.858316 Hartree
Epoch 12: -1.847396 Hartree
Epoch 13: -1.8391119999999999 Hartree
Epoch 14: -1.855492 Hartree
Epoch 15: -1.8314179999999998 Hartree
Epoch 16: -1.8660059999999998 Hartree
Epoch 17: -1.853266 Hartree
Epoch 18: -1.87206 Hartree
Epoch 19: -1.8629600000000002 Hartree
Epoch 20: -1.8583160000000003 Hartree
Final VQE Energy: -1.8583160000000003 Hartree
Time: 14.92269778251648

Reference

1. J. Tilly, H. Chen, S. Cao, et al. "The variational quantum eigensolver: a review of methods and best practices." Physics Reports, 2022, 986: 1-128.
2. P. O’Malley et al. "Scalable Quantum Simulation of Molecular Energies." Physics Review X, 6, 031007, 2016.

Hamiltonian simulation

Requirement: isqdeployer

isqdeployer is an all-in-one Python package. It serves as an Object-Oriented Interface for isq. You can easily install it using pip without the need for any additional installations or configurations.

pip install isqdeployer

Spin model

Wavefunction evolving with time

For a quantum system with Hamiltonian $H$, its wavefunction is evolving as $|\psi(t)\rangle = e^{-i\int dt H(t)} |\psi(0)\rangle$, where we have assume that $\hbar=1$. If we further have $H\ne H(t)$, then we have $|\psi(t)\rangle = e^{-i H t} |\psi(0)\rangle = \sum_ie^{-iE_i t}\langle i|\psi(0)\rangle|i\rangle$, where we have used the spectral decomposition of $H=\sum_i E_i |i\rangle\langle i|$. This usually needs us to obtain a diagonalization of $H$. Practically, a more efficient way [1] to simulate the time evolution is used if we have an efficient way to implement operators directly. Here we recall its major ideal.

For any Hamiltonian in the form of $H=\sum_k h_k$, we have $e^{-iHt}=( e^{-iH\tau} )^N$, where $\tau=t/N$ for some large $N$ and therefore $\tau\to0$. Then we have the following approximation:
$$e^{-i\sum_k h_k \tau}\approx \prod_k e^{-i h_k \tau}.$$
So the basic element for simulating the whole Hamiltonian is to simulate each term $e^{-i h_k \tau}$.

Basic gates in spin model

In spin model, the Hamiltonian can be represented as $$H= \sum _{ij;mn}\lambda^{mn}_{ij}\sigma_i^m \sigma_j^n,$$ where ${\{\sigma_i^m}\}=\{I,X,Y,Z\}$ at each spin site $m$. Correspondingly the basic gate we need to implement in the circuit is in the form of $e^{i \theta \sigma_i \sigma_j}$ for the pair of any two sites.

$e^{i a \sigma}$ gate implementation

We first consider the term $e^{i a Z}$ on one site. The corresponding implementation is given by

Similarly, we can directly implement $e^{i a X}$ and $e^{i a Y}$ by Rx and Ry. These terms represent the effect of the magnetic field in each direction.

$e^{i a \sigma_i^m \sigma_j^n}$ gate implementation

This is a two-spin interaction with action $\sigma_i$ on one site $m$ and $\sigma_j$ on the other site $n$. We first consider the term $e^{i a ZZ}$. This term represents an Ising type of interaction between two spins. Its implementation is given by:

For this gate, we can calculate the output state of each input state in $\{ |00\rangle,|01\rangle,|10\rangle,|11\rangle\}$ and see that it is equivalent to $e^{iaZZ}$. Based on this gate, we can easily implement other types. For implementing $e^{iaXZ}$, for instance, we can use the transformation $XZ =(U^+ZU)Z=U^+ZZU$, where $U=Ry(\frac{\pi}{2})$ acting on the first qubit (spin). Its implementation is shown by

By the same spirit, we can implement any $\sigma_i^m\sigma_j^n$ terms.

Example: Single spin in the magnetic field

Consider a system with Hamiltonian $H=X$, for input state $|\psi\rangle=|0\rangle$, the evolution of the wavefunction is give by

$$|\psi(t)\rangle= e^{-it}\frac{1}{\sqrt{2}}|+\rangle +e^{it}\frac{1}{\sqrt{2}}|-\rangle = \cos t |0\rangle - i \sin t |0\rangle .$$

We then simulate this process by quantum circuit.

from isqdeployer import Circuit 
from isqdeployer.circuitDeployer import PauliHamDeployer 
from isqdeployer.utils.pauliHamiltonian import PauliHamiltonian
import numpy as np 
import matplotlib.pyplot as plt 

ham = PauliHamiltonian(nq=1) # define a Hamiltonian
ham.setOneTerm(xi=1.0,p=[1]) # Hamiltonian = Pauli X, (0,1,2,3) represents (I,X,Y,Z)

cir = Circuit(
    nq=1, # only 1 qubit
    workDir="/home/user/Desktop/test", # see how it works by inspecting internal process 
    isInputArg=True, # use additional parameter (t) to increase the speed of the calculation
    ) # circuit, init state is |0>

t = cir.getInputArg(FI='F', id=0) # set the circuit to contains a parameter t

dep = PauliHamDeployer(circuit=cir) # use a deployer to implement gates
dep.expHt(
    h=ham,
    t=t,
    N=15, # Suzuki decomposition
    ) 
cir.setMeasurement([0])

Tlist = np.arange(0,3,0.1)# the time range we study
results = cir.runJob(paramList=[{'F':[t]} for t in Tlist ]) # batch submit all jobs

#---- plot 
prob0 = [res.get('0',0) for res in results]
prob1 = [res.get('1',0) for res in results]
plt.plot(Tlist,prob0,label="Prob(0)")
plt.plot(Tlist,prob1,label="Prob(1)")
plt.xlabel("Time (t)")
plt.ylabel('Probability')
plt.legend()
plt.show()

In the example above, we have set workDir="/home/user/Desktop/test". You can observe intermediate resources in that directory. The original isq file is named resource.isq. The results are calculated using an internal default simulator and are displayed below:

Example: Spin chain model

We now study a less trivial model whose Hamiltonian is given by
$$H = h\sum_{n}Z_n+\sum_{n}[J_{x}^{(n)}X_nX_{n+1}+J_{y}^{(n)}Y_nY_{n+1}+J_{z}^{(n)}Z_nZ_{n+1}].$$
The first term represents a uniform magnetic field on the system. The second term represents complex XYZ interactions. This model has very important applications in condensed matter physics.

In this demo, we set up a 3-site system. Parameters are set as $h=-1$, $J_x=-0.2$, $J_y=-0.3$, and $J_z=-0.1$. Here we do not consider the connection between the first and third sites. One interested in periodic boundary conditions can set additional interaction. In addition, in the beginning, we set the first spin in the state $|1\rangle$ while others are in $|0\rangle$.

from isqdeployer import Circuit 
from isqdeployer.circuitDeployer import PauliHamDeployer 
from isqdeployer.utils.pauliHamiltonian import PauliHamiltonian
import numpy as np 
import matplotlib.pyplot as plt 

# --- define the Hamiltonian ---
ham = PauliHamiltonian(nq=3) # define a Hamiltonian
# --- hz term
ham.setOneTerm(xi=-1.,p=[3,0,0])
ham.setOneTerm(xi=-1.,p=[0,3,0])
ham.setOneTerm(xi=-1.,p=[0,0,3])
# ---jx term 
ham.setOneTerm(xi=-0.2,p=[1,1,0])
ham.setOneTerm(xi=-0.2,p=[0,1,1])
# ---jy term 
ham.setOneTerm(xi=-0.3,p=[2,2,0])
ham.setOneTerm(xi=-0.3,p=[0,2,2])
# ---jz term 
ham.setOneTerm(xi=-0.1,p=[3,3,0])
ham.setOneTerm(xi=-0.1,p=[0,3,3])


cir = Circuit(nq=3,isInputArg=True,)  

t = cir.getInputArg(FI='F', id=0) # set the circuit to contains a parameter t

cir.X(0) # set the first spin up

dep = PauliHamDeployer(circuit=cir) # use a deployer to implement gates
dep.expHt( 
    h=ham, 
    t=t, 
    N=50, # Suzuki decomposition
    ) 
cir.setMeasurement([0,1,2])

Tlist = np.arange(0,10,0.1)# the time range we study
results = cir.runJob(paramList=[{'F':[t]} for t in Tlist ]) # batch submit all jobs


#---- plot 
def countProb(i,res):
    r = 0.0
    for k, v in res.items():r += (k[i] == "0")*v
    return r*2-1

Z0 = [ countProb(0,res) for res in results]
Z1 = [ countProb(1,res) for res in results]
Z2 = [ countProb(2,res) for res in results]

plt.plot(Tlist,Z0)
plt.plot(Tlist,Z1)
plt.plot(Tlist,Z2)

The result of the time evolution of each spin state is shown as follows:

As a comparison, an exact calculation by QuTip is shown in the figure. The precision can further increase by increasing $N$. To be more clear, we check $\langle\sigma_z^1(t)\rangle$ more carefully by the following figure.

As shown in the figure, for $N=2$ case, in a small time regime the result of the calculation matches the exact result very well. However, at large t regime, to obtain a better result we need a larger $N$.

Any other spin mode

One can easily extend the code to make a circuit to simulate an arbitrary spin model with more complicated interaction types. For any (2-local) spin system, we can represent it as $H=\sum_k \xi_k \sigma_{i(k)}^{m(k)} \sigma_{j(k)}^{n(k)}$. An it can be calcualted by isqdeployer.

Fermion-Hubbard model

The Hubbard model is one of the most important systems in condensed matter physics. Studying this model will reveal much exotic physics in strongly correlated electron systems such as superconductivity and quantum Hall effect. In this tutorial, we study the time evolution of a wave function in the system with such a Hamiltonian.

Hamiltonian

We start from a most general Hamiltonian given by $H=H_0+H_I$ where $$H_0=\sum_{ij}t_{ij}c_{i}^{\dagger}c_{j}$$ and
$$H_I=\sum_{ijkl}u_{ijkl}c_{i}^{\dagger}c_{j}^{\dagger}c_{k}c_{l}.$$ For a system with spins, the indices $i$ and $j$ should also be implicit. $H_0$ represents the single particle terms, such as onsite energy and hopping terms. $H_I$ contains many-body terms. For example, if we set $i=k$ and $j=l$, we obtain the Coulomb interaction $$H_U = U\sum_i n_{i\uparrow} n_{i\downarrow}.$$ In the most general case, $H_I$ can consider other interactions, such as Hund's interaction, and is written as $$H_I = H_U + H_J + ...$$ In the standard Hubbard, we consider $H_0$ and $H_U$. Nevertheless, the following discussion is also suited for other two-particle interactions.

Our task is to find an encoding method to transform this Hamiltonian into one with Pauli tensors only. By Jordan–Wigner transformation $c_{i}=\sigma_{i}^{-}\prod_{\alpha<i}\left(-Z_{\alpha}\right)$, the diagonal term of $H_0$ is given by $$t_{ii}c^\dagger_i c_i = \frac{t_{ii}}{2}\left(1+Z_{i}\right),$$ the hopping term (off-diagonal term) is given by

$$\begin{aligned} t_{ij}c^\dagger_{i}c_j+t_{ji}c^\dagger_jc_i= -\chi_{ij}X_{2-}^{i}Y_{2+}^{j}\mathcal{Z}_{ij}X_{2+}^{i}Y_{2-}^{j} -\xi_{ij}X_{2-}^{i}X_{2-}^{j}\mathcal{Z}_{ij}X_{2+}^{i}X_{2+}^{j}\\ -\xi_{ij}Y_{2+}^{i}Y_{2+}^{j}\mathcal{Z}_{ij}Y_{2-}^{i}Y_{2-}^{j} +\chi_{ij}Y_{2+}^{i}X_{2-}^{j}\mathcal{Z}_{ij}Y_{2-}^{i}X_{2+}^{j},\\ \end{aligned}$$
where $\chi_{ij}= i\frac{t_{ij}-t_{ji}}{2}\frac{\left(-1\right)^{j-i}}{2}$, $\xi_{ij}=\frac{t_{ij}+t_{ji}}{2}\frac{\left(-1\right)^{j-i}}{2}$, and $$\mathcal{Z}_{ij}= Z_i\otimes Z_{i+1}...\otimes Z_{j},$$ and $X_{2\pm} = R_x (\pm\frac{\pi}{2} )$, and $Y_{2\pm}=R_y(\pm\frac{\pi}{2})$. For real $t_{ij}$ case, we have $\chi_{ij}=0$. Finally the Coulomb interaction is given by $$H_U=\frac{U}{4}\left(1+Z_{j}+Z_{i}+Z_{i}Z_{j}\right).$$

To streamline configuration, we employ a specialized package, qsci, for deploying circuits in scientific calculations. One can install it by:

pip install qsci

With qsci, it's straightforward to deploy various quantum algorithms such as VQE, Hubbard model, and many others into circuits with different types of backends.

Example

In the following, we will provide a step-by-step example demonstrating coding with isqdeployer and qsci. Let's get started!

Let study a 4-site system with Hamiltonian given by $$H=t\sum_{i=1}^{3}(c^\dagger_{i\sigma}c_{i+1\sigma}+h.c.)-\mu\sum_{i\sigma} c^\dagger_{i\sigma } c_{i\sigma}+U\sum_{i}c^\dagger_{i\uparrow}c_{i\uparrow}c^\dagger_{i\downarrow}c_{i\downarrow},$$

where we have defined $i=4$ means $i=0$.We define the Hamiltonian of a Fermion system as follows:

from qsci.secondQuantization.hamiltonian.fermion import SpinfulHamiltonian
hamHubbard = SpinfulHamiltonian()

The hopping terms are acting like a ring. We set the parameters as $t=-1.0$, $\mu=1.0$, $U=2.0$ by:

t=-1.0
U=2.0
mu=1.0
for i in range(4): hamHubbard.setHopping(i=i,j=(i+1)%4,tij=t) # set hopping
for i in range(4): hamHubbard.setCoulombU(i=i,U=U) # set U
for i in range(4): hamHubbard.setOnSiteEnergy(i=i,m=-mu) # set mu

Since there is the spin degree of freedom, we need to use 2 qubits to represent a single site in the Hubbard model. The initial state, as shown in figure, is set by:

cir.X(0).X(1) 

Next, we need to choose the Jordan-Wigner (JW) transformation to encode the Hamiltonian into a circuit model:

from qsci.secondQuantization.qubitEncoding import JordanWignerEncoding
jw = JordanWignerEncoding() # use Jordan-Wigner encoding
ham = hamHubbard.exportPauliOperator(jw).exportPauliHamiltonian() # Pauli Hamiltonian

The following section is similar to previous examples, where the simulation of time evolution (controlled by the exp factor) is managed by isqdeployer. A complete runnable code is provided below:

from qsci.secondQuantization.hamiltonian.fermion import SpinfulHamiltonian
from qsci.secondQuantization.qubitEncoding import JordanWignerEncoding
from isqdeployer.circuit import Circuit
from isqdeployer.circuitDeployer.pauliHamiltonian import Deployer
import matplotlib.pyplot as plt
import numpy as np 
#---- parameters of Hubbard model ----
t=-1.0
U=2.0
mu=1.0
#-------------- set Hubbard model Hamiltonian -------------------
hamHubbard = SpinfulHamiltonian()
for i in range(4): hamHubbard.setHopping(i=i,j=(i+1)%4,tij=t) # set hopping
for i in range(4): hamHubbard.setCoulombU(i=i,U=U) # set U
for i in range(4): hamHubbard.setOnSiteEnergy(i=i,m=-mu) # set mu

#-------------- encode Hamiltonian into Pauli type
jw = JordanWignerEncoding() # use Jordan-Wigner encoding
ham = hamHubbard.exportPauliOperator(jw).exportPauliHamiltonian() # Pauli Hamiltonian
#-------------- run on circuit
cir = Circuit(nq=8,isInputArg=True) # 4 site Hubbard model map into 8 qubits
cir.X(0).X(1) #-- set initial state 
dpl = Deployer(circuit=cir) 
dpl.expHt(h=ham,t=cir.getInputArg('F',0),N=50) 
cir.setMeasurement(range(8))

#---- run job
tList = np.arange(0,5,0.1)
results = cir.runJob(paramList=[{'F':[t]} for t in tList])

#--- plot 
def N(res,i):
    r = 0 
    for k,v in res.items():
        r += (k[2*i] == "1")*v
    return r 
niList = [ [N(res,i) for res in results] for i in range(4) ]
for i in range(4):plt.plot( tList, niList[i], '-o', label=f"$N_{{{i}\\uparrow}}$" )
plt.legend()
plt.title("simulation with N=50")
plt.xlabel("t")
plt.ylabel("particle number")
plt.show()

The results with $N=50$ and $N=500$ are shown below. We also calculate the result by exact method (classical simulation). As we can see, for the $N=50$ case, the result at a small time region agrees with the exact results very well. At the large time region, the time slide $\tau=T/N$ becomes large, therefore the error increase. As one can see from the geometry of the system, sites 1 and 3 are symmetric, therefore we know that $\langle N_{1\sigma}\rangle=\langle N_{3\sigma}\rangle$. To obtain a better result, we can set $N=500$ as shown in the figure.

Reference

1. S. Lloyd. “Universal Quantum Simulators.” Science 273 (5278), 1996:1073–78.
2. J. Hubbard. "Electron Correlations in Narrow Energy Bands." Proc. R. Soc.Lond. A 276 (1365), 1963: 238–57.
3. P. Jordan and E. Wigner. "Ber Das Paulische Quivalenzverbot." Z. Physik 47 (9-10), 1928: 631–51.

Ended: Examples