CTLWE homomorphic encryption scheme

CTLWE.md

These are code snippets of the CTLWE holomorphic encryption scheme implemented on Optalysys' first silicon photonics optical computing demonstrator, described in more details in the Medium article https://medium.com/p/ee817ddc861c .

While the CTLWE scheme is, to the best of our knowledge, secure against classical and quantum attacks, this particular implementation was designed as a demonstration rather than a ready-to-use suite. It was used as a proof of principle that holomorphic encryption can be realized optically, and will form the basis of more efficient and fully secure optical encryption systems we are currently working on.

decrypt.cpp

v2D<double> CTLWE::decrypt(v3D<double> &cipher) {
    
    // dot product of the cipher and key
    v2D<double> s_dot_a = dot_product(key, cipher);

    // compute the decrypted message
    v2D<double> decrypted = cipher[k];
    for (int i=0; i<N1; i++) {
        for (int j=0; j<N2; j++) {
            
            // subtract the the coefficient of s_dot_a and round the result
            double x = floor(
                modulo_1(decrypted[i][j] - s_dot_a[i][j]) * M + 0.5
                ) / M;
            
            // convert 1 to 0 (they represent the same value in this scheme)
            decrypted[i][j] = x < 1 ? x : 0.; 
        }
    }

    return decrypted;
}

dot_product.cpp

template <typename T, typename U>
v2D<U> CTLWE::dot_product(v3D<T> &a1, v3D<U> &a2) {

    // initialize the result with the first product
    v2D<U> res = multinomial_multiplication(a1[0],a2[0]);

    // add the other products of multinomials 
    for (int l=1; l<a1.size(); l++) {
        v2D<U> plane = multinomial_multiplication(a1[l],a2[l]);
        add_multinomial(res, plane);
    }
    
    // return the result
    return res;
}

encrypt.cpp

v3D<double> CTLWE::encrypt(v2D<double> &message) {
    
    // choose the array a at random
    v3D<double> a; 
    for (int l=0; l<k; l++) {
        v2D<double> plane;
        for (int i=0; i<N1; i++) {
            v1D<double> row;
            for (int j=0; j<N2; j++) {
                row.push_back(u(gen));
            }
            plane.push_back(row);
        }
        a.push_back(plane);
    }

    // at this stage, a is a 3D array of numbers randomly chosen between 0 and 
    // 1 with a uniform probability density

    // compute b
    v2D<double> b = message;
    v2D<double> s_dot_a = dot_product(key, a); // this is the bottleneck
    add_multinomial(b, s_dot_a);

    // at this stage, b is equal to the message plus the dot product of a and the 
    // key

    // add Gaussian noise to b and cast the coefficients of a and b between 0. and 1.
    for (int i=0; i<N1; i++) {
        for (int j=0; j<N2; j++) {
            b[i][j] += d(gen);
            b[i][j] = modulo_1(b[i][j]);
            for (int l=0; l<k; l++) {
                a[l][i][j] = modulo_1(a[l][i][j]);
            }
        }
    }
    
    // append b to a
    a.push_back(b);
    
    // return the ciphertext
    return a;
}

mulMult.cpp

co* mulMult(co* m1, co* m2, int N1, int N2) {

    // 1. Dephasing
    
    // initialize the arrays
    co* m1_d = (co*) malloc(N1*N2*sizeof(co));
    co* m2_d = (co*) malloc(N1*N2*sizeof(co));
    
    // compute each of their coefficients
    for (int i=0; i<N1; i++) {
        for (int j=0; j<N2; j++) {
            
            // exponent giving the right phase
            co exponent = {0., (M_PI*i)/N1 + (M_PI*j)/N2};
            co exp_ = co_exp(&exponent);
            
            // multiply the coefficients (i,j) of m1 and m2 by exp_
            m1_d[i*N2+j] = co_mul(m1+i*N2+j, &exp_);
            m2_d[i*N2+j] = co_mul(m2+i*N2+j, &exp_);
        }
    }

    // 2. Fourier transforms
    co* F_m1_d = Fourier2(m1_d, N1, N2);
    co* F_m2_d = Fourier2(m2_d, N1, N2);

    // 3. Multiplication in Fourier space
    
    // complex array which will store the result
    co* F_res_d = (co*) malloc(N1*N2*sizeof(co));
    
    // loop over the coefficients
    for (int i=0; i<N1; i++) {
        for (int j=0; j<N2; j++) {
            // insert the product of the coefficients of the two multinomials
            F_res_d[i*N2+j] = co_mul(F_m1_d+i*N2+j, F_m2_d+i*N2+j);
        }
    }

    // 4. Inverse Fourier transform
    co* res_d = iFourier2(F_res_d, N1, N2);

    // 5. Inverse dephasing
    co* res = (co*) malloc(N1*N2*sizeof(co));
    for (int i=0; i<N1; i++) {
        for (int j=0; j<N2; j++) {
            
            // opposite of the exponent of step 1
            co exponent = {0., -(M_PI*i)/N1-(M_PI*j)/N2};
            co exp_ = co_exp(&exponent);
            
            // multiply the coefficients i,j of res_d by exp_
            res[i*N2+j] = co_mul(res_d+i*N2+j, &exp_);
        }
    }

    // free the temporary pointers
    free(res_d);
    free(F_res_d);
    free(F_m1_d);
    free(F_m2_d);
    free(m1_d);
    free(m2_d);

    // return the result
    return res;
}