В этой статье будет описано, как один чересчур поверивший в себя программист пытался взломать один из фундаментальных алгоритмов криптографии. Эта статья признана огородить других от подобных попыток, или наоборот, заинтересовать новых смельчаков для подобной авантюры. Сначала я опишу суть алгоритма на простом коде, затем перечислю методы и идеи, которыми я пытался его взломать.

Об алгоритме

Данный алгоритм лежит в основе популярных криптовалют, таких как Solana и Ton, а также используется для шифрования http трафика с помощью TLS. Его гипотетический взлом позволяет проводить транзакции от имени чужих крипто кошельков, что делает потенциального взломщика миллиардером. В основе алгоритма лежит эллиптическая кривая Twisted Edward Curve25519.

Про математику эллиптических кривых написаны целые книги, но в её основе лежат всего две операции - сложение двух точек и удвоение одной точки. Ниже приведена простейшая реализация этих операций

Q = 2**255 - 19  # Модуль конечного поля, в котором производятся вычисления
L = 2**252 + 27742317777372353535851937790883648493  # Общее количество точек на кривой
d = 37095705934669439343138083508754565189542113879843219016388785533085940283555


def inv(x):
    return pow(x, Q-2, Q) % Q


def add_point(pt1, pt2):
    x1, y1 = pt1
    x2, y2 = pt2
    x3 = (x1 * y2 + y1 * x2) * inv(1 + d * x1 * x2 * y1 * y2) % Q
    y3 = (y1 * y2 + x1 * x2) * inv(1 - d * x1 * x2 * y1 * y2) % Q
    return x3, y3


def double_point(pt):
    x, y = pt
    x3 = 2*x*y * inv(y**2 - x**2) % Q
    y3 = (x**2 + y**2) * inv(2 + x**2 - y**2) % Q
    return x3, y3


# базовая точка, эквивалент единицы в стандартной математике
bpt = (15112221349535400772501151409588531511454012693041857206046113283949847762202,
       46316835694926478169428394003475163141307993866256225615783033603165251855960)

На основе операций сложения и удвоения точек можно вывести дополнительную операцию - умножения точки на произвольное число через алгоритм double-and-add.

def scalarmult(pt, scalar):
    if scalar == 0:
        return 0, 1
    _ = double_point(scalarmult(pt, scalar >> 1))
    return add_point(_, pt) if scalar & 1 else _

Теперь самое главное, за счёт чего достигается невероятная крипто-стойкость алгоритма. Скаляр в данной операции может быть любым в диапазоне от 1 до L, при этом операция scalarmult вычисляется за O(log n). Приватным ключом здесь служит скаляр, а публичным - точка, вычисленная через scalarmult(bpt, scalar).

random_bytes = os.urandom(32)
num = int.from_bytes(random_bytes, byteorder='little')
private_key = num % L
public_key = scalarmult(bpt, private_key)

Восстановить скаляр из точки считается вычислительно невозможным для современных компьютеров. Наиболее эффективный из известных алгоритмов подбора ключа: "Baby-step giant-step". Этот алгоритм сначала вычисляет первые sqrt(n) точек и сохраняет их в хэш таблицу, после чего выполняет sqrt(n) операций add_point, каждый раз прибавляя к точке scalarmult(bpt, sqrt(n)). Для curve25519 ему требуется не менее 2^127 памяти и 2^128 операций, что на порядки превосходит современные возможности компьютеров.

Шифрование на основе алгоритма

В этом разделе я опишу практическую реализацию подписи сообщения с помощью ED25519.

В данном стандарте приватный ключ не используется в эллиптических кривых напрямую, вместо этого на основе приватного ключа вычисляется хэш sha512, который затем разделяется на два 32 байтовых ключа.

private_key = os.urandom(32)
key_hash = hashlib.sha512(private_key).digest()
a_bytes = key_hash[:32]
prefix = key_hash[32:]

Первый ключ преобразуется в int, после чего над ним выполняется несколько битовых операций для повышения криптостойкости. Первый бит скаляра устанавливается в 1, а последние три устанавливаются в ноль, что защищает от атак на малые подгруппы и от атак малого логарифма.

def bytes_to_clamped_scalar(s):
    a_unclamped = int.from_bytes(s, byteorder='little')
    AND_CLAMP = (1<<254) - 1 - 7
    OR_CLAMP = (1<<254)
    a_clamped = (a_unclamped & AND_CLAMP) | OR_CLAMP
    return a_clamped
a = bytes_to_clamped_scalar(a_bytes)

На основе второго ключа и подписываемого сообщения вычисляется одноразовое число r, которое защищает основной ключ от утечек. Сигнатурой подписываемого сообщения становится точка на основе скаляра r и скаляр S, представляющий собой произведение одноразового хэша с приватным ключом и суммы с числом r.

def Hint(m):
    h = hashlib.sha512(m).digest()
    return int.from_bytes(h, byteorder='little')

def to_bytes(pt):
    return pt[0].to_bytes(32, 'big') + pt[1].to_bytes(32, 'big')

public_key = scalarmult(bpt, a)
pk_bytes = to_bytes(public_key)
message = b'hello world'
r = Hint(prefix + message) % L
R = scalarmult(bpt, r)
R_bytes = to_bytes(R)
h = Hint(R_bytes + pk_bytes + message)
S = (r + h * a) % L
S_bytes = S.to_bytes(32, 'big')
signature = R_bytes + S_bytes

Валидация сигнатуры производится принимающей стороной с помощью вычисления точки r+h*a двумя разными способами и их сравнения.

def from_bytes(pt):
    return int.from_bytes(pt[:32], 'big'), int.from_bytes(pt[32:], 'big')
  
S_bytes = signature[64:]
R_bytes = signature[:64]
R = from_bytes(R_bytes)
A = from_bytes(pk_bytes)
S = int.from_bytes(S_bytes, 'big')
h = Hint(R_bytes + pk_bytes + message)
v1 = scalarmult(bpt, S)
v2 = add_point(R, scalarmult(A, h))
assert v1 == v2

Из трех скаляров у принимающей стороны доступен только h, другие две переменные он имеет в виде точек, восстановление из которых скаляра считается вычислительно невозможным. Гипотетически можно разделить скаляр на h, используя scalarmult(S, inv(h)), в результате получим (r/h + a), что также не позволяет однозначно вычислить a или r.

Для подделки сигнатуры необязательно иметь исходный ключ и prefix, достаточно знать параметр a. Ниже реализация фейковой сигнатуры

def fake_sign(m, a, pk):
    prefix = b'12345'
    r = Hint(prefix + m)
    R = Base.scalarmult(r)
    R_bytes = R.to_bytes()
    h = Hint(R_bytes + pk + m)
    S = r + h * a
    return R_bytes + scalar_to_bytes(S)

В результате задача взлома ED25519 сводится к восстановлению скаляра k на основе точки scalarmult(bpt, k).

Первые попытки взлома

Первое что я сделал это дополнил базовые операции над эллиптической кривой. На основе сложения, удвоения и умножения можно вывести вычитание, деление и отрицание.

def rev_point(pt):
    return (-pt[0]) % Q, pt[1]


def sub_point(pt1, pt2):
    pt2 = rev_point(pt2)
    return add_point(pt1, pt2)


def l_inv(val):
    return pow(val, L-2, L)


def div_point(pt, scalar):
    return scalarmult(pt, l_inv(scalar))

Отрицание точки работает так, что точка kG превращается в (L-k)G, где k это скаляр, L - размер конечного поля, а G - базовая точка bpt. Например, при L=13 и k=2, -k будет равно 11. Сложение с отрицанием точки эквивалентно вычитанию её оригинала. Деление в свою очередь эквивалентно умножению точки на инверсию скаляра, инверсия здесь производится по размеру поля L.

После дополнения операций я решил, что самым простым способом узнать скаляр k на основе точки kG будет уменьшение точки через математические операции до вычисленного заранее диапазона. Например, если вычислить первые 2^40 точек и довести kG в этот диапазон, сохраняя все шаги, то можно затем выполнить операции в обратном порядке и восстановить k.

Уменьшение скаляра через операции на кривой

Понимая, что все формулы на кривой работают по правилам модульной арифметики, я написал упрощенную версию всех операций.

def inv(val):
    return pow(val, L-2, L) % L


def rev(val):
    return (L - val) % L


def div(v1, v2=None, iv2=None):
    if not iv2:
        iv2 = inv(v2)
    return (v1 * iv2) % L


def add(v1, v2):
    return (v1 + v2) % L


def mult(v1, v2):
    return (v1 * v2) % L


def sub(v1, v2):
    return (v1 - v2) % L


def equal(v1, v2):
    return (v1 % L) == (v2 % L)

Я перебрал множество комбинаций всех доступных операций. Результат показал, что любое изменение числа без знания его свойств подвержено нормальному распределению. Чтобы один раз попасть в диапазон 2^192, нужно выполнить в среднем 2^64 операций. В ходе перебора я узнал несколько основных свойств кривой

• Деление действует на скаляр практически хаотично

• kG + rev(kG) = 0 (казалось бы очевидно)

• Деление на два чётного скаляра гарантированно уменьшает его в два раза

Третье свойство показалось мне самым интересным. С его помощью задача по восстановлению числа k из точки kG сводится к определению четности числа k из точки kG. Если число четное, то мы делим его на 2, а если нечётное, то вычитаем единицу и делим на два. Таким образом, можно быстро попасть в диапазон известных точек и восстановить приватный ключ с помощью операций в обратном порядке.

random_bytes = os.urandom(32)
num = int.from_bytes(random_bytes, byteorder='little')
first_num = num
steps = []
while num > 1:
    if num % 2 != 0:
        num -= 1
        steps.append('-')
    num = num // 2
    steps.append('/')

steps.reverse()
rec_num = 1
for step in steps:
    if step == '-':
        rec_num += 1
    else:
        rec_num *= 2
        
assert first_num == rec_num

Далее основной проблемой стало то, как определить четность точки. С этой проблемой я застрял на многие недели.

Методы определения четности скаляра

Первое, что я сделал, это определил рамки, какая точность считается приемлемой. Стандартная вероятность угадать четность точки 50%. Если повысить эту вероятность до 90%, то шанс что прогноз окажется верным 255 раз подряд будет ~2.15e-12. Это мало, но при достаточно мощном сервере уже позволяет сломать один целевой ключ за пару часов. Таким образом, задача сводится к тому, чтобы написать функцию оракула, которая угадывает четность точки с шансом не менее 90%.

Сначала я попытался найти закономерности в байтах самой точки, побитовый анализ показал, что их изменение не имеет прямых корреляций с четностью. После я опробовал нейросети на tensorflow.

import numpy as np
from tensorflow.keras.models import Sequential, load_model
from tensorflow.keras.layers import Dense, Dropout
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.callbacks import EarlyStopping
data = np.array(data)
np.random.shuffle(data)
features = data[:, :-1].astype(np.float64)
labels = data[:, -1].astype(np.float64)

max_chunk_value = (1 << 32) - 1
features = features / max_chunk_value
model = Sequential([
        Dense(512, input_dim=16, activation='relu'),
        Dropout(0.5),
        Dense(256, activation='relu'),
        Dropout(0.5),
        Dense(128, activation='relu'),
        Dropout(0.5),
        Dense(1, activation='sigmoid')
    ])
    model.compile(optimizer=Adam(learning_rate=0.001), loss='binary_crossentropy', metrics=['accuracy'])
history = model.fit(
    features,
    labels,
    validation_split=0.5,
    epochs=100000,
    batch_size=1024,
    verbose=1,
)

Нейросеть быстро обучалась на тестовой выборке, но на случайных данных вероятность всё равно стремилась к 50%, что подтвердило отсутствие явных закономерностей.

После этого я прочитал несколько книг, посвященных эллиптическим кривым, особенно интересной мне показалась эта. Также в одной иностранной статье я нашёл информацию об альтернативном способе деление точки на два, который использует вычисления на кривой Вейерштрасса и создает несколько корней, только один из которых лежит на кривой.

Я решил поискать закономерность между корнями при делении на два от четности скаляра. Для этого я написал реализацию этого алгоритма на python

def spec_half_element(pt):
    def is_square(n):
        return pow(n, (Q - 1) // 2, Q) == 1
    a = -1
    x, y = pt
    u = (a - d) * (1 + y) * inv(1 - y) % Q
    w = 2 * inv(x) % Q
    A = 2 * (a + d) % Q
    B = (a - d) ** 2 % Q
    Ap = -2 * A % Q
    Bp = (A ** 2 - 4 * B) % Q
    us = (u * 4) % Q
    ws = (w * 2) % Q
    wp = sqrt_mod(us, Q) # First halving (reversing psi2)
    if wp is None:
        return None  # Point is not divisible by 2
    up = (us - Ap - wp * ws) * inv(2) % Q
    if not is_square(up):
        up = Bp * inv(up) % Q
        wp = -wp % Q
    w = sqrt_mod(up, Q) # Second halving (reversing psi1)
    if w is None:
        return None  # Point is not divisible by 4
    u = (up - A - w * wp) * inv(2) % Q
    u2 = B * inv(u) % Q
    w2 = -w % Q
    x1 = 2 * inv(w) % Q  # Convert back to Edwards coordinates
    y1 = (-a + d + u) * inv(a - d + u) % Q
    x2 = 2 * inv(w2) % Q
    y2 = (-a + d + u2) * inv(a - d + u2) % Q
    return [(x1, y1), (x2, y2)]

Проверка показала, что первый и второй корни с равной вероятностью оказываются внутри кривой, а также что у них нету явного периода или корреляции с четностью скаляра. Из полезного я научился вычислять новый вид точек, которых нет на кривой, но при этом операции над ними, такие как удвоение, возвращают корректный результат, как если бы точка была на кривой. Также изучение литературы показало, что точки из кривой Эдвардса можно перемещать на кривые Вейерштрасса и Монтгомери, после чего производить вычисления уже по их формулам.

Завершив попытки определять четность скаляра, я решил попробовать принципиально новый подход - свести одну из рекурсивных формул в закрытую форму, которая может быть решена за полиномиальное время.

Удвоение точки на кривой Монтгомери

Изучая альтернативные кривые я нашёл интересную формулу на кривой Монтгомери. Эта формула представляла собой удвоение точки от одной координаты.

def to_montgomery(pt):
    x, y = pt
    u = (1 + y) * inv(1 - y) % Q
    v = (1 + y) * inv((1 - y) * x) % Q
    return u, v

def to_edward(mpt):
    u, v = mpt
    x = u * inv(v) % Q
    y = (u - 1) * inv(u + 1) % Q
    return x, y
  
def m_double_element_xonly(x):  # A = 486662
    numerator = (x ** 2 - 1) ** 2 % Q
    denominator = 4 * x * (x ** 2 + A * x + 1) % Q
    x3 = (numerator * inv(denominator)) % Q
    return x3

Данная формула позволяет удваивать точку и зависит при этом всего от одной переменной. Я сделал проверку на малой кривой со схожими свойствами и выяснил, что многократным удвоением базовой точки можно получить любую другую точку кроме нуля. Следовательно, если преобразовать эту формулу в закрытую форму, то можно будет её обратить и восстановить количество удвоений.

def double_recursive_form(x, k):
    for _ in range(k):
        x = ...

def double_closed_form(x, k):
    kx = ...

После этого я начал пытаться упростить эту формулу. Множество преобразований привели меня к такой форме.

def m_double_element_xonly(x):
    s = (x + inv(x) + 486662) * inv(4)
    x = (s + 14802493890*inv(s) - 243331) % Q
    return x

Здесь производится одна замена и все числа являются максимально малыми. Инверсию x можно представить как степень x^(Q-2), в обычной математике это является делением, но возможно в модульной у неё есть скрытые свойства. Тем не менее, несмотря на визуальную простоту, дальнейших преобразований этой формы я найти не смог.

Одной из интересных попыток был разворот формулы в обратную сторону, используя рекурсивные вычисления над s.

def half_element(x):
    k = ...
    s = (x + inv(x)) % Q
    for _ in range(k):
        v1 = m_sqrt(s**2 - 4)
        s = (s + v1 + m_sqrt(2*(s + v1) * (s+486662))) % Q
    return x

Данная формула работает напрямую с переменной S, избегая взаимодействия с x. Напрямую извлечь из неё x не получится, поскольку формула требует знания inv(x), но можно посчитать S для целевой точки и сравнивать с ней. Данная формула имеет 4 возможных решения в зависимости от знаков +-, и только одно из них есть на кривой, другие решения при повторном вызове не смогут извлечь корень и вернут ошибку. Полная формула изображена ниже

Основная проблема данной формулы, это наличие sqrt_mod, которые вычисляются также сложно, как и инверсии, что не позволяет превратить её из рекурсивной в закрытую форму. Сохранил её на случай новых открытий в модульной арифметике.

Примерно на этом месте я принял поражение и перестал исследовать curve25519.

Выводы

Я потратил на исследование многие месяцы, не все мои тесты попали в эту статью, но все они потерпели неудачу. Вероятно криптография эллиптических кривых будет безопасной до появления квантовых компьютеров, или возможно появится мощный AGI, который найдёт в ней изъян, который остальные не замечают.

Какой бы легкой на первый взгляд не казалась задача, если за ней стоят огромные деньги и тысячи ученых и математиков, то наверное лучше стоит обойти её стороной, и потратить время на что-то более полезное. Хотя как знать, может я остановился в шаге от мирового богатства.