복잡한 다항식 수식에서 계산기 유효 자릿수에 따른 approx() 오차

TI-nspire 기종에 대한 solve 질문글에 대해 답변을 하던 중 이상한 점을 발견하였습니다. 

https://allcalc.org/55823

평소에는 방정식의 해를 numeric 한 방식으로 solve 를 이용해 찾는 것보다, 

(소숫점을 없애서) exact 방식으로 참 값을 먼저 구하고 → 그 값에 대한 근사값 approx(참 값) 을 구하는 것이 더 정확했습니다. 

그런데 이번에는 반대로 오차가 커지는 겁니다. 

이 궁금증을 해결하기 위해 직접 분석을 진행해 보았습니다.

1. 분석 대상 수식

분석에 사용된 수식은 다음과 같습니다. cos, sin의 단위는 모두 degree입니다.

$$ \dfrac{120 \left( 150000 \left( \cos\left(\dfrac{7001}{5000}\right) \sqrt{2} - 2 \sin\left(\dfrac{217999}{5000}\right) \right) \cos\left(\dfrac{217999}{5000}\right) - 136342 \sin\left(\dfrac{217999}{5000}\right) \cos\left(\dfrac{7001}{5000}\right) - \left( 150000 \sin\left(\dfrac{7001}{5000}\right) \sin\left(\dfrac{217999}{5000}\right) - 68171 \right) \sqrt{2} \right)}{\sin\left(\dfrac{7001}{5000}\right) \left( 150000 \cos\left(\dfrac{7001}{5000}\right) \cos\left(\dfrac{217999}{5000}\right) - 150000 \sin\left(\dfrac{7001}{5000}\right) \sin\left(\dfrac{217999}{5000}\right) + 68171 \right)} $$

((120*(150000*(cos(((7001)/(5000)))*√(2)-2*sin(((217999)/(5000))))*cos(((217999)/(5000)))-136342*sin(((217999)/(5000)))*cos(((7001)/(5000)))-(150000*sin(((7001)/(5000)))*sin(((217999)/(5000)))-68171)*√(2)))/(sin(((7001)/(5000)))*(150000*cos(((7001)/(5000)))*cos(((217999)/(5000)))-150000*sin(((7001)/(5000)))*sin(((217999)/(5000)))+68171)))

2. 가장 정확한 기준값은? (고정밀도 계산)

오차를 측정하려면 가장 정확한 '참값'이 필요합니다. 일반적인 PC 계산 환경(64비트 float)의 한계를 넘어서기 위해, 파이썬의 mpmath 라이브러리를 사용하여 100자리의 정밀도로 기준값을 계산했습니다.

기준값 (100자리 정밀도):
73.04950705847862934420128091048894148771096960598761210206551516481655390211661403804814403351531886

3. 최종 정밀도 분석표

위 100자리 기준값을 바탕으로, 각기 다른 유효자릿수를 가진 가상의 계산기를 시뮬레이션하여 실제 오차를 측정한 결과입니다.

4. 왜 이런 오차가 발생할까?

컴퓨터의 숫자 저장 방식: 2진법의 한계

가장 근본적인 원인은 컴퓨터가 숫자를 2진법으로 저장하는 데 있습니다. 우리가 사용하는 10진수 소수 중 상당수는 2진수로 변환하면 무한소수가 되어, 정해진 비트(bit) 안에 완벽하게 담지 못하고 근사치로 저장됩니다. 이 작은 근사 오차가 계산 과정에서 계속 누적되어 최종 결과에 영향을 미칩니다.

누적 오차의 예시

(1 / 3) * 3 을 유효자릿수 4자리 계산기로 계산하는 상황을 가정해 봅시다.

1. 1 / 3 계산: 결과는 0.333333... 이지만, 4자리만 저장할 수 있으므로 0.3333으로 반올림됩니다. (첫 오차 발생)
2. 0.3333 * 3 계산: 결과는 0.9999가 됩니다. 참값인 1.0과 미세한 차이가 생깁니다.

복잡한 수식은 이런 과정이 수십, 수백 번 반복되는 것과 같으므로 작은 오차들이 모여 눈에 띄는 차이를 만들게 됩니다.

문제가 된 처음의 수식 역시 계산기 입장에서

   * 나눗셈 (Division): 12회
       * 분자와 분모의 각 sin, cos 함수 안에 있는 분수 계산 11회
       * 마지막에 분자를 분모로 나누는 계산 1회
   * 사인 (sin): 7회 
   * 코사인 (cos): 5회
   * 루트 (√): 2회 
   * 곱셈 (Multiplication): 15회
   * 뺄셈 (Subtraction): 5회
   * 덧셈 (Addition): 1회

를 모두 계산하는 과정에서 오차가 매번 발생하며 때로는 증폭되기도 하고, 또 누적되기도 하여 

최종적으로 solve의 정상 범주를 벗어난 x값이 구해진 것이라고 분석할 수 있겠습니다. 

부록: 분석에 사용된 전체 Python 코드

이 분석을 직접 재현해보고 싶으신 분들을 위해, 최종 분석에 사용된 전체 코드를 공유합니다. (mpmath 라이브러리 설치가 필요합니다: pip install mpmath)

import math
import sys
from mpmath import mp

# --- Simulation Function ---
def round_to_significant_digits(n, p):
    if n == 0 or not math.isfinite(n):
        return n
    try:
        factor = 10 ** (p - math.ceil(math.log10(abs(n))))
    except ValueError:
        return 0.0
    return round(n * factor) / factor

def calculate_simulation(precision=None):
    def apply_prec(value):
        if precision is not None:
            return round_to_significant_digits(value, precision)
        return value

    c_120, c_150000, c_7001, c_5000, c_2, c_217999, c_136342, c_68171 = \
        120.0, 150000.0, 7001.0, 5000.0, 2.0, 217999.0, 136342.0, 68171.0

    angle1_deg = apply_prec(c_7001 / c_5000)
    angle2_deg = apply_prec(c_217999 / c_5000)
    angle1_rad = apply_prec(math.radians(angle1_deg))
    angle2_rad = apply_prec(math.radians(angle2_deg))

    cos_a1, sin_a1, cos_a2, sin_a2, sqrt_2 = \
        apply_prec(math.cos(angle1_rad)), apply_prec(math.sin(angle1_rad)), \
        apply_prec(math.cos(angle2_rad)), apply_prec(math.sin(angle2_rad)), \
        apply_prec(math.sqrt(c_2))

    termA_inner_part1 = apply_prec(cos_a1 * sqrt_2)
    termA_inner_part2 = apply_prec(c_2 * sin_a2)
    termA_inner_sub = apply_prec(termA_inner_part1 - termA_inner_part2)
    termA = apply_prec(c_150000 * termA_inner_sub)
    mul_term1 = apply_prec(termA * cos_a2)
    termB = apply_prec(apply_prec(c_136342 * sin_a2) * cos_a1)
    termC_inner_sub = apply_prec(apply_prec(apply_prec(c_150000 * sin_a1) * sin_a2) - c_68171)
    termC = apply_prec(termC_inner_sub * sqrt_2)
    numerator = apply_prec(c_120 * apply_prec(apply_prec(mul_term1 - termB) - termC))

    den_term1 = sin_a1
    den_term2 = apply_prec(apply_prec(c_150000 * cos_a1) * cos_a2)
    den_term3 = apply_prec(apply_prec(c_150000 * sin_a1) * sin_a2)
    inner_den_final = apply_prec(apply_prec(den_term2 - den_term3) + c_68171)
    denominator = apply_prec(den_term1 * inner_den_final)

    return apply_prec(numerator / denominator) if denominator != 0 else float('inf')

# --- High-Precision Ground Truth Function ---
def get_high_precision_ground_truth():
    mp.dps = 100
    c_120, c_150000, c_7001, c_5000, c_2, c_217999, c_136342, c_68171 = \
        mp.mpf(120), mp.mpf(150000), mp.mpf(7001), mp.mpf(5000), mp.mpf(2), \
        mp.mpf(217999), mp.mpf(136342), mp.mpf(68171)

    angle1_rad = mp.radians(c_7001 / c_5000)
    angle2_rad = mp.radians(c_217999 / c_5000)

    cos_a1, sin_a1, cos_a2, sin_a2, sqrt_2 = \
        mp.cos(angle1_rad), mp.sin(angle1_rad), mp.cos(angle2_rad), mp.sin(angle2_rad), mp.sqrt(c_2)

    termA = c_150000 * (cos_a1 * sqrt_2 - c_2 * sin_a2)
    mul_term1 = termA * cos_a2
    termB = c_136342 * sin_a2 * cos_a1
    termC = (c_150000 * sin_a1 * sin_a2 - c_68171) * sqrt_2
    numerator = c_120 * (mul_term1 - termB - termC)

    inner_den_final = (c_150000 * cos_a1 * cos_a2) - (c_150000 * sin_a1 * sin_a2) + c_68171
    denominator = sin_a1 * inner_den_final

    return numerator / denominator

# --- Main Execution ---
if __name__ == "__main__":
    ground_truth = get_high_precision_ground_truth()

    print("### 최종 정밀도 분석표 (100자리 기준값 사용) ###")
    print("-" * 90)
    print(f"기준값 (mpmath 100자리): {ground_truth}")
    print("-" * 90)
    print(f"{ '유효자릿수':<12} | { '시뮬레이션 계산 결과':<45} | { '실제 오차 (절대값)':<30}")
    print("-" * 90)

    for p in range(6, 18):
        sim_result = calculate_simulation(precision=p)
        error = abs(mp.mpf(sim_result) - ground_truth)
        print(f"{p:<12} | {sim_result:<45} | {error}")

    print("-" * 90)

    standard_float_result = calculate_simulation()
    standard_float_error = abs(mp.mpf(standard_float_result) - ground_truth)
    print(f"참고: 표준 64비트 float 계산")
    print(f" - 계산 결과: {standard_float_result}")
    print(f" - 실제 오차: {standard_float_error}")
    print("-" * 90)