量子計算及量子計算的模擬

IT是一個繁榮的行業，寄託著無數人的夢想，充斥著無數的造夢神話。
IT是一個悲催的行業，層出不窮的新概念讓人應接不暇，幾乎只要有一天不學習，都可能讓你寢食不安。

量子計算機是一個炒的比較熱的概念，目前還處於上升期，感覺上已經到了爆發的邊緣，似乎隨時可以呼之欲出。
通常對於量子計算機的理解就是，因為量子計算機的儲存特徵，可以處理很大的資料，而不是像傳統計算機那樣只是處理1、0二進位制數，因此計算效率更高。從而有可能顛覆現有的計算機架構，甚至現有的所有的加密演算法，因為新的、高速的量子計算機的出現，都將因為可能會被快速的解密而失去效用。
這個概念對，但也不全對，並且可能丟失了很多重要的內容。所以這裡試圖更通俗的解釋一番。

傳統計算機

說量子計算之前，我們首先要看一下傳統的計算機是如何工作的：
體系結構、硬碟、記憶體、CPU啥的就不用說了，對於計算本身來說，這些體現不出來什麼不同。
我們要從CPU來解析，當前不管多麼複雜的計算機，計算的根本來自於兩個部件：
暫存器 ：用於儲存計算用的資料，及計算的結果，比如當前的64位元CPU,其實就代表暫存器是由64位元二進位制陣列成的。
邏輯閘 ：用於通過各種複雜拼接、組合，從而完成所需的計算。常見的邏輯有：與、或、非。
下面對應傳統的計算機，我們說說量子計算機的的原理。

測量

對於傳統計算機來講，測量就是讀取，讀取又是所有運算的基礎，功能沒啥好多說的，就是知道暫存器中儲存的數值是什麼。
量子計算機就複雜了，你應當聽說過薛定諤的貓，貓放在箱子裡面，箱子裡面的放射性元素是否衰變決定著貓的死活，但是在你開啟箱子之前，你無法確定，而在開啟箱子的瞬間，你實際也改變了箱子的狀態，也就是影響了貓的生死，這是量子物理重要的一個特徵。
這個特徵產生了以下幾個重要的概念，或者說重要的特徵，這些特徵決定了量子計算機的幾個重要功能。

疊加態

傳統的計算機，無論是暫存器還是記憶體，很重要的一個指標就是確定性或者說穩定性。在暫存器方面，其中的每一位(bit),可以確定是1或者是0，64位元計算機，暫存器儲存的數值最大是2的64次方（2^64）。
量子計算機採用疊加態在暫存器中儲存資料，因此同一個暫存器，可能儲存了多個資料，這導致量子暫存器可以儲存的數值，在同樣位數下，比傳統計算機儲存的資料更多。
當前的量子計算機，每一位(bit),通過疊加態可以儲存2對可以互相轉換的狀態，可以分別表示兩位二進位制數位。因為兩位二進位制數共有4個狀態，因此量子計算機暫存器指定位長下，可以最大表達4的N次方(4^N)的數位。
而同時因為疊加態及可互相轉換的特徵，實際上每個指定位長的暫存器，都可能儲存2^N個資料，而不是1個，這就是量子計算機的超強儲存能力（本項能力只是基於理論設想，在當前的各種量子機實現中，還沒有看到資料介紹實際的實現）。

概率

因為疊加態和測量也會對其中數值造成的影響，實際其中的值是由概率決定的。這在量子計算機的製造和演算法的研究中，都必須考慮到的問題。

量子密碼

因為不可測的特徵帶來的無法竊聽和不可克隆特徵，強大的量子計算能力雖然對傳統的密碼學是一個災難，但同時也會出現新的、更強大的加密演算法。從公開資料上看，我們國家所實現的星-地量子通訊的主要基礎就是來自於此。

向量計算和平行計算

因為上面說的疊加態的兩對狀態，每一個計算實際都轉換為了向量運算。現在計算機中最耗能的幾項運算：挖礦、圖形影象、深度學習等，都是轉換成向量和矩陣進行運算的。在傳統的計算中，因為二進位制的特徵，一般複雜度都是O(2^N),而在量子計算機中，因為這種特徵，每次是兩個資料參與運算，所以複雜度是O(N)。所以運算量不再呈現指數飆升，而是線性增加。
此外，因為疊加態的存在，我們在量子運算中，不必像傳統計算一樣同時只處理一個數值，而是同時處理疊加的多個數值，實現了真正的平行計算。而傳統計算中的並行，不過是把大的並行轉換成多個小的時間片再序列而已。這更進一步的提高了運算能力。當然這也需要更複雜的並行演算法來支援。

單量子位元門

如同傳統計算機一樣，量子計算機也是通過邏輯閘的運算來完成實際運算的。
同時因為上面說過的那些特徵，單量子位元門也有一個很重要的特徵：可逆性，我們知道，傳統計算機的邏輯閘是不可逆的，而單量子位元門的各種運算都可逆；另外一個重要特徵就是上面說過的：可並行。
這些常用的門中，包括以下幾個：
Hadamard：旋轉門
CNOT:受控非門，如果第一位置1，倒置第二位，否則保持不變。
Toffoli：控-控-非門，也叫CCNOT，如果前兩位置1，它將倒置第三位，否則所有位保持不變
noop：單位門，等於不做任何操作。
X門：求非變換，NOT門
Z門：相位移動操作
Y門：相當於上面兩個門的組合，Y=ZX

量子計算的模擬

目前的情況，除非是在相關單位工作，否則一般的開發人員尚無法親身體驗量子計算機。更談不上學習量子計算機的開發了。
但是通過上面的介紹你可以發現，除非只做一個拿來主義的使用者，否則這些顛覆性的特徵將帶來演算法上的重大改變，不及早的研究、積累，這真的可能成為程式設計師的一個“灰犀牛”。
除了在實際的量子計算機上實驗，目前也有很多軟體提供了量子計算的模擬能力，從而可以嘗試自己的演算法和實驗，達到學習的目的。
微軟甚至創立了一種新的語言叫“Q#”來應對將來的量子計算，相信應當也不錯，不過最近對於非開源的專案還是有些障礙，所以我們來嘗試另外一個工具庫DLIB。

環境搭建

首先要假設你已經有了基本的開發環境，比如g++/clang。Ubuntu要保證APT包管理的正常使用，macOS則預先安裝好Homebrew包管理。

在Ubuntu16上：

apt install libdlib-dev

在macOS：

brew install dlib

DEMO原始碼

本原始碼來自DLIB的Examples程式碼quantum_computing_ex.cpp，未做任何修改。
程式碼中使用上面介紹過的這幾個邏輯閘，實現了兩個最常用的基本演算法：Grover Search和Shor ECC校驗。這兩個演算法也是相當有名，網上一搜資料大把，我這半瓶水就不畫蛇添足了。

// The contents of this file are in the public domain. See LICENSE_FOR_EXAMPLE_PROGRAMS.txt
/*    This is an example illustrating the use of the quantum computing    simulation classes from the dlib C++ Library.    This example assumes you are familiar with quantum computing and    Grover's search algorithm and Shor's 9 bit error correcting code    in particular.  The example shows how to simulate both of these    algorithms.    The code to simulate Grover's algorithm is primarily here to show    you how to make custom quantum gate objects.  The Shor ECC example    is simpler and uses just the default gates that come with the    library.*/

#include <iostream>
#include <complex>
#include <ctime>
#include <dlib/quantum_computing.h>
#include <dlib/string.h>

using namespace std;
using namespace dlib;

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

// This declares a random number generator that we will be using below 
dlib::rand rnd;

// ----------------------------------------------------------------------------------------

void shor_encode (
    quantum_register& reg
)
/*!
    requires
        - reg.num_bits() == 1
    ensures
        - #reg.num_bits() == 9
        - #reg == the Shor error coding of the input register
!*/
{
    DLIB_CASSERT(reg.num_bits() == 1,"");

    quantum_register zeros;
    zeros.set_num_bits(8);
    reg.append(zeros);

    using namespace dlib::quantum_gates;
    const gate<1> h = hadamard();
    const gate<1> i = noop();

    // Note that the expression (h,i) represents the tensor product of the 1 qubit
    // h gate with the 1 qubit i gate and larger versions of this expression
    // represent even bigger tensor products.  So as you see below, we make gates
    // big enough to apply to our quantum register by listing out all the gates we
    // want to go into the tensor product and then we just apply the resulting gate
    // to the quantum register.

    // Now apply the gates that constitute Shor's encoding to the input register. 
    (cnot<3,0>(),i,i,i,i,i).apply_gate_to(reg);
    (cnot<6,0>(),i,i).apply_gate_to(reg);
    (h,i,i,h,i,i,h,i,i).apply_gate_to(reg);
    (cnot<1,0>(),i,cnot<1,0>(),i,cnot<1,0>(),i).apply_gate_to(reg);
    (cnot<2,0>(),cnot<2,0>(),cnot<2,0>()).apply_gate_to(reg);
}

// ----------------------------------------------------------------------------------------

void shor_decode (
    quantum_register& reg
)
/*!
    requires
        - reg.num_bits() == 9
    ensures
        - #reg.num_bits() == 1
        - #reg == the decoded qubit that was in the given input register
!*/
{
    DLIB_CASSERT(reg.num_bits() == 9,"");

    using namespace dlib::quantum_gates;
    const gate<1> h = hadamard();
    const gate<1> i = noop();

    // Now apply the gates that constitute Shor's decoding to the input register

    (cnot<2,0>(),cnot<2,0>(),cnot<2,0>()).apply_gate_to(reg);
    (cnot<1,0>(),i,cnot<1,0>(),i,cnot<1,0>(),i).apply_gate_to(reg);

    (toffoli<0,1,2>(),toffoli<0,1,2>(),toffoli<0,1,2>()).apply_gate_to(reg);

    (h,i,i,h,i,i,h,i,i).apply_gate_to(reg);

    (cnot<6,0>(),i,i).apply_gate_to(reg);
    (cnot<3,0>(),i,i,i,i,i).apply_gate_to(reg);
    (toffoli<0,3,6>(),i,i).apply_gate_to(reg);

    // Now that we have decoded the value we don't need the extra 8 bits any more so
    // remove them from the register.
    for (int i = 0; i < 8; ++i)
        reg.measure_and_remove_bit(0,rnd);
}

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

// This is the function we will use in Grover's search algorithm.  In this
// case the value we are searching for is 257.
bool is_key (unsignedlong n)
{
    return n == 257;
}

// ----------------------------------------------------------------------------------------

template <int bits>
class uf_gate;

namespace dlib {
template <int bits>
struct gate_traits<uf_gate<bits> >
{
    static const long num_bits = bits;
    static const long dims = dlib::qc_helpers::exp_2_n<num_bits>::value;
};}

template <int bits>
class uf_gate : public gate_exp<uf_gate<bits> >
{
    /*!
        This gate represents the black box function in Grover's search algorithm.
        That is, it is the gate defined as follows:
            Uf|x>|y> = |x>|y XOR is_key(x)>

        See the documentation for the gate_exp object for the details regarding
        the compute_state_element() andoperator() functions defined below.
    !*/
public:
    uf_gate() : gate_exp<uf_gate>(*this) {}

    static const long num_bits = gate_traits<uf_gate>::num_bits;
    static const long dims = gate_traits<uf_gate>::dims;

    const qc_scalar_type operator() (long r, long c) const{
        unsigned long output = c;
        // if the input control bit is set
        if (is_key(output>>1))
        {
            output = output^0x1;
        }

        if ((unsigned long)r == output)
            return 1;
        else
            return 0;
    }

    template <typename exp>
    qc_scalar_type compute_state_element (
        const matrix_exp<exp>& reg,
        long row_idx
    ) const{
        unsigned long output = row_idx;
        // if the input control bit is set
        if (is_key(output>>1))
        {
            output = output^0x1;
        }

        return reg(output);
    }
};

// ----------------------------------------------------------------------------------------

template <int bits>
class w_gate;

namespace dlib {
template <int bits>
struct gate_traits<w_gate<bits> >
{
    static const long num_bits = bits;
    static const long dims = dlib::qc_helpers::exp_2_n<num_bits>::value;
}; }

template <int bits>
class w_gate : public gate_exp<w_gate<bits> >
{
    /*!
        This is the W gate from the Grover algorithm
    !*/
public:

    w_gate() : gate_exp<w_gate>(*this) {}

    static const long num_bits = gate_traits<w_gate>::num_bits;
    static const long dims = gate_traits<w_gate>::dims;

    const qc_scalar_type operator() (long r, long c) const{
        qc_scalar_type res = 2.0/dims;
        if (r != c)
            return res;
        else
            return res - 1.0;
    }

    template <typename exp>
    qc_scalar_type compute_state_element (
        const matrix_exp<exp>& reg,
        long row_idx
    ) const{
        qc_scalar_type temp = sum(reg)*2.0/dims;
        // compute this value: temp = temp - reg(row_idx)*2.0/dims + reg(row_idx)*(2.0/dims - 1.0)
        temp = temp - reg(row_idx);

        return temp;
    }
};

// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------

int main()
{
    // seed the random number generator
    rnd.set_seed(cast_to_string(time(0)));

    // Pick out some of the gates we will be using below
    using namespace dlib::quantum_gates;
    const gate<1> h = quantum_gates::hadamard();
    const gate<1> z = quantum_gates::z();
    const gate<1> x = quantum_gates::x();
    const gate<1> i = quantum_gates::noop();

    quantum_register reg;

    // We will be doing the 12 qubit version of Grover's search algorithm.
    const int bits=12;
    reg.set_num_bits(bits);

    // set the quantum register to its initial state
    (i,i, i,i,i,i,i, i,i,i,i,x).apply_gate_to(reg);

    // Print out the starting bits
    cout << "starting bits: ";
    for (int i = reg.num_bits()-1; i >= 0; --i)
        cout << reg.probability_of_bit(i);
    cout << endl;

    // Now apply the Hadamard gate to all the input bits
    (h,h, h,h,h,h,h, h,h,h,h,h).apply_gate_to(reg);

    // Print out the status
    cout << "after Hadamard gate: ";
    for (int i = reg.num_bits()-1; i >= 0; --i)
        cout << reg.probability_of_bit(i) << " ";
    cout << endl;

    // Here we do the grover iteration
    for (int j = 0; j < 35; ++j)
    {
        (uf_gate<bits>()).apply_gate_to(reg);
        (w_gate<bits-1>(),i).apply_gate_to(reg);

        cout << j << " probability: bit 1 = "<< reg.probability_of_bit(1) << ", bit 9 = "<< reg.probability_of_bit(9) << endl;
    }

    cout << endl;

    // Print out the final probability of measuring a 1 for each of the bits
    for (int i = reg.num_bits()-1; i >= 1; --i)
        cout << "probability for bit "<< i << " = "<< reg.probability_of_bit(i) << endl;
    cout << endl;

    cout << "The value we want grover's search to find is 257 which means we should measure a bit pattern of 00100000001"<< endl;
    cout << "Measured bits: ";
    // finally, measure all the bits and print out what they are.
    for (int i = reg.num_bits()-1; i >= 1; --i)
        cout << reg.measure_bit(i,rnd);
    cout << endl;

    // Now let's test out the Shor 9 bit encoding
    cout << "nnnnNow let's try playing around with Shor's 9bit error correcting code"<< endl;

    // Reset the quantum register to contain a single bit
    reg.set_num_bits(1);
    // Set the state of this single qubit to some random mixture of the two computational bases
    reg.state_vector()(0) = qc_scalar_type(rnd.get_random_double(),rnd.get_random_double());
    reg.state_vector()(1) = qc_scalar_type(rnd.get_random_double(),rnd.get_random_double());
    // Make sure the state of the quantum register is a unit vector
    reg.state_vector() /= sqrt(sum(norm(reg.state_vector())));

    cout << "state: "<< trans(reg.state_vector());

    shor_encode(reg);
    cout << "x bit corruption on bit 8"<< endl;
    (x,i,i,i,i,i,i,i,i).apply_gate_to(reg); // mess up the high order bit
    shor_decode(reg); // try to decode the register

    cout << "state: "<< trans(reg.state_vector());

    shor_encode(reg);
    cout << "x bit corruption on bit 1"<< endl;
    (i,i,i,i,i,i,i,x,i).apply_gate_to(reg);
    shor_decode(reg);

    cout << "state: "<< trans(reg.state_vector());

    shor_encode(reg);
    cout << "z bit corruption on bit 8"<< endl;
    (z,i,i,i,i,i,i,i,i).apply_gate_to(reg);
    shor_decode(reg);

    cout << "state: "<< trans(reg.state_vector());

    cout << "nThe state of the input qubit survived all the corruptions in tact so the code works."<< endl;

}

編譯的方法如下：
Ubuntu:

g++ -std=c++11 -O3 -o quantum_computing_ex quantum_computing_ex.cpp -l dlib -l cblas

macOS:

g++ -std=c++11 -O3 -o quantum_computing_ex quantum_computing_ex.cpp (pkg-config --cflags --libs dlib-1) -l cblas