F4 Phantom

With F-4 Phantom II, we want to break the encryption! Please help us!!

‍‍‍nc 66.172.27.77 54979

We get a key-generation function, as seen below.

def gen_pubkey(e, p):
    assert gmpy.is_prime(p) != 0
    B = bin(p).strip('0b')
    k = random.randrange(len(B))
    k, l = min(k, len(B) - k), max(k, len(B) - k)
    assert k != l
    BB = B[:k] + str(int(B[k]) ^ 1) + B[k+1:l] + str(int(B[l]) ^ 1) + B[l+1:]
    q = gmpy.next_prime(int(BB, 2))
    assert p != q
    n = p*q
    key = RSA.construct((long(n), long(e)))
    pubkey = key.publickey().exportKey("PEM")
    return n, p, q, pubkey

So, , where . The density of primes (well, asymptotically) is , so we expect to be quite small.

We define . Now, we solve the following equation

Now, we know that the roots should be numerically close to and . So, we may simply call on them and check if they divide . Note that we need to do this for different hypotheses (different ).

import gmpy, base64
from Crypto.PublicKey import RSA

def solve(a):
   a2 = a ** 2/4
   L = gmpy.sqrt(a2 + N) # L > a/2
   C_1 = gmpy.next_prime(-a/2 + L)
   if N % C_1 == 0: return C_1
   C_2 = gmpy.next_prime(a/2 + L)
   if N % C_2 == 0: return C_2
   return False

# just some public key
pem_data1 = '''
MIGmMA0GCSqGSIb3DQEBAQUAA4GUADCBkAKBhyW0uQDhy7/eShVn6VQWQisC8sda
zl8xJmCe8eEPGNfDVMFgN7Ll4WdRLYE4T8dOvoMb99Cmjhym7Orb/qfP5q/BpTQy
5hr1w5RZVtYMqQ5I+M4qg7E8kkVhGJh/13bNwqEYvV6VNSCK+bgFWVGMelTQam31
Fiq6wsTLbIBCrLie52Cil1584QIEC8rPFw==
'''

pub =  RSA.importKey(base64.b64decode(pem_data1))
N = pub.n
m = (len(bin(N))+1)/2-2
print '[ ] Solving equations...'

for mm in range(m, m+3):
  for j in range(2, m/2+1):
    k = mm - j
    l = j
    retval = solve(2**(k-1) + 2**(l-1))
    if retval: 
      q = retval
      break
    retval = solve(2**(k-1) - 2**(l-1))
    if retval: 
      q = retval
      break

p = N / q
print '[+] Found p = {}... and q = {}...'.format(str(p)[:40], str(q)[:40])

λ python solve.py
[ ] Solving equations...
[+] Found p = 1381282812140256921334162381757305071089... and q = 1381282811302268925712750063033928508701...

Now, it turns out that . We cannot decrypt the message! If we look closely, we see that . Hmm… so, we can decrypt the message . And the message is the flag, which is constant. What if we sample two ciphertexts encrypting the same message but under different keys? Yes… so, basically, we have and . Then, we can combine it using CRT! This is easy!

import gmpy
from Crypto.Util.number import long_to_bytes


enc1 = 63188108518214820361083256140053967663112132356420859206347143811869148973386950682507343981284848159232100220605963292020722612854139075311063776086523677522160515093598202087755508512714885329251980275085813640578839753523295579661559881983237388178475574713319340090857313275483787001349560782781513895950569634984688753665
q1 = 33452043035349425454164058954054458228134102234436666511159820871022348004023966976017694615018111901825497223286811050478588079083387098695346723120647425399335947
p1 = 33452043035349425454164058954054813129854949698738692548175391185736387867969615080539316436404430573352896343366560167202569408949342999951142479436993639588965567
e1 = 3562731839

enc2 = 162331112890791758781057932826106636167735138703054666826574266304486608255768782351247144197937186145526374008317633308191215438756014724244242639042178681790803086016105986729819920057318114565244866632716401408212076926367974085730803964312385570202673351687846963322223962280999264618753873289802828995706526584379102468
q2 = 1381282812140256921334162381757305071089685391539884775199049048751523002134390007428038278188586333291047282664366863789424963612326970767160301759006158962675449
p2 = 1381282811302268925712750063033928508701820008572424411412024462643800411901779755548441592138469189655615818433739872652769585433967353091413641137354055362924329
e2 = 197840663

d1 = gmpy.invert(e1, p1-1)
d2 = gmpy.invert(e2, p2-1)

assert(gmpy.gcd(e1, p1-1) == 1)
assert(gmpy.gcd(e2, p2-1) == 1)

m = (pow(enc1, d1, p1) * p2 * gmpy.invert(p2, p1) + pow(enc2, d2, p2) * p1 * gmpy.invert(p1, p2)) % (p1*p2)
print long_to_bytes(m)

…prints:

********** great job! the flag is: ASIS{Still____We_Can_Solve_Bad_F4!}

Alright!

A fine OTP server

Connect to OTP generator server, and try to find one OTP.

nc 66.172.27.77 35156

We get the code for generating OTPs:

def gen_otps():
    template_phrase = 'Welcome, dear customer, the secret passphrase for today is: '

    OTP_1 = template_phrase + gen_passphrase(18)
    OTP_2 = template_phrase + gen_passphrase(18)

    otp_1 = bytes_to_long(OTP_1)
    otp_2 = bytes_to_long(OTP_2)

    nbit, e = 2048, 3
    privkey = RSA.generate(nbit, e = e)
    pubkey  = privkey.publickey().exportKey()
    n = getattr(privkey.key, 'n')

    r = otp_2 - otp_1
    if r < 0:
        r = -r
    IMP = n - r**(e**2)
    if IMP > 0:
    	c_1 = pow(otp_1, e, n)
    	c_2 = pow(otp_2, e, n)
    return pubkey, OTP_1[-18:], OTP_2[-18:], c_1, c_2

We note that . Because of this, the task is really trivial. We can simply compute the cubic root without taking any modulus into consideration. can solve this efficiently! This gives us:

ASIS{0f4ae19fefbb44b37f9012b561698d19}

Secured OTP server

Connect to OTP generator server, and try to find one OTP.
This is secure than first server 🙂

nc 66.172.33.77 12431

This is almost identical to the previous, but the OTP is chosen such that it will overflow the modulus when cubed.

    template_phrase = '*************** Welcome, dear customer, the secret passphrase for today is: '
    A = bytes_to_long(template_phrase + '00' * 18)

Note that we can write an OTP as , where . So, we have . The observant reader have noticed that if we remove , it will not wrap around . Now, we can mod out the terms containing , i.e., find . Now we are back in the scenario of previous challenge… so, we can use the method from before! Hence, . Finally, we get

ASIS{gj____Finally_y0u_have_found_This_is_Franklin-Reiter’s_attack_CongratZ_ZZzZ!_!!!}

DLP

You should solve a DLP challenge, but how? Of course , you don’t expect us to give you a regular and boring DLP problem!

nc 146.185.143.84 28416

The ciphertext is generated using the following function:

def encrypt(nbit, msg):
    msg = bytes_to_long(msg)
    p = getPrime(nbit)
    q = getPrime(nbit)
    n = p*q
    s = getPrime(4)
    enc = pow(n+1, msg, n**(s+1))
    return n, enc

We have that . We can write this as . Now, we can eliminate the ordo term by taking . Assuming that , we can determine by simple division. So, . In Python, we have that

Turns out this is the case, and, hence, we obtain

ASIS{Congratz_You_are_Discrete_Logarithm_Problem_Solver!!!}

unsecure ASIS sub-d

ASIS has many insecure sub-domains, but we think they are over HTTPS and attackers can’t leak the private data, what do you think?

So, we get a PCAP. The first thing we do is to extract data binwalk -e a.pcap. This generates a folder with all the certificates used.

Then, we look at the moduli.

#!/bin/bash
FILES=*.crt
for f in $FILES
do
  openssl x509 -inform der -in $f -noout -text -modulus 
done

This generates a file with all moduli. Let us try something simple! Common moduli! For each pair, we check if . If so, we have found a factor. Turns out two moduli have a common factor, so we can factor each of them and decrypt their traffic:

p1 = 146249784329547545035308340930254364245288876297216562424333141770088412298746469906286182066615379476873056564980833858661100965014105001127214232254190717336849507023311015581633824409415804327604469563409224081177802788427063672849867055266789932844073948974256061777120104371422363305077674127139401263621 
q1 = 136417036410264428599995771571898945930186573023163480671956484856375945728848790966971207515506078266840020356163911542099310863126768355608704677724047001480085295885211298435966986319962418547256435839380570361886915753122740558506761054514911316828252552919954185397609637064869903969124281568548845615791

p2 = 159072931658024851342797833315280546154939430450467231353206540935062751955081790412036356161220775514065486129401808837362613958280183385111112210138741783544387138997362535026057272682680165251507521692992632284864412443528183142162915484975972665950649788745756668511286191684172614506875951907023988325767 
q2 = 136417036410264428599995771571898945930186573023163480671956484856375945728848790966971207515506078266840020356163911542099310863126768355608704677724047001480085295885211298435966986319962418547256435839380570361886915753122740558506761054514911316828252552919954185397609637064869903969124281568548845615791

We can now generate two PEM-keys

d1 = gmpy.invert(e, (p1 - 1)*(q1 - 1))
key = RSA.construct((long(p1*q1), long(e), long(d1)))
f = open('privkey.pem','w')
f.write(key.exportKey('PEM'))
f.close()

Putting it into Wireshark, we obtain two images:

I totally agree.

Alice, Bob and Rob

We have developed a miniature of a crypto-system. Can you break it?
We only want to break it, don’t get so hard on our system!

This is McElice PKC. The ciphertexts are generated by splitting each byte in blocks of four bits. The following matrix is used as public key:

The ciphertext is generated as , which is a function from 4 bits to a byte. is an error (or pertubation) vector with only one bit set. This defines a map . So, each plaintext byte is two ciphertext bytes.

We can first create a set of codewords

P = numpy.matrix([[1, 1, 0, 0, 0, 1, 1, 0], [1, 1, 1, 1, 1, 1, 1, 1], [0, 1, 1, 0, 1, 1, 0, 0],[0, 1, 1, 1, 0, 0, 1, 0]])
image = {} # set of codewords
for i in range(0, 2**4):
    C = (numpy.array([int(b) for b in (bin(i))[2:].zfill(4)]) * P % 2).tolist()[0]
    image[int(''.join([str(c) for c in C]), 2)] = i

Then, go through each symbol in the ciphertext, flip all possible bits (corresponding to zeroing out ) and perform lookup in the set of codewords (compute the intersection between the Hamming ball of the ciperhext block and ).

f = open('flag.enc','r')
out = ''
for i in xrange(18730/2):
    blocks = f.read(2)
    j = ord(blocks[0])
    C1 = 0
    C2 = 0
    for i in range(0, 8):
        if j ^ 2**i in image:
            C1 = image[j ^ 2**i] << 4
    j = ord(blocks[1])
    for i in range(0, 8):
        if j ^ 2**i in image:
            C2 = image[j ^ 2**i]
    out += chr(C1+C2)
f=open('decrypted','w')
f.write(out)

Turns out it is a PNG:

Log in as admin to get the flag.

nc pki.hackable.software 1337

This was a fun challenge, because it required many steps to complete. The first thing we notice is that MD5 is being used:

def h(x):
    return int(hashlib.md5(x).hexdigest(), 16)

def makeMsg(name, n):
    return 'MSG = {n: ' + n + ', name: ' + name + '}'

def makeK(name, n):
    return 'K = {n: ' + n + ', name: ' + name + ', secret: ' + SECRET + '}'

def sign(name, n):
    k = h(makeK(name, n))
    r = pow(G, k, P) % Q
    s = (modinv(k, Q) * (h(makeMsg(name, n)) + PRIVATE * r)) % Q
    return (r*Q + s)

Naturally, we cannot register as . An immediate conclusion is that we need to some kind of collision attack. But how? We cannot exploit any collision in the account. What if we can create a collision in while keeping distinct?

from coll import Collider, md5pad, filter_disallow_binstrings
import os

def nullpad(payload):
    return payload + b'\x00' * (64 - len(payload) % 64)

def left_randpad(payload):
    random_string = os.urandom(64)
    return random_string[:len(payload) % 64] + payload

# nullpad the prefix so that it does not interfere
# with the two collision blocks
prefix = nullpad(b'K = {n: ')
suffix = b', name: groci, secret:'

collider = Collider()
collider.bincat(prefix)
collider.safe_diverge()
c1, c2 = collider.get_last_coll()
collider.bincat(suffix)
cols = collider.get_collisions()

# here are our collisions!
A = next(cols)
B = next(cols)

The above code creates a collision such that k = h(makeK(name, n)) generates the same value. The problem is that the server does not return the signature, but str(pow(sign(name, n), 65537, int(n.encode('hex'), 16))).

To get the signature, which is computed as , we need to perform some additional computations. Obviously, we could factor the modulus and compute and then obtain the signature as easy as … but factoring…

The probability that it is smooth enough is not very high. Also, the number is around digits so not a good candidate for msieve or other factoring software… so what then?

Well… what if and both are prime? Then, we can easily recover the signature as ! No need for time-consuming factoring! Embodied in Python, it could look like this:

import hashlib
from Crypto.Util.number import isPrime

def num2b(i):
    c = hex(i).strip('L')
    c = (len(c) % 2) * '0' + c[2:]
    return c.decode('hex')

# get the payloads
A = int(A[64:64*3].encode('hex'), 16)
B = int(B[64:64*3].encode('hex'), 16)

# PURE BRUTEFORCE, NOTHING FANCY HERE!
# append the same data to both payload until
# they resulting numbers are BOTH prime
mul = 2 ** 96
for i in range(1, 2**20, 2):
    if isPrime(AA + i) and isPrime(BB + i):
        print AA+i
        print BB+i
        break

As an example, we obtain

and

Putting this into action, we might do something like:

def connect(name, n):
   name_encoded = base64.b64encode(name)
   n_encoded = base64.b64encode(n)
   server = remote('pki.hackable.software', 1337)
   payload = 'register:' + name_encoded + ',' + n_encoded
   server.send(payload + '\n')
   return pow(int(server.recvline()), modinv(65537, n1-1), n1)

payload1 = nullpad('K = {n: ') + num2b(n1)
payload2 = nullpad('K = {n: ') + num2b(n2)

name = 'groci'
PORT = 1331

assert(h(makeK(name, payload1)) == h(makeK(name, payload2)))
assert(h(makeMsg(name, payload1)) != h(makeMsg(name, payload2)))

# Get first signature
sig = connect(name, payload1)
r1, s1 = sig / Q, sig % Q
z1 = h(makeMsg(name, payload1))

# Get second signature
connect(name, payload2)
sig = pow(int(server.recvline()), modinv(65537, n2-1), n2)
r2, s2 = sig / Q, sig % Q
z2 = h(makeMsg(name, payload2))

# Make sure we got a nonce re-use
assert(r1 == r2)

# OK, use standard nonce re-use technique...
delta = modinv(((s1 - s2) % Q), Q)
k = ( ((z1 - z2) % Q) * delta) % Q
r_inv = modinv(r1, Q)
PRIVATE = (((((s1 * k) % Q) - z1) % Q) * r_inv) % Q

# Now that we know the private key, 
# lets forge/sign the admin account!
name = 'admin'
n = 'snelhest'
k = h(makeK(name, n))

r = pow(G, k, P) % Q
s = (modinv(k, Q) * (h(makeMsg(name, n)) + PRIVATE * r)) % Q
sig = r*Q + s

# connect and login
name_encoded = base64.b64encode(name)
n_encoded = base64.b64encode(n)
server = remote('pki.hackable.software', 1337)
payload = 'login:' + name_encoded + ',' + n_encoded + ',' + base64.b64encode(num2b(sig))
server.send(payload + '\n')
print server.recvline()

We obtain the following signatures:

Now, with our forced nonce re-use, we can use standard techniques to obtain the private key and sign our account, which gives the flag
DrgnS{ThisFlagIsNotInterestingJustPasteItIntoTheScoreboard}!

Note to reader: I apologize for the excessive use of memes.

Integrity

Just a simple scheme.
nc 202.120.7.217 8221

The encrypt function in the scheme takes the username as input. It hashes the username with MD5, appends the name to the hash and encrypts with a secret key , i.e. . Then, the secret becomes . Notably, the first block contains the hash, but encrypted.

Recall how the decryption is defined:

So, what would happen if we input as username? Then, we have encrypted , but we want only . As mentioned before, the hash fits perfectly into a single block. So, by removing the ,  becomes the new  (which has no visible effect on the plaintext anymore!). Then, we have , which is all what we need.

The flag is flag{Easy_br0ken_scheme_cann0t_keep_y0ur_integrity}.

OTP1

I swear that the safest cryptosystem is used to encrypt the secret!

We start off analyzing the code. Seeing process(m, k) function, we note that this is actually something performed in with the mapping that, for instance, an integer is in binary, which corresponds to a polynomial . The code is doing in .

The keygen function repeatedly calls key = process(key, seed). The first value for key is random, but the remaining ones does not. seed remains the same. Define to be the key and the seed. Note that all elements are in . The first stream value is unknown.

So, we can compute the seed and key as and . The individual square roots exist and are unique.

def num2poly(num):
    poly = R(0)
    for i, v in enumerate(bin(num)[2:][::-1]):
        if (int(v)):
            poly += x ** i
    return poly

def poly2num(poly):
    bin = ''.join([str(i) for i in poly.list()])
    return int(bin[::-1], 2)

def gf2num(ele):
    return ele.polynomial().change_ring(ZZ)(2)

P = 0x10000000000000000000000000000000000000000000000000000000000000425L

fake_secret1 = "I_am_not_a_secret_so_you_know_me"
fake_secret2 = "feeddeadbeefcafefeeddeadbeefcafe"
secret = str2num(urandom(32))

R = PolynomialRing(GF(2), 'x')
x = R.gen()
GF2f = GF(2**256, name='a', modulus=num2poly(P))

f = open('ciphertext', 'r')
A = GF2f(num2poly(int(f.readline(), 16)))
B = GF2f(num2poly(int(f.readline(), 16)))
C = GF2f(num2poly(int(f.readline(), 16)))

b = GF2f(num2poly(str2num(fake_secret1)))
c = GF2f(num2poly(str2num(fake_secret2)))

# Retrieve partial key stream using known plaintexts
Y = B + b
Z = C + c

Q = (Z + Y**2)
K = (Y + Q).sqrt()

print 'flag{%s}' % hex(gf2num(A + K)).decode('hex')

This gives the flag flag{t0_B3_r4ndoM_en0Ugh_1s_nec3s5arY}.

OTP2

Well, maybe the previous one is too simple. So I designed the ultimate one to protect the top secret!

There are some key insights:

• The process1(m, k) function is basically the same as in previous challenge, but it computes the multiplication with the exception that elements are in this time.  We omitt the multiplication symbol from now on.
• The process2(m, k) function might look involved, but all that it does is to compute the matrix multplication between two matrices (with elements in ), i.e.,
• We start with matrices  and .
• Raising to a power yields has a closed form formula: .
• The nextrand(rand) function takes the integral value of , we call this and computes via a square-and-multiply type algorithm. In python, it would be

  def proc2(key):
      AN = A**gf2num(N)
      return key*AN+(AN+1)/(A+1)*B

Let us look at the nextrand(rand) function a little more. Let be the random value fed to the function. Once is computed, it returns

Define . Adding this to the above yields

.

So, . Note that given two elements of the key stream, all these elements are known. Once determined, we compute the (dicrete) to find . And once we have , we also have . Then, all secrets have been revealed!

From the plaintext, we can immediately get the key by XORing the first part of the plaintext with the corresponding part of the ciphertext. This gives .

R = PolynomialRing(GF(2), 'x')
x = R.gen()

GF2f = GF(2**128, name='a', modulus=num2poly(0x100000000000000000000000000000087))

A = GF2f(num2poly(0xc6a5777f4dc639d7d1a50d6521e79bfd))
B = GF2f(num2poly(0x2e18716441db24baf79ff92393735345))
G1 = GF2f(num2poly(G[1]))
G0 = GF2f(num2poly(0x2fe7878d67cdbb206a58dc100ad980ef))

U = B/(A+1)
Z = (G1+U)/(G0+U)
N = discrete_log(Z, A, K.order()-1)

We can then run the encryption (the default code) with the parameters fixed to obtain the flag flag{LCG1sN3ver5aFe!!}.

In this challenge, we are given a hash function, which essentially splits the input into chunks of ten bytes. It then appends six bytes of null data to it and XORs with the current state. The state is encrypted with the all-zero key and updated. This process is repeated until all blocks have been exhausted.

A simple case is to consider what happens when we hash ten bytes of null data.

If we can bring it back to the same state, we have a local collision. Here is an example of such a collision:

So, how did we create the above collision? Well, actually, it is not too complicated… first, note that we cannot control the last six bytes. Recall that . Let us reorder the function as follows:

If we force the trailing six bytes to null and then decrypt that block for different values on , since we can control . Equivalently, we can encrypt for different values of , where the trailing six bits will be the trailing six bits of .

Utilizing the birthday paradox, we can find a collision in the six bytes in time and space. This is done by putting about values of the encryption (or decryption) in a table. Then, we generate the decryptions (or encryptions, respectively) and look in the table.

In Python, it could look something like:

trailing_bytes_first = AES.new('\x00'*16).encrypt('\x00'*16)[-6:]

for i in range(0, 2**24):
    plaintext = os.urandom(10) + trailing_bytes_first
    data[AES.new('\x00'*16).encrypt(plaintext)[-6:]] = plaintext

for i in range(0, 2**24):
    ciphertext = os.urandom(10) + '\x00\x00\x00\x00\x00\x00'
    if AES.new('\x00'*16).decrypt(ciphertext)[-6:] in data:
        print [data[AES.new('\x00'*16).decrypt(ciphertext)[-6:]]], [ciphertext]

Two such values are

local_collision_blocks = ['!\x0e\x7f=\xc8\xb6$\xd4\x108\xfaY\xca4+.', '\xd9y\xb8\xa8RgG4\xc8U\x00\x00\x00\x00\x00\x00']

Finally, we generate a collision as follows

GIVEN = 'I love using sponges for crypto'
A = AES.new('\x00'*16).encrypt(local_collision_blocks[0])
B = AES.new('\x00'*16).decrypt(local_collision_blocks[1])
local_collision = '\x00'*10 + xorstring(local_collision_blocks[0][:10], AES.new('\x00'*16).encrypt('\x00'*16)) + xorstring(A,B)[:10] + xorstring(local_collision_blocks[1][:10], GIVEN[:10])

Today, I drew some replacement icons for some apps I use daily (i.e., Terminal, Wireshark, Hex Fiend and Hopper) since I did not like default ones. I release them under the GPL v.3.0 License, so feel free to use them in your OS 🙂

 

In action (with the dark grey and light grey terminal icons, respectively):

Here are some additional icons:

Consider a situation in which one part Alice wants to communicate with another part Bob over a discrete memoryless channel with crossover probability  while ensuring a low level of detection in the presence of a warden Willie, who observes the communication through another discrete memoryless channel with crossover probability . We assume that all parts are equally computationally bounded.

In this situation, we require that for non-zero . Let be a sequence of secret bits she wants to transmit.

NOTE: Alice and Bob may on forehand agree upon any strategy they desire, but the strategy is known also to Willie.

Communicating with a secret key

If Alice and Bob share a secret key , they may use it to pick a common encoding, i.e., an error-correcting code. This encoding in turn is used to encode the message, which of course must be able to correct the errors produced by the channel. Assuming that the correction capability is below the error rate of Willie’s channel, he cannot decode. Let bits be the length of the publicly transmitted sequence. An established result states that if bits are shared between Alice and Bob, they can transmit secret bits without Willie seeing them (with high probability).

Communicating without a secret key

Now, assume that Alice and Bob do not share a common key. Alice performs the following steps:

1. Picks a cryptographically secure random vector .
2. Computes the scalar product .
3. Sends the bits  over the channel.

This reduces to the problem of (LPN) and can be solved with the BKW algorithm (or similar) if is sufficiently low. In particular, if Bob receives

such sequences, or equivalently,

bits, he will be able to decode with high probability. Here we have exploited the piling-up lemma disregrading the fact that some bits in  are zero and does not contribute. For some probabilities  and natural number , the information is hidden from Willie. The information rate is determined as follows: , so

.

This bound can be improved upon by an increase in the number of parities.

Here is a more detailed paper on a similar topic. Another paper is here.

Suppose that we are given a system of modular equations

Trivially, this can be solved for unknown and using a simple Gaussian elimination step, i.e., writing the equations as

.

This is perfectly fine as long as the equations share a common modulus , but what about when the equations share unknowns but are defined under different moduli? Let us take a real-world scenario from the realm of cryptography.

Example: DSA using linear congruential generator (LCG)

The DSA (Digital Signature Algorithm) has two functions and . To sign a message (invoking ), the following steps are performed:

1. Let be a hash function and the message to be signed:
2. Generate a (random) ephemeral value where .
3. Compute . If , go to step 2.
4. Compute . If , go to step 2.
5. The signature is .

To verify as signature (invoking ), the below steps are performed:

1. If or  is not satisfied, reject the signature.
2. Compute .
3. Compute .
4. Compute .
5. Compute .
6. If , accept. Otherwise, reject.

In the case of using LCG as a psuedorandom-number generator for the values , two consecutively generated values (we assume one signature was generated right after the other) will be correlated as for some public parameters . Assuming that , we obtain equations

from the fourth step of .

Chinese Remainder Theorem (CRT)

A well-known theorem in number theory is the Chinese Remainder Theorem (commonly referred to with the acronym CRT), which deals with simple equations over different and pairwise-prime moduli. For instance,

which has the solution . In solving actual multivariate equations, we will hit a brick wall, as needed operations such as row reduction and modular inversion does not work.

The following code takes any moduli, which do not need to be relatively prime, and computes the CRT:

def crt_lll(m, n):
	p = n[0]
	q = n[1]
	a = m[0] % p
	b = m[1] % q
	L = p*q # can be tweaked

	A = matrix([
		[1, 1, 1/L, 0],
		[p, 0, 0, 0],
		[0, q, 0, 0],
		[0, 0, p*q/L, 0],
		[-a, -b, 0, L]
	])

	B = A.LLL()

	assert(B[-1, 0] == 0 and B[-1, 1] == 0)

	return B[-1, 2] * L


print crt_lll(
	[777 % 234, 777 % 567],
	[234, 567]
)

Thue’s lemma

Thue’s lemma states that given a modulus and an integer , there exists two integers such that

.

We can solve this with lattice reduction as follows:

def thues_lemma(a, mod):
    '''
            Find small elements x, y such that a * x = y
            and |x| < sqrt(mod) and |y| < sqrt(mod).
    '''
    R = Zmod(mod)
    A = matrix([
        [1,   a],
        [0, mod]
    ])
    A = A.LLL()
    return A[0][0], A[0][1]

# test
p = random_prime(2^512)
q = random_prime(2^512)
n = p * q
a = randint(0, 2**1024)
u, v = thues_lemma(a, n)
assert(a * u % n == v % n)

Lattice reduction

A lattice is a discrete linear space spanned by a set of basis vectors. For instance, and are two basis vectors that spans the space below. They are not unique in that sense, but evidently, they are the shortest basis vectors possible (such a basis can be found using the LLL algorithm). If the two column vectors and are basis vectors, then the corresponding lattice is denoted . Here, .

The problem of finding the shortest vector is called . Sloppily formulated: for a given lattice , find the shortest (in terms of Euclidan norm) non-zero vector. The answer to that question would be . Starting with different basis vectors, the problem will show be trickier.

A related problem is to find the closest vector, which is commonly called . In this scenario, we are given a vector in the vector space (but it might not be in ) and we want to find the closest vector in . There is a simple reduction from to by setting . There is also a reduction in the other direction, which involves extending the basis of to also include . This is called embedding.

Let us return to the example of DSA and different-moduli equations. The equations we got from before

can be formulated as basis vectors. In matrix form, we have

or (by embedding technique):

The idea is now to force the values corresponding to the guessed values to be close in the normed space. By adding additional constraints we may perform the following minimization:

where ,  and . Finding a closest approximation (preferably, using LLL/Babai) yields a solution to the DSA equations using LCG.

Example: Knapsack

The knapsack problem (sometimes subset-sum problem) is stated as follows. Given a weight and items/weights , find the sequence such that

It is not too hard to prove that this is NP-complete, but we omit the reduction here.

In a cryptographic setting, this can be used to encode data in the sequence . This is called the Merkle-Hellman public-key cryptosystem. It is easy to see that encryption actually is the encoding procedure mentioned above. However, the decryption procedure is a bit more involved; in fact it requires a special class of instances. If the weights can be transformed into a super-increasing sequence, the problem becomes trivial to retrieve the sequence .

Think of it this way. Assume that all weights sum to something smaller than the weight. Then, if the target sum (the ciphertext) is larger than (if not,  must be smaller assuming there is a unique solution), we know that . We can now remove and from the equation by solving for the remaining weights and a recalculated . This procedure can be repeated until all weights have been found (or ).

Merkle-Hellman provided a way to transform the super-increasing sequence into a hard one. We omit the details here, but the procedure can be found all over internet.

It turns out we can use lattice reduction for this problem. We create a basis matrix in the following way

and perform lattice reduction on it. Since the LLL/BKZ algorithm achieves a set of short basis, it will try to find something that makes the entries of the left-most column small. What this actually means is that it will find a sum of the rows that achieves a solution to the knapsack problem.

Of course, the solution must only contain values in . Depending on the instance, this may or may not be the case. So, how do we penalize the algorithm for choosing values outside the allow set?

A new approach*

Let us create a new basis matrix  in the following manner:

The algorithm performs the following steps:

1. Randomly (or deterministically) pick a guess on the number of non-zero values in the solution.
2. Update the matrix run LLL/BKZ on .
3. If a satisfiable reduced-basis vector is found, return it. Otherwise, goto 1.

It does not guarantee that an incorrect solution is penalized, but it increases the probability of it (reduces the set of ‘false’ basis vectors). We omit a formal proof, but think of it this way: Assume that is false solution and a reduced-basis vector of . In it also has to sum to the number of non-zero weights in the correct solution. Assume all false vectors appear randomly (they do not, but let us assume it!). Then for a correct guess of , the probability of the vector being is . If , this is , which is a noticeable reduction for most values of .

Here is an example of a CTF challenge which could be solved with the above technique.

* I do not know if this should attributed to someone.

Example: Ring-LWE

TBA

I have spent this week bathing in the ocean, catching a glimpse of sun and writing a library for post-quantum cryptography in Python. Check it out!

After fiddling with Project Euler, I started to wonder how to compute the sum of divisors in the most efficient way (albeit, most the problems in PE does not require any super-sophisticated solution). The trivial approach is doing division and adding up the divisors sequentially. This may be good enough for small , but as grows this approach gives an infeasible time complexity.

The next thing one may think of is to generate a list of prime numbers (using a sieve) and for each number try all the numbers in the list. This is faster than the trivial approach, but still requires . Since the sieve requires time , we do not need to do better than that. But can we match it? Turns out we can.

First, we need to establish a simple number theoretical fact: the is the divisor counting function which has the properties (and, thus, ) and .

The solution is quite simple; consider the first element in the prime set . That is . We split the range into even numbers and for even compute

.

Then for , compute , and so forth up to . The requires operations.

Running it for all powers and primes, we use in total operations. Now, we have . So the complexity is . This is a harmonic series and treating it as an integral, we have . Of course, we do only generate primes up to , so the final complexity is . So we are done 🙂

As always, let us look at a Python representation:

def sum_of_divisors(N):
    N += 1 # extend limit to inclusive
    primes = primetools.sieve(N)
    numbers = [1]*N
    for q in primes:
        tmp = {}
        mult = q
        while mult <= N:
            for i in range(mult, N, mult):
                if i in tmp: tmp[i] += mult
                else: tmp[i] = 1 + mult
            mult *= q
        for j in tmp: numbers[j] *= tmp[j]
    return numbers[2:]

Below is a plot of the timing of the code. Since the code stores data in memory, some undesired effects comes into play for largers values (due to the need of shuffling data back and forth) and is therefore not present in the graph.

Bounded by the sieving complexity, the algoritm cannot be any faster using the sieve. Can we do better? Note that the algorithm returns a list of sums. If we are really interested in the sum of sums, we may instead compute , which is and uses basically no memory. This can be written as

def sum_of_sum_of_divisors(N):
    div_sum = 0
    for i in range(1, N):
        div_sum += i * (N / i)
    return div_sum

How does it work? For each , it computes the cardinality of the set of numbers divisible by ,i.e., , which is exactly . This contributes with for every element, so the total contribution is . Doing it for all possible divisors, we end up with the sum of sum of divisors.

Turns out, we can actually do it in . Here is how:

def square_root_time_sum_of_sum_of_divisors(N):
    global dcached
    if N in dcached: return dcached[N]
    div_sum, i = 0, 1
    q = int(math.sqrt(N))
    while i <= q:
        div_sum += (i * (N / i))
        i += 1
    i = 1
    while i <= N / (q + 1):
        m = N / i
        k = N / (i + 1)
        div_sum += (i * (m * (m + 1) - k * (k + 1)) / 2)
        i += 1
    dcached[N] = div_sum
    return div_sum

Consider the below example. In the left table, we have the sum and in the right table where :

Clearly, it is the same table/matrix but transposed, so the sum of elements are the same. As we notice in the matrices above, their rows become sparser for increasing . Think about it! What if we were to compute only the bold numbers and adding them? Then we need only two function calls each for left and right. The point where we stop computing left and start to compute right can be chosen arbitrarily, but the optimal is . Then .

Example: Project Euler 23
This is a simple problem, where we are supposed to compute all numbers which cannot be expressed as a sum of two abundant numbers. Using the above code, we can solve the problem in about 100 ms using pypy:

def abundant_numbers(N):
    numbers = sum_of_divisors(N)
    abundants = []
    for i in range(0, N):
        if 2*i < numbers[i]: abundants.append(i)
    return abundants[1:]

N = 28123

abundants = abundant_numbers(N)
s = [0]*N

for x in abundants:
    for y in abundants:
        ls = x + y
        if ls < N: s[ls] = 1
        else: break

print sum(i for i in range(0, N) if not s[i])

Example: Project Euler 439
This problem is quite nice, because it requires a very well-thoughtout approach for it to be feasible solving it (it is also among the hardest problems on PE). Here is the question (notation slightly altered):

Let  be the sum of all divisors of . We define the function . For example,

You are given that and . Find .

So, running our algorithm would be fine for lower values, but if we would try to store values in memory, we will run into a wall very fast (this is about TBs assuming that each number is a 64-bit integer – note the modular reduction which can be done in every step). Clearly, we need to use the memoryless approach.

Let us first disregard the fact that only holds when . Then . The case of , we will have

OK, so what about ? We can estimate the how many numbers in the product for which it that holds . Take, . We have , so

In particular, . Let us consider the case . The (negative) contribution from this set

.

Running over all possible gcds (all numbers ), we get

The following code computes the smaller results in matter of milliseconds, but for , it takes too long (as we are aiming for the sub-minute mark):

cached = {}
dcached = {}
mod = 10**9

def S(N):
    global cached
    if N in cached: return cached[N]
    dv = square_root_time_sum_of_sum_of_divisors(N)
    S_sum, i = 0, 2
    while i <= N:
        # max recursion depth log2(10^11) ~ 33
        S_sum += (i * S(N / i)) % mod
        i += 1
    cached[N] = (dv * dv - S_sum) % mod
    return cached[N]

OK. Further optimization needed… note that

Using the same approach as before, we can derive that

where is the Möbius function, defined

We can compute this in time. The Möbius function can be pre-computed in like so:

def mobius(N):
    L = [1]*(N+1)
    i = 2
    while i <= N:
        # if i is not sq. free, then no multiple is
        if L[i] == 1 or L[i] == -1:
            L[i] = -L[i]
            sq = i * i
            j = i + i
            while j <= N:
                # we use 2 to mark non-primes
                if abs(L[j]) == 2:
                    L[j] = -L[j]
                else:
                    L[j] = -L[j]*2
                j += i
            j = sq
            while j <= N:
                L[j] = 0
                j += sq
        i += 1
    i = 2
    while i <= (N):
        if abs(L[i]) == 2:
            L[i] = L[i] / 2
        i += 1
    return L