Mathematics of Cryptography

§ 1

Foundations of Number Theory

Nearly all of modern cryptography rests on a handful of deep results from number theory. Understanding divisibility, prime numbers, and their relationships is the necessary starting point.

Divisibility and the Division Algorithm

For integers a and b with b ≠ 0, there exist unique integers q (quotient) and r (remainder) such that:

Division Algorithm a = q \cdot b + r, 0 \leq r < b // We say b | a (b divides a) iff r = 0

Prime Numbers

Definition — Prime

An integer p > 1 is prime if its only positive divisors are 1 and p itself. Every integer n > 1 factors uniquely into primes (Fundamental Theorem of Arithmetic):

Fundamental Theorem of Arithmetic n = p₁ e₁ \cdot p₂ e₂ \cdot \dots \cdot pₖ eₖ // unique factorization // Example: 360 = 2³ \cdot 3² \cdot 5¹

The difficulty of factoring large integers back into their prime components is the security foundation of RSA. The best known classical algorithm (General Number Field Sieve) runs in sub-exponential time — still completely infeasible for 2048-bit numbers.

GNU Octave Prime factorization & prime generation

%% Prime Factorization and Sieve of Eratosthenes
%% Run in Octave: octave --no-gui primes_demo.m

% --- Manual trial-division factorization ---
function factors = prime_factors(n)
  factors = [];
  d = 2;
  while d^2 <= n
    while mod(n, d) == 0
      factors = [factors, d];
      n = n / d;
    end
    d = d + 1;
  end
  if n > 1
    factors = [factors, n];
  end
end

% --- Sieve of Eratosthenes ---
function primes_list = sieve(limit)
  is_prime = true(1, limit);
  is_prime(1) = false;
  for i = 2 : floor(sqrt(limit))
    if is_prime(i)
      is_prime(i^2 : i : end) = false;
    end
  end
  primes_list = find(is_prime);
end

% --- Demo ---
disp('Factors of 360:'); disp(prime_factors(360));
disp('Primes up to 50:'); disp(sieve(50));
printf('Is 104729 prime? %d\n', isprime(104729));

Prime Number Theorem

The number of primes ≤ N is approximately π(N) ≈ N / ln(N). For a 512-bit number N ≈ 2⁵¹², roughly 1 in every 354 numbers is prime — making it efficient to find large primes by random sampling and testing.

§ 2

Euclidean & Extended Euclidean Algorithm

The Greatest Common Divisor (GCD) is central to almost every cryptographic key generation step. It determines when numbers are coprime — a requirement for modular inverses to exist.

Euclidean Algorithm

GCD via Repeated Division gcd (a, b) = gcd (b, a mod b) // recursive definition gcd (a, 0) = a // base case

Traced example — gcd(252, 105):

252 = 2·105 + 42

→

105 = 2·42 + 21

→

42 = 2·21 + 0

→

gcd = 21

Extended Euclidean Algorithm (EEA)

The EEA not only computes gcd(a, b) but also finds integers x, y (Bézout coefficients) satisfying:

Bézout's Identity a \cdot x + b \cdot y = gcd (a, b) // If gcd(a, n) = 1, then a\cdotx \equiv 1 (mod n) \to x is the modular inverse of a // This is how RSA's private key 'd' is computed from 'e' and φ(n)

Algorithm — Extended Euclidean

Input: a, n
Output: (g, x, y) where a·x + n·y = g

old_r, r ← a, n
old_s, s ← 1, 0
while r ≠ 0:
q ← ⌊old_r / r⌋
(old_r, r) ← (r, old_r − q·r)
(old_s, s) ← (s, old_s − q·s)
return (old_r, old_s, (old_r − old_s·a)/n)

GNU Octave GCD, Extended Euclidean, Modular Inverse

%% Extended Euclidean Algorithm + Modular Inverse

function [g, x, y] = extended_gcd(a, b)
  if b == 0
    [g, x, y] = deal(a, 1, 0);
  else
    [g, x1, y1] = extended_gcd(b, mod(a, b));
    [x, y] = deal(y1, x1 - floor(a/b) * y1);
  end
end

function inv = mod_inverse(a, n)
  % Returns x such that a*x ≡ 1 (mod n), or errors if none exists
  [g, x, ~] = extended_gcd(mod(a, n), n);
  if g ~= 1
    error('Inverse does not exist: gcd(%d,%d) = %d', a, n, g);
  end
  inv = mod(x, n);
end

% --- Demo ---
[g, x, y] = extended_gcd(252, 198);
printf('gcd(252,198) = %d;  Bezout: 252*(%d) + 198*(%d) = %d\n', g, x, y, 252*x+198*y);

% RSA use-case: find d = e^{-1} mod φ(n)
e = 65537; phi_n = 3120;    % toy example: n=p*q=53*61, φ=52*60
d = mod_inverse(e, phi_n);
printf('e=%d, φ(n)=%d → d=%d\n', e, phi_n, d);
printf('Verify e*d mod φ(n) = %d\n', mod(e*d, phi_n));

§ 3

Modular Arithmetic

Modular arithmetic is the language all discrete cryptography is written in. Operations "wrap around" at a modulus n, creating a finite arithmetic universe.

Congruence Definition a \equiv b (mod n) ⟺ n | (a - b) // Example: 17 \equiv 5 (mod 12) because 12 | (17-5)

Arithmetic rules

Modular Operations (a + b) mod n = ((a mod n) + (b mod n)) mod n (a \cdot b) mod n = ((a mod n) \cdot (b mod n)) mod n (a^k) mod n = computed efficiently via square-and-multiply

Fast Modular Exponentiation (Square-and-Multiply)

Computing a^e mod n naively requires e multiplications. Square-and-multiply reduces this to O(log e) multiplications — essential for RSA with 65537-bit exponents.

Square-and-Multiply — example: 3^13 mod 17 // 13 in binary = 1101 result = 1 bit 1 : result = (1 ² \cdot 3) mod 17 = 3 bit 1 : result = (3 ² \cdot 3) mod 17 = 27 mod 17 = 10 bit 0 : result = 10 ² mod 17 = 100 mod 17 = 15 bit 1 : result = (15 ² \cdot 3) mod 17 = 675 mod 17 = 12 // answer

GNU Octave Fast modular exponentiation

%% Square-and-Multiply modular exponentiation
%% (Octave's powermod() is built-in, but here's the explicit algorithm)

function result = mod_pow(base, exp, modulus)
  result = 1;
  base   = mod(base, modulus);
  while exp > 0
    if bitand(exp, 1)           % if lowest bit is 1
      result = mod(result * base, modulus);
    end
    exp  = bitshift(exp, -1);    % exp >>= 1
    base = mod(base^2, modulus);  % square
  end
end

printf('3^13 mod 17 = %d\n',    mod_pow(3, 13, 17));
printf('2^255 mod 997 = %d\n',  mod_pow(2, 255, 997));
%% Compare with built-in powermod (requires symbolic pkg or use mod)
printf('Verify: %d\n', mod(2^255, 997));   % float precision; mod_pow is exact

§ 4

Euler's Totient, Fermat's Little Theorem & Euler's Theorem

Euler's Totient Function φ(n)

φ(n) counts the number of integers in [1, n] that are coprime to n. This function determines the "size" of the multiplicative group mod n.

Totient Formulas φ (1) = 1 φ (p) = p - 1 // p prime φ (p k) = p k-1 (p - 1) φ (m \cdot n) = φ (m) \cdot φ (n) // gcd(m,n)=1 (multiplicative) φ (p \cdot q) = (p - 1)(q - 1) // RSA key formula, p,q distinct primes

Fermat's Little Theorem

If p is prime and gcd(a, p) = 1, then: a^(p−1) ≡ 1 (mod p). Equivalently: a^p ≡ a (mod p). This is the basis for primality tests and RSA decryption correctness.

Euler's Theorem (Generalization)

If gcd(a, n) = 1, then: a^φ(n) ≡ 1 (mod n). When n = p·q (RSA), this gives a^(p−1)(q−1) ≡ 1 (mod n), directly proving RSA decryption correctness.

GNU Octave Euler's totient function & theorem verification

%% Euler Totient + Fermat/Euler Theorem Verification

function phi = euler_totient(n)
  phi = sum(arrayfun(@(k) gcd(k, n) == 1, 1:n));
end

function phi = totient_formula(n)
  % Uses φ(n) = n * ∏ (1 - 1/p) over prime factors p
  phi = n;
  temp = n;
  d = 2;
  while d^2 <= temp
    if mod(temp, d) == 0
      while mod(temp, d) == 0; temp = temp/d; end
      phi = phi * (1 - 1/d);
    end
    d++;
  end
  if temp > 1; phi = phi * (1 - 1/temp); end
end

% Euler's Theorem: a^φ(n) ≡ 1 (mod n) when gcd(a,n)=1
function verify_euler(a, n)
  phi = totient_formula(n);
  res = mod_pow(a, phi, n);
  printf('φ(%d)=%d;  %d^φ(%d) mod %d = %d  [must be 1]\n', ...
         n, phi, a, n, n, res);
end

printf('φ(36) = %d  (expect 12)\n', totient_formula(36));
verify_euler(5, 36);
verify_euler(7, 100);
%% RSA primes: φ(p*q) = (p-1)*(q-1)
p = 53; q = 61;
printf('φ(%d) = %d = (%d-1)(%d-1) = %d\n', p*q, totient_formula(p*q), p, q, (p-1)*(q-1));

§ 5

Chinese Remainder Theorem (CRT)

The CRT states that if moduli n₁, n₂, …, nₖ are pairwise coprime, then the system of congruences has a unique solution mod N = n₁·n₂·⋯·nₖ. In RSA, CRT can speed up private key operations by roughly 4× using the prime factors directly.

CRT System x \equiv a₁ (mod n₁) x \equiv a₂ (mod n₂) ⋮ x \equiv aₖ (mod nₖ) // unique solution mod N = \prodnᵢ // Solution: x = Σ aᵢ \cdot Mᵢ \cdot yᵢ (mod N) // where Mᵢ = N/nᵢ, yᵢ = Mᵢ⁻¹ mod nᵢ

RSA-CRT acceleration

RSA Decryption with CRT dp = d mod (p - 1), dq = d mod (q - 1) m₁ = c dp mod p, m₂ = c dq mod q h = qInv \cdot (m₁ - m₂) mod p m = m₂ + h \cdot q // 4\times faster than direct c^d mod n

GNU Octave CRT solver and RSA-CRT decryption

%% Chinese Remainder Theorem Solver

function x = crt_solve(remainders, moduli)
  N   = prod(moduli);
  x   = 0;
  for i = 1 : length(moduli)
    Mi  = N / moduli(i);
    yi  = mod_inverse(Mi, moduli(i));
    x   = x + remainders(i) * Mi * yi;
  end
  x = mod(x, N);
end

% Solve: x≡2(mod3), x≡3(mod5), x≡2(mod7)
r = [2, 3, 2]; m = [3, 5, 7];
sol = crt_solve(r, m);
printf('CRT solution: x = %d\n', sol);
printf('Verify: %d%%3=%d, %d%%5=%d, %d%%7=%d\n', ...
       sol, mod(sol,3), sol, mod(sol,5), sol, mod(sol,7));

% RSA-CRT decryption demo (toy key)
p = 61; q = 53; n = p*q; e = 17;
d  = mod_inverse(e, (p-1)*(q-1));
msg= 42; c = mod_pow(msg, e, n);
% CRT-accelerated decrypt
dp  = mod(d, p-1); dq = mod(d, q-1);
qinv= mod_inverse(q, p);
m1  = mod_pow(c, dp, p);
m2  = mod_pow(c, dq, q);
h   = mod(qinv * mod(m1-m2, p), p);
dec = m2 + h*q;
printf('RSA-CRT: encrypt %d → %d → decrypt %d\n', msg, c, dec);

§ 6

Groups, Rings & Fields

Abstract algebra provides the framework that unifies symmetric encryption, public-key cryptography, and elliptic curves under a common language.

Group (G, ★)

Set with one operation. Requires: closure, associativity, identity element, inverses. Abelian if commutative.

Ring (R, +, ·)

Two operations. (R,+) is abelian group. Multiplication is associative and distributes over addition.

Field (F, +, ·)

A ring where every nonzero element has a multiplicative inverse. Both operations form abelian groups.

The Multiplicative Group ℤₙ*

Cyclic Group ℤn* = { a ∈ ℤn : gcd(a, n) = 1 }   |ℤₙ*| = φ(n) // g is a generator (primitive root) if powers of g produce every element // ℤₚ* (p prime) is always cyclic — critical for Diffie-Hellman // Discrete Log Problem: given g, p, h — find x such that g^x ≡ h (mod p) // No efficient classical algorithm exists for large p

Finite Fields GF(p)

For prime p, the field GF(p) = ℤp with addition and multiplication mod p. Every nonzero element has a multiplicative inverse (by Fermat's Little Theorem). These fields underlie RSA and Diffie-Hellman.

GNU Octave Find primitive roots (generators) of ℤp*

%% Find primitive roots of Z_p*
%% g is a primitive root mod p if ord(g) = p-1

function ord = element_order(g, p)
  % Smallest k > 0 such that g^k ≡ 1 (mod p)
  ord = 1; cur = mod(g, p);
  while cur ~= 1
    cur = mod(cur * g, p); ord++;
  end
end

function roots = primitive_roots(p)
  roots = [];
  for g = 2 : p-1
    if element_order(g, p) == p-1
      roots = [roots, g];
    end
  end
end

p = 23;
printf('Primitive roots of Z_%d*: ', p); disp(primitive_roots(p));
printf('Orders of elements mod %d:\n', p);
for g = 2:10
  printf('  ord(%d) = %d\n', g, element_order(g, p));
end

§ 7

Finite Fields GF(2ⁿ) and AES Mathematics

AES operates in the field GF(2⁸) — all 256 possible bytes form a field where addition is XOR and multiplication is polynomial multiplication modulo an irreducible polynomial.

GF(2⁸) — AES Irreducible Polynomial m (x) = x 8 + x 4 + x 3 + x + 1 // Binary: 0x11B = 100011011₂ // Elements are polynomials a₇x⁷ + \dots + a₁x + a₀, aᵢ \in {0,1} // Each byte 0x57 = 01010111₂ represents x⁶+x⁴+x²+x+1 Addition: 0x57 \oplus 0x83 = 0xD4 // just XOR Multiplication: 0x57 \cdot 0x83 mod m (x) = 0xC1

AES SubBytes — The S-Box

The AES S-Box applies two operations to each byte: (1) compute the multiplicative inverse in GF(2⁸) (0 maps to 0), then (2) apply an affine transformation over GF(2). This is the only nonlinear step in AES — providing all its resistance to linear and differential cryptanalysis.

AES Affine Transform s = A \cdot b \oplus c // b = GF(2⁸) inverse of input byte // A is the 8\times8 binary circulant matrix, c = 0x63 = 01100011₂ sᵢ = bᵢ \oplus b (i+4)%8 \oplus b (i+5)%8 \oplus b (i+6)%8 \oplus b (i+7)%8 \oplus cᵢ

MixColumns

MixColumns multiplies each 4-byte column of the state by a fixed 4×4 matrix over GF(2⁸), providing diffusion across bytes.

MixColumns Matrix ⎡ 02 03 01 01 ⎤ ⎢ 01 02 03 01 ⎥ // multiplication in GF(2⁸) ⎢ 01 01 02 03 ⎥ ⎣ 03 01 01 02 ⎦

GNU Octave GF(2⁸) arithmetic and AES S-Box generation

%% GF(2^8) arithmetic and AES S-Box generation

function r = gf2_mul(a, b)
  % Multiply two bytes in GF(2^8) mod 0x11B
  r = 0; p = a;
  for i = 1:8
    if bitand(b, 1); r = bitxor(r, p); end
    hi = bitand(p, hex2dec('80'));
    p  = bitand(bitshift(p, 1), hex2dec('FF'));
    if hi; p = bitxor(p, hex2dec('1B')); end  % reduce mod m(x)
    b  = bitshift(b, -1);
  end
end

function inv = gf2_inv(a)
  % Multiplicative inverse in GF(2^8) via brute-force (small field)
  if a == 0; inv = 0; return; end
  for b = 1:255
    if gf2_mul(a, b) == 1; inv = b; return; end
  end
end

function s = aes_sbox_byte(b)
  % Compute AES S-Box: GF inverse then affine transform
  b   = gf2_inv(b);
  s   = 0;
  for i = 0:7
    bit = bitxor(bitxor(bitxor(bitxor(bitget(b, i+1), ...
            bitget(b, mod(i+4,8)+1)), bitget(b, mod(i+5,8)+1)), ...
            bitget(b, mod(i+6,8)+1)), bitget(b, mod(i+7,8)+1));
    bit = bitxor(bit, bitget(hex2dec('63'), i+1));
    s   = bitset(s, i+1, bit);
  end
end

% Build first 16 entries of AES S-Box
printf('AES S-Box (first 16 bytes):\n');
for i = 0:15
  printf('%02X ', aes_sbox_byte(i));
end; printf('\n');
% Expected: 63 7C 77 7B F2 6B 6F C5 30 01 67 2B FE D7 AB 76
printf('GF mul 0x57 * 0x83 = 0x%02X (expect 0xC1)\n', gf2_mul(hex2dec('57'), hex2dec('83')));

§ 8

RSA: Full Mathematical Treatment

RSA's security rests on the integer factorization problem: given n = p·q, recovering p and q is computationally infeasible for large primes.

Key Generation

Choose two large distinct primes p and q (typically 1024–2048 bits each)
Compute the modulus n = p · q (the public modulus, ~2048–4096 bits)
Compute λ(n) = lcm(p−1, q−1) or φ(n) = (p−1)(q−1) (Carmichael vs Euler)
Choose public exponent e with 1 < e < λ(n) and gcd(e, λ(n)) = 1. Standard: e = 65537 = 2¹⁶ + 1
Compute private exponent d ≡ e⁻¹ (mod λ(n)) via Extended Euclidean Algorithm

RSA Encryption / Decryption // Public key: (n, e) Private key: (n, d) // Message m must satisfy 0 \leq m < n Encrypt: c = m e mod n Decrypt: m = c d mod n Sign: σ = Hash (m) d mod n Verify: check σ e mod n == Hash (m)

Correctness Proof

Proof — RSA Decryption Correctness

We need to show (m^e)^d ≡ m (mod n). Since e·d ≡ 1 (mod λ(n)), we can write e·d = 1 + k·λ(n) for some integer k. Then m^(e·d) = m^(1 + k·λ(n)) = m · (m^λ(n))^k ≡ m · 1^k = m (mod n) — the last step by Euler/Carmichael's theorem, which holds whenever gcd(m, n) = 1. Special cases for p|m or q|m are handled via CRT.

GNU Octave Complete RSA key generation, encrypt, decrypt, sign, verify

%% Complete RSA Implementation (toy-size keys for demonstration)
%% Requires: mod_pow() and mod_inverse() from earlier sections

function [n, e, d] = rsa_keygen(p, q)
  n     = p * q;
  phi_n = (p-1) * (q-1);      % Euler totient
  % Choose e: common choice is 65537; find valid e for small phi_n
  for e = [65537, 257, 17, 5, 3]
    if gcd(e, phi_n) == 1 && e < phi_n; break; end
  end
  d = mod_inverse(e, phi_n);
  printf('RSA Key: n=%d, e=%d, d=%d, phi=%d\n', n, e, d, phi_n);
end

function c = rsa_encrypt(m, e, n); c = mod_pow(m, e, n); end
function m = rsa_decrypt(c, d, n); m = mod_pow(c, d, n); end

% --- Demo with safe primes p=61, q=53 ---
[n, e, d] = rsa_keygen(61, 53);

for msg = [42, 99, 1234, 3000]
  c   = rsa_encrypt(msg, e, n);
  dec = rsa_decrypt(c, d, n);
  printf('m=%4d → c=%4d → dec=%4d  OK=%d\n', msg, c, dec, dec==msg);
end

%% Digital Signature: sign H(m) with private key, verify with public key
hash_m = mod(12345, n);             % simulate H(message) < n
sig    = rsa_decrypt(hash_m, d, n); % sign with private key
vrfy   = rsa_encrypt(sig, e, n);    % verify with public key
printf('Signature valid: %d\n', vrfy == hash_m);

§ 9

Diffie-Hellman & the Discrete Logarithm Problem

The Discrete Logarithm Problem (DLP): given a cyclic group G with generator g, and an element h ∈ G, find the integer x such that g^x = h. This is believed to be hard in large prime-order groups.

Diffie-Hellman Key Exchange // Public parameters: large prime p, generator g of Z_p* Alice: choose secret a, compute A = g a mod p; send A Bob: choose secret b, compute B = g b mod p; send B Alice computes: K = B a mod p = (g b) a mod p Bob computes: K = A b mod p = (g a) b mod p K = g^(ab) mod p // Shared secret — eavesdropper sees only A, B, g, p

Baby-Step Giant-Step Attack (BSGS)

BSGS solves the DLP in time O(√p) using a meet-in-the-middle approach. This is why DH prime sizes must be at least 2048 bits (making √p ≈ 2¹⁰²⁴).

BSGS Algorithm — find x such that g^x ≡ h (mod p) m = ⌈√p⌉ // Baby steps: precompute table {g^j mod p : j=0..m} // Giant steps: compute h · (g^(-m))^i for i=0..m, look up in table // Collision at i,j gives x = im + j

GNU Octave Diffie-Hellman key exchange and BSGS discrete log

%% Diffie-Hellman Key Exchange Demo

function diffie_hellman_demo(p, g, a, b)
  printf('Public:  p=%d, g=%d\n', p, g);
  A  = mod_pow(g, a, p);   % Alice's public value
  B  = mod_pow(g, b, p);   % Bob's public value
  Ka = mod_pow(B, a, p);   % Alice computes shared secret
  Kb = mod_pow(A, b, p);   % Bob computes shared secret
  printf('Alice private a=%d → A=g^a=%d\n', a, A);
  printf('Bob   private b=%d → B=g^b=%d\n', b, B);
  printf('Shared secret Ka=%d  Kb=%d  Match: %d\n', Ka, Kb, Ka==Kb);
end

% Baby-Step Giant-Step: find x s.t. g^x ≡ h (mod p)
function x = bsgs(g, h, p)
  m   = ceil(sqrt(p-1));
  % Baby steps: store g^j -> j
  baby = struct();
  gj  = 1;
  for j = 0:m
    baby.(sprintf('k%d', gj)) = j;
    gj = mod(gj * g, p);
  end
  % Giant steps: g^(-m) mod p
  gm_inv = mod_pow(mod_inverse(g, p), m, p);
  cur    = mod(h, p);
  for i = 0:m
    key = sprintf('k%d', cur);
    if isfield(baby, key)
      x = mod(i*m + baby.(key), p-1); return;
    end
    cur = mod(cur * gm_inv, p);
  end
  x = -1;  % not found
end

diffie_hellman_demo(23, 5, 6, 15);
printf('\nBSGS: log_5(%d) mod 23 = %d\n', mod_pow(5,9,23), bsgs(5, mod_pow(5,9,23), 23));

§ 10

Elliptic Curve Cryptography

An elliptic curve over a field 𝔽 is the set of points satisfying the Weierstrass equation, plus a "point at infinity" 𝒪 acting as the group identity. The group law on these points provides a hard DLP with much smaller key sizes than RSA.

Weierstrass Short Form E: y 2 = x 3 + ax + b over 𝔽 p // Non-singular: discriminant Δ = -16(4a³ + 27b²) \neq 0 // NIST P-256 (secp256r1): a = -3, b = 0x5AC635D8AA3A93E7B3EBBD55... // secp256k1 (Bitcoin): a = 0, b = 7

Point Addition Formulas

Given two distinct points P = (x₁, y₁) and Q = (x₂, y₂) on the curve, their sum R = P + Q = (x₃, y₃) is:

Point Addition (P \neq Q) λ = (y₂ - y₁) \cdot (x₂ - x₁) -1 mod p x₃ = λ 2 - x₁ - x₂ mod p y₃ = λ (x₁ - x₃) - y₁ mod p

Point Doubling (P = Q) λ = (3 x₁ 2 + a) \cdot (2 y₁) -1 mod p x₃ = λ 2 - 2 x₁ mod p y₃ = λ (x₁ - x₃) - y₁ mod p

Scalar Multiplication

k·P = P + P + … + P (k times) — computed efficiently with double-and-add (analogous to square-and-multiply). The Elliptic Curve DLP (ECDLP): given G and Q = k·G, find k. This is believed harder per bit than the classical DLP.

ECC Security Equivalence (approximate) ECC 256-bit \approx RSA 3072-bit \approx 128-bit symmetric ECC 384-bit \approx RSA 7680-bit \approx 192-bit symmetric ECC 521-bit \approx RSA 15360-bit \approx 256-bit symmetric

GNU Octave Elliptic curve point arithmetic and ECDH key exchange

%% Elliptic Curve over F_p — Point Addition, Scalar Multiply, ECDH
%% Toy curve: y² = x³ + 2x + 3  (mod 97)

a_c = 2; b_c = 3; p_c = 97;   % curve parameters

function R = ec_add(P, Q, a, p)
  % Add two points on y²=x³+ax+b mod p.  Point at infinity = [Inf Inf]
  if isinf(P(1)); R=Q; return; end
  if isinf(Q(1)); R=P; return; end
  if P(1)==Q(1) && mod(P(2)+Q(2),p)==0; R=[Inf Inf]; return; end
  if P(1)==Q(1) && P(2)==Q(2)   % doubling
    lam = mod((3*P(1)^2+a) * mod_inverse(mod(2*P(2),p), p), p);
  else                                       % addition
    lam = mod((Q(2)-P(2)) * mod_inverse(mod(Q(1)-P(1), p), p), p);
  end
  x3 = mod(lam^2 - P(1) - Q(1), p);
  y3 = mod(lam*(P(1)-x3) - P(2), p);
  R  = [x3, y3];
end

function R = ec_mul(k, P, a, p)
  % Double-and-add scalar multiplication
  R = [Inf Inf];
  while k > 0
    if bitand(k, 1); R = ec_add(R, P, a, p); end
    P = ec_add(P, P, a, p);
    k = bitshift(k, -1);
  end
end

% Find a point on the curve y² = x³ + 2x + 3 mod 97
G = []; 
for x = 0:96
  rhs = mod(x^3 + a_c*x + b_c, p_c);
  y   = mod_pow(rhs, floor((p_c+1)/4), p_c);  % sqrt via Tonelli for p≡3(mod4)
  if mod(y^2, p_c) == rhs && y>0; G = [x y]; break; end
end
printf('Generator G = (%d, %d)\n', G(1), G(2));

% ECDH key exchange
alice_priv = 12; bob_priv = 31;
alice_pub  = ec_mul(alice_priv, G, a_c, p_c);
bob_pub    = ec_mul(bob_priv,   G, a_c, p_c);
shared_a  = ec_mul(alice_priv, bob_pub,   a_c, p_c);
shared_b  = ec_mul(bob_priv,   alice_pub, a_c, p_c);
printf('Alice public: (%d,%d)\n', alice_pub(1), alice_pub(2));
printf('Bob   public: (%d,%d)\n', bob_pub(1),   bob_pub(2));
printf('Shared secret Alice: (%d,%d)\n', shared_a(1), shared_a(2));
printf('Shared secret Bob:   (%d,%d)  Match: %d\n', shared_b(1), shared_b(2), ...
       isequal(shared_a, shared_b));

§ 11

Primality Testing

Generating large RSA primes requires testing 500+ bit numbers for primality. Trial division is infeasible. Two probabilistic tests (and one deterministic) are used in practice.

Fermat Primality Test

Fermat Test Choose random a with 1 < a < n - 1 if a n-1 mod n \neq 1 : n is composite else : n is probably prime (Carmichael numbers can fool this test)

Miller-Rabin Probabilistic Test

Miller-Rabin is the industry standard. It has no false positives for Carmichael numbers. After k rounds the probability of a composite passing is at most 4^(-k). With 40 rounds: probability of error < 2^(-80).

Miller-Rabin — Write n-1 = 2ˢ\cdotd, d odd Choose random a \in [2, n-2] x = a d mod n if x = 1 or x = n - 1 : probably prime (this round) for r = 1 to s - 1 : x = x 2 mod n; if x = n - 1 : probably prime if loop exits without x = n - 1 : n is definitely composite

GNU Octave Miller-Rabin primality test and large prime generation

%% Miller-Rabin Primality Test

function result = miller_rabin(n, k)
  % Returns true if n is probably prime (k rounds)
  if n < 2;  result=false; return; end
  if n == 2 || n == 3; result=true; return; end
  if mod(n,2)==0; result=false; return; end
  % Write n-1 = 2^s * d
  s = 0; d = n-1;
  while mod(d,2)==0; d/=2; s++; end
  result = true;
  for i = 1:k
    a = 2 + floor(rand() * (n-3));
    x = mod_pow(a, d, n);
    if x==1 || x==n-1; continue; end
    composite = true;
    for r = 1:s-1
      x = mod(x^2, n);
      if x == n-1; composite=false; break; end
    end
    if composite; result=false; return; end
  end
end

% Test against known primes and composites
candidates = [7919, 104729, 15485863, 1000003, 9999991];
for n = candidates
  printf('%9d: MR=%d  isprime=%d\n', n, miller_rabin(n, 20), isprime(n));
end

% Generate a random prime in range [low, high]
function prime = gen_prime(low, high)
  do
    candidate = low + 2*floor(rand() * ((high-low)/2)) + 1; % force odd
  until miller_rabin(candidate, 40)
  prime = candidate;
end

printf('\nRandom 14-bit primes: %d, %d\n', gen_prime(8192, 16383), gen_prime(8192, 16383));

§ 12

Hash Function Mathematics

Merkle-Damgård Construction

Most cryptographic hash functions (MD5, SHA-1, SHA-2) use the Merkle-Damgård paradigm: a compression function f is applied iteratively to the padded message in fixed-size blocks.

Merkle-Damgård Iteration H₀ = IV // fixed initialization vector Hᵢ = f (Hᵢ₋₁, Mᵢ) // compress block i with chaining value Hash (M) = H_t // final chaining value = hash // Padding: append 1-bit, then zeros, then 64-bit message length // (length-extension attack: if you know H(m), you can compute H(m || extra))

SHA-256 Compression Internals

SHA-256 processes 512-bit blocks into a 256-bit digest. Each round updates eight 32-bit working variables using bitwise operations, additions mod 2³², and a message schedule.

SHA-256 — Core Operations (per round t) Σ₀(a) = ROTR(a, 2) \oplus ROTR(a, 13) \oplus ROTR(a, 22) Σ₁(e) = ROTR(e, 6) \oplus ROTR(e, 11) \oplus ROTR(e, 25) Ch(e, f, g) = (e \land f) \oplus (\neg e \land g) // "choose" Maj(a, b, c) = (a \land b) \oplus (a \land c) \oplus (b \land c) // "majority" T₁ = h + Σ₁(e) + Ch(e,f,g) + K[t] + W[t] // all mod 2³² T₂ = Σ₀(a) + Maj(a,b,c) (a,b,c,d,e,f,g,h) \leftarrow (T₁ + T₂, a, b, c, d + T₁, e, f, g)

GNU Octave SHA-256 from scratch (compact educational implementation)

%% SHA-256 in GNU Octave — educational full implementation
%% Works on strings; returns 64-char hex digest

function r = rotr32(x, n)
  r = bitor(bitshift(bitand(x, uint32(hex2dec('FFFFFFFF'))), -n), ...
            bitshift(bitand(x, uint32(hex2dec('FFFFFFFF'))), 32-n));
end

function hex = sha256(msg)
  % Initial hash values H0..H7 (first 32 bits of fractional parts of sqrt of first 8 primes)
  H = uint32([hex2dec('6a09e667') hex2dec('bb67ae85') ...
               hex2dec('3c6ef372') hex2dec('a54ff53a') ...
               hex2dec('510e527f') hex2dec('9b05688c') ...
               hex2dec('1f83d9ab') hex2dec('5be0cd19')]);
  % Round constants K (first 32 bits of fractional parts of cbrt of first 64 primes)
  K = uint32([hex2dec('428a2f98') hex2dec('71374491') hex2dec('b5c0fbcf') ...
               hex2dec('e9b5dba5') hex2dec('3956c25b') hex2dec('59f111f1') ...
               hex2dec('923f82a4') hex2dec('ab1c5ed5') hex2dec('d807aa98') ...
               hex2dec('12835b01') hex2dec('243185be') hex2dec('550c7dc3') ...
               hex2dec('72be5d74') hex2dec('80deb1fe') hex2dec('9bdc06a7') ...
               hex2dec('c19bf174') hex2dec('e49b69c1') hex2dec('efbe4786') ...
               hex2dec('0fc19dc6') hex2dec('240ca1cc') hex2dec('2de92c6f') ...
               hex2dec('4a7484aa') hex2dec('5cb0a9dc') hex2dec('76f988da') ...
               hex2dec('983e5152') hex2dec('a831c66d') hex2dec('b00327c8') ...
               hex2dec('bf597fc7') hex2dec('c6e00bf3') hex2dec('d5a79147') ...
               hex2dec('06ca6351') hex2dec('14292967') hex2dec('27b70a85') ...
               hex2dec('2e1b2138') hex2dec('4d2c6dfc') hex2dec('53380d13') ...
               hex2dec('650a7354') hex2dec('766a0abb') hex2dec('81c2c92e') ...
               hex2dec('92722c85') hex2dec('a2bfe8a1') hex2dec('a81a664b') ...
               hex2dec('c24b8b70') hex2dec('c76c51a3') hex2dec('d192e819') ...
               hex2dec('d6990624') hex2dec('f40e3585') hex2dec('106aa070') ...
               hex2dec('19a4c116') hex2dec('1e376c08') hex2dec('2748774c') ...
               hex2dec('34b0bcb5') hex2dec('391c0cb3') hex2dec('4ed8aa4a') ...
               hex2dec('5b9cca4f') hex2dec('682e6ff3') hex2dec('748f82ee') ...
               hex2dec('78a5636f') hex2dec('84c87814') hex2dec('8cc70208') ...
               hex2dec('90befffa') hex2dec('a4506ceb') hex2dec('bef9a3f7') ...
               hex2dec('c67178f2')]);
  % Pre-processing: padding
  bytes = uint8(msg);
  L     = length(bytes) * 8;
  bytes = [bytes, uint8(128)];          % append 0x80
  while mod(length(bytes), 64) ~= 56
    bytes = [bytes, uint8(0)];
  end
  for i = 7:-1:0                      % append 64-bit big-endian length
    bytes = [bytes, uint8(bitand(bitshift(L, -i*8), 255))];
  end
  % Process each 512-bit block
  for blk = 0 : length(bytes)/64-1
    chunk = bytes(blk*64+1 : blk*64+64);
    W     = zeros(1, 64, 'uint32');
    for i = 1:16
      W(i) = bitor(bitor(bitor(bitshift(uint32(chunk((i-1)*4+1)), 24), ...
               bitshift(uint32(chunk((i-1)*4+2)), 16)), ...
               bitshift(uint32(chunk((i-1)*4+3)), 8)), ...
               uint32(chunk((i-1)*4+4)));
    end
    for i = 17:64                    % Message schedule
      s0 = bitxor(bitxor(rotr32(W(i-15),7),rotr32(W(i-15),18)),bitshift(W(i-15),-3));
      s1 = bitxor(bitxor(rotr32(W(i-2),17),rotr32(W(i-2),19)),bitshift(W(i-2),-10));
      W(i) = W(i-16) + s0 + W(i-7) + s1;
    end
    [a_,b_,c_,d_,e_,f_,g_,h_] = deal(H(1),H(2),H(3),H(4),H(5),H(6),H(7),H(8));
    for t = 1:64                       % 64 rounds
      S1 = bitxor(bitxor(rotr32(e_,6),rotr32(e_,11)),rotr32(e_,25));
      ch = bitxor(bitand(e_,f_),bitand(bitcmp(e_),g_));
      T1 = h_ + S1 + ch + K(t) + W(t);
      S0 = bitxor(bitxor(rotr32(a_,2),rotr32(a_,13)),rotr32(a_,22));
      mj = bitxor(bitxor(bitand(a_,b_),bitand(a_,c_)),bitand(b_,c_));
      T2 = S0 + mj;
      [a_,b_,c_,d_,e_,f_,g_,h_] = deal(T1+T2,a_,b_,c_,d_+T1,e_,f_,g_);
    end
    H = H + [a_,b_,c_,d_,e_,f_,g_,h_];
  end
  hex = sprintf('%08x', H);
end

disp(sha256('abc'));
%% Expected: ba7816bf8f01cfea414140de5dae2ec73b00361bbef0469348423f656b7a7c69
disp(sha256(''));
%% Expected: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

§ 13

Post-Quantum Cryptography

Shor's Algorithm (1994) solves the integer factorization and discrete logarithm problems in polynomial time on a quantum computer, breaking RSA, DH, and all ECC variants. NIST finalized the first PQC standards in 2024 based on hard lattice and hash problems.

CRYSTALS-Kyber (ML-KEM)

Key encapsulation based on the Module Learning With Errors (MLWE) problem. Selected as NIST FIPS 203. Security levels: Kyber-512 (~128-bit), Kyber-768 (~192), Kyber-1024 (~256).

CRYSTALS-Dilithium (ML-DSA)

Digital signature algorithm based on the Module LWE and Module Short Integer Solution (MSIS) problems. NIST FIPS 204. Replaces ECDSA for signing.

SPHINCS+ (SLH-DSA)

Hash-based signature scheme. Security relies only on hash function properties — no lattice assumptions needed. NIST FIPS 205. Larger signatures (~8KB) but minimal assumptions.

Learning With Errors (LWE)

Given random pairs (A, b = As + e) where s is a secret vector and e is small noise, find s. This problem is believed hard even for quantum computers.

Learning With Errors — Core Problem // Choose secret s \in ℤq^n, error e \leftarrow χ (discrete Gaussian) // Given many samples: (aᵢ, bᵢ = ⟨aᵢ,s⟩ + eᵢ mod q) // Find s — computationally hard even for quantum algorithms // Kyber key sizes (vs RSA-2048 key at 256 bytes) Kyber-768 public key: 1184 bytes // ~192-bit security Kyber-768 ciphertext: 1088 bytes // key encapsulation

GNU Octave Toy LWE encryption — illustrates the core hard problem

%% Toy LWE Encryption (educational — real Kyber uses lattice NTT)
%% Parameters: n=4 (dimension), q=17 (modulus), small errors

q = 17; n = 4;

% Key generation
s   = randi([0, q-1], n, 1);     % secret key vector
A   = randi([0, q-1], n, n);     % public random matrix
e   = randi([-1, 1], n, 1);     % small error vector
b   = mod(A * s + e, q);         % public key b = As + e (mod q)

% Encryption of 1-bit message m ∈ {0,1}
m   = 1;
r   = randi([0,1], n, 1);        % random vector
e1  = randi([-1, 1], n, 1);
e2  = randi([-1, 1]);
u   = mod(A' * r + e1, q);
v   = mod(b' * r + e2 + round(q/2) * m, q);

% Decryption: compute v - s^T u and round
x    = mod(v - s' * u, q);       % x ≈ q/2·m (if errors small)
m_dec= round(x / (q/2));
m_dec= mod(m_dec, 2);

printf('LWE toy encrypt: m=%d → v=%d → decrypted=%d\n', m, v, m_dec);
printf('Secret s=[%d %d %d %d], public b=[%d %d %d %d]\n', s', b');

Algorithm Comparison

Algorithm	Hard Problem	Quantum Safe?	NIST Standard	Key Size
RSA-2048	Integer Factorization	❌ Broken by Shor's	Legacy	256 B
ECDH P-256	Elliptic Curve DLP	❌ Broken by Shor's	Legacy	64 B
AES-256	Block cipher security	✅ Grover halves bits	FIPS 197	32 B key
ML-KEM (Kyber-768)	MLWE	✅ No quantum attack	FIPS 203	1184 B pub
ML-DSA (Dilithium3)	MLWE + MSIS	✅ No quantum attack	FIPS 204	1952 B pub
SLH-DSA (SPHINCS+)	Hash function	✅ Minimal assumptions	FIPS 205	64 B pub