code reading HEAAN implementation in Julia

HEAAN implementation in Julia

HEAAN の実装を読む

Overview

test/api.test.jl > @testcase “encrypt()/decrypt()” 簡素化

• log_polynomial_length(n): 多項式の長さ, 長いほど安全だが遅い
• log_lo_modulus(q_L): 暗号文の長さ, 長いほどたくさん計算できるが遅い

    # 0. Params
    log_precision = 30
    log_cap = 100
    params = Params(log_polynomial_length=8, log_lo_modulus=300)
    # 1. KeyGen
    sk = SecretKey()
    pk = EncryptionKey(sk)
 
    # 2. Encryption
    z1 = [complex(0.3)]
    z2 = [complex(0.4)]
 
    c1, c2 = encrypt(pk, z1), encrypt(pk, z2)
 
    # 3. Operation
    c3 = add(c1, c2)
 
    # 4. Decryption
    z3 = decrypt(sk, c3)
    @assert is_approx(z1 .+ z2, c3, 23)

1. KeyGen

# 乱数生成器: メルセンヌ・ツイスタ
rng = MersenneTwister(123)
 
params = Params(log_polynomial_length=8, log_lo_modulus=300)
secret_key = SecretKey(rng, params)
enc_key = EncryptionKey(rng, secret_key)

Params(log_polynomial_length, log_lo_modulus)

HEAANスキームのパラメータ

src/params.jl(簡素化)

struct Params
    # in HEAAN: N
    log_polynomial_length::Int
    # in HEAAN: logQ, q_L in the paper
    log_lo_modulus::Int
    # in HEAAN: == logQ, P in the paper
    log_hi_modulus::Int
    # in HEAAN: sigma
    gaussian_noise_stddev::Float64
    # in HEAAN: h
    secret_key_length::Int
end

• gaussian_noise_stddev = 3.2
• secret_key_length = 64
• log_hi_modulus = log_lo_modulus(=300) + 50

SecretKey(params)

• 秘密鍵は 乱数 $Uni(\mathrm{0..}{2}^N)$から 一様にインデックスを secret_key_length = 64個を選びランダムビット
• [1 => true, 3 => false, 7 => true ] のようにスパースに保持

src/secret_key.jl(抜粋)

struct SecretKey
    params :: Params
    # [(pow, is -1?)]
    nonzero_entries :: Array{Pair{Int, Bool}, 1}
 
    function SecretKey(rng::AbstractRNG, params::Params)
        # 0 ~ 2^params.log_polynomial_length-1 の一様分布
        positions = sample_unique(
            rng, 0, 2^params.log_polynomial_length-1, params.secret_key_length)
        bits = rand(rng, Bool, params.secret_key_length)
        new(params, [pos => bit for (pos, bit) in zip(positions, bits)])
    end
end

EncryptionKey(sk)

src/public_keys.jl(抜粋)

struct EncryptionKey
    params :: Params
    key :: PublicKeyRNS
 
    function EncryptionKey(rng::AbstractRNG, secret_key::SecretKey)
        params = secret_key.params
        log_plen = params.log_polynomial_length
        plen = 2^log_plen
        log_modulus = params.log_lo_modulus + params.log_hi_modulus
 
        # rand_big_int(rng, 250, 256)
        # 0:2^log_modulusの乱数 を plen個 -> poly
        ax::Polynomial = rand_big_int(rng, log_modulus, plen)
        # ガウスノイズをかける
        bx::Polynomial = discrete_gaussian(rng, params.gaussian_noise_stddev, plen) - secret_key * ax
 
        plan = rns_plan(params)
        rax = to_rns_transformed(plan, ax, params.log_lo_modulus)
        rbx = to_rns_transformed(plan, bx, params.log_lo_modulus)
 
        new(params, PublicKeyRNS(rax, rbx))
    end
end

NTT(数論変換)

• Number Theoric Transform

NTTは，特殊な法を用いたFFT（高速フーリエ変換）である[13]．FFTとは，DFT（離散フーリエ変換）を高速に計算するアルゴリズムである．多項式環上での乗算，すなわち多項式同士の乗算には，次数をNとしてO(N2)の時間計算量を要する．多項式同士の乗算は，多項式の係数列同士の畳み込みとみなせる．係数列畳み込みの時間計算量も，多項式同士の乗算同様に$O(N^2)$である．一方，多項式の係数列をNTT形式に変換した場合，多項式同士の乗算（係数列の畳み込み）を$O(N\log N)$で得ることができる．RNS型CKKS方式においても，NTTを用いて乗算の計算コストを抑えている[7]．$R_q$上の多項式の係数ベクトル$a\in {\mathbb{Z}}_q^n$のNTT形式を$\tilde{a}$と表記する．$a$ を $q\equiv 1\mathrm{mod}2N$ を満たす素数$q$を法としてNTT形式$\tilde{a}$へ変換する関数を$\mathrm{NTT}(\mathrm{a};\mathrm{q})$と定義する．また，$\mathrm{NTT}(\mathrm{a};\mathrm{q})$の逆変換を$\mathrm{INTT}(\tilde{\mathrm{a}};\mathrm{q})$と定義する．

CKKS方式準同型暗号におけるRescale演算のGPU実装と演算性能評価 より

to_rns_transformed(…)

RNS(剰余数系)

有限個の整数からなる順序集合 $B=p_0,\cdots ,p_{k-1}$で，要素同士が互いに素である場合，$B$を基底と呼ぶ．
さらに，$P={\prod }_{i=0}^{k-1}{p}_i$ とおく．
$[\cdot {]}_B:{\mathbb{Z}}_p\mapsto {\prod }_{i=0}^{k-1}{Z}_{p_i}$ を，$a\mapsto [a{]}_B=([a{]}_{p_i}{)}_{0\le i<k}$ と置いた時，$[a{]}_B$を，$a\in {\mathbb{Z}}_P$ のRNS表現と呼ぶ．つまり，RNSを用いることで，$Z{\mathbb{Z}}_p$上の演算を，よりサイズの小さい $Z_{p_i}$に分割して計算を行うことができる．

互いに素な素数で割ったあまりを並べた表現

primes = [5, 3, 2] の時
29 = [29 % 5, 29 % 3, 29 % 2] = [4, 2, 1]

to_rns_transformed

function to_rns_transformed(
    plan::RNSPlan, 
    x::Polynomial{BinModuloInt{T, Q}}, 
    add_range::Int=0
) where {T, Q}
    # log_rangeを計算
    log_range = _log_modulus_mul(Q, add_range, length(x.coeffs))

    _ntt_forward(_to_rns(plan, x, log_range))
end

function _to_rns(
    plan::RNSPlan, 
    x::Polynomial{BinModuloInt{T, Q}}, 
    log_range::Int
) where {T, Q}
    # log_modulusより大きい最初の plan.max_bin_moduli
    np = min_nprimes(plan, log_range)
    # plan.max_bin_moduli[nprimes]
    log_range = max_log_modulus(plan, np)
 
    plen = length(x.coeffs)
    res = Array{UInt64}(undef, plen, np)
    for i in 1:plen
        res[i,:] .= to_rns_signed(plan, x.coeffs[i], np)
    end
    RNSPolynomial(plan, res, Q, log_range, x.modulus)
end

• to_tns_signed() は to_tns() をしてから符号逆転

2. Encode & Encryotion

encrypt(pk, z)

src/ciphertext.jl(抜粋)

function encrypt(
    pk::EncryptionKey, 
    z::Array{Complex{Float64}, 1},
    log_precision::Int, log_cap::Int
)
    m = encode(
        pk, z, 
        log_precision + params.log_hi_modulus,
        log_cap + params.log_hi_modulus
    )
    encrypt(pk, m)
end

encode(z)

$N/2$個の複素数の配列 $z$ を $\Delta$ でスケーリングしたあとパッキングを行う. パッキングした平文を $m$ とする

Array{Complex{Float64}, 1} -> PlainText

function encode(
    z::Array{Complex{Float64}, 1}, 
    log_precision::Int, 
    log_modulus::Int,
)
    slots = length(z)
    m = τ′(params, z, log_precision, log_modulus, minimize_range)
    Plaintext(params, m, log_modulus, log_precision, slots)
end

decode(m)

encodeの逆

PlainText -> Array{Complex{Float64}, 1}

function decode(m::Plaintext)
    z = τ(m)
    embed(m.params, z)
end

encrypt(pk, m)

function encrypt(rng::AbstractRNG, pk::EncryptionKey, m::Plaintext)
    params = plain.params
    plen = 2 ^ params.log_polynomial_length
    log_modulus = plain.log_cap
 
    # random bool polyonmial
    vx = sample_ZO(plen)
 
    # vx --[to_rns]-> vx_rns
    plan = rns_plan(params)
    vx_rns = to_rns_transformed(plan, vx, params.log_lo_modulus + params.log_hi_modulus)
 
    # vx_rns * rax
    rax = from_rns_transformed(vx_rns * key.key.rax, log_modulus)
    noise = discrete_gaussian(rng, params.gaussian_noise_stddev, plen)
    # ax = rax: Poly + noise: Poly
    ax = (rax + noise)
 
    # bx = plain: Poly + (vx_rns: RNS * rax: Poly) + noise: Poly
    # vx_rns * rbx
    rbx = from_rns_transformed(vx_rns * key.key.rbx, log_modulus)
    noise = discrete_gaussian(rng, params.gaussian_noise_stddev, plen)
    # bx = rbx: Poly + plain: Poly + noise: Poly
    bx = (plain.polynomial + rbx + noise)
 
    # ax, bx: Polynomial{Z_450, 256}
    # 450 = log_lo + log_hi
 
    # ax.map(i -> right_shift_rounded(i, params.log_hi_modulus))
    ax = broadcast_into_polynomial(right_shift_rounded, ax, params.log_hi_modulus)
    bx = broadcast_into_polynomial(right_shift_rounded, bx, params.log_hi_modulus)
 
    Ciphertext(ax, bx)
end

3. Operation

add(c1, c2)

function add(c1::Ciphertext, c2::Ciphertext)
    Ciphertext(c1.ax + c2.ax, c1.bx + c2.bx)
end

• 同型に演算ができる
• log_cap, log_precisionは変化なし

mul(mk, c1, c2)

function mul(
    key::MultiplicationKey, 
    c1::Ciphertext, 
    c2::Ciphertext
)
    params = cipher1.params
    plan = rns_plan(params)
 
    log_cap = cipher1.log_cap
 
    ra1 = to_rns_transformed(plan, cipher1.ax, log_cap)
    rb1 = to_rns_transformed(plan, cipher1.bx, log_cap)
    ra2 = to_rns_transformed(plan, cipher2.ax, log_cap)
    rb2 = to_rns_transformed(plan, cipher2.bx, log_cap)
 
    axax = from_rns_transformed(ra1 * ra2, log_cap)
    bxbx = from_rns_transformed(rb1 * rb2, log_cap)
 
    ra1 = ra1 + rb1
    ra2 = ra2 + rb2
 
    axbx = from_rns_transformed(ra1 * ra2, log_cap)
 
    key = key.key
 
    raa = to_rns_transformed(plan, axax, params.log_lo_modulus + params.log_hi_modulus)
 
    #=
    In the paper, we want to find `(P^(-1) * x * y) mod q`,
    where `P = 2^log_hi_modulus`, `Q = 2^log_lo_modulus`, `x` is a polynomial modulo `q <= Q`,
    and `y` is a polynomial modulo `P * Q`.
 
    So, we perform the multiplication using the full range `q * Q * P`
    (times the polynomial length to prevent overflow during NTT in RNS),
    then drop the high `Q` bits (during the conversion from RNS to big integer),
    then shift by `P`.
    =#
 
    log_modulus = cipher1.log_cap + params.log_hi_modulus
    ax = broadcast_into_polynomial(
        right_shift_rounded,
        from_rns_transformed(raa * key.rax, log_modulus), params.log_hi_modulus)
    bx = broadcast_into_polynomial(
        right_shift_rounded,
        from_rns_transformed(raa * key.rbx, log_modulus), params.log_hi_modulus)
 
    ax = ax + axbx - bxbx - axax
    bx = bx + bxbx
 
    Ciphertext(
        params,
        ax,
        bx,
        cipher1.log_cap,
        cipher1.log_precision + cipher2.log_precision,
        cipher1.slots)
end

3.1 Decryption Structure of HEAAN

既存のHE方式の多くは、ZtやZt[X]/(ΦM(X))のようなモジュロ空間上で演算を行い、同相計算後の再帰メッセージのLSBを暗号化した暗号文を計算することを目的としている。 例えば、BGV型暗号[5,25,33]の場合、平文は暗号文のモジュラスの最下位ビットに配置され、つまり、暗号化されたメッセージの秘密情報khに対する暗号化情報は、〈c,sk〉=m+te(modq)の形式の復号化構造を持つ。ｍ1,m2の暗号化の乗算は、ｍ1m2の一部のＬＳＢ（すなわち、[ｍ1m2]t）を保持し、そのＭＳＢ（すなわち、ｂｍ1m2/te）は誤りによって破壊される。一方、FV型方式[4,22,2]では、暗号文の最上位ビットにメッセージを配置し、その復号構造が〈c,sk〉=bq/te-m+e(modq)を満たすようにしています。しかし、結果として得られるメッセージのＭＳＢは、〈ci,sk〉=q-Ii+bq/te-m+eiの同型乗算の間にも破壊され、それぞれがメッセージの左位置に追加のエラーＩｉｉを含む。

我々の目標は、暗号化されたデータに対して近似演算を行うこと、あるいは等価的に、同型演算の後に得られるメッセージのMSBを計算することである。主な考え方は、入力メッセージの重要な数字に続く暗号化ノイズを付加することである。より正確には、我々のスキームは、いくつかの小さなRoreに対して、〈c,sk〉=m+e(modq)の形式の復号化構造を持っている。この暗号化誤差は、スキームの安全性を保証するために挿入しているが、近似計算の際に発生する誤差と考える。つまり、復号化アルゴリズムの出力は、高精度で元のメッセージの近似値として扱われます。平文のサイズは、同型演算のための暗号文弾性率に比べて十分に小さく、算術演算の結果が暗号文弾性率よりも小さくなるようにする。

より慎重に検討しなければならない問題がいくつかあります。暗号化されていない近似計算では、連続して演算を行うと小さな誤差が爆発する可能性があるため、計算結果とアルゴリズムの正確な値との近さを考慮することは非常に重要です。また、メッセージのサイズの管理も課題の一つです。 メッセージの丸めを行わずに多倍深さLの回路を計算すると、出力値のビットサイズはLとともに指数関数的に大きくなります。このようなナイーブな方法は、暗号文のモジュラスが巨大化してしまうため、実用には適していません。この問題を解決するために、中間値を基底値で分割する新しい手法を提案します。これにより、7
メッセージのサイズをほぼ同じにして、必要な暗号文のモジュラスを深さLの中で線形にする方法です。この方法により、メッセージの大きさをほぼ同じに保ち、必要な暗号文のモジュラスを深さLで線形にすることができます。

3.3 Leveled Homomorphic Encryption Scheme for Approximate Arithmetic
この節の目的は、近似アリスメティックのためのレベル化されたHEスキームを構築することである。

$p>0$, モジュラス $q_0$, $q_{\ell }=p^{\ell }\cdot q_0$ for $0<l\le L$
セキュリティパラメータ $\lambda$については、$M=M(\lambda ,q_L)$ を選択します。
レベル $0\le \ell \le L$ の場合、レベル$l$の暗号文は固定整数 $k$ に対して $R_{q_{\ell }}^k$ のベクトルである。我々のスキームは、ノイズ推定のための定数BcleanとBmult(`)を持つ5つのアルゴリズム(KeyGen,Enc,Dec,Add,Mult)から構成されている。

• $\mathrm{KeyGen}(1^{\lambda })$: $\mathrm{sk}$ を作成

• ${\mathrm{Enc}}_{\mathrm{pk}}(\mathrm{m})$ : $<c,\mathrm{sk}>=\mathrm{m}+\mathrm{e}$

• ${\mathrm{Dec}}_{\mathrm{sk}}(\mathrm{c})$: レベル $\ell$ の $c$ から $m^{\prime }\leftarrow <C,\mathrm{sk}>(\mathrm{mod}{\mathrm{q}}_{\mathrm{\ell }})$

• $\mathrm{Add}({\mathrm{c}}_1,{\mathrm{c}}_2)$: $m_1,m_2$ を計算(ノイズは合計)

• ${\mathrm{Mul}}_{\mathrm{evk}}({\mathrm{c}}_1,{\mathrm{c}}_2)$: $<c_{\mathrm{mult}},\mathrm{sk}>=<{\mathrm{c}}_1,\mathrm{sk}>\cdot <{\mathrm{c}}_2,\mathrm{sk}>+{\mathrm{e}}_{\mathrm{mult}}(\mathrm{mod}{\mathrm{q}}_{\mathrm{\ell }})$

• ${\mathrm{RS}}_{\mathrm{\ell }\to {\mathrm{\ell }}^{\prime }}(\mathrm{c})$: $\ell$ から より低い ${\ell }^{\prime }$ への変換. $c^{\prime }\leftarrow \lfloor \frac{q_{\ell }^{\prime }}{q_{\ell }}c\rceil$

$\nu$: レベル $0\le \ell \le L$ 上限
$B\in \mathbb{R}$: ノイズの上限

各パラメータの推移

異なるレベル間の演算

• 低い方に合わせる $RS:m\to \frac{q_{\ell }^{\prime }}{q_{\ell }}m$
• 単純な modswitch を行う

異なるレベルの暗号文の同相演算.暗号文sc,c′ofm,m′が異なるレベルに属している場合、同相演算を行う前に、より大きなレベルの暗号文をより小さなレベルの暗号文に近づけるべきである。単純なモジュール化とRS手続きの2つの方法があります。この論文では，特に断りのない限り，異なるレベルの暗号文に対して計算を行う前に，単純なモジュラス削減を行い，より小さなモジュラスへと変換することにしている．

• $DG({\sigma }^2)$: 分散 ${\sigma }^2$ の離散ガウス分布: ${\mathbb{Z}}^N$
• $ZO(\rho )$: $\frac{\rho }{2}$ で1, $\frac{\rho }{2}$ で-1, $1-\rho$ で0
• ${\mathrm{Enc}}_{\mathrm{pk}}(\mathrm{m})$:
  $v\to \mathrm{ZO}(0.5)$
  $e_1,e_2\to \mathrm{DG}({\sigma }^2)$
  $\leftarrow v\cdot \mathrm{pk}+(\mathrm{m}+{\mathrm{e}}_0,{\mathrm{e}}_1)\mathrm{mod}{\mathrm{q}}_{\mathrm{L}}$
• ${\mathrm{Dec}}_{\mathrm{sk}}(\mathrm{c})$: for $c=(b,a),\leftarrow b+a\cdot s(\mathrm{mod}{\mathrm{q}}_{\mathrm{\ell }})$
• ${\mathrm{Mult}}_{\mathrm{evk}}({\mathrm{c}}_1,{\mathrm{c}}_2)$:
  $c_1=(b_1,a_1),c_2=(b_2,a_2)$
  $\mathrm{let}({\mathrm{d}}_0,{\mathrm{d}}_1,{\mathrm{d}}_2)=({\mathrm{b}}_1{\mathrm{b}}_2,{\mathrm{a}}_1{\mathrm{b}}_2+{\mathrm{a}}_2{\mathrm{b}}_1,{\mathrm{a}}_1{\mathrm{a}}_2)(\mathrm{mod}{\mathrm{q}}_{\mathrm{\ell }})$
  $\leftarrow (d_0,d_1)+\lfloor P^{-1}\cdot d_2\cdot \mathrm{evk}\rceil (\mathrm{mod}{\mathrm{q}}_{\mathrm{\ell }})$
  ($P$ は log_hi_modulus: Int のこと)

ロジスティック関数の近似

$f(x)=\frac{1}{1+e^{-x}}$

s