#!/usr/bin/ruby
# coding: utf-8

##################################################################################
#
# Marcin Woźniak
# s434812
#
##################################################################################

require 'openssl'
require 'securerandom'
require 'prime'
require 'thread'
require 'thwait'

#################################################################################
# Funkcja nwd(a,b)
#
# Oblicza nwd podanych liczb np
# nwd(10,14) => 2
#
#################################################################################
def nwd(a, b)
  if a == 0
    return false
  end
  b == 0 ? a : nwd(b, a.modulo(b))
end

#################################################################################
# Funkcja extended_euklides(a,b)
#
# Oblicza rozszerzony algorytm euklidesa
# zwracajac u,v
# extended_euklides(10,14) => [3, -2]
#
##################################################################################
def extended_euklides(a, b)
  return 1, 0 if b == 0

  q, r = a.divmod b
  s, t = extended_euklides(b, r)

  return t, s - q * t
end

#################################################################################
# Funkcja random_gen_Zn(k,n)
#
# Oblicza losowy element z zbioru Z_n
# random_gen_Zn(1,10) => 1
#
#################################################################################
def random_gen_Zn(k,n)

  if n == 0
    n = 2 ** k
  end

  if k == 1
    max = 1
  else
    kb = k.to_s(2)
    minimum = []
    minimum << 1
    k = kb.length - 1

    while (k != 0) do
      j = SecureRandom.random_number(2)
      minimum << j
      k = k - 1
    end
    min = minimum.join.to_i(2)
    max = n - 1
    if min < max
      return SecureRandom.random_number(min..max)
    end
  end
end

#################################################################################
# Funkcja reciprocal_Phi_p(n,p)
#
# Oblicza odwrotnosc w grupie Phi(n)
# reciprocal_Phi_p(10,7) => 5
#
# Funkcja oblicza 10^(-1) mod 7 = 5
# Dowód: 5 * 10 mod 7 = 1
#        a^(-1) * a mod p = e = 1
#################################################################################
def reciprocal_Phi_p(n,p)
  u = extended_euklides(n,p)[0]
  v = extended_euklides(n,p)[1]

  if (u * n % p) == 1
    if u < 0
      return u + p
    end
    return u
  else
    if v < 0
      return v + p
    end
    return v
  end
end

#################################################################################
# Funkcja betterExponentiation(x,k,n)
#
# Oblicza potęgi x^k mod n
# betterExponentiation(8,2,30) => 4
#
#################################################################################
def betterExponentiation(x,k,n)
  if n == 0
    return false
  end

  if x == 0
    return 0
  end

  if x < n && x > 0
    b = k.to_s(2).reverse
    l = b.count "[0-1]"
    y = 1
    i = l - 1

    while i >= 0
      y = y**2 % n
      if b[i]=="1"
        y = y * x % n
      end
      i = i - 1
    end
    return y
  end
end

#################################################################################
# Funkcja remSqEuler(a,b)
#
# Sprawdzenie czy element a jest
# reszta kwadratowa w Z_p
# remSqEuler(3,13) => true
#
#################################################################################
def remSqEuler(a,p)
  ans = betterExponentiation(a,(p-1)/2,p)
  if ans == 1 && primalityTest(p)
    return true
  else
    return false
  end
end

#################################################################################
# Funkcja squareRootFp(a,b)
#
#################################################################################
def squareRootFp(p,b)
  if p % 4 == 3 && remSqEuler(p,b) == true
      a = betterExponentiation(b, (p+1)/4, p)
      return a
  end
end

#################################################################################
# Funkcja primalityTest(n)
#
# Test pierwszości
# primalityTest(13) => true
#
#################################################################################
def primalityTest(n)
  if n == 1
    return false
  end

  if n == 2 || n == 3
    return true
  end

  counter = 10
  while (counter != 0) do
    b = SecureRandom.random_number(2..n-2) # Tez dziala n-1
      if betterExponentiation(b,n-1,n) != 1
        return false
      end
      counter = counter - 1
  end
  return true
end

#################################################################################
# Funkcja specyficPrimaryNumber
#
# Generuje potrzebne do pierwszego
# ElGamala
##################################################################################
def specyficPrimaryNumber
  p = 0
  q = 0

  qThread = Thread.new {
  while true
    q = SecureRandom.random_number(2 ** 256)
    if primalityTest(q)
      break
    end
  end
  }

  qThread.join

  while true do
    q = SecureRandom.random_number(2 ** 256)
    p = 2 * q + 1
    if primalityTest(p)
      return p,q
    end
  end
end

#################################################################################
# Funkcja generator(a,b)
#
# Generuje generator dla podanych liczb
# do RSA i ElGamala
#
#################################################################################
def generator(p,q)
  while true
    g = SecureRandom.random_number(2..p-2)
    if betterExponentiation(g,q,p) == 1
      next
    else
      return g
    end
  end
end

#################################################################################
# Funkcja generate(n)
#
# Generator liczby losowej o podanej
# n-bitowej liczbie
#
#################################################################################
def generate(n)
  return `openssl prime -generate -bits '#{n}'`.gsub(/\n$/, '').to_i
end

#### MODULE 2 ####

#################################################################################
# Funkcja returnRownanie(a,b,p)
#
# Zwraca ładniejszą formę równania
#
#################################################################################
def returnRownanie(a,b,p)
  puts
  puts "Równanie krzywej jest równe: " + "Y^2 = X^3+" + a.inspect + "X+" + b.inspect + " mod "+ p.inspect
  puts
end

#################################################################################
# Funkcja delta(a,b,p)
#
# Oblicza deltę dla krzywej
# eliptycznej o podanych parametrach
# a,b nad ciałem F_p
#
#################################################################################
def delta(a,b,p)
  d = ((4 * betterExponentiation(a,3,p) % p) + (27 * betterExponentiation(b,2,p) % p)) % p
  return d
end

#################################################################################
# Funkcja rownanieKrzywej(a,b,p,x)
#
# Oblicza wartość fx dla podanej
# krzywej eliptycznej o podanych
# parametrach a,b,p,x
#
#################################################################################
def rownanieKrzywej(a,b,p,x)
  fx = ((betterExponentiation(x,3,p) + (a * x) % p + b % p) % p) % p
  return fx
end

#################################################################################
# Funkcja generatorKrzywej(p)
#
# Generuje krzywą o podanej liczbie
# pierwszej p.
#
#################################################################################
def generatorKrzywej(p)
  a = 0
  b = 0

  while true
    if primalityTest(p) &&  (p % 4 == 3)
      threads = []

      threads << Thread.new {
        a = SecureRandom.random_number(1..p-1)
      }

      threads << Thread.new {
        b = SecureRandom.random_number(1..p-1)
      }

      ThreadsWait.all_waits(*threads)
      if delta(a,b,p) != 0
        #returnRownanie(a,b,p)
        return a,b
      end
    else
      puts "Liczba nie spełnia wymagań"
      break
    end
  end
end

#################################################################################
# Funkcja punktNaKrzywej(a,b,p)
#
# Zwraca losowy punkt na krzywej
# eliptycznej o podanych wartościach
# a,p,p.
#
#################################################################################
def punktNaKrzywej(a,b,p)
  if (delta(a,b,p) != 0) &&  (p % 4 == 3)
    while true
      x = SecureRandom.random_number(0..p-1)
      fx = rownanieKrzywej(a,b,p,x)
      if remSqEuler(fx,p)
        y = betterExponentiation(fx,((p+1)/4),p)
        return x,y
      end
    end
  end
end

#################################################################################
# Funkcja czyPunktNalezyDoKrzywej(a,b,p,x,y)
#
# Sprawdza czy podany punkt nalezy
# do krzywej eliptycznej.
#
#################################################################################
def czyPunktNalezyDoKrzywej(a,b,p,x,y)
  fx = rownanieKrzywej(a,b,p,x)

  if fx == betterExponentiation(y,2,p)
    return true
  else
    return false
  end
end

#################################################################################
# Funkcja punktPrzeciwny(x,y)
#
# Zwraca przeciwny punkt dla (x,y)
#
#################################################################################
def punktPrzeciwny(x,y)
  return x,-y
end

#################################################################################
# Funkcja sumaPunktow(a,b,p,x1,y1,x2,y2)
#
# Oblicza sumę punktów na krzywej
# eliptycznej dla podanych a,b,p,x1,y1,x2,y2
#
#################################################################################
def sumaPunktow(a,b,p,x1,y1,x2,y2)
  # 0 - element neutrany --> P + 0 = P
  if (x1 == "e" && y1 == "e" )
    return x2,y2
  elsif (x2 == "e" && y2 == "e")
    return x1,y1
  end
  # P + Q = R
  if (x1 != x2)
    lambda = (((y2 - y1) % p) * reciprocal_Phi_p((x2 - x1),p)) % p
    x3 = (betterExponentiation(lambda,2,p) - (x1 % p) - (x2 % p)) % p
    y3 = (lambda * (x1 - x3) - y1) % p
    return x3,y3
  end
  # P + -Q = 0 DZIALA
  if (x1 == x2) && (y1 == -y2)
    #puts "0 - el.neutralny"
    return "e","e"
  end
  # P + P = 2P DZIALA
  if (x1 == x2) && (y1 == y2)
    lambda = (((3 * betterExponentiation(x1,2,p) % p + a) % p) * reciprocal_Phi_p(2 * y1,p)) % p
    x3 = (betterExponentiation(lambda,2,p) - (x1 % p) - (x2 % p)) % p
    y3 = (lambda * (x1 - x3) - y1) % p
    return x3,y3
  end
  # P + -Q = 0 DZIALA
  if (x1 == x2) && (y1 == y2%p)
    #puts "0 - el.neutralny"
    return "e","e"
  end
end

#### MODULE 3 ####


#################################################################################
# Funkcje mysqrt(n)
#
# Oblicza wartość pierwiastka z dużych
# liczb
#
# Źródło: https://stackoverflow.com/questions/8226087/how-do-i-get-math-sqrt-to-return-a-bignum-and-not-a-float
#################################################################################
def mysqrt(x)
  return 0 if x==0
  m=x
  p=x
  loop do
    r=(m+p/m)/2
    return m if m<=r
    m=r
  end
end

#################################################################################
# Funkcja liczenieOrd(p)
#
# Oblicza wartość ord dla podanych p.
#################################################################################
def liczenieOrd(p)
  ord = (p + 1 - (2 * mysqrt(p))) % p
  return ord
end

#################################################################################
# Zadanie 1
# Funkcja generowanieKluczyElGamalKrzywaEliptyczna(k)
#
# Generatoruje klucz publiczny
# oraz prywatny.
#################################################################################
def generowanieKluczyElGamalKrzywaEliptyczna(k)
  while true
  p = generate(k)
    if (primalityTest(p)) && (p % 4 == 3)
      krzywa = generatorKrzywej(p)
      a = krzywa[0].to_i
      b = krzywa[1].to_i

      punktyNaKrzywej = Array.new

      punktP = punktNaKrzywej(a,b,p)
      ord = liczenieOrd(p)

      while true
        x = SecureRandom.random_number(1..ord)
        if x < ord
          punktQ = wielokrotnoscPunktu(a,b,p,x,punktP[0],punktP[1])
          pubKey = [a,b,p,punktP[0],punktP[1],punktQ[0],punktQ[1]]
          privKey = [a,b,p,punktP[0],punktP[1],punktQ[0],punktQ[1],x]
          return a,b,p,punktP[0],punktP[1],punktQ[0],punktQ[1],x
        end
      end
    end
  end
end

#################################################################################
# Zadanie 2
# Funkcja wielokrotnoscPunktu(a,b,p,n,x,y)
#
# Oblicza wielokrotność punktów np.
# 2n = n + n
#################################################################################
def wielokrotnoscPunktu(a,b,p,n,x,y)
  punktQ = [x,y]
  punktR = ["e","e"]

  while n > 0
    if n % 2 == 1
      punktR = sumaPunktow(a,b,p,punktR[0],punktR[1],punktQ[0],punktQ[1])
      n = n - 1
    end
    punktQ = sumaPunktow(a,b,p,punktQ[0],punktQ[1],punktQ[0],punktQ[1])
    n = n / 2
  end
  return punktR
end

#################################################################################
# Zadanie 3
# Funkcja algorytmKodowania(a,b,p,m,n,u)
#
#################################################################################
def algorytmKodowania(a,b,p,m,u)
  n = m + SecureRandom.random_number(0..1000000)
  if (m < n) && (p > n*u)
    for i in 1..u
      x = (m * u % p) + (i % p)
      fx = rownanieKrzywej(a,b,p,x)
      if remSqEuler(fx,p)
        y = betterExponentiation(fx,((p+1)/4),p)
      end
    end
  else
    puts "Nieprawidołowe dane"
  end
  puts "Punkt na krzywej to #{[x,y].inspect}"
  return [x,y]
end

#################################################################################
# Zadanie 4
# Funkcja szyfrowanieElGamala(a,b,p,m,u,px,py,qx,qy)
#
# Zwraca szyfrogramy.
#################################################################################
def szyfrowanieElGamala(a,b,p,u,px,py,qx,qy,pmx,pmy)
  y = SecureRandom.random_number(0..liczenieOrd(p))
  puts "y = #{y}"
  c1 = wielokrotnoscPunktu(a,b,p,y,px,py)
  yq = wielokrotnoscPunktu(a,b,p,y,qx,qy)
  c2 = sumaPunktow(a,b,p,pmx,pmy,yq[0],yq[1])
  puts "Ciphers: C1=#{c1}, C2=#{c2}"
  return c1,c2
end

#################################################################################
# Zadanie 5
# Funkcja deKodowanieElGamala(a,b,p,c1x,c1y,c2x,c2y,x)
#
# Zwraca odszyfrowany punkt początkowy na którym była wiadomość.
#################################################################################
def deKododwanieElGamala(a,b,p,c1x,c1y,c2x,c2y,x)
  xc1 = wielokrotnoscPunktu(a,b,p,x,c1x,c1y)
  pmd = sumaPunktow(a,b,p,c2x,c2y,xc1[0],-xc1[1])
  #puts "-------------------DEBUG MODE-------------------\na = #{a}\nb = #{b}\np = #{p}\nc1x = #{c1x} \nc1y = #{c1y} \nc2x = #{c2x} \nc2y = #{c2y} \nx = #{x}\nXC1 = #{xc1.inspect}\npm = #{pmd.inspect}\n-------------------"
  return pmd
end

#################################################################################
# Zadanie 6
# Funkcja algorytmDeSzyfrowania(x,y,u)
#
# Zwraca odszyfrowaną wiadomość.
#################################################################################
def algorytmDeSzyfrowania(x,y,u)
  m = (x - 1) / u
  return m
end

########################### Module 4 ############################################
# Podstawowe operacje na Galois Field GF(2^8)
# Źródła:
# * https://people.scs.carleton.ca/~maheshwa/courses/4109/Seminar11/The_Advanced_Encryption_Standard_AES_.pdf
# * https://cs465.internet.byu.edu/static/lectures/w19/AES.pdf
# * https://en.wikipedia.org/wiki/Finite_field_arithmetic
# * https://swarm.cs.pub.ro/~mbarbulescu/cripto/Understanding%20Cryptography%20by%20Christof%20Paar%20.pdf
# * http://www.cs.man.ac.uk/~banach/COMP61411.Info/CourseSlides/Wk2.2.FinField.pdf
#################################################################################

#################################################################################
# Zadanie 1
# Funkcja suma(a,b) wykorzystujac liczby hex
#################################################################################
def suma(a,b)
  binA = a.to_i(16).to_s(2)
  binB = b.to_i(16).to_s(2)
  return (binA.to_i(2) ^ binB.to_i(2)).to_s(16)
end

#################################################################################
# Zadanie 2
# Funkcja xtime(a) wykorzystujac liczby hex
#################################################################################
def xtime(a)
  binA = a.to_i(16).to_s(2)
  const = "1B"

  dl = binA.length
  while dl != 8
    binA = "0" + binA
    dl = dl + 1
  end

  if binA[0].to_i == 1
    binA[0] = ''
    return suma((binA.to_i(2) << 1 ).to_s(16), const)
  elsif binA[0].to_i == 0
    return (binA.to_i(2) << 1 ).to_s(16)
  end
end

#################################################################################
# Zadanie 3
# Funkcja iloczyn(a,b) wykorzystujac liczby hex
# {53} • {CA} = {01}
# {57} • {13} = {fe}
#################################################################################
def iloczyn(a,b)
  solve = "0"
  binA = a.to_i(16).to_s(2)

  len = binA.length - 1
  binA.split('').each { |a|
    if a == "1"
      tmp = b
      counter = len
      while counter != 0
        tmp = xtime(tmp)
        counter = counter - 1
      end
      solve = suma(tmp.to_i(16).to_s(16), solve)
    end
    len = len - 1
  }
  return solve
end

#################################################################################
# Zadanie 4
# Funkcja odwrotnosc(a) wykorzystujac liczby hex
#################################################################################
def odwrotnosc(a)
  b = a
  i = 13
  const = 1
  while i > 0
    if i.to_s(2).to_i & const.to_s(2).to_i == 1
      tmp = b
    else
      tmp = a
    end
    b = iloczyn(b,tmp)
    i = i - 1
  end
  return b
end