From 371730107438b45b76beec54dcb31699dabff515 Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Thu, 20 Jun 2019 14:12:07 +0800 Subject: [PATCH 1/3] update Shor algorithm --- examples/order_finding.jl | 282 -------------------------------------- src/Mod.jl | 105 ++++++++++++++ src/QuAlgorithmZoo.jl | 1 + src/number_theory.jl | 101 ++++++++++++++ src/shor.jl | 80 +++++++++++ test/shor.jl | 65 +++++++++ 6 files changed, 352 insertions(+), 282 deletions(-) delete mode 100644 examples/order_finding.jl create mode 100644 src/Mod.jl create mode 100644 src/number_theory.jl create mode 100644 src/shor.jl create mode 100644 test/shor.jl diff --git a/examples/order_finding.jl b/examples/order_finding.jl deleted file mode 100644 index b2628f6..0000000 --- a/examples/order_finding.jl +++ /dev/null @@ -1,282 +0,0 @@ -using Yao, YaoBlocks, BitBasis, LuxurySparse, LinearAlgebra - -function YaoBlocks.cunmat(nbit::Int, cbits::NTuple{C, Int}, cvals::NTuple{C, Int}, U0::Adjoint, locs::NTuple{M, Int}) where {C, M} - YaoBlocks.cunmat(nbit, cbits, cvals, copy(U0), locs) -end - -"""x^Nz%N""" -function powermod(x::Int, k::Int, N::Int) - rem = 1 - for i=1:k - rem = mod(rem*x, N) - end - rem -end - -Z_star(N::Int) = filter(i->gcd(i, N)==1, 0:N-1) -Eulerφ(N) = length(Z_star(N)) - -"""obtain `s` and `r` from `ϕ` that satisfies `|s/r - ϕ| ≦ 1/2r²`""" -continued_fraction(ϕ, niter::Int) = niter==0 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) -continued_fraction(ϕ::Rational, niter::Int) = niter==0 || ϕ.den==1 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) - -""" -Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group. -""" -function mod_inverse(x::Int, N::Int) - for i=1:N - (x*i)%N == 1 && return i - end - throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!")) -end - -is_order(r, x, N) = powermod(x, r, N) == 1 - -"""get the order, the classical approach.""" -function get_order(::Val{:classical}, x::Int, N::Int) - findfirst(r->is_order(r, x, N), 1:N) -end - -function rand_primeto(L) - while true - x = rand(2:L-1) - d = gcd(x, L) - if d == 1 - return x - end - end -end - -function shor(L, ver=Val(:quantum); maxiter=100) - L%2 == 0 && return 2 - # some classical method to accelerate the solution finding - for i in 1:maxiter - x = rand_primeto(L) - r = get_order(ver, x, L) - # if `x^(r/2)` is non-trivil, go on. - # Here, we do not condsier `powermod(x, r÷2, L) == 1`, since in this case the order should be `r/2` - if r%2 == 0 && powermod(x, r÷2, L) != L-1 - f1, f2 = gcd(powermod(x, r÷2, L)-1, L), gcd(powermod(x, r÷2, L)+1, L) - if f1!=1 - return f1 - elseif f2!=1 - return f2 - else - error("Algorithm Fail!") - end - end - end -end - -""" - Mod{N} <: PrimitiveBlock{N} - -calculates `mod(a*x, L)`, notice `gcd(a, L)` should be 1. -""" -struct Mod{N} <: PrimitiveBlock{N} - a::Int - L::Int - function Mod{N}(a, L) where N - @assert gcd(a, L) == 1 && L<=1<= m.L ? i+1 : mod(i*m.a, m.L)+1 - for j in 1:B - @inbounds nstate[_i,j] = reg.state[i+1,j] - end - end - reg.state = nstate - reg -end - -function Yao.mat(::Type{T}, m::Mod{N}) where {T, N} - perm = Vector{Int}(undef, 1<= m.L ? i+1 : mod(i*m.a, m.L)+1] = i+1 - end - PermMatrix(perm, ones(T, 1< (b&mask, b>>k) -end - -function Yao.apply!(reg::ArrayReg{B}, m::KMod{N, K}) where {B, N, K} - YaoBlocks._check_size(reg, m) - nstate = zero(reg.state) - - reader = bint2_reader(Int, K) - for b in basis(reg) - k, i = reader(b) - _i = i >= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L) - _b = k + _i<= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L) - _b = k + _i< c; nshots=nshots) - reader = bint2_reader(Int, K) - for r in res - k, i = reader(r) - # get s/r - ϕ = bfloat(k) # - ϕ == 0 && continue - - order = order_from_float(ϕ, x, L) - if order === nothing - continue - else - return order - end - end - return nothing -end - -function order_from_float(ϕ, x, L) - k = 1 - rnum = continued_fraction(ϕ, k) - while rnum.den < L - r = rnum.den - @show r - if is_order(r, x, L) - return r - end - k += 1 - rnum = continued_fraction(ϕ, k) - end - return nothing -end - -using Test -function check_Euler_theorem(N::Int) - Z = Z_star(N) - Nz = length(Z) # Eulerφ - for x in Z - @test powermod(x,Nz,N) == 1 # the order is a devisor of Eulerφ - end -end - -@testset "Euler" begin - check_Euler_theorem(150) -end - -@testset "Mod" begin - @test_throws AssertionError Mod{4}(4,10) - @test_throws AssertionError Mod{2}(3,10) - m = Mod{4}(3,10) - @test mat(m) ≈ applymatrix(m) - @test isunitary(m) - @test isunitary(mat(m)) - @test m' == Mod{4}(7,10) -end - -@testset "KMod" begin - @test_throws AssertionError KMod{6, 2}(4,10) - @test_throws AssertionError KMod{4, 2}(3,10) - m = KMod{6, 2}(3,10) - @test mat(m) ≈ applymatrix(m) - @test isunitary(m) - @test isunitary(mat(m)) - @test m' == KMod{6, 2}(7,10) -end - -using Random -@testset "shor_classical" begin - Random.seed!(129) - L = 35 - f = shor(L, Val(:classical)) - @test f == 5 || f == 7 - - L = 25 - f = shor(L, Val(:classical)) - @test_broken f == 5 - - L = 7*11 - f = shor(L, Val(:classical)) - @test f == 7 || f == 11 - - L = 14 - f = shor(L, Val(:classical)) - @test f == 2 || f == 7 -end - -using Random -@testset "shor quantum" begin - Random.seed!(129) - L = 15 - f = shor(L, Val(:quantum)) - @test f == 5 || f == 3 -end diff --git a/src/Mod.jl b/src/Mod.jl new file mode 100644 index 0000000..cdd07df --- /dev/null +++ b/src/Mod.jl @@ -0,0 +1,105 @@ +# TODO +# compile Mod and KMod to elementary gates. + +export Mod, KMod + +""" + Mod{N} <: PrimitiveBlock{N} + +calculates `mod(a*x, L)`, notice `gcd(a, L)` should be 1. +""" +struct Mod{N} <: PrimitiveBlock{N} + a::Int + L::Int + function Mod{N}(a, L) where N + @assert gcd(a, L) == 1 && L<=1<= m.L ? i+1 : mod(i*m.a, m.L)+1 + for j in 1:B + @inbounds nstate[_i,j] = reg.state[i+1,j] + end + end + reg.state = nstate + reg +end + +function Yao.mat(::Type{T}, m::Mod{N}) where {T, N} + perm = Vector{Int}(undef, 1<= m.L ? i+1 : mod(i*m.a, m.L)+1] = i+1 + end + PermMatrix(perm, ones(T, 1< (b&mask, b>>k) +end + +function Yao.apply!(reg::ArrayReg{B}, m::KMod{N, K}) where {B, N, K} + YaoBlocks._check_size(reg, m) + nstate = zero(reg.state) + + reader = bint2_reader(Int, K) + for b in basis(reg) + k, i = reader(b) + _i = i >= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L) + _b = k + _i<= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L) + _b = k + _i< Vector + +returns the Z* group elements of `N`, i.e. {x | gcd(x, N) == 1} +""" +Z_star(N::Int) = filter(i->gcd(i, N)==1, 0:N-1) +Eulerφ(N) = length(Z_star(N)) + +""" + continued_fraction(ϕ, niter::Int) -> Rational + +obtain `s` and `r` from `ϕ` that satisfies `|s/r - ϕ| ≦ 1/2r²` +""" +continued_fraction(ϕ, niter::Int) = niter==0 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) +continued_fraction(ϕ::Rational, niter::Int) = niter==0 || ϕ.den==1 ? floor(Int, ϕ) : floor(Int, ϕ) + 1//continued_fraction(1/mod(ϕ, 1), niter-1) + +""" + mod_inverse(x::Int, N::Int) -> Int + +Return `y` that `(x*y)%N == 1`, notice the `(x*y)%N` operations in Z* forms a group and this is the definition of inverse. +""" +function mod_inverse(x::Int, N::Int) + for i=1:N + (x*i)%N == 1 && return i + end + throw(ArgumentError("Can not find the inverse, $x is probably not in Z*($N)!")) +end + +""" + is_order(r, x, N) -> Bool + +Returns true if `r` is the order of `x`, i.e. `r` satisfies `x^r % N == 1`. +""" +is_order(r, x, N) = powermod(x, r, N) == 1 + +""" + find_order(x::Int, N::Int) -> Int + +Find the order of `x` by brute force search. +""" +function find_order(x::Int, N::Int) + findfirst(r->is_order(r, x, N), 1:N) +end + +""" + rand_primeto(N::Int) -> Int + +Returns a random number `2 ≦ x < N` that is prime to `N`. +""" +function rand_primeto(N::Int) + while true + x = rand(2:N-1) + d = gcd(x, N) + if d == 1 + return x + end + end +end + +""" + order_from_float(ϕ, x, L) -> Int + +Estimate the order of `x` to `L`, `r`, from a floating point number `ϕ ∼ s/r` using the continued fraction method. +""" +function order_from_float(ϕ, x, L) + k = 1 + rnum = continued_fraction(ϕ, k) + while rnum.den < L + r = rnum.den + if is_order(r, x, L) + return r + end + k += 1 + rnum = continued_fraction(ϕ, k) + end + return nothing +end + +""" + factor_a_power_b(N::Int) -> (Int, Int) or nothing + +Factorize `N` into the power form `a^b`. +""" +function factor_a_power_b(N::Int) + y = log2(N) + for b = 2:ceil(Int, y) + x = 2^(y/b) + u1 = floor(Int, x) + u1^b == N && return (u1, b) + (u1+1)^b == N && return (u1+1, b) + end +end + +end diff --git a/src/shor.jl b/src/shor.jl new file mode 100644 index 0000000..4c0c9f5 --- /dev/null +++ b/src/shor.jl @@ -0,0 +1,80 @@ +include("number_theory.jl") +include("Mod.jl") + +export shor, order_finding_circuit, get_order + +""" + shor(L::Int, ver=Val(:quantum); maxtry=100) + +factorize an integer `L`, `ver` can be either `Val(:quantum)` or `Val(:classical)`. +""" +function shor(L::Int, ver=Val(:quantum); maxtry=100) + L%2 == 0 && return 2 + + # find solutions like `a^b` + res = factor_a_power_b(L) + res !== nothing && return res[1] + + for i in 1:maxtry + x = NumberTheory.rand_primeto(L) + r = get_order(ver, x, L) + # if `x^(r/2)` is non-trivil, go on. + # Here, we do not condsier `powermod(x, r÷2, L) == 1`, since in this case the order should be `r/2` + if r%2 == 0 && powermod(x, r÷2, L) != L-1 + f1, f2 = gcd(powermod(x, r÷2, L)-1, L), gcd(powermod(x, r÷2, L)+1, L) + if f1!=1 + return f1 + elseif f2!=1 + return f2 + else + error("Algorithm Fail!") + end + end + end +end + +"""estimate the required size of the output register.""" +estimate_ncbit(nbit::Int, ϵ::Real) = 2*nbit + 1 + ceil(Int,log2(2+1/2ϵ)) + +""" + order_finding_circuit(x::Int, L::Int; nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25)) -> AbstractBlock + +Returns the circuit for finding the order of `x` to `L`, +feeding input `|1>⊗|0>` will get the resulting quantum register with the desired "phase" information. +""" +function order_finding_circuit(x::Int, L::Int; nbit::Int=bit_length(L-1), ncbit::Int=estimate_ncbit(nbit, 0.25)) + N = nbit+ncbit + chain(N, repeat(N, H, 1:ncbit),KMod{N, ncbit}(x, L), concentrate(N, QFTBlock{ncbit}()', 1:ncbit)) +end + +""" + find_order([ver], x::Int, N::Int; nshots=10) -> Union{Int, Nothing} + +Get the order of `x`, `ver` can be `Val(:classical)` (default) or `Val(:quantum)`, +when using the quantum approach, we can set key word arguments `nshot`, +`nbit` (size of input data register) and `ncbit` (size of control register, or output register). +""" +get_order(::Val{:classical}, x::Int, N::Int; kwargs...) = NumberTheory.find_order(x, N) +function get_order(::Val{:quantum}, x::Int, N::Int; nshots=10, kwargs...) + c = order_finding_circuit(x, N; kwargs...) + n = nqubits_data(c[2]) + ncbit = nqubits_control(c[2]) + reg = join(product_state(n, 1), zero_state(ncbit)) + + res = measure(copy(reg) |> c; nshots=nshots) + reader = bint2_reader(Int, ncbit) + for r in res + k, i = reader(r) + # get s/r + ϕ = bfloat(k) # + ϕ == 0 && continue + + order = NumberTheory.order_from_float(ϕ, x, N) + if order === nothing + continue + else + return order + end + end + return nothing +end diff --git a/test/shor.jl b/test/shor.jl new file mode 100644 index 0000000..b909d26 --- /dev/null +++ b/test/shor.jl @@ -0,0 +1,65 @@ +using Test, QuAlgorithmZoo, Yao +using Random +using QuAlgorithmZoo.NumberTheory + +"""Euler theorem states that the order is a devisor of Eulerφ (or the size of Z* group of `N`)""" +function check_Euler_theorem(N::Int) + Z = Z_star(N) + Nz = length(Z) # Eulerφ + for x in Z + @test powermod(x,Nz,N) == 1 # the order is a devisor of Eulerφ + end +end + +@testset "Euler" begin + check_Euler_theorem(150) +end + +@testset "Mod" begin + @test_throws AssertionError Mod{4}(4,10) + @test_throws AssertionError Mod{2}(3,10) + m = Mod{4}(3,10) + @test mat(m) ≈ applymatrix(m) + @test isunitary(m) + @test isunitary(mat(m)) + @test m' == Mod{4}(7,10) +end + +@testset "KMod" begin + @test_throws AssertionError KMod{6, 2}(4,10) + @test_throws AssertionError KMod{4, 2}(3,10) + m = KMod{6, 2}(3,10) + @test mat(m) ≈ applymatrix(m) + @test isunitary(m) + @test isunitary(mat(m)) + @test m' == KMod{6, 2}(7,10) +end + +@testset "shor_classical" begin + Random.seed!(129) + L = 35 + f = shor(L, Val(:classical)) + @test f == 5 || f == 7 + + L = 25 + f = shor(L, Val(:classical)) + @test f == 5 + + L = 7*11 + f = shor(L, Val(:classical)) + @test f == 7 || f == 11 + + L = 14 + f = shor(L, Val(:classical)) + @test f == 2 || f == 7 + + @test factor_a_power_b(25) == (5, 2) + @test factor_a_power_b(15) == nothing +end + +@testset "shor quantum" begin + Random.seed!(129) + L = 15 + f = shor(L, Val(:quantum)) + @test f == 5 || f == 3 +end From 9593b69e268a102462c839bb8aa59096dd275f3b Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Thu, 20 Jun 2019 16:56:22 +0800 Subject: [PATCH 2/3] fix test --- src/shor.jl | 2 +- test/runtests.jl | 4 ++++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/shor.jl b/src/shor.jl index 4c0c9f5..c64fe89 100644 --- a/src/shor.jl +++ b/src/shor.jl @@ -12,7 +12,7 @@ function shor(L::Int, ver=Val(:quantum); maxtry=100) L%2 == 0 && return 2 # find solutions like `a^b` - res = factor_a_power_b(L) + res = NumberTheory.factor_a_power_b(L) res !== nothing && return res[1] for i in 1:maxtry diff --git a/test/runtests.jl b/test/runtests.jl index 331bbd6..b540221 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -49,3 +49,7 @@ end @testset "Diff" begin include("Diff.jl") end + +@testset "Shore" begin + include("shor.jl") +end From f152c0da31ae8eacabca031f751637746ffed71e Mon Sep 17 00:00:00 2001 From: GiggleLiu Date: Thu, 20 Jun 2019 21:44:27 +0800 Subject: [PATCH 3/3] add shor to readme --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index 1a052ba..433f5de 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,7 @@ Disclaimer: **this package is still under development and needs further polish.* - [x] Hadamard Test - [x] State Overlap Algorithms - [x] Quantum SVD +- [x] Shor ## License