MODULAR ARITHMETIC

크립토핵 사이트에서 크립토 공부를 시작하며 문제를 풀기 전에 강의 course 부터 쭉 볼 생각이다.

가장 기초 강의인 MODULAR ARITHMETIC을 시작한다.

기초 수학인 것 같으니 간단한 설명과 해결 코드만 남기면서 지나가겠다.

 

첫번째는 Greatest Common Divisor로 그냥 최대공약수 GCD이다. 

GCD를 유클리드 알고리즘으로 구해보라는 문제이다.

def gcd(a,b):
    if a==0 or b==0:
        return a+b
    elif a>b:
        return gcd(a%b, b)
    else:
        return gcd(b%a, a)
    

a = 66528
b = 52920

print(gcd(a,b))

그냥 이렇게 짜면 된다.

아주 간단!

 

Extended Euclidean algorithm은  a*u+b*v = gcd(a, b) 를 만족하는 u와 v를 효과적으로 찾아주는 알고리즘이다.

유클리드 알고리즘으로 gcd를 찾아나가는 과정을 거꾸로 거슬러 올라가면 된다.

a=bq+r일 때 b=rq2+r2 이런식으로 찾는 과정에서 r=a-bq이고, r2=b-rq2=b-(a-bq)q2이므로 r2도 a*u+b*v 형태로 나타낼 수 있다.

따라서 r1, r2, ..., rn을 전부 a*u+b*v 형태로 표현 가능하다.

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0  # Base case: gcd(a, 0) = a, so we can return a, 1, 0
    
    # Recursive case: 
    gcd, x1, y1 = extended_gcd(b, a % b)
    # g=bx1+ry1=bx1+ay1-bqy1
    
    x = y1
    
    # y = x' - (a // b) * y'
    y = x1 - (a // b) * y1
    
    return gcd, x, y

# Example usage:
a = 26513
b = 32321

gcd, x, y = extended_gcd(a, b)

print(f"gcd({a}, {b}) = {gcd}")
print(f"Equation: {a}({x}) + {b}({y}) = {gcd}")

재귀함수를 이용해서 간단하게 코드를 작성할 수 있다.

 

그 다음은 그냥 mod에 대한 설명이다.

11 ≡ x mod 6
8146798528947 ≡ y mod 17

x, y 중에 더 작은걸 구하라는데 걍 %로 나머지 계산을 하면 된다.

 

그 다음은 페르마 소정리이다.

p가 소수일 때 a^(p-1)을 p로 나눈 나머지가 1이라는 기본적인 정리이다.

여기부터는 그냥 기본적인 정수론적 내용들에 대한 설명이 쭉 있다.

 

그 다음은 잉여역수에 대한 내용이다. 

g*d≡1(mod p) 를 만족하는 d를 찾아야 하는데 그냥 g*d+p*?=1이면 되니까 아까 본 유클리드 확장 알고리즘을 쓰면 된다.

 

그 다음은 이차잉여에 대한 내용이다.

a가 p에 대한 이차잉여이면, 제곱해서 모듈러로 a가 나오는 숫자가 존재한다는 뜻이다.

이와 관련해서 르장드르 심볼이라는 것이 있는데 

(a/p)=1이면 a는 p에 대한 이차잉여이고 -1이면 이차비잉여이다. 0이면 a가 p의 배수이다.

페르마 소정리에 의해 a^(p-1)이 1이기 때문에 a^((p-1)/2)이 르장드르 심볼 (a/p)와 p에 대해 합동이다.

즉, a^((p-1)/2) 이 1이면, a=x^2인 x값이 존재하고 -1이면 존재하지 않는다.

p = 101524035174539890485408575671085261788758965189060164484385690801466167356667036677932998889725476582421738788500738738503134356158197247473850273565349249573867251280253564698939768700489401960767007716413932851838937641880157263936985954881657889497583485535527613578457628399173971810541670838543309159139

ints = [25081841204695904475894082974192007718642931811040324543182130088804239047149283334700530600468528298920930150221871666297194395061462592781551275161695411167049544771049769000895119729307495913024360169904315078028798025169985966732789207320203861858234048872508633514498384390497048416012928086480326832803, 45471765180330439060504647480621449634904192839383897212809808339619841633826534856109999027962620381874878086991125854247108359699799913776917227058286090426484548349388138935504299609200377899052716663351188664096302672712078508601311725863678223874157861163196340391008634419348573975841578359355931590555, 17364140182001694956465593533200623738590196990236340894554145562517924989208719245429557645254953527658049246737589538280332010533027062477684237933221198639948938784244510469138826808187365678322547992099715229218615475923754896960363138890331502811292427146595752813297603265829581292183917027983351121325, 14388109104985808487337749876058284426747816961971581447380608277949200244660381570568531129775053684256071819837294436069133592772543582735985855506250660938574234958754211349215293281645205354069970790155237033436065434572020652955666855773232074749487007626050323967496732359278657193580493324467258802863, 4379499308310772821004090447650785095356643590411706358119239166662089428685562719233435615196994728767593223519226235062647670077854687031681041462632566890129595506430188602238753450337691441293042716909901692570971955078924699306873191983953501093343423248482960643055943413031768521782634679536276233318, 85256449776780591202928235662805033201684571648990042997557084658000067050672130152734911919581661523957075992761662315262685030115255938352540032297113615687815976039390537716707854569980516690246592112936796917504034711418465442893323439490171095447109457355598873230115172636184525449905022174536414781771, 50576597458517451578431293746926099486388286246142012476814190030935689430726042810458344828563913001012415702876199708216875020997112089693759638454900092580746638631062117961876611545851157613835724635005253792316142379239047654392970415343694657580353333217547079551304961116837545648785312490665576832987, 96868738830341112368094632337476840272563704408573054404213766500407517251810212494515862176356916912627172280446141202661640191237336568731069327906100896178776245311689857997012187599140875912026589672629935267844696976980890380730867520071059572350667913710344648377601017758188404474812654737363275994871, 4881261656846638800623549662943393234361061827128610120046315649707078244180313661063004390750821317096754282796876479695558644108492317407662131441224257537276274962372021273583478509416358764706098471849536036184924640593888902859441388472856822541452041181244337124767666161645827145408781917658423571721, 18237936726367556664171427575475596460727369368246286138804284742124256700367133250078608537129877968287885457417957868580553371999414227484737603688992620953200143688061024092623556471053006464123205133894607923801371986027458274343737860395496260538663183193877539815179246700525865152165600985105257601565]

for i in ints :
	if (pow(i, (p-1)//2, p) == 1) :
		print(i, 'is qr')
		print(pow(i, (p+1)//4, p))

문제를 해결하기 위해서는 이렇게 코드를 짜면 된다.

그냥 (p-1)/2만큼 pow 함수로 거듭제곱을 해보면 된다.

그리고 만약 이차잉여이면, 제곱근 값을 출력해주면 된다.

(p+1)/4번 거듭제곱한 값이 제곱근인 이유는 a^(p+1)이 a제곱이니까 루트a는 a^(p+1)/4이다. 

문제에서 p가 4k+3꼴이라고 했으니까 가능하다.

그런데 p가 너무 큰데 실행속도가 이상하게 빠르다.

1억번 연산하는데 1초가 걸리는데 pow로 거듭제곱을 일일이 하면 p-1번 해야 하는데 시간이 너무 많이 걸려야 정상인데 이상해서 찾아보았다.

pow 함수에서 3번째 인자까지 사용해 모듈러 지수화를 하면 이진 분할 기법이라는 알고리즘으로 빠르게 계산이 가능하다.

지수를 이진수로 변환하고 각 이진수 자릿수에 대해 제곱을 반복하며 이진 표현에서 1이 있는 위치만 계산을 하면 

우리가 지수로 입력한 p-1의 log값만큼만 계산을 반복할 수 있다.

 

그 다음은 Tonelli-Shanks 알고리즘이다. 이 알고리즘은 제곱해서 모듈러p에 대해  n이 되는 값을 시간복잡도 log의 제곱만에 찾아낼 수 있는 알고리즘이다.

방법과 원리에 대한 설명은 수식을 다 적기 너무 복잡해서

https://sean9892.tistory.com/27

Introduction and Proof of Tonelli-Shanks Algorithm

이 링크를 참조하면 될 것 같다.

처음에는 너무 복잡해서 이해하기 힘들긴 한데 한 4,5번 읽어보면 이해 된다.

def tonelli_shanks(n,p):
    def isQR(x,p):
        return pow(x,(p-1)//2,p)==1
    
    if not isQR(n,p):
        return -1

    Q,S=p-1,0
    while Q%2==0:
        S+=1
        Q//=2

    z=None
    for x in range(2,p):
        if not isQR(x,p):
            z=x
            break

    M,c,t,R=S,pow(z,Q,p),pow(n,Q,p),pow(n,(Q+1)//2,p)

    while True:
        if t==0:
            return 0
        elif t==1:
            return R
        k=t*t%p
        ii=None
        for i in range(1,M):
            if k==1:
                ii=i
                break
            k*=k
            k%=p
        b=pow(c,2**(M-i-1),p)%p
        M=ii%p
        c=b*b%p
        t=t*c%p
        R=R*b%p

a = 8479994658316772151941616510097127087554541274812435112009425778595495359700244470400642403747058566807127814165396640215844192327900454116257979487432016769329970767046735091249898678088061634796559556704959846424131820416048436501387617211770124292793308079214153179977624440438616958575058361193975686620046439877308339989295604537867493683872778843921771307305602776398786978353866231661453376056771972069776398999013769588936194859344941268223184197231368887060609212875507518936172060702209557124430477137421847130682601666968691651447236917018634902407704797328509461854842432015009878011354022108661461024768
p = 30531851861994333252675935111487950694414332763909083514133769861350960895076504687261369815735742549428789138300843082086550059082835141454526618160634109969195486322015775943030060449557090064811940139431735209185996454739163555910726493597222646855506445602953689527405362207926990442391705014604777038685880527537489845359101552442292804398472642356609304810680731556542002301547846635101455995732584071355903010856718680732337369128498655255277003643669031694516851390505923416710601212618443109844041514942401969629158975457079026906304328749039997262960301209158175920051890620947063936347307238412281568760161

print(tonelli_shanks(a,p))

문제 해결 코드는 이렇게 작성하면 된다.

 

그 다음은 중국인의 나머지 정리 CRT에 대한 내용이다.

m1~mn은 서로소인 정수이고 이 n개의 합동식을 만족하는 유일한 정수해가 존재한다는 것이 CRT의 내용이다.

그래서 그 유일한 해를 찾는 방법은 다음과 같다.

저 x 시그마값을 m1~mn으로 머릿속으로 나눠보면 a1~an이 나온다는 것을 알아차릴 수 있다.

def extended_gcd(a, b):
    """ 확장 유클리드 알고리즘을 사용하여 gcd(a, b)와 a의 b에 대한 역원을 찾는다. """
    if b == 0:
        return a, 1, 0
    gcd, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return gcd, x, y

def mod_inverse(a, m):
    """ 모듈로 m에 대한 a의 역원을 계산한다. """
    gcd, x, _ = extended_gcd(a, m)
    if gcd == 1:
        return x % m
    else:
        return None  # 역원이 존재하지 않음

def chinese_remainder_theorem(a, m):
    """ 중국인의 나머지 정리를 사용하여 연립 합동식의 해를 찾는다. """
    n = len(a)
    M = 1
    x = 0
    
    # 전체 모듈로 M을 계산
    for mi in m:
        M *= mi
    
    # 각 합동식에 대해 계산
    for i in range(n):
        Mi = M // m[i]
        yi = mod_inverse(Mi, m[i])
        if yi is None:
            raise ValueError("역원이 존재하지 않습니다. m_i가 서로소인지 확인하세요.")
        x += a[i] * Mi * yi
    
    # 최종 해를 모듈로 M으로 나타냄
    x %= M
    return x

# 입력: 각 합동식의 a값과 m값
a = [2, 3, 1]  # x ≡ 2 (mod 3), x ≡ 3 (mod 4), x ≡ 1 (mod 5)
m = [3, 4, 5]  # 모듈로 값

# CRT를 사용하여 해를 찾기
result = chinese_remainder_theorem(a, m)
print("합동식을 만족하는 x의 값은:", result)

CRT로 n개의 합동식을 만족하는 해를 찾는 코드이다.

이 코드에 의하면, 해를 찾는데 걸리는 시간은 nlogm이다. n은 합동식의 개수이고 m은 m1~mn중의 최대값이다.

유클리드 확장 알고리즘으로 역원을 찾는 과정이 logm만큼 걸리기 때문이다.

그런데 CRT로 해를 찾으면 m값이 커지거나 n이 많아지면 시간이 오래 걸리게 되는데, 찾다보니까 더 효율적인 알고리즘이 있다고 해서 소개하겠다.

 

Garner's Algorithm이라는 것인데 효율적으로 큰 수의 모듈로 연산을 처리할 수 있는 방법이다.

 

이 알고리즘을 구현한 파이썬 코드이다.

def extended_gcd(a, b):
    """확장 유클리드 알고리즘을 사용하여 gcd(a, b)와 a의 b에 대한 역원을 찾습니다."""
    if b == 0:
        return a, 1, 0
    gcd, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return gcd, x, y

def mod_inverse(a, m):
    """모듈로 m에 대한 a의 역원을 계산합니다."""
    gcd, x, _ = extended_gcd(a, m)
    if gcd != 1:
        return None  # 모듈로 역원이 존재하지 않음
    else:
        return x % m

def garner_algorithm(a, m):
    """Garner의 알고리즘을 사용하여 연립 합동식의 해를 찾습니다."""
    n = len(a)
    u = [0] * n
    u[0] = a[0]
    M = [1] * n  # M_i 누적 곱을 저장
    x = u[0]
    
    for i in range(1, n):
        M[i] = M[i-1] * m[i-1]
        inv = mod_inverse(M[i], m[i])
        
        if inv is None:
            raise ValueError("모듈로 역원이 존재하지 않습니다.")
        
        u[i] = ((a[i] - x) * inv) % m[i]
        x += u[i] * M[i]

    return x

# 예제
a = [2, 3, 5]
m = [5, 11, 17]
print("result x:", garner_algorithm(a, m))

 

그 다음은 encrypt 코드와 결과물을 보고 원본 문자열 flag를 찾으라는 문제이다.

from random import randint

a = 288260533169915
p = 1007621497415251

FLAG = b'crypto{????????????????????}'


def encrypt_flag(flag):
    ciphertext = []
    plaintext = ''.join([bin(i)[2:].zfill(8) for i in flag])
    print(plaintext)
    for b in plaintext:
        e = randint(1, p)
        n = pow(a, e, p)
        if b == '1':
            ciphertext.append(n)
        else:
            n = -n % p
            ciphertext.append(n)
    return ciphertext


print(encrypt_flag(FLAG))

encrypt 코드이다.

plaintext는 flag의 문자를 하나씩 i에 담고, bin(i)[2:]로 이진수로 만들고 맨앞에 0b부분을 제외한다.

zfill로 문자열의 길이를 8로 맞춰주고 이 과정을 전부 다 한 다음에 join으로 하나의 문자열로 합친다.

그리고 b가 1이냐 0이냐에 따라 다른 방식으로 계산한 n값을 ciphertext에 추가해서 출력한다.

원본 문자열 flag를 얻어내려면 plaintext를 알아내서 8글자씩 자른 다음에 문자열로 복원하면 된다.

plaintext를 알아내려면 출력값을 보고 b값들을 한글자씩 알아내면 되는데 b가 1인지 0인지 판단만 하면 된다.

a = 288260533169915
p = 1007621497415251

print(pow(a,(p-1)//2,p))

이렇게 실행해보면 1이 나온다. 

즉, a는 p에 대한 이차잉여이기 때문에 a^e도 이차잉여이다.

그리고 -a^e는 이차비잉여이기 때문에 ciphertext에 있는 값이 이차잉여이면 b=1이고 이차비잉여이면 b=0임을 알 수 있다.

from Crypto.Util.number import *

a = 288260533169915
p = 1007621497415251

print(pow(a,(p-1)//2,p))

output = [67594220461269, 501237540280788, 718316769824518, 296304224247167, 48290626940198, 30829701196032, 521453693392074, 840985324383794, 770420008897119, 745131486581197, 729163531979577, 334563813238599, 289746215495432, 538664937794468, 894085795317163, 983410189487558, 863330928724430, 996272871140947, 352175210511707, 306237700811584, 631393408838583, 589243747914057, 538776819034934, 365364592128161, 454970171810424, 986711310037393, 657756453404881, 388329936724352, 90991447679370, 714742162831112, 62293519842555, 653941126489711, 448552658212336, 970169071154259, 339472870407614, 406225588145372, 205721593331090, 926225022409823, 904451547059845, 789074084078342, 886420071481685, 796827329208633, 433047156347276, 21271315846750, 719248860593631, 534059295222748, 879864647580512, 918055794962142, 635545050939893, 319549343320339, 93008646178282, 926080110625306, 385476640825005, 483740420173050, 866208659796189, 883359067574584, 913405110264883, 898864873510337, 208598541987988, 23412800024088, 911541450703474, 57446699305445, 513296484586451, 180356843554043, 756391301483653, 823695939808936, 452898981558365, 383286682802447, 381394258915860, 385482809649632, 357950424436020, 212891024562585, 906036654538589, 706766032862393, 500658491083279, 134746243085697, 240386541491998, 850341345692155, 826490944132718, 329513332018620, 41046816597282, 396581286424992, 488863267297267, 92023040998362, 529684488438507, 925328511390026, 524897846090435, 413156582909097, 840524616502482, 325719016994120, 402494835113608, 145033960690364, 43932113323388, 683561775499473, 434510534220939, 92584300328516, 763767269974656, 289837041593468, 11468527450938, 628247946152943, 8844724571683, 813851806959975, 72001988637120, 875394575395153, 70667866716476, 75304931994100, 226809172374264, 767059176444181, 45462007920789, 472607315695803, 325973946551448, 64200767729194, 534886246409921, 950408390792175, 492288777130394, 226746605380806, 944479111810431, 776057001143579, 658971626589122, 231918349590349, 699710172246548, 122457405264610, 643115611310737, 999072890586878, 203230862786955, 348112034218733, 240143417330886, 927148962961842, 661569511006072, 190334725550806, 763365444730995, 516228913786395, 846501182194443, 741210200995504, 511935604454925, 687689993302203, 631038090127480, 961606522916414, 138550017953034, 932105540686829, 215285284639233, 772628158955819, 496858298527292, 730971468815108, 896733219370353, 967083685727881, 607660822695530, 650953466617730, 133773994258132, 623283311953090, 436380836970128, 237114930094468, 115451711811481, 674593269112948, 140400921371770, 659335660634071, 536749311958781, 854645598266824, 303305169095255, 91430489108219, 573739385205188, 400604977158702, 728593782212529, 807432219147040, 893541884126828, 183964371201281, 422680633277230, 218817645778789, 313025293025224, 657253930848472, 747562211812373, 83456701182914, 470417289614736, 641146659305859, 468130225316006, 46960547227850, 875638267674897, 662661765336441, 186533085001285, 743250648436106, 451414956181714, 527954145201673, 922589993405001, 242119479617901, 865476357142231, 988987578447349, 430198555146088, 477890180119931, 844464003254807, 503374203275928, 775374254241792, 346653210679737, 789242808338116, 48503976498612, 604300186163323, 475930096252359, 860836853339514, 994513691290102, 591343659366796, 944852018048514, 82396968629164, 152776642436549, 916070996204621, 305574094667054, 981194179562189, 126174175810273, 55636640522694, 44670495393401, 74724541586529, 988608465654705, 870533906709633, 374564052429787, 486493568142979, 469485372072295, 221153171135022, 289713227465073, 952450431038075, 107298466441025, 938262809228861, 253919870663003, 835790485199226, 655456538877798, 595464842927075, 191621819564547]

plaintext = ""

for i in output:
    if pow(i, (p-1)//2, p) == 1:
        plaintext += '1'
    else:
        plaintext += '0'

plaintext = int(plaintext,2)
print(long_to_bytes(plaintext))

이렇게 flag를 찾을 수 있다.

이차잉여인지 판별해서 plaintext를 복구하고 이를 long_to_bytes 함수를 이용해 원본 문자열 flag로 바꿀 수 있다.

long_to_bytes는 알아서 바이트(8비트)로 나눈뒤에 아스키 문자로 출력해준다.

 

그 다음은 modular binomials로 두 합동식에서 N=p*q일 때 p, q를 찾아야 한다.

여기서 N, c1, c2, e1, e2가 주어져 있다.

이렇게 해주면 결국 q의 배수를 뽑아낼 수 있으니까 이 값과 N의 최대공약수를 구하면 q가 나온다.

마찬가지로 p도 구할 수 있다.

N = 14905562257842714057932724129575002825405393502650869767115942606408600343380327866258982402447992564988466588305174271674657844352454543958847568190372446723549627752274442789184236490768272313187410077124234699854724907039770193680822495470532218905083459730998003622926152590597710213127952141056029516116785229504645179830037937222022291571738973603920664929150436463632305664687903244972880062028301085749434688159905768052041207513149370212313943117665914802379158613359049957688563885391972151218676545972118494969247440489763431359679770422939441710783575668679693678435669541781490217731619224470152467768073
e1 = 12886657667389660800780796462970504910193928992888518978200029826975978624718627799215564700096007849924866627154987365059524315097631111242449314835868137
e2 = 12110586673991788415780355139635579057920926864887110308343229256046868242179445444897790171351302575188607117081580121488253540215781625598048021161675697
c1 = 14010729418703228234352465883041270611113735889838753433295478495763409056136734155612156934673988344882629541204985909650433819205298939877837314145082403528055884752079219150739849992921393509593620449489882380176216648401057401569934043087087362272303101549800941212057354903559653373299153430753882035233354304783275982332995766778499425529570008008029401325668301144188970480975565215953953985078281395545902102245755862663621187438677596628109967066418993851632543137353041712721919291521767262678140115188735994447949166616101182806820741928292882642234238450207472914232596747755261325098225968268926580993051
c2 = 14386997138637978860748278986945098648507142864584111124202580365103793165811666987664851210230009375267398957979494066880296418013345006977654742303441030008490816239306394492168516278328851513359596253775965916326353050138738183351643338294802012193721879700283088378587949921991198231956871429805847767716137817313612304833733918657887480468724409753522369325138502059408241232155633806496752350562284794715321835226991147547651155287812485862794935695241612676255374480132722940682140395725089329445356434489384831036205387293760789976615210310436732813848937666608611803196199865435145094486231635966885932646519

def extended_gcd(a, b):
    if b == 0:
        return a, 1, 0
    gcd, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return gcd, x, y

def mod_inverse(a, m):
    gcd, x, _ = extended_gcd(a, m)
    if gcd != 1:
        return None
    else:
        return x % m

a1 = 2
b1 = 3
a2 = 5
b2 = 7

forq = pow(a2, -e1*e2, N) * pow(c2, e1, N) - pow(a1, -e1*e2, N) * pow(c1, e2, N)
q,_,_ = extended_gcd(forq, N)
forp = pow(b2, -e1*e2, N) * pow(c2, e1, N) - pow(b1, -e1*e2, N) * pow(c1, e2, N)
p,_,_ = extended_gcd(forp, N)
print(p, q)

이렇게 p, q를 구할 수 있다.