falcon_parameterized.cry

module falcon_parameterized where

parameter
    type _k : #
    type _n : #
    type constraint (fin _k, _n == 2^^_k, _k > 0)
    type pkbytelen : #
    type sbytelen : #
    type constraint (fin sbytelen)
    type coeffs_number_of_bits : #
    sigma : Float64
    type Q = Rational
    type C = (Float64,Float64)
    Omega_phi: [_n]C
    type constraint (8 * sbytelen >= 328)


type slen = 8 * sbytelen - 328
type q = 3 * 2^^12 + 1

type Byte = [8]Bit
type ZZ = Integer

while : {a} (a -> Bit) -> (a -> a) -> a -> a
while condition body initial_state =
  if(condition initial_state) then while condition body (body initial_state)
  else initial_state


dowhile : {a} (a -> Bit) -> (a -> a) -> a -> a
dowhile condition body initial_state =
  if(condition next_state) then while condition body next_state else next_state
  where next_state = body initial_state

CMul : (C, C) -> C
CMul(x, y) = (x.0*y.0-x.1*y.1, x.0*y.1+x.1*y.0)

CAdd : (C, C) -> C 
CAdd(x, y) = (x.0 + x.1, y.0 + y.1)

CAddList : {n} (fin n) => [n]C -> C 
CAddList(l) = sums ! 0
    where sums = [zero] # [CAdd(el,sums') | el <- l 
                                          | sums' <- sums
                          ]

Cinv : C -> C
Cinv((a,b)) = (a/.denominator, -(b/.denominator)) where 
    denominator = a^^2 + b^^2

Conjugate : C -> C 
Conjugate((x, y)) = (x, -y)

PolyConj : {k, n} (ispoweroftwo k n) => Poly n C -> Poly n C
PolyConj(f) = map Conjugate f

Normalize_FFT : {k, n} (ispoweroftwo k n) => (Float64, Poly n C) -> Poly n C
Normalize_FFT(_q, F) = dot(F, [Cinv(_q,0) | i <- [0 .. (n-1)]])

type constraint ispoweroftwo a b = (fin a, b == 2^^a)

type constraint positive_power a = (fin a, a > 0, (2*(2^^(a-1)*a)+2^^a == 2^^a*(1+a)), 2 ^^ a == 2 * 2 ^^ (a - 1))

// helper function for recursive functions. Will be removed
// when Cryptol is updated.
resize : {m,n,a} (fin m, fin n, Zero a) => [m]a -> [n]a
resize xs = take`{n} (xs # repeat`{inf} zero)



FFT : {k, n} (ispoweroftwo k n) => Poly n C -> FFT n
FFT(x) | k == 0 => x
       | k >  0 => resize(join([result0,result1])) where
            even = [x@(2*i  ) | i <- [0 .. (2^^(k-1)-1)]]
            odd =  [x@(2*i+1) | i <- [0 .. (2^^(k-1)-1)]]
            left = FFT`{k-1}(even)
            right = FFT`{k-1}(odd)
            X = join([left,right])
            result0 = [X@i + CMul(Omega_phi@(`n*i),(X@(i+`n/2))) | i <- [0 .. (n/2-1)]]
            result1 = [X@i - CMul(Omega_phi@(`n*i),(X@(i+`n/2))) | i <- [0 .. (n/2-1)]]

FFTInv  : {k, n} (ispoweroftwo k n) => (FFT n) -> (Poly n C)
FFTInv(x) = zero


FFT' : {k, n} (ispoweroftwo k n) => Poly n ZZ -> FFT n
FFT'(x) = FFT`{k}(map IntToCmplx x : (Poly n C))

FFT'' : {k, n} (ispoweroftwo k n) => Poly n (Z q) -> FFT n
FFT''(x) = FFT'`{k}(map fromZ x)

// Number Theoretic Transform
NTT : {k, n} (ispoweroftwo k n) => Poly n (Z q) -> Poly n (Z q)
NTT(f) = zero

NTTInv : {k, n} (ispoweroftwo k n) => Poly n (Z q) -> Poly n (Z q)
NTTInv(f) = zero

ModInv : Z q -> Z q
ModInv(f) =  1/. f

phi : {k, n} (ispoweroftwo k n) => Poly (n+1) (Z q)
phi = [1]#zero#[1]

IntToCmplx : ZZ -> C
IntToCmplx(x) = ((fromInteger x),zero)

dot : {n} (fin n) => (FFT n, FFT n) -> FFT n
dot(f,g) = [CMul(fi,gi) | fi <- f | gi <- g]

HadamardDivision : {n} (fin n) => (FFT n, FFT n) -> FFT n
HadamardDivision(f,g) = dot(f, (map Cinv g))

PolyMul : {k, n} (ispoweroftwo k n) => ((Poly n C), (Poly n C)) -> (Poly n C)
PolyMul(f,g) = FFTInv`{k}(dot(FFT`{k}(f), FFT`{k}(g)))

PolyInv : {k, n} (ispoweroftwo k n) => (Poly n C) -> (Poly n C)
PolyInv(f) = FFTInv`{k}(map Cinv (FFT`{k}(f)))

PolyMulInZ : {k, n} (ispoweroftwo k n) => ((Poly n ZZ), (Poly n ZZ)) -> (Poly n ZZ)
PolyMulInZ(f, g) = zero // TODO!

PolyDivInZ : {k, n} (ispoweroftwo k n) => ((Poly n ZZ), (Poly n ZZ)) -> (Poly n Q)
PolyDivInZ(f, g) = zero // TODO!

innerProduct : {k, n} (ispoweroftwo k n) => ((Poly n C), (Poly n C)) -> C 
innerProduct(f, g) = zero // TODO!

norm_sq : {k, n} (ispoweroftwo k n) => (Poly n C) -> Float64
norm_sq(f) = zero // TODO!

splitfft : {k, n} (ispoweroftwo k n, n > 2) => (FFT n) -> [2](FFT (n/2))
splitfft FFTf = [resize(FFTf0), resize(FFTf1)] where
    FFTf0 = [FFTf@(2*i  ) | i <- [0 .. (2^^(k-1)-1)]]
    FFTf1 = [FFTf@(2*i+1) | i <- [0 .. (2^^(k-1)-1)]]

mergefft : {k, n} (ispoweroftwo k n, n>=2) => [2](FFT (n/2)) -> (FFT n)
mergefft ([f0, f1]) = resize FFTf where
    FFTf = join[[f0@i,f1@i] | i <- [0 .. (2^^(k-1)-1)]]

type Poly degree a = [degree]a

type FFT degree = Poly degree C
type NTT degree = Poly degree (Z q)

// public key types
type publicKey = Poly _n (Z q)
type encodedPublicKey = [14 * _n /^ 8 + 1]Byte

// private key types
type privateKey = ([2][2](FFT _n), falconTree _k)
type encodedPrivateKey = [2*((_n+1)*coeffs_number_of_bits) + ((_n+1)*8)]Bit

type falconTree k = [2^^k*(1+k)]C

getValue : {k, n} (ispoweroftwo k n) => falconTree k -> Poly n C
getValue T = take`{n} T

// getChildren: {k, n} (ispoweroftwo k n, positive_power k) => falconTree k -> (falconTree (k-1), falconTree (k-1))
// getChildren T = children where
//     children = split`{2,2^^(k-1)*k} (drop`{n} T)

get_leftchild : {k, n} (ispoweroftwo k n, k > 0) => falconTree k -> falconTree (k-1)
get_leftchild T = left where
    children = resize (drop`{n} T): [2*2^^(k-1)*k]C
    [left, right] = split`{2,2^^(k-1)*k} children

get_rightchild : {k, n} (ispoweroftwo k n, k > 0) => falconTree k -> falconTree (k-1)
get_rightchild T = right where
    children = resize (drop`{n} T): [2*2^^(k-1)*k]C
    [left, right] = split`{2,2^^(k-1)*k} children
    
newTree : {k, n} (ispoweroftwo k n, k >= 1) => (Poly n C, falconTree (k-1), falconTree (k-1)) -> falconTree k
newTree(value, leftchild, rightchild) = T where
    T = zero

// The leaves contain real numbers
// For consistency their type is C
newLeaf : C -> falconTree 0
newLeaf sigma_leaf = [sigma_leaf]

get_leaves : {k, n} (ispoweroftwo k n) => falconTree k -> [2^^k]C 
get_leaves T | k == 0 => T
             | k >  0 => resize leaves where
                left_child = get_leftchild`{k} T : falconTree (k-1)
                right_child = get_rightchild`{k} T : falconTree (k-1)
                left_leaves  = get_leaves`{k-1} left_child : [2^^(k-1)]C
                right_leaves  = get_leaves`{k-1} right_child : [2^^(k-1)]C
                leaves = left_leaves # right_leaves : [2*2^^(k-1)]C

set_leaves : {k, n} (ispoweroftwo k n) => (falconTree k, [2^^k]C) -> falconTree k
set_leaves(T, leaves) | k == 0 => leaves
                      | k >  0 => newTree`{k}(value, newleftchild, newrightchild) where 
                        value = getValue`{k} T
                        leftchild = get_leftchild`{k} T 
                        rightchild = get_rightchild`{k} T
                        [leftleaves, rightleaves] = split`{2,2^^(k-1)} (resize leaves)
                        newleftchild = set_leaves`{k-1}(leftchild, leftleaves)
                        newrightchild = set_leaves`{k-1}(rightchild, rightleaves)

Compress : Poly _n ZZ -> [slen]
Compress s = zero

Decompress : [slen] -> Poly _n ZZ
Decompress str = zero


import Primitive::Keyless::Hash::SHAKE::SHAKE256
HashToPoint : {q', k, n, len} (q' <= 2^^16, q' >= 1, ispoweroftwo k n, fin len) => [len] -> Poly n (Z q')
HashToPoint str = [c i | i <- [0..(n-1)]] where
    k = 2^^16 / `q
    c i = fromInteger(t i % `q)
    t i = zero // TODO!

// phi and q are fixed so we do not parameterize this algorithm
// f and g cannot be sampled in Cryptol so we take them as input
Keygen : (Poly _n ZZ, Poly _n ZZ) -> (privateKey, publicKey)
Keygen (f', g') = (sk, pk) where
    (f, g, F, G) = NTRUGen(f', g'): (Poly _n ZZ, Poly _n ZZ, Poly _n ZZ, Poly _n ZZ)
    B = [[map IntToCmplx g,map IntToCmplx (-f)],[map IntToCmplx G,map IntToCmplx (-F)]] : [2][2](Poly _n C)
    Bhat = [[FFT'`{_k}(g),FFT'`{_k}(-f)],[FFT'`{_k}(G),FFT'`{_k}(-F)]] : [2][2](FFT _n)
    GG = B * Bhat
    T = ffLDL`{_k}(GG)
    leaves = get_leaves`{_k} T
    new_leaves = [(sigma/.r_leaf,0) | (r_leaf, _) <- leaves]
    T' = set_leaves`{_k}(T, new_leaves)
    sk = (Bhat, T')
    h = NTTInv`{_k}(NTT`{_k}(map fromInteger g)*(map ModInv (NTT`{_k}(map fromInteger f))))
    pk = h

N : {k, n} (ispoweroftwo k n, k > 0) => (Poly n ZZ) -> (Poly (2^^(k-1)) ZZ)
N(f) = zero

xSquared : {k, n} (ispoweroftwo k n) => (Poly n ZZ) -> (Poly (2^^(k+1)) ZZ)
xSquared(F) = resize (join[[c,0] | c <- F])

minusx : {k, n} (ispoweroftwo k n) => (Poly n ZZ) -> (Poly n ZZ)
minusx(g) = resize (join[ [g@i, -(g@(i+1))] | i <- [0, 2 .. (n-1)]])

// phi and q are fixed so we do not parameterize this algorithm
// f and g cannot be sampled in Cryptol so we take them as input
// We assume they pass conditions in lines 7 and 9
NTRUGen : (Poly _n ZZ, Poly _n ZZ) -> (Poly _n ZZ, Poly _n ZZ, Poly _n ZZ, Poly _n ZZ)
NTRUGen(f, g) = (f, g, F, G) where
    (F, G) = NTRUSolve`{_k}(f, g)

NTRUSolve : {k, n} (ispoweroftwo k n) => (Poly n ZZ, Poly n ZZ) -> (Poly n ZZ, Poly n ZZ)
NTRUSolve(f, g) | k == 0 => ([F], [G]) where
                                (gcd, u, v) = eGCD(f@0, g@0)
                                (F, G) = (u*`q, v*`q)
                | k >  0 => (_F, _G) where 
                                f' = N`{k}(f)
                                g' = N`{k}(g)
                                (F', G') = NTRUSolve`{k-1}(f',g')
                                F = PolyMulInZ`{k}(xSquared`{k-1}(F'), minusx`{k}(g))
                                G = PolyMulInZ`{k}(xSquared`{k-1}(G'), minusx`{k}(f)) 
                                (_F, _G) = Reduce`{k}(f,g,F,G)

eGCD : (ZZ, ZZ) -> (ZZ, ZZ, ZZ)
eGCD(a, b) =
    if a == 0
    then (b, 0, 1)
    else (g, t - (b / a) * s, s) where
        (g, s, t) = eGCD((b%a), a)

Reduce : {k, n} (ispoweroftwo k n) => (Poly n ZZ, Poly n ZZ, Poly n ZZ, Poly n ZZ) -> (Poly n ZZ, Poly n ZZ)
Reduce(f,g,F,G) = (_F, _G) where 
    numerator = PolyMulInZ`{k}(F, adjoint`{k}(f)) + PolyMulInZ`{k}(G, adjoint`{k}(g))
    denominator = PolyMulInZ`{k}(f, adjoint`{k}(f)) + PolyMulInZ`{k}(g, adjoint`{k}(g))
    _k = map roundToEven (PolyDivInZ`{k}(numerator, denominator)): Poly n ZZ        
    F' = F - PolyMulInZ`{k}(_k,f)
    G' = G - PolyMulInZ`{k}(_k, g)
    (_F, _G, _) = Reduce'`{k}(_k,f,g,F',G')

adjoint : {k, n} (ispoweroftwo k n) => (Poly n ZZ) -> (Poly n ZZ)
adjoint(f) = reverse f

Reduce' : {k, n} (ispoweroftwo k n) => (Poly n ZZ, Poly n ZZ, Poly n ZZ, Poly n ZZ, Poly n ZZ) -> (Poly n ZZ, Poly n ZZ, Poly n ZZ)
Reduce'(_k,f,g,F,G) =
    if _k == zero
    then (F, G, _k)
    else Reduce'`{k}(_k', f, g, F', G') where 
        numerator = PolyMulInZ`{k}(F, adjoint`{k}(f)) + PolyMulInZ`{k}(G, adjoint`{k}(g))
        denominator = PolyMulInZ`{k}(f, adjoint`{k}(f)) + PolyMulInZ`{k}(g, adjoint`{k}(g))
        _k' = map roundToEven (PolyDivInZ`{k}(numerator, denominator)): Poly n ZZ
        F' = F - PolyMulInZ`{k}(_k',f)
        G' = G - PolyMulInZ`{k}(_k', g)

one : {k, n} (ispoweroftwo k n) => FFT n
one = [(1.0,0.0) | i <- [0 .. (n-1)]]

LDL : {k, n} (ispoweroftwo k n) => [2][2](FFT n) -> (([2][2](FFT n)), ([2][2](FFT n)))
LDL(G) = (L, D) where 
    D00 = G@0@0
    L10 = HadamardDivision(G@1@0,G@0@0): FFT n
    D11 = G@1@1 - dot(L10,dot(PolyConj`{k}(L10),G@0@0))
    L = [[one`{k}, zero],[L10, one`{k}]]
    D = [[D00, zero],[zero, D11]]

ffLDL : {k, n} (ispoweroftwo k n, n >= 2) => ([2][2](FFT n)) -> (falconTree k)
ffLDL(G) | k == 1 => T where
            (L, D) = LDL`{k}(G) : ([2][2](FFT n), [2][2](FFT n))
            left_child = D@0@0 : FFT 2
            left_leaf = newLeaf(left_child@0) // Revise this
            right_child = D@1@1 : FFT 2
            right_leaf = newLeaf(right_child@0) // Revise this
            T = newTree`{k}(L@1@0, left_leaf, right_leaf)
         | k >  1 => T where
            (L, D) = LDL`{k}(G) : ([2][2](FFT n), [2][2](FFT n))
            [d00, d01] = splitfft`{k}(D@0@0)
            [d10, d11] = splitfft`{k}(D@1@1)
            G0 = [[resize d00,resize d01],[resize d01, resize d00]] // Revise this
            G1 = [[resize d10,resize d11],[resize d11, resize d10]]
            Tleftchild = ffLDL`{k-1}(G0)
            Trightchild = ffLDL`{k-1}(G1)
            T = newTree`{k}(L@1@0, Tleftchild, Trightchild)


type signature = ([320], [slen])
Sign : {len} (fin len) => ([len], privateKey, Integer, [320]) -> signature
Sign(m, sk, beta, r) = zero where 
    (Bhat, T) = sk
    c = HashToPoint`{q,_k}(r#m)
    [[FFTg,FFT_f],[FFTG,FFT_F]] = Bhat // Bhat stores them in FFT representation (some are negated)
    t0 = Normalize_FFT`{_k}(`q, dot(FFT''`{_k,_n}(c),FFT_F))
    t1 = Normalize_FFT`{_k}(`q, dot(FFT''`{_k,_n}(c),FFT_F))
    tt = [t0, t1] : [2](FFT _n)
    s_init : [slen]
    s_init = zero
    is_bot : [slen] -> Bit
    is_bot(_s) = (_s == zero)
    body' : [slen] -> [slen]
    body'(_s) = zero where
        ss_init = zero : FFT _n
        norm_greater_than_beta : FFT _n -> Bit
        norm_greater_than_beta(_ss) = norm_sq`{_k}(_ss) > (fromInteger beta)
        body'' : FFT _n -> FFT _n
        body''(_ss) = zero where 
            zz = ffSampling`{_k}(tt, T) : [2](FFT _n)
            // _ss' = [ dot`{_n}((tt_ - zz_),Bhat_) | tt_ <- tt | zz_ <- zz | Bhat_ <- Bhat]
        _ss_final = dowhile norm_greater_than_beta body'' ss_init
    s_final = dowhile is_bot body' s_init



ffSampling : {k, n} (ispoweroftwo k n) => ([2](FFT n), falconTree k) -> [2](FFT n)
ffSampling(t, T) = z where
    [t0, t1] = t
    z = [z0, z1]
    z0 = zero
    z1 = zero

RCDT = [
    3024686241123004913666,
    1564742784480091954050,
    636254429462080897535,
    199560484645026482916,
    47667343854657281903,
    8595902006365044063,
    1163297957344668388,
    117656387352093658,
    8867391802663976,
    496969357462633,
    20680885154299,
    638331848991,
    14602316184,
    247426747,
    3104126,
    28824,
    198,
    1,
    0
] : [19][72]

BaseSampler : [72] -> ZZ
BaseSampler(u) = z0 where 
    z0 = zero

import Float
ApproxExp : (Float64, Float64) -> [64]
ApproxExp(x, ccs) = y'' where 
    C = [
            0x00000004741183A3,
            0x00000036548CFC06,
            0x0000024FDCBF140A,
            0x0000171D939DE045,
            0x0000D00CF58F6F84, 
            0x000680681CF796E3, 
            0x002D82D8305B0FEA, 
            0x011111110E066FD0,
            0x0555555555070F00, 
            0x155555555581FF00, 
            0x400000000002B400, 
            0x7FFFFFFFFFFF4800,
            0x8000000000000000
        ]
    y = C@0
    z = fromInteger (floor(2^^63 * x))
    y' = foldl (\yy -> \u -> C@u - (z*y)>>63) y [1 .. 12]
    z' = fromInteger (floor(2^^63 * ccs))
    y'' = (z'*y')>>63

BerExp : (Float64, Float64) -> Bit
BerExp(x, ccs) = w < 0 where 
    ln2 = 0.69314718056
    s = x /. ln2
    r = floor(x - s * ln2)
    s' = min(s, 63)
    // TODO: Finish this
    w = zero

SamplerZ : (Float64, Float64) -> ZZ
SamplerZ(m, sigma') = z where 
    z = zero

Verify : {len} (fin len) => ([len], signature, publicKey, Integer) -> Bit
Verify(m, sig, pk, beta) = norm^^2 <= beta where 
    norm = zero