Paillier.m

# define some functions

# modular exponentiation
# source: https://gist.github.com/ttezel/4635562
function result = modexp (x, y, n)
    %anything raised to 0th power = 1 so return 1
    if (y == 0)
        result = 1;
        return;
    end

    %recurse
    z = modexp(x, floor(y/2), n);

    %if even square the result
    if (mod(y, 2) == 0)
        result = mod(z*z, n);
        return;
    else
        %odd so square the result & multiply by itself
        result = mod(x*z*z, n);
        return;
    end
end

# extended euclidean algorithm
# source: http://stanoyevitch.net/CRYPT/CRYPT/EuclidAlgExt.m
function OutVec = extended_euclidean_algorithm(a,b)

  aa = max(a,b); bb = min(a,b);
  U = [aa 1 0]; V = [bb 0 1];
  
  while V(1) > 0
    W = U - floor(U(1)/V(1))*V;
    U = V; V = W;
  end
  
  d = U(1); x = U(2); y = U(3);
  
  OutVec = [d x y];
endfunction

# inverse modulation
function retval = inverse_of(n,p)
  gcd = extended_euclidean_algorithm(n,p)(1)
  y   = extended_euclidean_algorithm(n,p)(2)
  x   = extended_euclidean_algorithm(n,p)(3)
  
  retval = mod(x,p)
  
endfunction

# L function
function retval = L(x,n)
  retval = floor((x-1)/n);
endfunction

# define primes
p=17
q=19

# out data to be encrypted
m=10

if (p==q)
  disp ("P and Q cannot be the same")
endif

n = p*q

gLambda = lcm(p-1,q-1)

g=randi([20,150])
if (gcd(g,n*n) == 1)
  disp("g and n*n are relatively prime (coprime)")
else
  disp("g and n*n are not relatively prime, generate a new g")
endif
  
r=randi([20,150])

l = floor( (modexp(g, gLambda, n^2)-1) / n )
gMu = inverse_of(l, n)

k1 = modexp(g, m, n^2)
k2 = modexp(r, n, n^2)

cipher = mod( (k1 * k2), (n*n) )

l = floor( (modexp(cipher, gLambda, n^2)-1) / n )

mess = mod( (l * gMu), n)

m1=2

k3 = modexp(g, m1, n*n)

cipher2 = mod( (k3 * k2), (n*n) ) 

ciphertotal = mod( (cipher * cipher2), (n*n) )

l = floor( (modexp(ciphertotal, gLambda, n^2)-1) / n )

mess2 = mod( (l * gMu), n )