{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GCD / Euclidean algorithm\n",
    "Explaination: https://brilliant.org/wiki/extended-euclidean-algorithm/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Returns the greatest common divisor of a and b.\n",
    "def gcd(a, b):\n",
    "    if b > a:\n",
    "        return gcd(b, a)\n",
    "    if a % b == 0:\n",
    "        return b\n",
    "    return gcd(b, a % b)\n",
    "\n",
    "gcd(35,30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Back substitution / Extended euclidean algorithm\n",
    "Explaination: https://brilliant.org/wiki/extended-euclidean-algorithm/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ax + by = gcd(a,b)    ==>    35*1 + 10*-3 = gcd(35,10)\n",
      "\n",
      "gcd(35,10) = 5\n",
      "x = 1, y = -3\n"
     ]
    }
   ],
   "source": [
    "def gcdExtended(a, b):  \n",
    "    # Base Case  \n",
    "    if a == 0 :   \n",
    "        return b,0,1\n",
    "\n",
    "    gcd,x1,y1 = gcdExtended(b%a, a)\n",
    "\n",
    "    # Update x and y using results of recursive call  \n",
    "    x = y1 - (b//a) * x1  \n",
    "    y = x1  \n",
    "\n",
    "    return gcd,x,y \n",
    "\n",
    "a, b = 35,10\n",
    "g, x, y = gcdExtended(a, b)\n",
    "print(f\"ax + by = gcd(a,b)    ==>    {a}*{x} + {b}*{y} = gcd({a},{b})\\n\")\n",
    "\n",
    "print(f\"gcd({a},{b}) = {g}\")\n",
    "print(f\"x = {x}, y = {y}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modular inverse, A^-1\n",
    "Explanation: https://www.khanacademy.org/computing/computer-science/cryptography/modarithmetic/a/modular-inverses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3 * 5 ≡ 1 (mod 7)\n",
      "Modular inverse for 3 (mod 7) is 5\n",
      "a^-1 = 5\n"
     ]
    }
   ],
   "source": [
    "# Naive method of finding modular inverse for A (mod C) \n",
    "def inverse(a, c):\n",
    "    for b in range(c):\n",
    "        if (a*b)%c == 1: return b\n",
    "    print(\"NO INVERSE FOUND\")\n",
    "\n",
    "a,c = 3,7\n",
    "b = inverse(a,c)\n",
    "\n",
    "print(f\"{a} * {b} ≡ 1 (mod {c})\")\n",
    "print(f\"Modular inverse for {a} (mod {c}) is {b}\")\n",
    "print(f\"a^-1 = {b}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chineese remainder theorem\n",
    "Source: https://rosettacode.org/wiki/Chinese_remainder_theorem#Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "23\n"
     ]
    }
   ],
   "source": [
    "from functools import reduce\n",
    "def crt(n, a):\n",
    "    sum = 0\n",
    "    prod = reduce(lambda a, b: a*b, n)\n",
    "    for n_i, a_i in zip(n, a):\n",
    "        p = prod // n_i\n",
    "        sum += a_i * mul_inv(p, n_i) * p\n",
    "    return sum % prod\n",
    "\n",
    "def mul_inv(a, b):\n",
    "    b0 = b\n",
    "    x0, x1 = 0, 1\n",
    "    if b == 1: return 1\n",
    "    while a > 1:\n",
    "        q = a // b\n",
    "        a, b = b, a%b\n",
    "        x0, x1 = x1 - q * x0, x0\n",
    "    if x1 < 0: x1 += b0\n",
    "    return x1\n",
    "\n",
    "n = [3, 5, 7]\n",
    "a = [2, 3, 2]\n",
    "print(crt(n, a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Euler function φ and Z*n\n",
    "From slides 2020-4135-l07 - Number Theory for Public Key Cryptography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "φ(81) = 54 since the 54 items in Z*81 are each relatively prime to 81.\n",
      "Z*81 = [1, 2, 4, 5, 7, 8, 10, 11, 13, 14, 16, 17, 19, 20, 22, 23, 25, 26, 28, 29, 31, 32, 34, 35, 37, 38, 40, 41, 43, 44, 46, 47, 49, 50, 52, 53, 55, 56, 58, 59, 61, 62, 64, 65, 67, 68, 70, 71, 73, 74, 76, 77, 79, 80]\n"
     ]
    }
   ],
   "source": [
    "# Relatively prime means that gcd equals 1\n",
    "def is_relative_prime(a,b):\n",
    "    if b == 0 or a == 0: return False\n",
    "    return gcd(a,b) == 1\n",
    "\n",
    "# The set of positive integers less than n and relatively prime \n",
    "# to n form the reduced residue class Z∗n.\n",
    "def z_star(n):\n",
    "    r = []\n",
    "    for i in range(1,n):\n",
    "        if is_relative_prime(i,n):\n",
    "            r.append(i)\n",
    "    return r\n",
    "\n",
    "n = 81\n",
    "ar = z_star(n)\n",
    "e = len(ar)\n",
    "print(f\"φ({n}) = {e} since the {e} items in Z*{n} are each relatively prime to {n}.\")\n",
    "print(f\"Z*{n} = {ar}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing for primality\n",
    "From slides 2020-4135-l07 - Number Theory for Public Key Cryptography"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "137 True\n"
     ]
    }
   ],
   "source": [
    "from math import sqrt\n",
    "# Inefficient, but accurate method\n",
    "def is_prime(n):\n",
    "    if n % 2 == 0 and n > 2: \n",
    "        return False\n",
    "    return all(n % i for i in range(3, int(sqrt(n)) + 1, 2))\n",
    "\n",
    "print(137, is_prime(137))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1105 is probable prime\n"
     ]
    }
   ],
   "source": [
    "# Fermat primality test\n",
    "\n",
    "# Inputs: \n",
    "#    n: a value to test for primality;\n",
    "#    k: a parameter that determines the number of \n",
    "#       times to test for primality\n",
    "from random import randint\n",
    "def fermat_prime(n, k):\n",
    "    for x in range(k):\n",
    "        a = randint(2,n-2)\n",
    "        # Fermat’s theorem says that if a number p is prime then \n",
    "        # a^p−1 mod p = 1 for all a with gcd(a,p) = 1.\n",
    "        if not (a**(n-1))%n == 1:\n",
    "            return \"composite\"\n",
    "    return \"probable prime\"\n",
    "\n",
    "n = 1105 # value to test for primality;\n",
    "k = 1 # number of times to test for primality\n",
    "\n",
    "print(f\"{n} is {fermat_prime(n,k)}\")\n",
    "\n",
    "# First few Carmichael numbers are: 561, 1105, 1729, 2465 \n",
    "# Will output probable prime for every a with gcd(a, n) = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Miller-Rabin\n",
    "From https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test  \n",
    "This test is not fooled by Charmichael numbers.\n",
    "\n",
    "To test `n`, you choose different `a` values, and test if \n",
    "\n",
    "> $a^{d*2^i} \\equiv 1 \\pmod{n}$ for $i = (0, 1, 2, ... r)$\n",
    "\n",
    "If this happens, look at the congruence on the calculation just before you got `1`. If that value is not -1 (aka n-1), the value is composite aka not prime."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "168 = 21*2^3 = 168\n",
      "65536 = 1*2^16 = 65536\n",
      "\n",
      "16127 is probably prime.\n",
      "561 is composite.\n",
      "2465 is composite.\n",
      "65537 is probably prime.\n"
     ]
    }
   ],
   "source": [
    "# Miller-Rabin primality test\n",
    "# This test is not fooled by Charmichael numbers.\n",
    "\n",
    "from math import log2, floor\n",
    "from random import randint\n",
    "\n",
    "def decompose_even_number(n):\n",
    "    '''Find the values r and d such that n == d * 2**r.'''\n",
    "    r = 0\n",
    "    d = n\n",
    "    for i in range(1, floor(log2(n)) + 1):\n",
    "        if (d/2).is_integer():\n",
    "            d = d/2\n",
    "            r += 1\n",
    "        else:\n",
    "            break\n",
    "    return (r, int(d))\n",
    "\n",
    "verbose = False  # Print during execution\n",
    "\n",
    "def miller_rabin_prime(n):\n",
    "    ''':param n: should be odd, and n>3. No primes can be even, for obvious reasons.'''\n",
    "    accuracy = 3  # How many different `a` to test\n",
    "    r, d = decompose_even_number(n-1)  # n is odd, thus n-1 is even\n",
    "    verbose and print(f\"Testing {n} a maximum of {accuracy} times...\")\n",
    "    \n",
    "    for _ in range(accuracy):\n",
    "        try:\n",
    "            a = randint(2, n-2)  # The \"witness\".  TODO: don't reuse a witness.\n",
    "            verbose and print(f\"\\tTesting {n} with a={a}\")\n",
    "\n",
    "            x = (a**d) % n\n",
    "            if x == 1 or x == n-1:\n",
    "                verbose and print(f\"\\t({a}**{d})%{n} was 1 or -1 already!\")\n",
    "                continue # This witness was immediately 1 or -1 (aka n-1), so it can't help us any more.\n",
    "            for squaring in range(r-1):\n",
    "                x = (x**2) % n\n",
    "                verbose and print(f\"\\tx={x}\")\n",
    "                if x == n-1:\n",
    "                    raise Exception(\"Try next witness\")  # This witness can't help any more.\n",
    "            return \"composite\" # The witness didnt uphold the criteria for a prime, so n is composite!\n",
    "        except Exception as tryNext:\n",
    "            continue\n",
    "    return \"probably prime\"\n",
    "    \n",
    "s, d = decompose_even_number(168)\n",
    "print(f\"168 = {d}*2^{s} = {d*(2**s)}\")\n",
    "s, d = decompose_even_number(65536)\n",
    "print(f\"65536 = {d}*2^{s} = {d*(2**s)}\")\n",
    "print()\n",
    "\n",
    "for value in [16127, 561, 2465, 65537]:\n",
    "    print(f\"{value} is {miller_rabin_prime(value)}.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integer factorisation\n",
    "Given an integer, find its prime factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The factors of 641361 are [3, 7, 7, 4363], verified\n",
      "Manual test: True\n"
     ]
    }
   ],
   "source": [
    "# Brute force solution\n",
    "def prime_factors(n):\n",
    "    i = 2\n",
    "    factors = []\n",
    "    while i * i <= n:\n",
    "        if n % i:\n",
    "            i += 1\n",
    "        else:\n",
    "            n //= i\n",
    "            factors.append(i)\n",
    "    if n > 1:\n",
    "        factors.append(n)\n",
    "    return factors\n",
    "\n",
    "# Check if the numbers in array a is the factors of n\n",
    "def factor_test(a, n):\n",
    "    k = 1\n",
    "    for i in factors:\n",
    "        k=k*i\n",
    "    return k==n\n",
    "\n",
    "n = 641361\n",
    "factors = prime_factors(n)\n",
    "print(f\"The factors of {n} are {factors}, {'verified' if factor_test(factors, n) else 'ERROR'}\")\n",
    "\n",
    "print(\"Manual test:\", factor_test([3, 7, 7, 4363], 641361))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discrete logarithm problem\n",
    "Given a prime `p` and an integer `y` with `0 < y < p`, find `x` such that \n",
    "```\n",
    "y = g^x mod p\n",
    "```\n",
    "\n",
    "No polynomial algorithims exists, but \"Baby-step giant-step\" is a reasonable one. \n",
    "\n",
    "Sources:\n",
    "- https://stackoverflow.com/a/1832648/5976426\n",
    "- https://en.wikipedia.org/wiki/Baby-step_giant-step\n",
    "- https://gist.github.com/0xTowel/b4e7233fc86d8bb49698e4f1318a5a73"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = 2570\n",
      "y = g^x mod p    ==>    34 = 7324^2570 mod 4363\n",
      "Test: True\n"
     ]
    }
   ],
   "source": [
    "from math import ceil, sqrt\n",
    "\n",
    "def bsgs(g, y, p):\n",
    "    '''\n",
    "    Solve for x in y = g^x mod p given a prime p.\n",
    "    If p is not prime, you shouldn't use BSGS anyway.\n",
    "    '''\n",
    "    N = ceil(sqrt(p - 1))  # phi(p) is p-1 if p is prime\n",
    "\n",
    "    # Store hashmap of g^{1...m} (mod p). Baby step.\n",
    "    tbl = {pow(g, i, p): i for i in range(N)}\n",
    "\n",
    "    # Precompute via Fermat's Little Theorem\n",
    "    c = pow(g, N * (p - 2), p)\n",
    "\n",
    "    # Search for an equivalence in the table. Giant step.\n",
    "    for j in range(N):\n",
    "        h = (y * pow(c, j, p)) % p\n",
    "        if h in tbl:\n",
    "            return j * N + tbl[h]\n",
    "\n",
    "    # Solution not found\n",
    "    return None\n",
    "\n",
    "def disc_log_test(y, g, x, p):\n",
    "    return y == (g**x) % p\n",
    "\n",
    "#y = 355407489\n",
    "#g = 7894352216\n",
    "#p = 604604729 # Must be prime\n",
    "\n",
    "# y = g^x mod p\n",
    "y = 34\n",
    "g = 7324\n",
    "p = 4363 # Must be prime\n",
    "\n",
    "x = bsgs(g, y, p)\n",
    "\n",
    "print(f\"x = {x}\")\n",
    "\n",
    "print(f\"y = g^x mod p    ==>    {y} = {g}^{x} mod {p}\")\n",
    "\n",
    "# Don't run the test with big numbers!\n",
    "if x<10000: print(f\"Test: {disc_log_test(y,g,x,p)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Block ciphers\n",
    "\n",
    "Notation:\n",
    "\n",
    "- `P`: Plaintext block (length n bits) \n",
    "- `C`: Ciphertext block (length n bits) \n",
    "- `K` : Key (length k bits)\n",
    "- `C = E(P, K)`: Encryption function \n",
    "- `P = D(C, K)`: Decryption function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
      "54\n"
     ]
    }
   ],
   "source": [
    "# Discarded code\n",
    "\n",
    "def z(n):\n",
    "    return [i for i in range(n)]\n",
    "print(z(10))\n",
    "\n",
    "n = 81\n",
    "zn = z(n)\n",
    "rel_primes = list(filter(lambda x: is_relative_prime(x,n), zn))\n",
    "print(len(rel_primes))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}