import math

def fx(x):
    if x == 0:
        x = 1e-10  # Avoid division by zero  #如果x能为0，注释掉这行##############
        pass
    return math.sin(x)/x   #把函数改成题干的形式###################

# 龙贝格方法 积分
def Romberg(a, b, err):
    table = [[],[],[],[]]
    t00 = (b-a)*(fx(a)+fx(b))/2
    table[0].append(t00)
    print(f"t0(0)={t00}")
    t01 = t10 = t11 = t20 = t21 = t30 = t31 = 0
    k = 1
    while True:
        tmp = 0
        for i in range(1, 2**(k-1)+1):
            tmp += fx(a + (b-a)*(2*i-1)/(2**k))
        t01 = t00/2+(b-a)*tmp/(2**k)
        print(f"t0({k})={t01}")
        table[0].append(t01)
        
        if k>1:
            t11 = (4*t01-t00)/3
            print(f"t1({k-1})={t11}")
            table[1].append(t11)
            if k>2:
                t21 = (16*t11-t10)/15
                print(f"t2({k-2})={t21}")
                table[2].append(t21)
                if k>3:
                    t31 = (64*t21-t20)/63
                    print(f"t3({k-3})={t31}")
                    table[3].append(t31)
                    if abs(t31-t30) < err:
                        break
                    t30 = t31
                else:
                    t30 = (64*t21-t20)/63
                    print(f"t3(0)={t30}")
                    table[3].append(t30)
                t20 = t21
            else:
                t20 = (16*t11-t10)/15
                print(f"t2(0)={t20}")
                table[2].append(t20)
            t10 = t11
        else:
            t10 = (4*t01-t00)/3
            print(f"t1(0)={t10}")
            table[1].append(t10)
        t00 = t01
        k += 1
    print("Romberg table:")
    for i in range(len(table)):
        print(f"t{i}: {table[i]}")
    return t31, k

# 把范围ab与精度换成题干的############
if __name__ == "__main__":
    a = 0
    b = 1
    err = 0.5e-6    #把精度要求改成题干的##################
    result, k = Romberg(a, b, err)
    print(f"Result: {result}, k={k}")