notebooks/TTM4135.ipynb

680 lines
18 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"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": 16,
"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\n",
"\n",
"> $\\phi(16) = \\phi(2^{4}) = 16*\\left(1-\\frac{1}{2}\\right)$"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"φ(21) = 12 since the 12 items in Z*21 are each relatively prime to 21.\n",
"Z*21 = [1, 2, 4, 5, 8, 10, 11, 13, 16, 17, 19, 20]\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 Zn.\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 = 21\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": [
"## Eulers Totient\n",
"- https://www.dcode.fr/euler-totient\n",
"- https://www.ibnus.io/eulers-totient-function-using-python/\n",
"\n",
"In general:\n",
"\n",
"> $\\varphi(n) = n \\prod_{p \\mid n} \\left( 1 - \\frac{1}{p} \\right)$\n",
"\n",
"Examples:\n",
"> $\\varphi(16) = \\varphi(2^{4}) = 16*\\left(1-\\frac{1}{2}\\right) = 8$\n",
"\n",
"> $\\varphi(1125) = \\varphi(33555) = 1125\\left(1-\\frac{1}{3}\\right)\\left(1-\\frac{1}{5}\\right) = 600$\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"𝜑(1125) = 600\n"
]
}
],
"source": [
"# Euler's Totient Function\n",
"# using Euler's product formula\n",
"def totient(n) :\n",
" result = n # Initialize result as n\n",
"\n",
" # Consider all prime factors\n",
" # of n and for every prime\n",
" # factor p, multiply result with (1 1/p)\n",
" p = 2\n",
" while(p * p <= n) :\n",
"\n",
" # Check if p is a prime factor.\n",
" if (n % p == 0) :\n",
"\n",
" # If yes, then update n and result\n",
" while (n % p == 0) :\n",
" n = n // p\n",
" result = result * (1.0 - (1.0 / (float) (p)))\n",
" p = p + 1\n",
"\n",
" # If n has a prime factor\n",
" # greater than sqrt(n)\n",
" # (There can be at-most one\n",
" # such prime factor)\n",
" if (n > 1) :\n",
" result = result * (1.0 - (1.0 / (float)(n)))\n",
"\n",
" return (int)(result)\n",
"\n",
"n = 1125\n",
"print(f\"𝜑({n}) = {totient(n)}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finding a generator of Z*p\n",
"- https://crypto.stanford.edu/pbc/notes/numbertheory/gen.html\n",
"- Lecture 2, page 19\n",
"\n",
"A generator of $Z^p$ is an element of order $p 1$\n",
"\n",
"To find a generator of Zp we can choose a value g and test it as follows:\n",
"1. compute all the distinct prime factors of p 1 and call them f1,f2,...,fr\n",
"2. then g is a generator as long as $g^{\\frac{p1}{fi}} \\neq 1 \\mod(p)$ for $i = 1,2,,...,r$"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Is 3 a generator for Z*4? True\n"
]
}
],
"source": [
"# Finding a generator g of Zp\n",
"def is_generator(p, g):\n",
" pf = prime_factors(p-1)\n",
" for f in pf:\n",
" if (g**((p-1)/f))%p == 1:\n",
" return False\n",
" return True\n",
"\n",
"g = 3 # Generator\n",
"p = 4 # Z*p\n",
"\n",
"print(f\"Is {g} a generator for Z*{p}? {is_generator(p, g)}\")"
]
},
{
"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": 7,
"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": 8,
"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 _ in range(k):\n",
" a = randint(2,n-2)\n",
" # Fermats theorem says that if a number p is prime then \n",
" # a^p1 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": 9,
"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": 10,
"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",
"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": 11,
"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": 12,
"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
}