image source: Atomium Public Domain
In about a decade (or so we are told) we may reach a tipping point in the world of cryptography, as a practical quantum computer will become a reality. Personally I think it will take longer than that, perhaps even a lot longer. Often the people who anticipate quantum computers in the shorter term are the very same people that are looking for funding to do quantum research. A classic conflict of interest, but let’s face it no-one these days is going to fund research that may not pay off in their lifetime!
But for cryptographers the advent of a working quantum computer will be a game changer. Basically most current methods for public key cryptography will be broken. RSA, Elliptic Curve Crypto, all rendered useless overnight.
Fortunately awareness of this problem has been around for a while, and already there are well known alternatives that can be deployed, going right back to the McEliece scheme from the 1970s.
Unfortunately these alternatives are often unwieldy, and their strength is often based on not-very-well studied or understood hard problems. In comparison the strength of RSA rested fairly securely on the bed-rock of the hardness of the integer factorisation problem, which has been extensively studied for centuries. However for a quantum computer factorisation is easy.
Its hard to let go of all that beautiful number theory, but for researchers certainly the time has come to move on. I suspect that the majority of students who will be starting Ph.D. studies in cryptography in the future, will be focused almost entirely on what are referred to as post-quantum methods.
Getting started on post-quantum cryptography is hard. The maths is quite different, and lots of methods have been proposed, many of which are already broken. And remember that now we can use (non-existent) quantum computers in our cryptanalysis. (We may not have a quantum computer, but fortunately we can simulate them!) But through the fog a few potential winners are beginning to emerge.
Its certainly a different world. Key sizes tend to be a lot bigger. But implementations are often surprisingly efficient as big number arithmetic may no longer be required. Public key cryptography isn’t too hard, but digital signature can be very tricky. Protocols may randomly fail. But there are new possibilities for multi-linear maps, homomorphic encryption and all kinds of goodies that were hard to achieve using the classic tools.
One method which seems to be gaining some traction is based on the lattice-based Ring-LWE (Learning With Errors) problem, which is believed to be hard. This can be used as the basis for public key cryptography. Many recent research papers have focused on it.
So we implemented a method for public key encryption using Ring-LWE– see below. Even though its not optimised, its impressively fast. The program is quite short – less than 400 lines of C in one self-contained file. Deliberately introduced “errors” are an important part of the algorithm – and sometimes the an unfortunate accumulation of errors cause the decryption to fail. Parameters are chosen carefully to minimise this failure rate, while optimising the security. There is still some uncertainty as to what parameter sizes are required to reach any given level of security. Key sizes are not as massive as we had imagined – but run the program and see for yourself!
The source code provided here is placed in the public domain as a service to those who would like to experiment with post-quantum crypto. No excuse now for not “getting started”. Enjoy!
/* Demo of LWE Public Key Encryption
For full description see for example http://eprint.iacr.org/2015/410 and its references
M.Scott 11/11/2015
gcc lwe.c -o lwe.exe
*/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
/* Choose parameters */
/* {256,7681,11.31} */
#define PRIME 7681 // prime modulus (=1 mod 2.2^LGN)
#define LGN 8 // polynomial degree is 2^LGN
#define SB 11
/* {512,12289,12.18}
#define PRIME 12289
#define LGN 9
#define SB 11
*/
#define DEGREE (1<<LGN)
#define SCALE (1<<SB) // Scaling for probability matrix
/* some number theory functions borrowed from MIRACL */
int spmd(int x,int n,int m)
{ /* returns x^n mod m */
int r,sx;
x%=m;
r=0;
if (x==0) return r;
r=1;
if (n==0) return r;
sx=x;
for (;;)
{
if (n%2!=0) r=(r*sx)%m;
n/=2;
if (n==0) return r;
sx=(sx*sx)%m;
}
}
int invers(int x,int y)
{ /* returns inverse of x mod y */
int r,s,q,t,p,pos;
if (y!=0) x%=y;
r=1;
s=0;
p=y;
pos=1;
while (p!=0)
{ /* main euclidean loop */
q=x/p;
t=r+s*q; r=s; s=t;
t=x-p*q; x=p; p=t;
pos=!pos;
}
if (!pos) r=y-r;
return r;
}
int sqrmp(int x,int m)
{ /* square root mod a small prime =1 mod 8 by Shanks method *
* returns 0 if root does not exist or m not prime */
int z,y,v,w,t,q,i,e,n,r,pp;
x%=m;
if (x==0) return 0;
if (x==1) return 1;
if (spmd(x,(m-1)/2,m)!=1) return 0; /* Legendre symbol not 1 */
q=m-1;
e=0;
while (q%2==0)
{
q/=2;
e++;
}
if (e==0) return 0; /* even m */
for (r=2;;r++)
{ /* find suitable z */
z=spmd(r,q,m);
if (z==1) continue;
t=z;
pp=0;
for (i=1;i<e;i++)
{ /* check for composite m */
if (t==(m-1)) pp=1;
t=(t*t)%m;
if (t==1 && !pp) return 0;
}
if (t==(m-1)) break;
if (!pp) return 0; /* m is not prime */
}
y=z;
r=e;
v=spmd(x,(q+1)/2,m);
w=spmd(x,q,m);
while (w!=1)
{
t=w;
for (n=0;t!=1;n++) t=(t*t)%m;
if (n>=r) return 0;
y=spmd(y,(1<<(r-n-1)),m);
v=(v*y)%m;
y=(y*y)%m;
w=(w*y)%m;
r=n;
}
return v;
}
/* LWE code */
/* precompute 2^n -th root of unity and its powers */
void init_ntt(int16_t *roots)
{
int logn=LGN;
int q=PRIME;
int j,proot=q-1;
int newn=DEGREE;
for (j=1;j<logn;j++) proot=sqrmp(proot,q);
roots[0]=proot; /* build table of powers */
for (j=1;j<newn;j++) roots[j]=(roots[j-1]*proot)%q;
}
/* number theory transform */
void ntt(int16_t *roots,int16_t *data)
{ /* roots precomputed */
int m,j,k,istep,i,ii,jj,w,temp;
int logn=LGN;
int q=PRIME;
int newn=DEGREE;
int mmax=newn;
for (k=0;k<logn;k++) {
istep=mmax;
mmax>>=1;
ii=newn;
jj=newn/istep;
ii-=jj;
for (m=0;m<mmax;m++) {
w=roots[ii-1];
ii-=jj;
for (i=m;i<newn;i+=istep) {
j=i+mmax;
temp=q+data[i]-data[j];
data[i]+=data[j]-q;
data[i]+=(data[i]>>16)&q;
data[j]=(w*temp)%q;
}
}
}
}
/* inverse number theory transform */
void intt(int16_t *roots,int16_t *data)
{ /* roots precomputed */
int m,j,k,i,istep,ii,jj,w,temp;
int logn=LGN;
int q=PRIME;
int newn=DEGREE;
int mmax=1;
for (k=0;k<logn;k++) {
istep=(mmax<<1);
ii=0;
jj=newn/istep;
ii+=jj;
for (m=0;m<mmax;m++) {
w=roots[ii-1];
ii+=jj;
for (i=m;i<newn;i+=istep) {
j=i+mmax;
temp=(w*data[j])%q;
data[j]=data[i]-temp;
data[j]+=(data[j]>>16)&q; // if (data[j]<0) data[j]+=q;
data[i]+= temp-q;
data[i]+=(data[i]>>16)&q; // avoid branches
}
}
mmax=istep;
}
}
/* create error polynomial according to normal distribution */
void sampler(uint8_t *P,int16_t *x)
{
int i,j,k,q=PRIME;
int newn=DEGREE;
for (i=0;i<newn;i++)
{ /* generate random index into probability matrix */
k=rand();
j=P[k%SCALE];
x[i]=j;
if (j>0)
{ /* toss a coin to decide if negative */
if ((k>>SB)%2==0) x[i]=q-j;
}
}
}
/* generated from this dist.c utility */
/*
#include <stdio.h>
#include <math.h>
// choose power of 2 as maximum such that truncated outputs can be stored in bytes
#define SCALE 2048
int main()
{
double x,y,s=0;
double c=11.31; // set standard deviation (or 12.18)
double delta=c/sqrt(2*M_PI);
delta=2*delta*delta;
y=1/c;
printf("%lf\n",SCALE*y);
s=y;
for (int i=1;i<32;i++)
{
x=i;
y=exp(-(x*x)/delta);
y/=c;
printf("%lf\n",SCALE*y);
s+=(y+y);
}
printf("Sum (should be 1)= %lf\n",s);
return 0;
}
*/
/* standard deviation = 11.31 */
#if LGN==8
int lim[]={182,177,164,145,122,98,75,54,38,25,16,9,5,3,1,1,0}; // 182+2*(177+164....) = 2048
#endif
/* standard deviation = 12.18 */
#if LGN==9
int lim[]={168,165,154,139,120,99,78,60,43,30,20,13,8,5,3,1,1,1,0} ;
#endif
int main()
{
int i,j,k,u,s,t,lo,hi,inv,err,q=PRIME;
int16_t m[128],ma[128];
int16_t r1[DEGREE],r2[DEGREE],a[DEGREE],p[DEGREE],md[DEGREE],e1[DEGREE],e2[DEGREE],e3[DEGREE],c1[DEGREE],c2[DEGREE],roots[DEGREE];
uint8_t P[SCALE];
int newn=DEGREE;
srand(1); /* initialise random number generator */
/* Initialise probability matrix */
for (j=0;j<lim[0];j++)
P[j]=0;
lo=lim[0];
for (i=1;lim[i]!=0;i++)
{
hi=lo+2*lim[i];
for (j=lo;j<hi;j++) P[j]=i;
lo=hi;
}
if (hi!=SCALE)
{
printf("hi= %d\n",hi);
printf("Problem with generation of probability matrix\n");
return 0;
}
/* pre-calculate powers of roots of unity */
init_ntt(roots);
/* http://eprint.iacr.org/2015/410 section 2.2 */
/* Key Generation */
for (i=0;i<newn;i++)
a[i]=rand()%q;
sampler(P,r1); /* create random error polynomials */
sampler(P,r2);
ntt(roots,r1); /* convert to "frequency" domain... */
ntt(roots,r2);
for (i=0;i<newn;i++)
p[i]=(r1[i]+(int)(q-a[i])*r2[i])%q; /* ...so that polynomial multiplication is cheap */
printf("Public key= [");
for (i=0;i<newn-1;i++) printf("%d, ",a[i]);
printf("%d]\n",a[newn-1]);
printf("[");
for (i=0;i<newn-1;i++) printf("%d, ",p[i]);
printf("%d]\n",p[newn-1]);
printf("\nPrivate key= [");
for (i=0;i<newn-1;i++) printf("%d, ",r2[i]);
printf("%d]\n",r2[newn-1]);
/* Encryption */
for (i=0;i<128;i++)
{ /* Generate an AES key */
if (rand()%2==0) m[i]=0;
else m[i]=1;
}
printf("\nAES key= \n");
for (i=0;i<128;i++) printf("%d",m[i]);
printf("\n");
/* Encode message to withstand noise - see https://www.sha.rub.de/media/sh/veroeffentlichungen/2013/08/14/lwe_encrypt.pdf */
u=newn/128;
for (i=0;i<128;i++)
{
for (j=0;j<u;j++)
md[u*i+j]=m[i]*(q/2);
}
sampler(P,e1); /* create random noise error polynomials */
sampler(P,e2);
sampler(P,e3);
ntt(roots,e1);
ntt(roots,e2);
for (i=0;i<newn;i++)
c1[i]=(e2[i]+(int)a[i]*e1[i])%q;
for (i=0;i<newn;i++)
e3[i]=(e3[i]+md[i])%q;
ntt(roots,e3);
for (i=0;i<newn;i++)
c2[i]=(e3[i]+(int)p[i]*e1[i])%q;
printf("\nCiphertext= [");
for (i=0;i<newn-1;i++) printf("%d, ",c1[i]);
printf("%d]\n",c1[newn-1]);
printf("[");
for (i=0;i<newn-1;i++) printf("%d, ",c2[i]);
printf("%d]\n",c2[newn-1]);
/* Decryption */
for (i=0;i<newn;i++)
md[i]=(c2[i]+(int)r2[i]*c1[i])%q;
inv=invers(newn,q);
intt(roots,md);
for (i=0;i<newn;i++)
md[i]=(md[i]*inv)%q;
/* Decode */
for (i=0;i<128;i++)
{ /* recover message from noise using threshold function */
s=0;
for (j=0;j<u;j++)
{
t=md[u*i+j]-(q/2);
if (t>0) s+=t;
if (t<0) s-=t;
}
if (s<(u*q)/4) ma[i]=1;
else ma[i]=0;
}
printf("\nDecrypted AES key= \n");
for (i=0;i<128;i++) printf("%d",ma[i]);
printf("\n");
err=0;
for (i=0;i<128;i++)
if (m[i]!=ma[i]) err++;
if (!err) printf("Decryption Success\n");
else printf("Decryption Failed \n");
return 0;
}
