ch6_4shor.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "PF3uTLhqxyEc",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 15882,
     "status": "ok",
     "timestamp": 1727772202484,
     "user": {
      "displayName": "Eun-Jin Im",
      "userId": "12874128355928418134"
     },
     "user_tz": -540
    },
    "id": "PF3uTLhqxyEc",
    "outputId": "8890dd5a-d46d-44af-f2f5-3f068d267745"
   },
   "outputs": [],
   "source": [
    "!pip install qiskit qiskit-ibm-runtime pylatexenc qiskit_aer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af8c299d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiskit\n",
    "qiskit.version.get_version_info()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf6c8b1f-c1e0-46d5-bc17-8541b4121610",
   "metadata": {
    "executionInfo": {
     "elapsed": 828,
     "status": "ok",
     "timestamp": 1727772264367,
     "user": {
      "displayName": "Eun-Jin Im",
      "userId": "12874128355928418134"
     },
     "user_tz": -540
    },
    "id": "bf6c8b1f-c1e0-46d5-bc17-8541b4121610",
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister, transpile\n",
    "from qiskit.visualization import plot_histogram, array_to_latex\n",
    "from qiskit.circuit.library import QFT,PhaseEstimation\n",
    "\n",
    "import numpy as np\n",
    "from numpy.random import randint\n",
    "from math import gcd,lcm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2eb96d97",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 15 \n",
    "a = 7\n",
    "\n",
    "xvals = range(N)\n",
    "yvals = [1]\n",
    "y = 1\n",
    "for x in range(1,N):\n",
    "    y = np.mod(a*y, N)\n",
    "    yvals += [y] \n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(xvals, yvals, linewidth=1, linestyle='solid', marker='o')\n",
    "ax.set(xlabel='$x$', ylabel=f'${a}^x$ mod ${N}$')\n",
    "ax.set_xticks(xvals,labels=xvals)\n",
    "ax.set_yticks(yvals,labels=yvals)\n",
    "\n",
    "try: \n",
    "    r = yvals[1:].index(1) + 1\n",
    "    plt.annotate('', xy=(0,1), xytext=(r,1),\n",
    "                 arrowprops=dict(arrowstyle='<->'))\n",
    "    plt.annotate(f'$r={r}$', xy=(r/3,1.5))\n",
    "except ValueError:\n",
    "    print(f'An order not found, {a} may be larger than {N} or {a} and {N} has a common factor {gcd(a,N)}.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70187c4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure(figsize=(15,1), dpi=100)\n",
    "ax = fig.add_subplot(111)\n",
    "plt.subplots_adjust()\n",
    "ax.spines[\"top\"].set_color(\"None\")\n",
    "ax.spines[\"right\"].set_color(\"None\")\n",
    "ax.spines[\"left\"].set_color(\"None\")\n",
    "ax.set_yticks([],labels=\"\")\n",
    "\n",
    "x = {0.0:'0', 1.0:'1'}\n",
    "for denom in range(7,-1,-1):\n",
    "    for num in range(1,denom):\n",
    "        x[float(num)/denom] = f'{num}/{denom}'\n",
    "x_sorted = sorted(x.items(), key=lambda x:x[0])\n",
    "x_ticks = [x[0] for x in x_sorted]\n",
    "x_labels = [x[1] for x in x_sorted]\n",
    "ax.set_xticks(x_ticks,labels=x_labels)\n",
    "mark = 108/128\n",
    "ax.annotate(\"108/128\",(mark,0),(mark-0.025,0.7), arrowprops={\"arrowstyle\":\"->\", \"color\":\"0.0\"})\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af23bdf6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from fractions import Fraction\n",
    "# 분모가 7 이하인 유리수 중에서 108/128 과 가장 가까운 유리수의 분자와 분모를 찾는다.\n",
    "f = Fraction(108/128).limit_denominator(7)\n",
    "r = f.denominator\n",
    "s = f.numerator\n",
    "print(f\"estimated phase = {s}/{r}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tanNhg1VG-O2",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 357,
     "status": "ok",
     "timestamp": 1727772568674,
     "user": {
      "displayName": "Eun-Jin Im",
      "userId": "12874128355928418134"
     },
     "user_tz": -540
    },
    "id": "tanNhg1VG-O2",
    "outputId": "161b73ea-b288-4d14-aa4c-8e93776cb501"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import randint\n",
    "from math import gcd,lcm\n",
    "# factoring N=21\n",
    "N =  21  # number to be factored\n",
    "n =   5   # number of qubits representing N\n",
    "t =  2*n + 1 # number of evaluation qubits\n",
    "a =  randint(2,N)\n",
    "g = gcd(a,N)\n",
    "print(f\"a = {a}, gcd({a},{N}) = {g}\") # if g != 1, g is a factor of N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f31620f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# myQPE() returns phase estimation circuit \n",
    "from qiskit.circuit.library import QFT\n",
    "\n",
    "def myQPE(t, unitary):\n",
    "    n = unitary.num_qubits\n",
    "    rx = QuantumRegister(t,'qx')\n",
    "    ry = QuantumRegister(n,'qy')    \n",
    "    qc = QuantumCircuit(rx, ry)\n",
    "    \n",
    "    # QuantumCircuit 객체 unitary 를 Gate 객체 gate 로 만들고 이름을 \"U\"로 정한다.\n",
    "    gate = unitary.to_gate()\n",
    "    gate.name = \"U\"\n",
    "\n",
    "    # 제어-U 게이트 c_gate 를 만들고 QPE 회로를 생성한다.\n",
    "    c_gate = gate.control(1)\n",
    "    \n",
    "    qc.h(rx)\n",
    "    for i in range(t):\n",
    "        for j in range(2**i):\n",
    "            qc.append(c_gate, rx[i:i+1]+ry[:])\n",
    "\n",
    "    qc.append(QFT(t,inverse=True),rx)\n",
    "    qc.name = 'myQPE'\n",
    "    return qc\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "WHTwxKUPycYD",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "executionInfo": {
     "elapsed": 1783051,
     "status": "ok",
     "timestamp": 1727774363744,
     "user": {
      "displayName": "Eun-Jin Im",
      "userId": "12874128355928418134"
     },
     "user_tz": -540
    },
    "id": "WHTwxKUPycYD",
    "outputId": "155059c6-7def-4928-f1bd-c49a44ea2dac"
   },
   "outputs": [],
   "source": [
    "# generate U matrix of (a*i) mod N\n",
    "two_power_n = 1 << n # 2^n\n",
    "U = np.eye(two_power_n) # Identity matrix\n",
    "for j in range(N):\n",
    "    U[j][j] = 0\n",
    "    U[(a*j)%N][j] = 1\n",
    "# U 행렬을 구현하는 양자 회로 qcU 를 생성\n",
    "qcU = QuantumCircuit(n)\n",
    "qcU.unitary(U, range(n))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7bb960c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from qiskit.visualization import array_to_latex\n",
    "array_to_latex(U, prefix=\"U = \", max_size=32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17a7c986",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a order-finding circuit using myQPE\n",
    "qx = QuantumRegister(t,'x')\n",
    "qy = QuantumRegister(n,'y')\n",
    "c = ClassicalRegister(t,'c')\n",
    "qc_of = QuantumCircuit(qx, qy, c)\n",
    "qc_of.x(qy[0])                 # input state = |1>\n",
    "qc_of.append(myQPE(t,qcU),qx[:]+qy[:])\n",
    "qc_of.measure(qx,c)\n",
    "qc_of.draw(output='mpl',fold=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc86c755",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Define two matrices\n",
    "A = np.array([[1, 2], [3, 4]])\n",
    "B = np.array([[5, 6], [7, 8]])\n",
    "print('A =\\n', A)\n",
    "print('B =\\n', B)\n",
    "# element-wise multiplication\n",
    "print('A * B =\\n', A*B)\n",
    "# Perform matrix multiplication using @ operator\n",
    "print('A @ B =\\n', A@B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae8e471e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate a order-finding circuit\n",
    "\n",
    "qx = QuantumRegister(t,'x')\n",
    "qy = QuantumRegister(n,'y')\n",
    "c = ClassicalRegister(t,'c')\n",
    "qc_of = QuantumCircuit(qx, qy, c)\n",
    "qc_of.x(qy[0])                 # input state = |1>\n",
    "\n",
    "qc_of.h(qx)\n",
    "\n",
    "two_power_n = 1 << n # 2^n\n",
    "u_power = np.eye(two_power_n) # Identity matrix\n",
    "for j in range(N):\n",
    "    u_power[j][j] = 0\n",
    "    u_power[(a*j)%N][j] = 1\n",
    "\n",
    "for i in range(t):\n",
    "    qc_u_power = QuantumCircuit(n)\n",
    "    qc_u_power.unitary(u_power, range(n))\n",
    "    gate_u_power = qc_u_power.to_gate()\n",
    "    gate_u_power.name = f'[*{a}mod{N}]^2^{i}'\n",
    "    c_gate_u_power = gate_u_power.control(1)    \n",
    "    qc_of.append(c_gate_u_power, [qx[i]]+qy[:])\n",
    "    u_power = u_power @ u_power\n",
    "\n",
    "    \n",
    "qc_of.append(QFT(t,inverse=True),qx)\n",
    "\n",
    "qc_of.measure(qx,c)\n",
    "qc_of.draw(output='mpl',fold=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "141139a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
    "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
    "from qiskit_aer import AerSimulator\n",
    "aer_sim = AerSimulator()\n",
    "pm = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)\n",
    "isa_circuit = pm.run(qc_of)\n",
    "sampler = Sampler(mode=aer_sim)\n",
    "job = sampler.run([isa_circuit], shots=1000)\n",
    "result = job.result()\n",
    "count = result[0].data.c.get_counts()\n",
    "print(count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5db4afec",
   "metadata": {},
   "outputs": [],
   "source": [
    "from fractions import Fraction\n",
    "r_pre = 1<<t  # r_pre == 2^t\n",
    "estimated = dict()\n",
    "for i in range(N):\n",
    "    estimated[i] = 0\n",
    "print('measured value    estimated order')\n",
    "for measured in count :\n",
    "    s_pre = int(measured,2) # convert a binary string to an integer    \n",
    "    f = Fraction(s_pre/r_pre).limit_denominator(N-1)\n",
    "    r = f.denominator\n",
    "    s = f.numerator\n",
    "    print(f'{measured} -> {s_pre}/{r_pre} -> {s}/{r} -> {r}' )\n",
    "    estimated[r] += count[measured]\n",
    "print('order : frequency')\n",
    "for i in range(N):\n",
    "    print(f'  {i}  :  {estimated[i]}') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10aba3bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_order(sampler, isa_circuit, t):\n",
    "    job = sampler.run([isa_circuit], shots=1)\n",
    "    result = job.result()\n",
    "    count = result[0].data.c.get_counts()\n",
    "    phase = list(count.keys())\n",
    "    s_pre = int(phase[0],2) # convert a binary string to an integer\n",
    "    r_pre = 1<<t  # r_pre == 2^t\n",
    "    f = Fraction(s_pre/r_pre).limit_denominator(N-1)\n",
    "    return  f.denominator\n",
    "\n",
    "def confirm_order(a, r, N):\n",
    "    return (a**r % N == 1)\n",
    "\n",
    "r = estimate_order(sampler, isa_circuit, t)\n",
    "if confirm_order(a, r, N):\n",
    "    found = True\n",
    "else : \n",
    "    r1 = r\n",
    "    r2 = estimate_order(sampler, isa_circuit, t)\n",
    "    if confirm_order(a, r2, N):\n",
    "        r = r2\n",
    "        found = True\n",
    "    else :\n",
    "        r = lcm(r1, r2)\n",
    "        if r >= N :\n",
    "            found = False\n",
    "        elif confirm_order(a, r, N):\n",
    "            found = True\n",
    "        else :\n",
    "            found = False\n",
    "    print(f\"r1={r1} r2={r2} estimated r={r}\")\n",
    "\n",
    "if (found) :\n",
    "    print(f\"estimated value {r} is an order of {a} for {N}\")\n",
    "else :\n",
    "    print(f\"{r} is NOT an order of {a} for {N} because {a}^{r} % {N} = {a**r%N}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "rGYqLVobG-O4",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 338,
     "status": "ok",
     "timestamp": 1727785167425,
     "user": {
      "displayName": "Eun-Jin Im",
      "userId": "12874128355928418134"
     },
     "user_tz": -540
    },
    "id": "rGYqLVobG-O4",
    "outputId": "d418c6ec-033e-458d-ea8b-05f1afd5e1ed"
   },
   "outputs": [],
   "source": [
    "if r % 2 == 0 : # r is even\n",
    "    a_power_r_over2 = a**int(r/2)\n",
    "    p = gcd(a_power_r_over2+1,N)\n",
    "    print (f\"p = gcd( {a}^({r}/2) + 1 (= {a_power_r_over2+1}), {N} ) = {p}\")\n",
    "    if (p == N)  :\n",
    "        print (f\"{a}^({r}/2) + 1 (= {a_power_r_over2+1}) is a multiple of {N}.\")\n",
    "        print (f\" a = {a} is not lucky. try again.\")\n",
    "    elif (p == 1) :\n",
    "        print (f\" r = {r} is not an order of {a}. try again.\")\n",
    "    else :\n",
    "        q = gcd(a_power_r_over2-1,N)\n",
    "        print (f\"q = gcd( {a}^({r}/2) - 1 (= {a_power_r_over2-1}), {N} ) = {q}\")\n",
    "        print (f\"{p} and {q} are factors of {N}\")\n",
    "else :\n",
    "    print (f\"r = {r} is an odd number. a = {a} is not lucky. try again.\")"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}