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
 
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<<N
-        new{N}(a, L)
-    end
-end
-
-function Yao.apply!(reg::ArrayReg{B}, m::Mod{N}) where {B, N}
-    YaoBlocks._check_size(reg, m)
-    nstate = zero(reg.state)
-    for i in basis(reg)
-        _i = i >= 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<<N)
-    for i in basis(N)
-        @inbounds perm[i >= m.L ? i+1 : mod(i*m.a, m.L)+1] = i+1
-    end
-    PermMatrix(perm, ones(T, 1<<N))
-end
-
-Base.adjoint(m::Mod{N}) where N = Mod{N}(mod_inverse(m.a, m.L), m.L)
-Yao.print_block(io::IO, m::Mod{N}) where N = print(io, "Mod{$N}: $(m.a)*x % $(m.L)")
-
-Yao.isunitary(::Mod) = true
-# Yao.ishermitian(::Mod) = true  # this is not true for L = 1.
-# Yao.isreflexive(::Mod) = true  # this is not true for L = 1.
-
-"""
-    KMod{N, K} <: PrimitiveBlock{N}
-
-The first `K` qubits are exponent `k`, and the rest `N-K` are base `a`,
-it calculates `mod(a^k*x, L)`, notice `gcd(a, L)` should be 1.
-"""
-struct KMod{N, K} <: PrimitiveBlock{N}
-    a::Int
-    L::Int
-    function KMod{N, K}(a, L) where {N, K}
-        @assert gcd(a, L) == 1 && L<=1<<(N-K)
-        new{N, K}(a, L)
-    end
-end
-
-nqubits_data(m::KMod{N, K}) where {N, K} = N-K
-nqubits_control(m::KMod{N, K}) where {N, K} = K
-
-function bint2_reader(T, k::Int)
-    mask = bmask(T, 1:k)
-    return b -> (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<<K + 1
-        for j in 1:B
-            @inbounds nstate[_b,j] = reg.state[b+1,j]
-        end
-    end
-    reg.state = nstate
-    reg
-end
-
-function Yao.mat(::Type{T}, m::KMod{N, K}) where {T, N, K}
-    perm = Vector{Int}(undef, 1<<N)
-    reader = bint2_reader(Int, K)
-    for b in basis(N)
-        k, i = reader(b)
-        _i = i >= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L)
-        _b = k + _i<<K + 1
-        @inbounds perm[_b] = b+1
-    end
-    PermMatrix(perm, ones(T, 1<<N))
-end
-
-Base.adjoint(m::KMod{N, K}) where {N, K} = KMod{N, K}(mod_inverse(m.a, m.L), m.L)
-Yao.print_block(io::IO, m::KMod{N, K}) where {N, K} = print(io, "Mod{$N, $K}: $(m.a)^k*x % $(m.L)")
-
-Yao.isunitary(::KMod) = true
-# Yao.ishermitian(::Mod) = true  # this is not true for L = 1.
-# Yao.isreflexive(::Mod) = true  # this is not true for L = 1.
-
-estimate_K(nbit::Int, ϵ::Real) = 2*nbit + 1 + ceil(Int,log2(2+1/2ϵ))
-
-using QuAlgorithmZoo: QFTBlock
-function order_finding_circuit(x::Int, L::Int; nbit::Int=bit_length(L-1), K::Int=estimate_K(nbit, 0.25))
-    N = nbit+K
-    chain(N, repeat(N, H, 1:K), KMod{N, K}(x, L), concentrate(N, QFTBlock{K}()', 1:K))
-end
-
-function shor(L::Int; nshots=10)
-    x = rand_primeto(L)
-end
-
-function get_order(::Val{:quantum}, x::Int, L::Int; nshots=10)
-    c = order_finding_circuit(x, L)
-    n = nqubits_data(c[2])
-    K = nqubits_control(c[2])
-    reg = join(product_state(n, 1), zero_state(K))
-
-    res = measure(copy(reg) |> 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/Diff.jl b/src/Diff.jl
index f237ee3..f54a895 100644
--- a/src/Diff.jl
+++ b/src/Diff.jl
@@ -1,4 +1,4 @@
-export Rotor, generator, AbstractDiff, BPDiff, QDiff, backward!, gradient, CPhaseGate, DiffBlock
+export Rotor, generator, Diff, backward!, gradient, CPhaseGate, DiffBlock
 import Yao: expect, content, chcontent, mat, apply!
 using StatsBase
 
@@ -13,81 +13,33 @@ Return the generator of rotation block.
 """
 generator(rot::RotationGate) = rot.block
 generator(rot::PutBlock{N, C, GT}) where {N, C, GT<:RotationGate} = PutBlock{N}(generator(rot|>content), rot |> occupied_locs)
-generator(c::CphaseGate{N}) where N = ControlBlock(N, c.ctrol_locs, ctrl_config, control(2,1,2=>Z), c.locs)
-
-abstract type AbstractDiff{GT, N, T} <: TagBlock{GT, N} end
-Base.adjoint(df::AbstractDiff) = Daggered(df)
-
-istraitkeeper(::AbstractDiff) = Val(true)
+generator(c::CphaseGate{N}) where N = ControlBlock{N}(c.ctrl_locs, c.ctrl_config, Z, c.locs)
 
 #################### The Basic Diff #################
 """
-    QDiff{GT, N} <: AbstractDiff{GT, N, T}
-    QDiff(block) -> QDiff
+    Diff{GT, N, T} <: TagBlock{GT, N}
+    Diff(block) -> Diff
 
 Mark a block as quantum differentiable.
 """
-mutable struct QDiff{GT, N, T} <: AbstractDiff{GT, N, T}
-    block::GT
-    grad::T
-    QDiff(block::DiffBlock{N, T}) where {N, T} = new{typeof(block), N, T}(block, T(0))
-end
-content(cb::QDiff) = cb.block
-chcontent(cb::QDiff, blk::DiffBlock) = QDiff(blk)
-
-@forward QDiff.block apply!
-mat(::Type{T}, df::QDiff) where T = mat(T, df.block)
-Base.adjoint(df::QDiff) = QDiff(content(df)')
-
-function YaoBlocks.print_annotation(io::IO, df::QDiff)
-    printstyled(io, "[̂∂] "; bold=true, color=:yellow)
-end
-
-#################### The Back Propagation Diff #################
-"""
-    BPDiff{GT, N, T, PT} <: AbstractDiff{GT, N, Complex{T}}
-    BPDiff(block, [grad]) -> BPDiff
-
-Mark a block as differentiable, here `GT`, `PT` is gate type, parameter type.
-
-Warning:
-    please don't use the `adjoint` after `BPDiff`! `adjoint` is reserved for special purpose! (back propagation)
-"""
-mutable struct BPDiff{GT, N, T} <: AbstractDiff{GT, N, T}
+mutable struct Diff{GT, N, T} <: TagBlock{GT, N}
     block::GT
     grad::T
-    input::AbstractRegister
-    BPDiff(block::AbstractBlock{N}, grad::T) where {N, T} = new{typeof(block), N, T}(block, grad)
+    Diff(block::DiffBlock{N, T}) where {N, T} = new{typeof(block), N, T}(block, T(0))
 end
-BPDiff(block::AbstractBlock) = BPDiff(block, zeros(parameters_eltype(block), nparameters(block)))
-BPDiff(block::DiffBlock{N, T}) where {N, T} = BPDiff(block, T(0))
+content(cb::Diff) = cb.block
+chcontent(cb::Diff, blk::DiffBlock) = Diff(blk)
 
-content(cb::BPDiff) = cb.block
-chcontent(cb::BPDiff, blk::AbstractBlock) = BPDiff(blk)
+istraitkeeper(::Diff) = Val(true)
 
-mat(::Type{T}, df::BPDiff) where T = mat(T, df.block)
-function apply!(reg::AbstractRegister, df::BPDiff)
-    if isdefined(df, :input)
-        copyto!(df.input, reg)
-    else
-        df.input = copy(reg)
-    end
-    apply!(reg, content(df))
-    reg
-end
-
-function apply!(δ::AbstractRegister, adf::Daggered{<:BPDiff{<:Rotor}})
-    df = adf |> content
-    apply!(δ, content(df)')
-    df.grad = -statevec(df.input |> generator(content(df)))' * statevec(δ) |> imag
-    δ
-end
+@forward Diff.block apply!
+mat(::Type{T}, df::Diff) where T = mat(T, df.block)
+Base.adjoint(df::Diff) = chcontent(df, content(df)')
 
-function YaoBlocks.print_annotation(io::IO, df::BPDiff)
+function YaoBlocks.print_annotation(io::IO, df::Diff)
     printstyled(io, "[∂] "; bold=true, color=:yellow)
 end
 
-
 #### interface #####
 export autodiff, numdiff, opdiff, StatFunctional, statdiff, as_weights
 
@@ -103,7 +55,7 @@ autodiff(mode::Symbol) = block->autodiff(mode, block)
 autodiff(mode::Symbol, block::AbstractBlock) = autodiff(Val(mode), block)
 
 # for BP
-autodiff(::Val{:BP}, block::DiffBlock) = BPDiff(block)
+autodiff(::Val{:BP}, block::DiffBlock) = Diff(block)
 autodiff(::Val{:BP}, block::AbstractBlock) = block
 # Sequential, Roller and ChainBlock can propagate.
 function autodiff(mode::Val{:BP}, blk::Union{ChainBlock, Sequential})
@@ -111,16 +63,16 @@ function autodiff(mode::Val{:BP}, blk::Union{ChainBlock, Sequential})
 end
 
 # for QC
-autodiff(::Val{:QC}, block::Union{RotationGate, CphaseGate}) = QDiff(block)
+autodiff(::Val{:QC}, block::Union{RotationGate, CphaseGate}) = Diff(block)
 # escape control blocks.
 autodiff(::Val{:QC}, block::ControlBlock) = block
 
 function autodiff(mode::Val{:QC}, blk::AbstractBlock)
     blks = subblocks(blk)
     isempty(blks) ? blk : chsubblocks(blk, autodiff.(mode, blks))
- end
+end
 
-@inline function _perturb(func, gate::AbstractDiff{<:DiffBlock}, δ::Real)
+@inline function _perturb(func, gate::Diff{<:DiffBlock}, δ::Real)
     dispatch!(-, gate, (δ,))
     r1 = func()
     dispatch!(+, gate, (2δ,))
@@ -129,7 +81,7 @@ function autodiff(mode::Val{:QC}, blk::AbstractBlock)
     r1, r2
 end
 
-@inline function _perturb(func, gate::AbstractDiff{<:Rotor}, δ::Real)  # for put
+@inline function _perturb(func, gate::Diff{<:Rotor}, δ::Real)  # for put
     dispatch!(-, gate, (δ,))
     r1 = func()
     dispatch!(+, gate, (2δ,))
@@ -139,21 +91,21 @@ end
 end
 
 """
-    numdiff(loss, diffblock::AbstractDiff; δ::Real=1e-2)
+    numdiff(loss, diffblock::Diff; δ::Real=1e-2)
 
 Numeric differentiation.
 """
-@inline function numdiff(loss, diffblock::AbstractDiff; δ::Real=1e-2)
+@inline function numdiff(loss, diffblock::Diff; δ::Real=1e-2)
     r1, r2 = _perturb(loss, diffblock, δ)
     diffblock.grad = (r2 - r1)/2δ
 end
 
 """
-    opdiff(psifunc, diffblock::AbstractDiff, op::AbstractBlock)
+    opdiff(psifunc, diffblock::Diff, op::AbstractBlock)
 
 Operator differentiation.
 """
-@inline function opdiff(psifunc, diffblock::AbstractDiff, op::AbstractBlock)
+@inline function opdiff(psifunc, diffblock::Diff, op::AbstractBlock)
     r1, r2 = _perturb(()->expect(op, psifunc()) |> real, diffblock, π/2)
     diffblock.grad = (r2 - r1)/2
 end
@@ -201,45 +153,68 @@ expect(stat::StatFunctional{1, <:Function}, xs::AbstractVector) = mean(stat.data
 Base.ndims(stat::StatFunctional{N}) where N = N
 
 """
-    statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{<:Any, <:AbstractArray}; initial::AbstractVector=probfunc())
-    statdiff(samplefunc, diffblock::AbstractDiff, stat::StatFunctional{<:Any, <:Function}; initial::AbstractVector=samplefunc())
+    statdiff(probfunc, diffblock::Diff, stat::StatFunctional{<:Any, <:AbstractArray}; initial::AbstractVector=probfunc())
+    statdiff(samplefunc, diffblock::Diff, stat::StatFunctional{<:Any, <:Function}; initial::AbstractVector=samplefunc())
 
 Differentiation for statistic functionals.
 """
-@inline function statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{2}; initial::AbstractVector=probfunc())
+@inline function statdiff(probfunc, diffblock::Diff, stat::StatFunctional{2}; initial::AbstractVector=probfunc())
     r1, r2 = _perturb(()->expect(stat, probfunc(), initial), diffblock, π/2)
     diffblock.grad = (r2 - r1)*ndims(stat)/2
 end
-@inline function statdiff(probfunc, diffblock::AbstractDiff, stat::StatFunctional{1})
+@inline function statdiff(probfunc, diffblock::Diff, stat::StatFunctional{1})
     r1, r2 = _perturb(()->expect(stat, probfunc()), diffblock, π/2)
     diffblock.grad = (r2 - r1)*ndims(stat)/2
 end
 
 """
-    backward!(δ::AbstractRegister, circuit::AbstractBlock) -> AbstractRegister
+    backward!(state, circuit::AbstractBlock) -> AbstractRegister
 
 back propagate and calculate the gradient ∂f/∂θ = 2*Re(∂f/∂ψ*⋅∂ψ*/∂θ), given ∂f/∂ψ*.
+`state` is a pair of output_register => the corresponding adjoint.
 
 Note:
 Here, the input circuit should be a matrix block, otherwise the back propagate may not apply (like Measure operations).
 """
-backward!(δ::AbstractRegister, circuit::AbstractBlock) = apply!(δ, circuit')
+function backward!(state, block::AbstractBlock)
+    out, outδ = state
+    adjblock = block'
+    backward_params!((out, outδ), block)
+    in = apply!(out, adjblock)
+    inδ = apply!(outδ, adjblock)
+    return (in, inδ)
+end
+
+function backward!(state, circuit::Union{ChainBlock, Concentrator})
+    for blk in Base.Iterators.reverse(subblocks(circuit))
+        state = backward!(state, blk)
+    end
+    return state
+end
+
+backward!(state, block::Measure) = throw(MethodError(backward!, (state, block)))
+
+backward_params!(state, block::AbstractBlock) = nothing
+function backward_params!(state, block::Diff{<:DiffBlock})
+    in, outδ = state
+    Σ = generator(content(block))
+    block.grad = -statevec(in |> Σ)' * statevec(outδ) |> imag
+    in |> Σ
+    nothing
+end
 
 """
     gradient(circuit::AbstractBlock, mode::Symbol=:ANY) -> Vector
 
-collect all gradients in a circuit, mode can be :BP/:QC/:ANY, they will collect `grad` from BPDiff/QDiff/AbstractDiff respectively.
+collect all gradients in a circuit, mode can be :BP/:QC/:ANY, they will collect `grad` from Diff respectively.
 """
-gradient(circuit::AbstractBlock, mode::Symbol=:ANY) = gradient!(circuit, parameters_eltype(circuit)[], mode)
+gradient(circuit::AbstractBlock) = gradient!(circuit, parameters_eltype(circuit)[])
 
-gradient!(circuit::AbstractBlock, grad, mode::Symbol) = gradient!(circuit, grad, Val(mode))
-function gradient!(circuit::AbstractBlock, grad, mode::Val)
+function gradient!(circuit::AbstractBlock, grad)
     for block in subblocks(circuit)
-        gradient!(block, grad, mode)
+        gradient!(block, grad)
     end
     grad
 end
 
-gradient!(circuit::BPDiff, grad, mode::Val{:BP}) = append!(grad, circuit.grad)
-gradient!(circuit::QDiff, grad, mode::Val{:QC}) = push!(grad, circuit.grad)
-gradient!(circuit::AbstractDiff, grad, mode::Val{:ANY}) = append!(grad, circuit.grad)
+gradient!(circuit::Diff, grad) = append!(grad, circuit.grad)
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<<N
+        new{N}(a, L)
+    end
+end
+
+function Yao.apply!(reg::ArrayReg{B}, m::Mod{N}) where {B, N}
+    YaoBlocks._check_size(reg, m)
+    nstate = zero(reg.state)
+    for i in basis(reg)
+        _i = i >= 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<<N)
+    for i in basis(N)
+        @inbounds perm[i >= m.L ? i+1 : mod(i*m.a, m.L)+1] = i+1
+    end
+    PermMatrix(perm, ones(T, 1<<N))
+end
+
+Base.adjoint(m::Mod{N}) where N = Mod{N}(NumberTheory.mod_inverse(m.a, m.L), m.L)
+Yao.print_block(io::IO, m::Mod{N}) where N = print(io, "Mod{$N}: $(m.a)*x % $(m.L)")
+
+Yao.isunitary(::Mod) = true
+# Yao.ishermitian(::Mod) = true  # this is not true for L = 1.
+# Yao.isreflexive(::Mod) = true  # this is not true for L = 1.
+
+"""
+    KMod{N, K} <: PrimitiveBlock{N}
+
+The first `K` qubits are exponent `k`, and the rest `N-K` are base `a`,
+it calculates `mod(a^k*x, L)`, notice `gcd(a, L)` should be 1.
+"""
+struct KMod{N, K} <: PrimitiveBlock{N}
+    a::Int
+    L::Int
+    function KMod{N, K}(a, L) where {N, K}
+        @assert gcd(a, L) == 1 && L<=1<<(N-K)
+        new{N, K}(a, L)
+    end
+end
+
+nqubits_data(m::KMod{N, K}) where {N, K} = N-K
+nqubits_control(m::KMod{N, K}) where {N, K} = K
+
+function bint2_reader(T, k::Int)
+    mask = bmask(T, 1:k)
+    return b -> (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<<K + 1
+        for j in 1:B
+            @inbounds nstate[_b,j] = reg.state[b+1,j]
+        end
+    end
+    reg.state = nstate
+    reg
+end
+
+function Yao.mat(::Type{T}, m::KMod{N, K}) where {T, N, K}
+    perm = Vector{Int}(undef, 1<<N)
+    reader = bint2_reader(Int, K)
+    for b in basis(N)
+        k, i = reader(b)
+        _i = i >= m.L ? i : mod(i*powermod(m.a, k, m.L), m.L)
+        _b = k + _i<<K + 1
+        @inbounds perm[_b] = b+1
+    end
+    PermMatrix(perm, ones(T, 1<<N))
+end
+
+Base.adjoint(m::KMod{N, K}) where {N, K} = KMod{N, K}(NumberTheory.mod_inverse(m.a, m.L), m.L)
+Yao.print_block(io::IO, m::KMod{N, K}) where {N, K} = print(io, "Mod{$N, $K}: $(m.a)^k*x % $(m.L)")
+
+Yao.isunitary(::KMod) = true
+# Yao.ishermitian(::Mod) = true  # this is not true for L = 1.
+# Yao.isreflexive(::Mod) = true  # this is not true for L = 1.
diff --git a/src/QCBM.jl b/src/QCBM.jl
index 9e4b694..51bd947 100644
--- a/src/QCBM.jl
+++ b/src/QCBM.jl
@@ -12,7 +12,7 @@ struct QCBM{BT<:AbstractBlock, KT<:AbstractKernel} <: QCOptProblem
     dbs
 end
 function QCBM(circuit::AbstractBlock, kernel::AbstractKernel, ptrain::Vector)
-    QCBM(circuit, kernel, ptrain, collect_blocks(AbstractDiff, circuit))
+    QCBM(circuit, kernel, ptrain, collect_blocks(Diff, circuit))
 end
 
 # INTERFACES
@@ -40,12 +40,12 @@ psi(qcbm) = zero_state(qcbm.circuit |> nqubits) |> qcbm.circuit
 """generated probability distribution"""
 probs(qcbm::QCBM) = qcbm |> psi |> probs
 """
-    mmdgrad(qcbm::QCBM, db::AbstractDiff; p0::Vector) -> Float64
+    mmdgrad(qcbm::QCBM, db::Diff; p0::Vector) -> Float64
 
 gradient of MMD two sample test loss, `db` must be contained in qcbm.
 `p0` is current probability distribution.
 """
-function mmdgrad(qcbm::QCBM, db::AbstractDiff; p0::Vector)
+function mmdgrad(qcbm::QCBM, db::Diff; p0::Vector)
     statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)), initial=p0 |> as_weights) -
         2*statdiff(()->probs(qcbm) |> as_weights, db, StatFunctional(kmat(qcbm.kernel)*qcbm.ptrain))
 end
diff --git a/src/QCOptProblem.jl b/src/QCOptProblem.jl
index 9e85749..ba25c1a 100644
--- a/src/QCOptProblem.jl
+++ b/src/QCOptProblem.jl
@@ -35,7 +35,7 @@ function gradient end
 
 collection of all differentiable units.
 """
-diff_blocks(qop::QCOptProblem) = collect_blocks(AbstractDiff, qop |> circuit)
+diff_blocks(qop::QCOptProblem) = collect_blocks(Diff, qop |> circuit)
 
 """
     num_gradient(qop::QCOptProblem) -> Vector
diff --git a/src/QSVD.jl b/src/QSVD.jl
index ab6b0b3..1a8d828 100644
--- a/src/QSVD.jl
+++ b/src/QSVD.jl
@@ -19,7 +19,7 @@ function train_qsvd!(reg, circuit_a::AbstractBlock{Na}, circuit_b::AbstractBlock
     obs = -mapreduce(i->put(nbit, i=>Z), +, (1:Na..., Na+Nc+1:Na+Nb...))
     params = parameters(c)
     for i = 1:maxiter
-        grad = opdiff.(() -> copy(reg) |> c, collect_blocks(AbstractDiff, c), Ref(obs))
+        grad = opdiff.(() -> copy(reg) |> c, collect_blocks(Diff, c), Ref(obs))
         QuAlgorithmZoo.update!(params, grad, optimizer)
         println("Iter $i, Loss = $(Na+expect(obs, copy(reg) |> c))")
         dispatch!(c, params)
diff --git a/src/QuAlgorithmZoo.jl b/src/QuAlgorithmZoo.jl
index 6726b42..5f04589 100644
--- a/src/QuAlgorithmZoo.jl
+++ b/src/QuAlgorithmZoo.jl
@@ -32,6 +32,7 @@ include("hamiltonian_solvers.jl")
 include("HadamardTest.jl")
 include("lin_diffEq_HHL.jl")
 include("QSVD.jl")
+include("shor.jl")
 
 
 end # module
diff --git a/src/QuGAN.jl b/src/QuGAN.jl
index 70b5917..26d3049 100644
--- a/src/QuGAN.jl
+++ b/src/QuGAN.jl
@@ -23,8 +23,8 @@ struct QuGAN{N} <: QCOptProblem
         N = nqubits(target)
         c = Sequence([gen, addbits!(1), dis])
         witness_op = put(N+1, (N+1)=>P0)
-        gdiffs = collect_blocks(AbstractDiff, gen)
-        ddiffs = collect_blocks(AbstractDiff, dis)
+        gdiffs = collect_blocks(Diff, gen)
+        ddiffs = collect_blocks(Diff, dis)
         new{N}(target, gen, dis, zero_state(N), witness_op, c, gdiffs, ddiffs)
     end
 end
diff --git a/src/hamiltonian_solvers.jl b/src/hamiltonian_solvers.jl
index 4fee6f2..d344692 100644
--- a/src/hamiltonian_solvers.jl
+++ b/src/hamiltonian_solvers.jl
@@ -51,7 +51,7 @@ variational quantum eigensolver, faithful simulation with optimizer Adam(lr=0.01
 """
 function vqe_solve!(circuit::AbstractBlock{N}, hamiltonian::AbstractBlock; niter::Int=100) where N
     optimizer = Adam(lr=0.01)
-    dbs = collect_blocks(AbstractDiff, circuit)
+    dbs = collect_blocks(Diff, circuit)
     params = parameters(circuit)
     for i = 1:niter
         grad = opdiff.(()->zero_state(N) |> circuit, dbs, Ref(hamiltonian))
diff --git a/src/number_theory.jl b/src/number_theory.jl
new file mode 100644
index 0000000..355e281
--- /dev/null
+++ b/src/number_theory.jl
@@ -0,0 +1,101 @@
+export NumberTheory
+
+module NumberTheory
+
+export Z_star, Eulerφ, continued_fraction, mod_inverse, rand_primeto, factor_a_power_b
+export is_order, order_from_float, find_order
+
+"""
+    Z_star(N::Int) -> 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..c64fe89
--- /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 = NumberTheory.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/CircuitBuild.jl b/test/CircuitBuild.jl
index d64687c..4e078e2 100644
--- a/test/CircuitBuild.jl
+++ b/test/CircuitBuild.jl
@@ -46,14 +46,14 @@ end
     reg = rand_state(4)
     dispatch!(c, randn(nparameters(c)))
 
-    dbs = collect_blocks(BPDiff, c)
+    dbs = collect_blocks(Diff, c)
     op = kron(4, 1=>Z, 2=>X)
     loss1z() = expect(op, copy(reg) |> c)  # return loss please
 
     # back propagation
     ψ = copy(reg) |> c
     δ = copy(ψ) |> op
-    backward!(δ, c)
+    backward!((ψ, δ), c)
     bd = gradient(c)
 
     # get num gradient
diff --git a/test/Diff.jl b/test/Diff.jl
index 29dc559..164fcbb 100644
--- a/test/Diff.jl
+++ b/test/Diff.jl
@@ -5,22 +5,23 @@ using LinearAlgebra, Test, Random
 @testset "BP diff" begin
     reg = rand_state(4)
     block = put(4, 2=>rot(X, 0.3))
-    df = BPDiff(block)
+    df = Diff(block)
     @test df.grad == 0
     @test nqubits(df) == 4
 
-    df2 = BPDiff(rot(CNOT, 0.3))
+    df2 = Diff(rot(CNOT, 0.3))
     @test df2.grad == 0
     @test nqubits(df2) == 2
+    @test_throws MethodError backward!((reg, reg), Measure(4))
 end
 
 @testset "Qi diff" begin
     reg = rand_state(4)
-    df2 = QDiff(rot(CNOT, 0.3))
+    df2 = Diff(rot(CNOT, 0.3))
     @test df2.grad == 0
     @test nqubits(df2) == 2
 
-    @test df2' isa QDiff
+    @test df2' isa Diff
     @test mat(df2) == mat(df2')'
 end
 
@@ -76,19 +77,20 @@ end
 
     circuit = chain(4, repeat(4, H, 1:4), put(4, 3=>Rz(0.5)) |> autodiff(:BP), control(4, 2, 1=>shift(0.4)) |> autodiff(:BP), control(2, 1=>X), put(4, 4=>Ry(0.2)) |> autodiff(:BP))
     op = put(4, 3=>Y)
-    θ = [0.1, 0.2, 0.3]
+    θ = [0.9, 0.2, 0.3]
     dispatch!(circuit, θ)
     loss! = loss_expect!(circuit, op)
     ψ0 = rand_state(4)
     ψ = copy(ψ0) |> circuit
 
     # get gradient
-    δ = ψ |> op
-    backward!(δ, circuit)
+    δ = copy(ψ) |> op
+    in, inδ = backward!((ψ, δ), circuit)
+    @test in ≈ ψ0
     g1 = gradient(circuit)
 
     g2 = zero(θ)
-    η = 0.01
+    η = 1e-5
     for i in 1:length(θ)
         θ1 = copy(θ)
         θ2 = copy(θ)
@@ -96,7 +98,7 @@ end
         θ2[i] += 0.5η
         g2[i] = (loss!(copy(ψ0), θ2) - loss!(copy(ψ0), θ1))/η |> real
     end
-    g3 = opdiff.(() -> copy(ψ0) |> circuit, collect_blocks(BPDiff, circuit), Ref(op))
+    g3 = opdiff.(() -> copy(ψ0) |> circuit, collect_blocks(Diff, circuit), Ref(op))
     @test isapprox.(g1, g2, atol=1e-5) |> all
     @test isapprox.(g2, g3, atol=1e-5) |> all
 end
@@ -106,8 +108,8 @@ end
     @test generator(Rx(0.1)) == X
     circuit = chain(put(4, 1=>Rx(0.1)), control(4, 2, 1=>Ry(0.3)))
     c2 = circuit |> autodiff(:BP)
-    @test c2[1] isa BPDiff
-    @test !(c2[2] isa BPDiff)
+    @test c2[1] isa Diff
+    @test !(c2[2] isa Diff)
 end
 
 @testset "numdiff & opdiff" begin
@@ -125,13 +127,12 @@ end
 
     reg = rand_state(4)
     c = chain(put(4, 1=>Rx(0.5)), control(4, 1, 2=>Ry(0.5)), control(4, 1, 2=>shift(0.3)),  kron(4, 2=>Rz(0.3), 3=>Rx(0.7))) |> autodiff(:QC)
-    dbs = collect_blocks(QDiff, c)
+    dbs = collect_blocks(Diff, c)
     loss1z() = expect(kron(4, 1=>Z, 2=>X), copy(reg) |> c) |> real  # return loss please
     nd = numdiff.(loss1z, dbs)
     ed = opdiff.(()->copy(reg) |> c, dbs, Ref(kron(4, 1=>Z, 2=>X)))
     gd = gradient(c)
-    @test gradient(c, :QC) == gd
-    @test gradient(c, :BP) == []
+    @test gradient(c) == gd
     @test isapprox(nd, ed, atol=1e-4)
     @test ed == gd
 end
@@ -146,7 +147,7 @@ end
     prs = [1=>2, 2=>3, 3=>1]
     c = ibm_diff_circuit(nbit, 2, prs) |> autodiff(:QC)
     dispatch!(c, :random)
-    dbs = collect_blocks(AbstractDiff,c)
+    dbs = collect_blocks(Diff,c)
 
     p0 = zero_state(nbit) |> c |> probs
     sample0 = measure(zero_state(nbit) |> c; nshots=5000)
@@ -162,7 +163,7 @@ end
     V = StatFunctional(h)
     c = ibm_diff_circuit(nbit, 2, prs) |> autodiff(:QC)
     dispatch!(c, :random)
-    dbs = collect_blocks(AbstractDiff, c)
+    dbs = collect_blocks(Diff, c)
 
     p0 = zero_state(nbit) |> c |> probs |> as_weights
     loss0 = expect(V, p0 |> as_weights)
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
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