affine_cipher(s,a,b):=block([tmp], 
if not(equal(gcd(a,26),1)) then 
return("Bad parameters: a is not coprime with 26") 
else 
(tmp:mod(a*(map(cint,charlist(s))-97)+b,26), 
simplode(map(ascii,tmp+65))) 
)$ 

affine_decipher(s,a,b):=block([tmp], 
tmp:mod(inv_mod(a,26)*(map(cint,charlist(s))-65-b),26), 
simplode(map(ascii,tmp+97)) 
)$

text_to_numbers(s):=map(cint,charlist(s))-96$

base27(list):=block( sum(reverse(list)[j]*27^(j-1),j,1,length(list)))$

dec2base(n,b):=block([remainders,tmp], 
remainders:[n], 
do (tmp:rest(remainders,1), 
remainders:flatten(push(divide(remainders[1],b),tmp)), 
if is(first(remainders)<b) then return (remainders) 
)
)$

load(distrib)$

Elgamal_parameters(n):=block([p,tmp,roots], 
p:next_prime(floor(random_normal(2^n+12345,100,100)[random(100)])), 
tmp:ifactors(totient(p)), 
m:random(p-100), 
roots:sublist(makelist(i,i,m,m+99),lambda([x],zn_primroot_p(x,p,tmp))), 
[p,first(roots)] 
)$

Elgamal_publickey(p,r,N):=power_mod(r,N,p)$

Elgamal_cipher(s,A,k,p,r):=block([M,K], 
M:base27(text_to_numbers(s)), 
K:power_mod(A,k,p), 
M1:power_mod(r,k,p), 
M2:mod(K*M,p), 
[M1,M2] 
)$

Elgamal_decipher(M,N,p):=block([K,tmp], 
K:power_mod(M[1],N,p), 
tmp:power_mod(power_mod(K,-1,p)*M[2],1,p), 
simplode(map(ascii,dec2base(tmp,27)+96))
)$

Hill_cipher(s,K):=block([tmp1,tmp2,tmp3,tmp4,tmp5,L1,L2,m], 
m:length(K), 
if is(equal(second(divide(slength(s),m)),0)) then ( 
tmp4:map(cint,charlist(s))-97, 
L2:length(tmp4), 
tmp5:makelist(create_list(tmp4[j],j,(i-1)*m+1,m*i),i,1,L2/m), 
simplode(map(ascii, mod(flatten( makelist(list_matrix_entries(tmp5[j].transpose(K)),j,1,length(tmp5))) ,26) +65)) ) 
else
( tmp1:smake(m-second(divide(slength(s),m)),ascii(120)), 
tmp2:map(cint,charlist(concat(s,tmp1)))-97, 
L1:length(tmp2), 
tmp3:makelist(create_list(tmp2[j],j,(i-1)*m+1,m*i),i,1,L1/m), 
simplode(map(ascii, mod(flatten( makelist(list_matrix_entries(tmp3[j].transpose(K)),j,1,length(tmp3))) ,26) +65)) ) 
)$


zn_linsyst(listec,listvar,n):= block([c,tmp,divs,L,i0,j0,cinv,A,B,a,b,d,q,mprime,u,x,y], 
c:coefmatrix(listec,listvar), 
tmp:list_matrix_entries(c), 
for i:1 thru 4 do 
divs[i]:gcd(tmp[i],n), 
L:makelist(divs[i],i,1,4), 
if not(member(1,L)) then return( "There are no solutions, no coefficient is coprime with the modulus" ) 
else ( 
if is(equal(L[1],1)) then (i0:0,j0:0), 
if is(equal(L[2],1)) then (i0:0,j0:1), 
if is(equal(L[3],1)) then (i0:1,j0:0), 
if is(equal(L[4],1)) then (i0:1,j0:1), 
cinv:inv_mod(c[i0+1,j0+1],n), 
A:c[mod(i0+1,2)+1,mod(j0+1,2)+1]- cinv*c[mod(i0+1,2)+1,j0+1]*c[i0+1,mod(j0+1,2)+1],
B:rhs(listec[mod(i0+1,2)+1])- inv_mod(c[i0+1,j0+1],n)*rhs(listec[i0+1])*c[mod(i0+1,2)+1,j0+1], 
a:mod(A,n), 
b:mod(B,n), 
d:gcd(a,n), 
if elementp(d,divisors(b)) then ( q:b/d, mprime:n/d, u:first(gcdex(a,n)), 
x:mod(mod(q*u,n),mprime) ) 
else return( "There are no solutions: in ax=b mod n, gcd(a,n) does not divide b" ), 
[ listvar[mod(j0+1,2)+1] =x, listvar[j0+1]=mod(cinv*(rhs(listec[i0+1])-c[i0+1,mod(j0+1,2)+1]*x),n) ] ) 
)$

matrix_mod_inv(A,n):=block([det,mdet], 
det:determinant(A), 
if equal(gcd(det,n),1) then (mdet:inv_mod(det,n), 
mod(mdet*adjoint(A),n)) 
else return("The matrix is not invertible mod n") 
)$