فهرست مقاله
برای پذیرفتن یک تابع عددی، خروجیِ بدون خطا و چند رقم آشنا کافی نیست. معادلهٔ درجه‌دو مثال خوبی است: فرمولی که روی کاغذ درست است، می‌تواند ریشهٔ کوچک را در محاسبهٔ ممیز شناور از دست بدهد. در این درس، قرارداد محدود و قابل‌آزمونی برای دو ریشهٔ حقیقی می‌سازیم؛ جواب دستی، مرجع مستقل و خطای نسبی را کنار هم می‌بینیم. هیچ مدار یا مدل واقعی اجرا نشده و نتیجه، تأیید یک سامانهٔ مهندسی نیست.
۱. قرارداد مسئله پیش از انتخاب تلورانس
معادلهٔ ما a*u²+b*u+c=۰ است. u و ضرایب در این تمرین بی‌بعدند؛ مثلاً اگر ولتاژ را بر ولتاژ مرجع تقسیم کرده باشیم، u عدد نرمال‌شده خواهد بود. کد واحد فیزیکی را حدس نمی‌زند و این نرمال‌سازی باید بیرون تابع درست تعریف شود. برای نمونهٔ دستی، ضرایب یک، منفی پنج و شش دو ریشهٔ دو و سه دارند؛ این جواب ساده، ترتیب و علامت پایه را کنترل می‌کند.
ورودی فقط int یا float دقیق است و ابتدا به float تبدیل می‌شود. مقدار غیرصفر هر ضریب باید از نظر قدرمطلق بین ده به توان منفی صد و ده به توان صد باشد. صفر برای b و c مجاز است، ولی a صفر نیست. این حدود، قرارداد آموزشی برای پرهیز از دامنه‌های افراطیِ واسطه‌هاست؛ نه دامنهٔ تمام چندجمله‌ای‌ها و نه تضمین دقت برای هر مسئلهٔ ممکن.
دو ریشه به ترتیب عددی برگردانده می‌شوند، نه ترتیب بزرگیِ قدرمطلق. ریشهٔ تکراری دو بار در خروجی می‌آید. معادلهٔ خطی و ریشه‌های مختلط پشتیبانی نمی‌شوند و خطای صریح دارند. برای حفظ دقت همهٔ ارقام یک int بسیار بزرگ یا معنای دقیقِ رشتهٔ اعشاری نیز ادعایی نداریم؛ مسئلهٔ محاسبه، ضرایبی است که واقعاً پس از تبدیل در float نمایش داده شده‌اند.
۲. چرا ریشهٔ کوچک در تفریق گم می‌شود؟
در فرمول مستقیم، صورت یکی از ریشه‌ها می‌تواند تفاضل دو عدد تقریباً برابر باشد. مثال تألیفیِ ما a و c برابر یک و b برابر ده به توان شانزده است. یکی از ریشه‌ها نزدیک منفی ده به توان شانزده و دیگری نزدیک منفی ده به توان منفی شانزده است. در اجرای مستقیمِ ثبت‌شده، صورت ریشهٔ کوچک صفر شده و خروجیِ آن صفر است؛ اندازهٔ ظاهراً بزرگِ ریشهٔ دیگر این خطا را نشان نمی‌دهد.
بحث حذف عددی Goldberg توضیح می‌دهد که بازنویسی بعضی عبارت‌ها می‌تواند از این افت جلوگیری کند؛ محاسبهٔ ممیز نیز مسئلهٔ جداگانه‌ای دارد. برای شکل معادلهٔ ما، q را از منفی یک‌دومِ مجموع b و رادیکالِ ممیز با علامتِ b می‌سازیم. سپس یک ریشه q/a و دیگری c/q است؛ برای ریشهٔ کوچک دوباره همان تفریق حساس را انجام نمی‌دهیم.
یادداشت اولیهٔ Kahan دربارهٔ درجه‌دو انتخاب علامت و دشواری ریشه‌های نزدیک را بررسی می‌کند. شکل نوشتن معادله در آن یادداشت A*z²−۲B*z+C است و با نام‌گذاری ما فرق دارد. کد این درس پیاده‌سازی تألیفی و محدودِ خودش است؛ برنامهٔ MATLAB آن منبع را منتقل یا اجرا نکرده‌ایم و ادعای دقت تمام بیت‌های آن را به این مثال نسبت نمی‌دهیم.
۳. کد کامل و دو مرز دقت مستقل
علامت ممیز D=b²−۴ac را با Fraction ساخته‌شده از همان floatهای ورودی تعیین می‌کنیم؛ بنابراین در این مرحله، حاصل‌ضرب و تفاضل کسری دقیق‌اند. D مثبت سپس یک بار برای sqrt به float تبدیل می‌شود. این انتخاب ترکیبی، هزینهٔ حساب گویا دارد و آزمون سرعت نیست. محاسبهٔ q و تقسیم‌های ریشه همچنان ممیز شناورند؛ دقیق بودن علامت D به معنی دقیق بودنِ همهٔ مراحل نیست.
"""Original bounded real quadratic-root lesson; exact represented-float D sign."""
from decimal import Decimal, localcontext
from fractions import Fraction
import itertools
import json
import math
import platform
def _coefficient(value):
if type(value) not in (int, float):
raise TypeError('coefficient must be exact int or float')
try:
value = float(value)
except OverflowError as error:
raise ValueError('coefficient not representable') from error
if not math.isfinite(value) or (value != 0 and not 1e-100 <= abs(value) <= 1e100):
raise ValueError('nonzero coefficient magnitude must be within authored [1e-100,1e100]')
return value
def real_quadratic_roots(a, b, c):
"""ax²+bx+c=0, dimensionless coefficients/unknown, real roots only."""
a, b, c = map(_coefficient, (a, b, c))
if a == 0:
raise ValueError('a=0 is not a quadratic in this API')
af, bf, cf = map(Fraction.from_float, (a, b, c))
exact_d = bf * bf - 4 * af * cf
if exact_d < 0:
raise ValueError('no real roots for the represented coefficients')
if exact_d == 0:
root = -b / (2 * a)
return root, root
discriminant = float(exact_d)
if not math.isfinite(discriminant) or discriminant <= 0:
raise ValueError('positive discriminant outside representable range')
q = -0.5 * (b + math.copysign(math.sqrt(discriminant), b))
if q == 0 or not math.isfinite(q):
raise ValueError('q outside representable nonzero range')
roots = q / a, c / q
if not all(math.isfinite(x) for x in roots):
raise ValueError('root not finite')
return tuple(sorted(roots))
def _decimal_reference(a, b, c):
"""Independent high-precision direct formula for the actual float inputs."""
with localcontext() as context:
context.prec = 120
da, db, dc = (Decimal.from_float(float(v)) for v in (a, b, c))
discriminant = db * db - 4 * da * dc
if discriminant < 0:
raise ValueError('negative reference discriminant')
root_d = discriminant.sqrt()
return tuple(sorted(((-db - root_d) / (2 * da), (-db + root_d) / (2 * da))))
def _naive(a, b, c):
discriminant = b * b - 4 * a * c
root_d = math.sqrt(discriminant)
return tuple(sorted(((-b - root_d) / (2 * a), (-b + root_d) / (2 * a))))
def _relative_decimal_error(actual, expected):
with localcontext() as context:
context.prec = 120
return abs((Decimal.from_float(actual) - expected) / expected)
def _raises(expected, function, *args):
try:
function(*args)
except expected:
return True
return False
def run_example():
main = (1.0, 1e16, 1.0)
stable = real_quadratic_roots(*main)
naive = _naive(*main)
reference = _decimal_reference(*main)
stable_error = _relative_decimal_error(stable[1], reference[1])
naive_error = _relative_decimal_error(naive[1], reference[1])
root_fraction = Fraction.from_float(stable[1])
exact_residual = root_fraction ** 2 + Fraction(10 ** 16) * root_fraction + 1
flipped_b = real_quadratic_roots(1, -1e16, 1)
known_results = []
for r1, r2 in itertools.combinations_with_replacement((-7.0, -2.0, 0.0, 0.5, 3.0), 2):
for scale in (0.5, 1.0, 2.0, 10.0):
coefficients = (scale, -scale * (r1 + r2), scale * r1 * r2)
observed = real_quadratic_roots(*coefficients)
oracle = _decimal_reference(*coefficients)
known_results.append(all(math.isclose(x, float(y), rel_tol=1e-14, abs_tol=1e-14)
and math.isclose(x, expected, rel_tol=1e-14, abs_tol=1e-14)
for x, y, expected in zip(observed, oracle, sorted((r1, r2)))))
c_below = math.nextafter(1.0, 0.0)
near_double = real_quadratic_roots(1, 2, c_below)
near_oracle = _decimal_reference(1, 2, c_below)
invalid_types, invalid_numbers = [], []
for index in range(3):
for value in (True, '1', 1 + 0j, None):
values = [1, -5, 6]
values[index] = value
invalid_types.append(_raises(TypeError, real_quadratic_roots, *values))
for value in (float('nan'), float('inf'), -float('inf'), 1e101, 1e-101, 10 ** 400):
values = [1, -5, 6]
values[index] = value
invalid_numbers.append(_raises(ValueError, real_quadratic_roots, *values))
checks = {
'hand_dimensionless_roots_two_and_three': real_quadratic_roots(1, -5, 6) == (2.0, 3.0),
'direct_formula_cancellation_loses_small_root': naive[1] == 0.0 and stable[1] != 0.0
and math.isclose(stable[1], -1e-16, rel_tol=1e-15),
'small_root_relative_error_compared_with_independent_decimal120': stable_error < Decimal('1e-15') and naive_error == 1,
'large_root_matches_independent_decimal_reference': _relative_decimal_error(stable[0], reference[0]) < Decimal('1e-15'),
'vieta_sum_and_product_for_cancellation_fixture': math.isclose(math.fsum(stable), -1e16)
and math.isclose(stable[0] * stable[1], 1, rel_tol=1e-15),
'exact_fraction_residual_is_small_but_not_used_as_only_oracle': 0 < abs(exact_residual) < Fraction(1, 10 ** 15)
and _relative_decimal_error(stable[1], reference[1]) < Decimal('1e-15'),
'copysign_branch_handles_negative_b': math.isclose(flipped_b[0], 1e-16, rel_tol=1e-15)
and flipped_b[1] == 1e16,
'repeated_nonzero_root_returned_twice': real_quadratic_roots(1, -2, 1) == (1.0, 1.0),
'zero_discriminant_and_b_zero_no_division_by_q_zero': real_quadratic_roots(1, 0, 0) == (0.0, 0.0),
'zero_constant_keeps_one_zero_root': real_quadratic_roots(2, -4, 0) == (0.0, 2.0),
'opposite_sign_roots_b_zero': real_quadratic_roots(1, 0, -4) == (-2.0, 2.0),
'degenerate_linear_case_and_negative_discriminant_rejected': _raises(ValueError, real_quadratic_roots, 0, 1, 1)
and _raises(ValueError, real_quadratic_roots, 1, 0, 1),
'invalid_types_bool_and_nonfinite_outside_bounds_rejected': all(invalid_types) and all(invalid_numbers),
'common_positive_coefficient_scaling_keeps_hand_roots': all(real_quadratic_roots(s, -5 * s, 6 * s) == (2.0, 3.0)
for s in (0.5, 2, 10)),
'common_negative_scaling_keeps_hand_roots': real_quadratic_roots(-1, 5, -6) == (2.0, 3.0),
'known_root_family_and_decimal_reference_sixty_cases': len(known_results) == 60 and all(known_results),
'near_repeated_represented_coefficients_checked_with_decimal_reference': near_double[0] < near_double[1]
and all(math.isclose(x, float(y), rel_tol=1e-14, abs_tol=1e-14) for x, y in zip(near_double, near_oracle))
and _raises(ValueError, real_quadratic_roots, 1, 2, math.nextafter(1.0, math.inf)),
'inclusive_coefficient_bounds_and_numeric_root_sort': real_quadratic_roots(1e-100, 0, -1e-100) == (-1.0, 1.0)
and real_quadratic_roots(1e100, 0, -1e100) == (-1.0, 1.0)
and stable[0] <= stable[1] and flipped_b[0] <= flipped_b[1],
}
assert len(checks) == 18 and all(type(v) is bool and v for v in checks.values())
return {'runtime': {'python': platform.python_version(), 'libraries': ['stdlib math/Fraction/Decimal only; no MATLAB/model or circuit execution']},
'contract': {'equation': 'a*u²+b*u+c=0, dimensionless normalized u and coefficients',
'inputs': 'exact int/float converted to float first; nonzero coefficient abs in[1e-100,1e100]; a!=0',
'discriminant': 'exact rational sign for actual represented floats, then rounded once to float',
'roots': 'real only, sorted numerically, repeated root counted twice; q formula and c/q',
'oracle': 'independent direct formula Decimal.from_float, precision120 for selected fixtures'},
'metrics': {'coefficients': main, 'stable_roots': stable, 'naive_roots': naive,
'decimal_reference_roots': [str(x) for x in reference],
'stable_small_root_relative_error': str(stable_error), 'naive_small_root_relative_error': str(naive_error),
'exact_residual_at_stable_small_root': float(exact_residual),
'known_root_reference_cases': len(known_results)},
'checks': checks,
'limitations': ['Authored hybrid exact-discriminant/q implementation, not a verbatim Kahan algorithm, universal correctly-rounded solver or speed benchmark.',
'Decimal reference applies to actual converted float coefficients, not exact user decimal spelling or arbitrary original large ints.',
'Bounded real coefficients only; linear/complex roots, interval uncertainty and physical unit inference are excluded.',
'No formula removes sensitivity near repeated roots or uncertainty in coefficients; small residual alone is not a general forward-error guarantee.']}
if __name__ == '__main__':
print(json.dumps(run_example(), ensure_ascii=False, indent=2))
برای D دقیقاً صفر، ریشهٔ تکراری از منفی b تقسیم بر ۲a برگردانده می‌شود؛ حالت b و c صفر نیز به تقسیم بر q صفر نمی‌رسد. bool، NaN، بی‌نهایت و ورودی نامعتبر رد می‌شوند. اگر D منفی باشد، تابع آن را با تلورانس دلخواه صفر نمی‌کند. در نزدیکی ریشهٔ تکراری، تفاوت کوچک ضرایب می‌تواند معنای مسئله را عوض کند؛ نباید با گردکردنِ پنهان، پاسخ حقیقیِ ساختگی ایجاد کنیم.
مرجع جدا از Decimal.from_float و localcontext با دقت ۱۲۰ رقم استفاده می‌کند و فرمول مستقیم را در آن دقت اجرا می‌کند. این سازنده، float واقعی را می‌گیرد؛ با ساختن Decimal از املای اعشاریِ دلخواه یکسان نیست. Fraction.from_float نیز همین تمایزِ مقدار نمایش‌داده‌شده را دارد. مقایسهٔ دو روش باید دربارهٔ یک مسئلهٔ ورودی باشد.
۴. معیار پذیرش با جواب مستقل و محدودیت باقیمانده
کد با Python 3.14.5 اجرا شده است. در مثال نامتوازن، روش q ریشهٔ کوچک منفی ده به توان منفی شانزده می‌دهد؛ خطای نسبی نسبت به مرجع Decimal حدود ۲٫۰۹ ضرب‌در ده به توان منفی هفده است. روش مستقیم، صفر داده و خطای نسبی آن یک است. این عددها نتیجهٔ همین ورودی و اجرا هستند؛ تضمین عمومیِ خطای همهٔ ریشه‌ها یا همهٔ دستگاه‌ها نیستند.
{
"runtime": {
"python": "3.14.5",
"libraries": [
"stdlib math/Fraction/Decimal only; no MATLAB/model or circuit execution"
]
},
"contract": {
"equation": "a*u²+b*u+c=0, dimensionless normalized u and coefficients",
"inputs": "exact int/float converted to float first; nonzero coefficient abs in[1e-100,1e100]; a!=0",
"discriminant": "exact rational sign for actual represented floats, then rounded once to float",
"roots": "real only, sorted numerically, repeated root counted twice; q formula and c/q",
"oracle": "independent direct formula Decimal.from_float, precision120 for selected fixtures"
},
"metrics": {
"coefficients": [
1.0,
1e+16,
1.0
],
"stable_roots": [
-1e+16,
-1e-16
],
"naive_roots": [
-1e+16,
0.0
],
"decimal_reference_roots": [
"-9999999999999999.9999999999999998999999999999999999999999999999989999999999999999999999999999999800000000000000000000000",
"-1.0000000000000000000000000000000100000000000000000000000000000002000000000000000000000000E-16"
],
"stable_small_root_relative_error": "2.09022132759654064381588850591530545415254830202726184056657811459382045899767159861841663303579297096490036287856144772E-17",
"naive_small_root_relative_error": "1",
"exact_residual_at_stable_small_root": 2.0902213275965406e-17,
"known_root_reference_cases": 60
},
"checks": {
"hand_dimensionless_roots_two_and_three": true,
"direct_formula_cancellation_loses_small_root": true,
"small_root_relative_error_compared_with_independent_decimal120": true,
"large_root_matches_independent_decimal_reference": true,
"vieta_sum_and_product_for_cancellation_fixture": true,
"exact_fraction_residual_is_small_but_not_used_as_only_oracle": true,
"copysign_branch_handles_negative_b": true,
"repeated_nonzero_root_returned_twice": true,
"zero_discriminant_and_b_zero_no_division_by_q_zero": true,
"zero_constant_keeps_one_zero_root": true,
"opposite_sign_roots_b_zero": true,
"degenerate_linear_case_and_negative_discriminant_rejected": true,
"invalid_types_bool_and_nonfinite_outside_bounds_rejected": true,
"common_positive_coefficient_scaling_keeps_hand_roots": true,
"common_negative_scaling_keeps_hand_roots": true,
"known_root_family_and_decimal_reference_sixty_cases": true,
"near_repeated_represented_coefficients_checked_with_decimal_reference": true,
"inclusive_coefficient_bounds_and_numeric_root_sort": true
},
"limitations": [
"Authored hybrid exact-discriminant/q implementation, not a verbatim Kahan algorithm, universal correctly-rounded solver or speed benchmark.",
"Decimal reference applies to actual converted float coefficients, not exact user decimal spelling or arbitrary original large ints.",
"Bounded real coefficients only; linear/complex roots, interval uncertainty and physical unit inference are excluded.",
"No formula removes sensitivity near repeated roots or uncertainty in coefficients; small residual alone is not a general forward-error guarantee."
]
}برای کنترل مستقل، شصت حالت از خانوادهٔ ریشه‌های معلوم می‌سازیم و هم با جواب ساخته‌شده و هم با مرجع Decimal مقایسه می‌کنیم. همهٔ این حالت‌ها یک گروه از هجده گروه نام‌دارند. تغییر علامت b، صفر بودن c، دو ریشهٔ مخالف، ریشهٔ تکراری، ضرب مشترک ضرایب در عدد مثبت یا منفی و ترتیب خروجی، گروه‌های معنادار دیگری هستند. تکرار اجرای همان برنامه، گروه تازه حساب نمی‌شود.
باقیماندهٔ چندجمله‌ای در ریشهٔ کوچک را نیز با حساب کسریِ دقیقِ همان مقدار خروجی می‌سنجیم؛ حدود ۲٫۰۹ ضرب‌در ده به توان منفی هفده است. این بررسی تنها مرجع ما نیست. باقیماندهٔ کوچک به‌تنهایی در مسئلهٔ بدشرط، خطای کوچکِ مکان ریشه را ثابت نمی‌کند؛ ارزیابی باقیمانده با همان حساب کم‌دقت هم ممکن است خطای دیگری پنهان کند. جوابِ مستقل و خطای ریشه را جدا نگه می‌داریم.
برای آزمون نزدیک به ریشهٔ تکراری، c را یک float مجاورِ کمتر از یک و سپس مجاورِ بیشتر از یک می‌گیریم، با a یک و b دو. اولی دو ریشهٔ حقیقی متمایز دارد و دومی در قرارداد ما رد می‌شود. این آزمایش حساسیت ضرایب را حذف نمی‌کند. قرارداد bisection مدار معیار توقف متفاوتی دارد و آزمون expm1 ضریب فیلتر نمونهٔ دیگری از انتخاب عبارت عددی مناسب است.
منابع در ۲۴ مهر ۱۴۰۵ خوانده شدند؛ یادداشت Kahan تاریخ بیستم نوامبر ۲۰۰۴ و وضعیت متنِ در حال کار دارد. مقالهٔ اصلی Goldberg مربوط به مارس ۱۹۹۱ است. نسخهٔ مستندات پایتون 3.14.8 از نسخهٔ اجرای کد جدا ثبت شده است؛ این درس گزارش مقالهٔ تازه یا ادعای برتری سرعت حل‌گر نیست.
۵. چهار پرسش با پاسخ کامل
پرسش اول: صفرِ روش مستقیم برای ریشهٔ کوچکِ مثال پذیرفتنی است؟ پاسخ: خیر؛ مرجع مستقل نزدیک منفی ده به توان منفی شانزده است. صفر در این مقایسه خطای نسبی یک دارد؛ روش q در همین اجرا افت ناشی از آن تفریق را دور می‌زند.
پرسش دوم: برای a یک، b و c صفر، c/q محاسبه کنیم؟ پاسخ: خیر؛ ممیز دقیق صفر است و شاخهٔ ریشهٔ تکراری، دو صفر برمی‌گرداند. حالت ویژه را پیش از تقسیم نامعتبر جدا می‌کنیم.
پرسش سوم: چرا مرجع Decimal را از رشتهٔ اعداد دلخواه نمی‌سازیم؟ پاسخ: چون قرارداد فعلی دربارهٔ floatهای تبدیل‌شده است. املای اعشاریِ دقیق می‌تواند چندجمله‌ای دیگری تعریف کند؛ مرجع باید با ورودی واقعیِ تابع منطبق باشد.
پرسش چهارم: باقیماندهٔ کوچک، همهٔ دقت را ثابت می‌کند؟ پاسخ: خیر؛ حساسیت ریشه به ضرایب و خطای خودِ ارزیابی باقیمانده مهم‌اند. این درس آن را همراه جواب مستقل و خطای نسبیِ ریشه کنترل می‌کند، نه به‌عنوان تضمین عمومیِ دقت.
منابع و مطالعهٔ بیشتر
- W. Kahan؛ Qdrtcs، انتخاب علامت و حساسیت ریشه‌های نزدیک
- D. Goldberg؛ حذف عددی در محاسبات ممیز شناور
- Python 3.14؛ Decimal.from_float و localcontext
- Python 3.14؛ Fraction.from_float
بررسی منابع: ۲۴ مهر ۱۴۰۵. تحلیل‌ها و پیشنهادهای اجرایی، جمع‌بندی نویسنده‌اند.
