做题随记

一道浮点位数很高的题目。。

题目

from decimal import *
getcontext().prec = 300 #给定精确度为300位

def pi():#生成三百位数的pi,小数点后299位
    lasts, t, s, n, na, d, da = 0, Decimal(3), 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n + na, na + 8
        d, da = d + da, da + 32
        t = (t * n) / d
        s += t
    return s

def sin(x):
    x = Decimal(x) % pi()
    p, factor = 0, x
    for n in range(10000):#sinx的泰勒展开
        p += factor
        factor *= - (x ** 2) / ((2 * n + 2) * (2 * n + 3))
    return p

flag = int.from_bytes(open('flag.txt', 'rb').read(), byteorder='big')
assert(flag < 2 ** 500)
print(sin(flag))
# 0.0244671923862258387777329311304061726117292338164472101939417405767980930400528701470775326212479053586551379126483335032568358369085951301715103743247208229702075426052832758137520953914256336461357368653604728080035204833015365117280775229337006339154544638317547072102876405128588259546493517773166
#实际上这个结果是小数点后301位。。

分析

pi()函数用一个不知道是什么的算法生成了300位（小数点后299位）的Π，sin(x)函数用泰勒展开到第9999轮来计算sinx的结果（反正对300位的结果来说绰绰有余了），然后给出了sin(flag)的结果，是小数点后301位的。之前强网杯的ezmath也是用高精度的浮点数，并且也用了e^x的泰勒展开。

尝试1

既然加密是用sinx的泰勒展开，最初想到用arcsin的展开式去求回去

def double_factorial(n):
    return __import__("functools").reduce(lambda x, y: x * y, range(n, 0, -2))

def arcsin(x):
    x = Decimal(x) % pi()
    p = x
    for i in range(1,3000):
        f=x
        p+=double_factorial(2*i-1)*(f**(2*i+1))/((2*i+1)*double_factorial(2*i))
    return p
#0.0244696342317182955380167878859997878648849860507225686198324843880327826384968171815366236530238829279631100698854243103679784022302933149263781311340548493687767140548481414565508236688055024007324091848224931595200235667864586900877342501435623376419548532684797634337279346936509642695145952284863

这个结果sin还原回去最后一位和结果不同，猜一位就能得到正确的余数

尝试2

牛顿法也能求根，吴师傅爆了3个小时才爆出来。。

from decimal import *
import multiprocessing

getcontext().prec = 300


def sin(x):
    p, factor = 0, x
    for n in range(10000):  # 感觉像泰勒
        p += factor
        factor *= - (x ** 2) / ((2 * n + 2) * (2 * n + 3))
    return p


l = Decimal(0)
r = Decimal(0.1)

des = Decimal(
    '0.0244671923862258387777329311304061726117292338164472101939417405767980930400528701470775326212479053586551379126'
    '483335032568358369085951301715103743247208229702075426052832758137520953914256336461357368653604728080035204833015'
    '365117280775229337006339154544638317547072102876405128588259546493517773166')
while True:
    mid2 = (l + r) / 2
    mid1 = (l + mid2) / 2
    mid3 = (mid2 + r) / 2
    p = multiprocessing.Pool(3)
    res = [p.apply_async(sin, args=(mid1,)), p.apply_async(sin, args=(mid2,)), p.apply_async(sin, args=(mid3,))]
    p.close()
    p.join()
    values = [i.get() for i in res]
    if des < values[0]:
        r = mid1
    elif values[0] < des < values[1]:
        l = mid1
        r = mid2
    elif values[1] < des < values[2]:
        l = mid2
        r = mid3
    elif values[2] < des:
        l = mid3
    else:
        break
    print('l=', end='')
    print(l)
    print('r=', end='')
    print(r)
    print('sin(x)=', end='')
    print(values[1])
print(l)
print(r)
print(mid1)
print(mid2)
print(mid3)
print(values[0])
print(values[1])
print(values[2])

#0.0244696342317182955380167878859997878648849860507225686198324843880327826384968171815366236530238829279631100698854243103679784022302933149263781311340548493687767140548481414565508236688055024007324091848224931595200235667864586900877342501435623376419548532684797634337279346936509642695145952284864

这个值用sin反展开回去和结果是一样的，nice，所以现在得到了flag和pi求余的结果了。

接下来是怎么求flag： \[ ∵n\%pi=x \\∴n=x+m*pi \] 因为n是整数，所以可以爆破，但是flag的范围太大了，感觉不太现实。

还有一个想法：利用n是整数这一点，x和pi都有小数部分，可以转换成同余方程： \[ 例如当精度是2的时候，3.14*m+0.02=n位整数\\ ∴14*m\%100=98 ...类推到300位也是如此 \]

但是最后算出来的m大小远超2**500，所以没能做出来。

对于这一题还有个疑问，题目给的结果小数点后有301位的，但是pi小数点后只有299位，这样用原始的pi是不可能出解的（试过将pi也增加到301位，但也没出）。