Floating-Point Precision — যখন ছোট্ট Error জমে বিপর্যয় হয়

আগে এটা বুঝি

একটা লাইন চালান, যেকোনো ভাষায়:

>>> 0.1 + 0.2
0.30000000000000004

0.30000000000000004। শূন্যের পর 4। এটা কি bug? কোনো compiler-এর ভুল? একটা অদ্ভুত Python-নির্দিষ্ট সমস্যা?

গত লেসনের পর এই প্রশ্নের উত্তর আপনার কাছে ইতিমধ্যেই আছে — আপনি জানেন 0.1-কে IEEE 754-এ store করলে আসলে যা জমা থাকে তা ঠিক 0.1 নয়। এই লেসনে আমরা সেই জ্ঞানটা সম্পূর্ণ করব — দেখাব ঠিক কেন যোগফলটা 0.30000000000000004 হয়, বিট-বাই-বিট।

কিন্তু এটা শুধু একটা কৌতূহল নয়। আসল প্রশ্ন হলো: এই ছোট্ট error যখন হাজার-লক্ষ-কোটিবার arithmetic operation-এর মধ্য দিয়ে যায়, তখন কী হয়? কখনো তারা একে অপরকে প্রায় বাতিল করে দেয় — নিরাপদ। কখনো একটা একক subtraction-ই নির্ভুলতা সম্পূর্ণ ধ্বংস করে দেয় — catastrophic। আর কখনো তারা নীরবে জমতে থাকে — কয়েক সেকেন্ডের ঘড়ির টিক, একটার পর একটা, একটানা চলতে থাকা কোনো সিস্টেমে — যতক্ষণ না একটা ক্ষেপণাস্ত্র প্রতিরক্ষা ব্যবস্থা তার লক্ষ্য হারিয়ে ফেলে।

এই লেসনের কাজ: এই তিনটা পরিস্থিতির মধ্যে পার্থক্য করতে শেখা, আর যেখানে ক্ষতি হয় সেখানে তা ঠেকানোর হাতিয়ার শেখা। “কতটা ছোট” আর “কতটা নিরীহ” এক জিনিস নয় — আর সেই পার্থক্যটাই আজকের মূল পাঠ।

মূল ধারণা

`0.1 + 0.2 ≠ 0.3` — সম্পূর্ণ প্রমাণ

গত লেসনে আমরা 0.1f (single precision)-এর bit pattern হাতে বের করেছিলাম। এখানে double precision (Python-সহ বেশিরভাগ ভাষার default float) ব্যবহার করি, কারণ 0.1 + 0.2 ≠ 0.3 উদাহরণটা double-এই সবচেয়ে বেশি দেখা যায়।

একই প্রক্রিয়া (normalize → বাইনারি ভগ্নাংশ বের করা → ৫২ বিটে round করা) double precision-এ প্রয়োগ করলে এই তিনটা exact stored মান পাওয়া যায়:

Expression	Store হওয়া exact decimal মান
`0.1`	`0.1000000000000000055511151231257827021181583404541015625`
`0.2`	`0.200000000000000011102230246251565404236316680908203125`
`0.1 + 0.2` (computed)	`0.3000000000000000444089209850062616169452667236328125`
`0.3` (stored directly)	`0.299999999999999988897769753748434595763683319091796875`

লক্ষ্য করুন — 0.2-এর exact মান ঠিক 0.1-এর exact মানের দ্বিগুণ। এটা কাকতালীয় নয়: × 2 মানে শুধু exponent-এ 1 যোগ করা (গত লেসনের bias trick), mantissa অপরিবর্তিত থাকে — তাই দ্বিগুণ করাটা সবসময় exact, কোনো নতুন rounding ছাড়াই।

0.1 + 0.2 যখন সত্যিকারের arithmetic দিয়ে যোগ হয় (দুটো exact stored মান যোগ করে, তারপর নিকটতম representable double-এ round করে), ফলাফল হয় 0.3-এর exact stored মান থেকে ঠিক এক ULP বেশি — পার্থক্য ঠিক 2^{-54} \approx 5.551 \times 10^{-17}।

0.30000000000000004-এর যে repr আমরা দেখি, সেটা এই 0.3000000000000000444089...-এর shortest round-tripping decimal string — display-এর জন্য বাছা প্রতিনিধি, raw মান নয়।

Error-এর ভাষা — absolute বনাম relative

একটা computed মান \hat{x}, প্রকৃত মান x-এর তুলনায়:

$\text{absolute error} = |\hat{x} - x|$

$\text{relative error} = \frac{|\hat{x} - x|}{|x|}$

Floating-point-এর precision relative, absolute নয় — এটাই exponent trick-এর সরাসরি ফলাফল। 1.0-এর কাছে দুই representable মানের ব্যবধান (ULP) 2^{-52} (double), কিন্তু 1,000,000,000.0-এর কাছে ব্যবধান অনেক বড় — magnitude-এর সমানুপাতিক।

Machine epsilon

$\varepsilon = \text{১.০ আর তার ঠিক পরের representable float-এর ব্যবধান}$

Precision	`ε`	আনুমানিক দশমিক
single (`float32`)	`2^{-23}`	`1.1920929 \times 10^{-7}`
double (`float64`)	`2^{-52}`	`2.220446049250313 \times 10^{-16}`

ε একটা universal ধ্রুবক নয় — এটা 1.0-এর কাছের precision। যেকোনো সংখ্যা x-এর কাছে ULP আনুমানিক |x| \times \varepsilon (যদি x কোনো subnormal range-এ না থাকে)। এই সূত্রটাই পরে “epsilon দিয়ে তুলনা” অংশে গুরুত্বপূর্ণ হয়ে উঠবে।

কতটা বদলায় সেটা concrete সংখ্যায় দেখা যাক — double precision-এ, বিভিন্ন magnitude-এ পরপর দুই representable মানের প্রকৃত ব্যবধান (exact binade হিসাব থেকে, শুধু আনুমানিক নয়):

| Magnitude (|x|) | Binade exponent | ULP (exact) | |---|---:|---| | \approx 1 | 0 | 2^{-52} \approx 2.22 \times 10^{-16} | | \approx 10^3 | 9 | 2^{-43} \approx 1.14 \times 10^{-13} | | \approx 10^6 | 19 | 2^{-33} \approx 1.16 \times 10^{-10} | | \approx 10^9 | 29 | 2^{-23} \approx 1.19 \times 10^{-7} | | \approx 10^{15} | 49 | 2^{-3} = 0.125 | | \approx 10^{20} | 66 | 2^{14} = 16,384 |

শেষ সারিটা লক্ষ্য করুন — 10^{20}-এর কাছে দুই representable double-এর ব্যবধান ১৬ হাজারের বেশি। এই magnitude-এ কোনো হিসাবের ফলাফলে 1.0 (বা তার কম) নির্ভুলতা আশা করাটাই ভুল — সেই নির্ভুলতা double precision-এ শারীরিকভাবে সম্ভবই না, algorithm যত ভালোই হোক।

একটা single arithmetic operation-এর error bound: IEEE 754 গ্যারান্টি দেয় প্রতিটা মৌলিক operation (+, -, ×, ÷, sqrt) correctly rounded — মানে ফলাফল ঠিক সেই representable মান যা প্রকৃত গাণিতিক ফলাফলের সবচেয়ে কাছে (round-to-nearest-even নিয়মে)। তার মানে একটা একক operation-এর relative error সবসময় \varepsilon/2-এর মধ্যে সীমাবদ্ধ।

সমস্যাটা একক operation-এ নয় — সমস্যা তখন হয় যখন অনেক operation একসাথে হয়, আর তাদের error একে অপরের সাথে interact করে।

এই “correctly rounded” গ্যারান্টিটাই লেসন ৮-এর round-to-nearest-even নিয়মের সরাসরি ফলাফল, আর এটা IEEE 754-এর একটা অসাধারণ শক্তি — প্রতিটা একক +, -, ×, ÷ operation-এর ফলাফল সম্পূর্ণ predictable, প্ল্যাটফর্ম-নিরপেক্ষ। যা predictable নয়, তা হলো বহু operation-এর সমষ্টিগত প্রভাব — আর সেটাই এই লেসনের বাকি অংশের বিষয়।

ভেতরে কী ঘটছে

যেখানে error জমে, বাতিল হয়, বা বিস্ফোরিত হয়

১. Catastrophic cancellation — তত্ত্ব

যখন দুইটা প্রায়-সমান, বড় floating-point সংখ্যা বিয়োগ করা হয়, ফলাফলের absolute error মোটামুটি একই থাকে (বিয়োগ operation নিজে correctly rounded), কিন্তু ফলাফলের magnitude হঠাৎ অনেক ছোট হয়ে যায় — কারণ leading (গুরুত্বপূর্ণ) digit-গুলো একে অপরকে বাতিল করে দেয়।

$\text{relative error} = \frac{\text{absolute error (অপরিবর্তিত)}}{\text{ফলাফলের magnitude (হঠাৎ ছোট)}}$

Absolute একই, denominator ছোট — relative error বিস্ফোরিত হয়। যে low-order বিটগুলো আগে “নগণ্য noise” ছিল, বিয়োগের পর সেগুলোই ফলাফলের প্রধান অংশ হয়ে যায়।

এটা কোনো নতুন error তৈরি করে না — বিয়োগ operation নিজে সঠিক। সমস্যাটা হলো এটা upstream থেকে আসা আগে থেকেই বিদ্যমান error-কে উন্মোচিত করে, magnify করে দেয়।

a = 1.234567 8[9]     ← শেষ digit-টা (বন্ধনীতে) আগে থেকেই সন্দেহজনক
b = 1.234567 3[2]     ← এটাও

a-এর leading digits "1.234567" আর b-এর "1.234567" — অভিন্ন,
তাই বিয়োগে বাতিল:

  1.234567 8[9]
− 1.234567 3[2]
──────────────────
  0.000000 5[7]      ← শুধু এই অংশটাই বাকি — যা আগে ছিল
                        নগণ্য শেষ digit, এখন এটাই *পুরো* উত্তর

বিয়োগের আগে leading digit-গুলো 'মূল্যবান তথ্য' মনে হচ্ছিল — বিয়োগের পর সেগুলো বাতিল, আর যা বাকি থাকল সেটা প্রায় পুরোটাই আগে থেকে জমে থাকা rounding noise।

২. Naive summation — জমতে থাকা error

n-টা সংখ্যা পরপর যোগ করলে (s = s + x_i, লুপে), প্রতিটা যোগে সর্বোচ্চ \varepsilon/2 relative error যোগ হতে পারে। Worst case-এ এই error-গুলো জমতে জমতে বাড়ে:

$\text{total relative error} = O(n \cdot \varepsilon)$

(Goldberg-এর পেপার আর Higham-এর টেক্সটবই — দুটোই এই bound-এর প্রামাণ্য উৎস, দুটোই আরও কঠোর, order-নির্ভর সংস্করণ দেয়।)

n ছোট হলে এটা অদৃশ্য। কিন্তু n = 10^7, 10^9 — বৈজ্ঞানিক simulation, financial aggregation, বা machine learning training-এ স্বাভাবিক সংখ্যা — এই bound বাস্তব হয়ে ওঠে।

একটা ছোট, সহজে যাচাইযোগ্য দৃষ্টান্ত — এমনকি মাত্র ১০ বার:

>>> 0.1 * 10 == 1.0
True
>>> sum([0.1] * 10) == 1.0
False
>>> sum([0.1] * 10)
0.9999999999999999

আশ্চর্যজনক — একটা single multiplication (0.1 \times 10) সঠিক 1.0-এ round হয় (কারণ সেই একটা operation-এর error half-ULP-এর নিচে), কিন্তু দশবার repeated addition-এ প্রতিটা ধাপের ছোট rounding জমে গিয়ে চূড়ান্ত ফলাফল 1.0-এর এক ULP নিচে থামে। একই গাণিতিক প্রশ্ন, ভিন্ন algorithm, ভিন্ন উত্তর — এটাই floating-point arithmetic-এ non-associativity-র সরাসরি ফল ((a+b)+c \ne a+(b+c) সবসময় সত্য নয়)।

৩. Kahan summation — compensation দিয়ে error পুনরুদ্ধার

William Kahan (হ্যাঁ, IEEE 754-এর সেই একই স্থপতি) ১৯৬৫ সালে একটা elegant কৌশল প্রকাশ করেন — প্রতিটা addition-এ যে বিট “হারিয়ে যায়” (rounding-এর কারণে), সেটা আলাদাভাবে track করে পরের addition-এ ফিরিয়ে দেওয়া:

sum = 0.0
c   = 0.0          # কতটা "হারিয়েছি" তার running হিসাব

প্রতিটা x-এর জন্য:
    y = x - c      # হারানো অংশ পুষিয়ে নিয়ে x সংশোধন করা
    t = sum + y    # স্বাভাবিক addition (এখানেও rounding হবে)
    c = (t - sum) - y   # ঠিক কতটা হারালাম, সেটা বের করে রাখা
    sum = t

মূল কৌশলটা c = (t - sum) - y লাইনে — t - sum (গাণিতিকভাবে) বলে দেয় addition-এ আসলে sum-এ কতটুকু যোগ হলো; সেটা থেকে y (যা যোগ করতে চেয়েছিলাম) বাদ দিলে বাকি থাকে ঠিক ততটুকু যা rounding-এ হারিয়ে গেছে। পরের iteration-এ y = x - c সেই হারানো অংশ ফিরিয়ে নিয়ে আসে।

ফলাফল: error bound O(n\varepsilon) থেকে নেমে আসে প্রায় O(\varepsilon) — n-এর উপর নির্ভরই করে না। হাজার হোক বা কোটি সংখ্যা, error প্রায় স্থির থাকে।

Kahan summation-ও অভ্রান্ত নয়

Kahan summation দুর্দান্ত কিন্তু সর্বশক্তিমান নয় — কিছু বিশেষভাবে সাজানো, চরম magnitude-ব্যবধানের ইনপুটে (যেমন খুব বড় আর খুব ছোট সংখ্যা পর্যায়ক্রমে, বিশেষ ধরনের tie-তে আটকে গেলে) এর সুবিধা কমে যেতে পারে। আরও শক্তিশালী সংস্করণ আছে — Neumaier summation (Kahan-Babuška variant), যেটা |sum| বনাম |y|-এর তুলনা করে আরও সাবধানে compensation হিসাব করে। ব্যবহারিক নিয়ম: n বড় হলে আর precision গুরুত্বপূর্ণ হলে, শুধু sum += x না লিখে অন্তত Kahan (আরও ভালো, Neumaier) ব্যবহার করুন — কিন্তু “একদম নিখুঁত” এর গ্যারান্টি নয়, “অনেক বেশি নির্ভরযোগ্য”-র গ্যারান্টি।

Neumaier summation-এর মূল পার্থক্য — Kahan-এ ধরে নেওয়া হয় sum সবসময় y-এর চেয়ে বড় (তাই y = x - c করলেই যথেষ্ট)। কিন্তু বাস্তবে মাঝে মাঝে নতুন x আগের sum-এর চেয়ে বড় হতে পারে (বিশেষত শুরুর দিকে যখন sum ছোট) — তখন Kahan-এর অনুমান ভেঙে পড়ে। Neumaier sum আর x-এর মধ্যে কোনটা বড় সেটা প্রতি ধাপে পরীক্ষা করে, compensation-টা সঠিক দিক থেকে হিসাব করে:

def neumaier_sum(values):
    total = 0.0
    c = 0.0                     # compensation
    for x in values:
        t = total + x
        if abs(total) >= abs(x):
            c += (total - t) + x    # low-order bits of x হারিয়েছে
        else:
            c += (x - t) + total    # low-order bits of total হারিয়েছে
        total = t
    return total + c            # শেষে compensation যোগ করা

লক্ষ্য করুন Kahan-এর সাথে একটা গুরুত্বপূর্ণ পার্থক্য — এখানে c লুপের ভেতরে y-তে ফিরিয়ে দেওয়া হয় না, বরং শেষে একবারে total + c করে যোগ করা হয়।

৪. Patriot missile (১৯৯১) — জমতে থাকা error-এর একটা প্রকৃত বিপর্যয়

এই ঘটনাটা প্রায়ই ভুলভাবে বলা হয় “একটা float == তুলনার bug”। প্রকৃত mechanism সম্পূর্ণ ভিন্ন — আর সেটাই এই লেসনের সবচেয়ে গুরুত্বপূর্ণ ইতিহাস।

পটভূমি: ২৫ ফেব্রুয়ারি, ১৯৯১, উপসাগরীয় যুদ্ধ। সৌদি আরবের দাহরানে একটা US Patriot air-defense battery একটা আগত Scud ক্ষেপণাস্ত্র ট্র্যাক করতে ব্যর্থ হয় — এটা একটা US army ব্যারাকে আঘাত করে, ২৮ জন সেনা নিহত হন।

System-এর ঘড়ি: Patriot-এর অভ্যন্তরীণ system clock সময় গোনে ০.১ সেকেন্ড একক-এ, একটা integer counter দিয়ে। এই integer counter-কে সেকেন্ডে রূপান্তর করতে হতো (target-এর ভবিষ্যৎ অবস্থান হিসাব করার জন্য), যার জন্য counter-কে 0.1-এর একটা ২৪-বিট fixed-point approximation দিয়ে গুণ করা হতো (গত লেসনের fixed-point লেসন মনে করুন — এটা IEEE 754 float নয়, একটা সাধারণ ২৪-বিট truncated binary fraction)।

সমস্যাটা কোথায়: 0.1-এর বাইনারি expansion অসীম (লেসন ৭-এ প্রমাণিত — 0.0(0011) repeating)। ২৪ বিটে chop (truncate) করলে — round নয়, শুধু কেটে ফেলা — একটা ছোট কিন্তু নির্দিষ্ট error থেকে যায়:

$\text{chopping error} \approx 0.000000095 \text{ সেকেন্ড প্রতি } 0.1\text{-সেকেন্ড tick-এ}$

কেন জমল: Patriot battery-গুলো নিয়মিত reboot করার কথা ছিল (যা error-কে ছোট রাখত), কিন্তু যুদ্ধের সময় দাহরানের battery-টা প্রায় ১০০ ঘণ্টা একটানা চালু ছিল, reboot ছাড়াই।

$100 \text{ ঘণ্টা} = 360,000 \text{ সেকেন্ড} = 3,600,000 \text{ টা } 0.1\text{-সেকেন্ড tick}$

$\text{মোট accumulated error} \approx 3,600,000 \times 0.000000095 \approx 0.342 \text{ সেকেন্ড}$

Scud missile-এর গতি প্রায় 1,676 m/s (Mach ৫)। এই সময়ের error-কে দূরত্বে রূপান্তর করি:

$0.342 \text{ s} \times 1,676 \text{ m/s} \approx 573 \text{ মিটার}$

Radar একটা “range gate” ব্যবহার করে — target ঠিক কোথায় থাকার কথা তার চারপাশে একটা সরু search window, noise বাদ দিতে। ৫৭০+ মিটার ভুল predicted অবস্থান radar-এর range gate-এর বাইরে ফেলে দিল — system Scud-কে কখনো একটা বৈধ target হিসেবে চিনতেই পারল না, তাই কোনো interceptor ছোঁড়া হয়নি।

কেন এটা 'float comparison bug' নয়

কোথাও কোনো if a == b তুলনা এই ব্যর্থতার কারণ ছিল না। এটা ছিল সময়ের সাথে জমতে থাকা একটা truncation error, একটা fixed-point (binary, না IEEE 754 float) representation-এ — ঠিক যে ধরনের জমতে থাকা error এই লেসনের “naive summation” অংশে আমরা তাত্ত্বিকভাবে দেখেছি, কিন্তু এখানে প্রতিটা “addition” ছিল একটা ঘড়ির টিক, আর জমার সময় ছিল ১০০ ঘণ্টা।

শিক্ষাটা তাই সাধারণীকৃত: যেকোনো systematic (একদিকে পক্ষপাতযুক্ত) error, যথেষ্ট সময় বা যথেষ্ট repetition পেলে, বিপর্যয়কর হয়ে উঠতে পারে — এমনকি প্রতিটা একক error অকল্পনীয়ভাবে ছোট হলেও। এটাই এই দুর্ঘটনার আসল শিক্ষা, শুধু “float সাবধানে ব্যবহার করুন” নয়।

উদাহরণ

সম্পূর্ণ hand-worked উদাহরণ — Catastrophic Cancellation

Quadratic formula-র চেয়ে ভালো classic উদাহরণ নেই (Goldberg-এর পেপারেই এই উদাহরণ আছে)। সমীকরণ:

$x^2 - 100000x + 1 = 0$

অর্থাৎ a=1, b=-100000, c=1। প্রকৃত (গাণিতিকভাবে সঠিক) মূল দুটো আনুমানিক x_1 \approx 0.00001 (ছোট) আর x_2 \approx 99999.99999 (বড়)।

আমরা এটা float32-এ (single precision) হিসাব করব, যেখানে ~৭ দশমিক digit নির্ভরযোগ্য — গত লেসনের precision সূত্র মনে করুন।

ধাপ ১ — b² হিসাব। 100000² = 10,000,000,000 (= 10^{10})। এই সংখ্যা কি float32-এ exact? 10^{10} = 1024 \times 9,765,625, আর 9,765,625 \lt 2^{24} — হ্যাঁ, exactly representable।

ধাপ ২ — discriminant। `disc = b^2 - 4ac = 10,000,000,000

4 = 9,999,999,996। কিন্তু এই magnitude-এ (\sim 10^10, 2^33আর2^34-এর মাঝে) float32-এর ULP 2^10 = 2^10 = 1024। -4 সংশোধনটা এই ULP-এর অর্ধেকেরও (512`) অনেক ছোট — rounding-এই হারিয়ে যায়:

$\text{disc}_{f32} = 10,000,000,000.0 \quad \text{(অর্থাৎ } b^2 \text{-এরই সমান!)}$

এখানেই প্রথম precision loss ঘটে গেছে — sqrt নেওয়ার আগেই।

ধাপ ৩ — sqrt। \sqrt{10^{10}} = 100,000 — exactly (যেহেতু disc_f32 নিজেই exactly 10^{10})।

$\sqrt{\text{disc}}_{f32} = 100,000.0$

ধাপ ৪ — naive formula (cancellation-এর শিকার):

$x_1 = \frac{-b - \sqrt{\text{disc}}}{2a} = \frac{100,000.0 - 100,000.0}{2} = \frac{0.0}{2} = 0.0$

সম্পূর্ণ ভুল — প্রকৃত উত্তর \approx 0.00001, কিন্তু আমরা পেলাম ঠিক শূন্য। Relative error: 100\%।

ধাপ ৫ — দ্বিতীয় মূল (কোনো cancellation নেই):

$x_2 = \frac{-b + \sqrt{\text{disc}}}{2a} = \frac{100,000.0 + 100,000.0}{2} = 100,000.0$

এটা প্রকৃত x_2 \approx 99,999.99999-এর খুব কাছাকাছি (relative error \sim 10^{-10}, চমৎকার) — কারণ এখানে যোগ হচ্ছে, দুইটা সমমানের সংখ্যার বিয়োগ নয়, তাই আগের ধাপের ছোট precision loss (disc-এ) magnify হয়নি।

ধাপ ৬ — সমাধান: Vieta’s formula ব্যবহার করা। দুই মূলের গুণফল সবসময় x_1 \times x_2 = c/a। যেহেতু x_2 নির্ভুল, x_1 কে ভাগ দিয়ে বের করা যায়, বিয়োগ ছাড়াই:

$x_1 = \frac{c}{a \cdot x_2} = \frac{1}{1 \times 100,000.0} = 0.00001$

এটাই প্রকৃত মূলের (\approx 0.0000100000...1) সাথে float32 precision পর্যন্ত হুবহু মেলে — cancellation সম্পূর্ণ এড়ানো গেল, শুধু একই তথ্য ভিন্নভাবে সাজিয়ে।

মূল শিক্ষা — গাণিতিকভাবে সমতুল্য ≠ numerically সমতুল্য

(-b - \sqrt{disc})/(2a) আর c / (a \cdot x_2) — গণিতের ক্লাসে এই দুইটা ঠিক একই সংখ্যা বোঝায় (algebraically সমতুল্য, Vieta’s formula-র সরাসরি ফল)। কিন্তু floating-point arithmetic-এ একটা প্রায় ধ্বংস হয়ে গেল (0.0), আরেকটা নির্ভুল রইল (0.00001)।

এটাই এই পুরো লেসনের কেন্দ্রীয় দাবি: floating-point-এ কীভাবে আপনি একটা হিসাব সাজান সেটা কী হিসাব করছেন তার সমান গুরুত্বপূর্ণ। পরের section-এ আপনি নিজেই এটা numpy.float32 দিয়ে চালিয়ে verify করবেন।

নিজে চালিয়ে দেখুন

EXPERIMENT

0.1 + 0.2 বিট-বাই-বিট verify করা

Python 3· ১০ মিনিট

from decimal import Decimal

a, b, c = 0.1, 0.2, 0.3

print("0.1 exact  :", Decimal(a))
print("0.2 exact  :", Decimal(b))
print("0.1+0.2    :", Decimal(a + b))
print("0.3 exact  :", Decimal(c))
print()
print("0.1 + 0.2 == 0.3 :", (a + b) == c)

diff = a + b - c
print("পার্থক্য       :", Decimal(diff))
print("পার্থক্য / ULP :", diff / 2**-54)   # আশা করি প্রায় 1.0

# non-associativity প্রমাণ
print()
print("0.1 * 10 == 1.0        :", 0.1 * 10 == 1.0)
print("sum([0.1]*10) == 1.0   :", sum([0.1] * 10) == 1.0)
print("sum([0.1]*10)          :", sum([0.1] * 10))

প্রত্যাশিত output:

0.1 exact  : 0.1000000000000000055511151231257827021181583404541015625
0.2 exact  : 0.200000000000000011102230246251565404236316680908203125
0.1+0.2    : 0.3000000000000000444089209850062616169452667236328125
0.3 exact  : 0.299999999999999988897769753748434595763683319091796875

0.1 + 0.2 == 0.3 : False
পার্থক্য       : 5.55111512312578270211815834045410156250E-17
পার্থক্য / ULP : 1.0

0.1 * 10 == 1.0        : True
sum([0.1]*10) == 1.0   : False
sum([0.1]*10)          : 0.9999999999999999

Decimal(x) কোনো round করে না — এটা x-এ যা আসলে store আছে তার হুবহু decimal প্রতিরূপ দেখায়। এটাই এই module-এর মূল থিসিসের চূড়ান্ত প্রমাণ: print(0.1) যা দেখায় তা display-এর সময় round করা রূপ, Decimal(0.1) যা দেখায় তা মেমরিতে আসলে যা আছে।

এটা কী প্রমাণ করে

0.1+0.2 আর 0.3-এর মধ্যে পার্থক্য ঠিক এক ULP (2⁻⁵⁴) — এটা এই লেসনের 'concept' অংশে হাতে করা হিসাবের সাথে হুবহু মিলে যায়, শুধু বিশ্বাসের ভিত্তিতে নয়।

EXPERIMENT

Catastrophic cancellation নিজে দেখুন — quadratic formula

Python 3 (numpy) অথবা C· ১৫ মিনিট

import numpy as np

def naive_roots(a, b, c):
    a, b, c = np.float32(a), np.float32(b), np.float32(c)
    disc = b * b - np.float32(4) * a * c
    sq = np.sqrt(disc)
    x1 = (-b - sq) / (np.float32(2) * a)
    x2 = (-b + sq) / (np.float32(2) * a)
    return x1, x2

def stable_roots(a, b, c):
    a, b, c = np.float32(a), np.float32(b), np.float32(c)
    disc = b * b - np.float32(4) * a * c
    sq = np.sqrt(disc)
    # যে মূলে কোনো cancellation নেই সেটা সরাসরি হিসাব করি
    x2 = (-b + sq) / (np.float32(2) * a) if b \< 0 else (-b - sq) / (np.float32(2) * a)
    x1 = c / (a * x2)                    # Vieta's formula — division, subtraction না
    return x1, x2

a, b, c = 1.0, -100000.0, 1.0

x1_n, x2_n = naive_roots(a, b, c)
x1_s, x2_s = stable_roots(a, b, c)

true_x1 = 1.00000000001e-05   # উচ্চ-precision reference মান

print(f"naive : x1 = {x1_n!r}, x2 = {x2_n!r}")
print(f"stable: x1 = {x1_s!r}, x2 = {x2_s!r}")
print(f"true  : x1 ≈ {true_x1}")
print(f"naive-এর relative error  : {abs(x1_n - true_x1) / true_x1:.2%}")
print(f"stable-এর relative error : {abs(x1_s - true_x1) / true_x1:.2e}")

প্রত্যাশিত output:

naive : x1 = 0.0, x2 = 100000.0
stable: x1 = 1.0000000116860974e-05, x2 = 100000.0
true  : x1 ≈ 1.00000000001e-05
naive-এর relative error  : 100.00%
stable-এর relative error : 1.17e-06

naive-এর x1 হুবহু 0.0 — এই লেসনের “example” অংশের হাতের হিসাবের সাথে নিখুঁত মিল। stable-এর x1 প্রকৃত মানের float32 precision-এর মধ্যেই (~1e-6 relative error, ঠিক যেমন ~7 দশমিক digit নির্ভরযোগ্যতা থেকে প্রত্যাশিত)।

C-তে একই পরীক্ষা, float (single precision) দিয়ে:

#include <stdio.h>
#include <math.h>

int main(void) {
    float a = 1.0f, b = -100000.0f, c = 1.0f;
    float disc = b*b - 4.0f*a*c;
    float sq = sqrtf(disc);

    float x1_naive  = (-b - sq) / (2.0f * a);
    float x2         = (-b + sq) / (2.0f * a);
    float x1_stable = c / (a * x2);

    printf("disc          = %.1f\n", disc);
    printf("sqrt(disc)    = %.1f\n", sq);
    printf("naive x1      = %.10f\n", x1_naive);
    printf("stable x1     = %.10f\n", x1_stable);
    return 0;
}

gcc -O2 -o cancel cancel.c -lm && ./cancel

এটা কী প্রমাণ করে

একই গাণিতিক সূত্র, দুই ভিন্ন algebraic রূপে লিখলে, float32-এ সম্পূর্ণ ভিন্ন নির্ভুলতা দেয় — এটাই এই লেসনের 'example' অংশের হাতে করা হিসাবের সরাসরি machine-verified সংস্করণ।

নিজে বানান

BUILD IT

Kahan Summation — accumulated error পুনরুদ্ধার করুন

Python এবং C · ●●●○○

একটা naive_sum() লিখুন — সাধারণ loop-এ যোগ করা
একটা kahan_sum() লিখুন — compensation variable c সহ
দুইটাকে একটা বড়, চ্যালেঞ্জিং input-এ (অনেক ছোট সংখ্যা + কিছু বড় সংখ্যা মিশিয়ে) তুলনা করুন
পার্থক্য measure করুন একটা high-precision reference-এর (Python decimal.Decimal, বা math.fsum) সাথে
দেখুন n বাড়ালে naive_sum-এর error কীভাবে বাড়ে, kahan_sum-এর error কীভাবে প্রায় স্থির থাকে

def naive_sum(values):
    total = 0.0
    for x in values:
        total += x
    return total


def kahan_sum(values):
    total = 0.0
    c = 0.0                      # compensation — হারানো precision-এর হিসাব
    for x in values:
        y = x - c
        t = total + y
        c = (t - total) - y
        total = t
    return total


import random
from decimal import Decimal, getcontext
getcontext().prec = 50

def reference_sum(values):
    """উচ্চ-precision decimal দিয়ে 'সত্যিকারের' যোগফল"""
    return sum(Decimal(x) for x in values)


random.seed(42)
# ইচ্ছাকৃতভাবে কঠিন input — অনেক ছোট মান, মাঝে মাঝে বড়
data = [random.uniform(0, 1) for _ in range(1_000_000)]

n_result = naive_sum(data)
k_result = kahan_sum(data)
ref = reference_sum(data)

print(f"naive  : {n_result!r}")
print(f"kahan  : {k_result!r}")
print(f"reference (decimal): {float(ref)!r}")
print(f"naive absolute error : {abs(Decimal(n_result) - ref)}")
print(f"kahan absolute error : {abs(Decimal(k_result) - ref)}")

চালালে দেখবেন — n = 10^6-এর মতো বড় sample-এ naive_sum-এর error সাধারণত kahan_sum-এর error-এর চেয়ে কয়েক অর্ডার অফ ম্যাগনিটিউড বড় হয়, ঠিক যেমন এই লেসনের “hood” অংশের O(n \varepsilon) বনাম O(\varepsilon) তত্ত্ব বলে। নির্দিষ্ট digit-গুলো random.seed, platform, আর Python-এর সংস্করণ ভেদে সামান্য বদলাতে পারে — নিজে চালিয়ে সংখ্যাগুলো দেখুন।

C সংস্করণ:

#include <stdio.h>

double naive_sum(const double *xs, size_t n) {
    double total = 0.0;
    for (size_t i = 0; i \< n; i++) total += xs[i];
    return total;
}

double kahan_sum(const double *xs, size_t n) {
    double total = 0.0, c = 0.0;
    for (size_t i = 0; i \< n; i++) {
        double y = xs[i] - c;
        double t = total + y;
        c = (t - total) - y;
        total = t;
    }
    return total;
}

Compiler optimization Kahan summation ভাঙতে পারে

-ffast-math-এর মতো aggressive optimization flag কম্পাইলার-কে (t - total) - y-কে algebraically সরল করে ফেলার অনুমতি দিতে পারে (গাণিতিকভাবে এটা প্রায়ই 0-এর সমান!) — যা পুরো compensation কৌশলটাই নষ্ট করে দেয়, কারণ Kahan summation-এর পুরো কাজই নির্ভর করে সেই “গাণিতিকভাবে শূন্য কিন্তু floating-point-এ শূন্য নয়” পার্থক্যের উপর। এই কারণেই numerical library-তে Kahan summation-এর কাছাকাছি কোডে -ffast-math বন্ধ রাখা বা volatile/compiler-specific pragma ব্যবহার করা হয়।

নিজে বাড়ান:

math.fsum() (Python built-in, Shewchuk algorithm ব্যবহার করে, Kahan-এর চেয়েও নির্ভুল) দিয়ে তুলনা করুন
Input-টাকে ইচ্ছাকৃতভাবে “খারাপ” বানান — একটা বিশাল সংখ্যা (1e10) আর তারপর লক্ষ লক্ষ ছোট সংখ্যা (1e-5) মিশিয়ে — naive_sum কতটা খারাপ হয় দেখুন
n-এর বিভিন্ন মানে (10^3, 10^5, 10^7) error measure করে log-log plot করুন — naive-এর error n-এর সাথে বাড়ে কি না যাচাই করুন
Kahan-এর বদলে Neumaier summation (Kahan-Babuška variant) implement করুন — |total| \ge |y| কি না চেক করে compensation term আলাদাভাবে হিসাব করে

বাস্তব সিস্টেমে

যেখানে এই precision সমস্যা সত্যিই ঘটেছে

Vancouver Stock Exchange index (১৯৮২-৮৩)। ১৯৮২ সালের জানুয়ারিতে VSE একটা নতুন সূচক (index) চালু করে, 1000.000 থেকে শুরু। প্রতিটা লেনদেনের পর সূচক পুনর্গণনা করা হতো, কিন্তু সিস্টেমটা ফলাফল round না করে truncate (কেটে ফেলা) করত — তিন দশমিক স্থানের পর বাকিটা ফেলে দিত, প্রতিদিন হাজার হাজার বার। প্রতিটা truncation-এর পক্ষপাত সবসময় একদিকে (নিচের দিকে, কারণ truncation ধনাত্মক সংখ্যাকে সবসময় ছোট করে) — ঠিক যেমন এই লেসনের “systematic bias” নীতি বলে। প্রায় ২২ মাস পর, নভেম্বর ১৯৮৩-এ, সূচক নেমে দাঁড়ায় মাত্র 524.811-এ, যদিও প্রকৃত বাজারমূল্য বেড়েছিল। সংশোধন করে (সঠিক rounding দিয়ে পুরো recalculation) একটা সপ্তাহান্তে সূচক লাফিয়ে ওঠে 1098.892-এ — একরাতে প্রায় দ্বিগুণ, শুধু ২২ মাসের truncation error মুছে ফেলার ফলে।

NumPy-র sum() — pairwise summation। এই লেসনের naive summation তত্ত্ব (O(n\varepsilon) error) এতটাই বাস্তব সমস্যা যে NumPy ২০১৩ সাল (সংস্করণ ১.৯) থেকে তার default sum()-এ সাধারণ left-to-right loop-এর বদলে pairwise (cascade) summation ব্যবহার করে — array-টাকে recursively অর্ধেক-অর্ধেক ভাগ করে, প্রতিটা অর্ধেকের যোগফল আলাদাভাবে বের করে, তারপর সেই দুইটা যোগফল যোগ করে। এতে error bound নেমে আসে O(\log n \cdot \varepsilon)-এ — Kahan summation-এর মতো O(\varepsilon) না হলেও, naive O(n\varepsilon)-এর চেয়ে বহুগুণ ভালো, আর Kahan-এর মতো প্রতি element-এ বাড়তি operation-এর খরচও নেই। এটাই দেখায় কেন “শুধু sum += x লিখলাম” আর “NumPy-র sum() ব্যবহার করলাম” — দুটো ভিন্ন array-তে ভিন্ন নির্ভুলতার ফলাফল দিতে পারে, এমনকি গাণিতিকভাবে একই কাজ করলেও।

Excel 2007-এর display বাগ। Microsoft-স্বীকৃত একটা বাগ — কিছু নির্দিষ্ট floating-point ফলাফল, যাদের প্রকৃত মান 65535-এর খুব কাছাকাছি ছিল, ভুলভাবে 100000 হিসেবে প্রদর্শিত হতো (প্রকৃত internal মান ঠিকই ছিল, শুধু binary-থেকে-decimal display conversion routine-এ একটা edge-case bug ছিল)। এটা arithmetic-এর ভুল ছিল না, কিন্তু floating-point representation-এর একটা নির্দিষ্ট bit pattern-এর জন্য display logic-এর একটা ভুল যা সরাসরি IEEE 754-এর জটিলতা থেকেই এসেছিল।

Spreadsheet-এ chained percentage rounding। একটা সাধারণ, এখনো প্রতিদিন ঘটে চলা সমস্যা — কর, discount, বা সুদ পরপর কয়েকবার হিসাব করলে (প্রতিবার একটা intermediate round করে), চূড়ান্ত ফলাফল “সঠিক” hand-calculation-এর সাথে কয়েক পয়সা মেলে না। এই কারণেই serious accounting software intermediate round এড়িয়ে শেষ ধাপে একবারই round করে — বা পুরোপুরি decimal/integer arithmetic ব্যবহার করে।

Python-এর decimal.Decimal — সঠিক সমাধান। যখন binary floating-point-এর 0.1-জাতীয় সমস্যা এড়ানো জরুরি (টাকার হিসাব, আইনি/নিয়ন্ত্রক compliance প্রয়োজন এমন হিসাব), Python-এর decimal module base-10 arithmetic দেয়:

from decimal import Decimal
Decimal('0.1') + Decimal('0.2') == Decimal('0.3')   # True!

Decimal('0.1') exactly 0.1 store করে, কারণ radix এখানে 10, 2 নয় — লেসন ৭-এর মূল সত্যের প্রত্যক্ষ প্রয়োগ (হর-এ 5 থাকলে base-2-এ সসীম নয়, কিন্তু base-10-এ সসীম)। খরচ: software-এ বাস্তবায়িত decimal arithmetic hardware-এর native binary float-এর চেয়ে উল্লেখযোগ্যভাবে ধীর।

Stripe, Shopify-জাতীয় payment API — integer cents। একটা আরও সহজ, দ্রুত বিকল্প — টাকা কখনো ভগ্নাংশ হিসেবে না রেখে smallest currency unit-এ (USD-তে cents) একটা সাধারণ integer হিসেবে রাখা। $19.99 মানে integer 1999। কোনো floating-point জড়িতই নয়, তাই কোনো rounding সমস্যাও নেই।

IEEE 754-2008-এর decimal floating-point। কম পরিচিত হলেও IEEE 754-2008 এবং পরবর্তী সংস্করণে decimal32/decimal64/decimal128 নামে base-10 floating-point format-ও সংজ্ঞায়িত আছে — কিছু IBM POWER hardware আর COBOL-ভিত্তিক mainframe financial system-এ native support আছে, যেখানে regulatory প্রয়োজনে decimal নির্ভুলতা বাধ্যতামূলক।

GPU/parallel reduction-এর non-determinism। আধুনিক GPU-তে বড় array-র sum (reduction) প্রায়ই বহু thread সমান্তরালে আংশিক যোগফল বের করে, শেষে একত্র করে। যেহেতু floating-point addition associative নয় ((a+b)+c \ne a+(b+c) সবসময় নয়), thread scheduling-এর ক্রম সামান্য বদলালে চূড়ান্ত sum-ও সামান্য বদলাতে পারে — একই কোড, একই input, তবু ভিন্ন run-এ সামান্য ভিন্ন ফলাফল। Machine learning training reproducibility নিয়ে গবেষণায় এটা একটা সক্রিয় বিষয়।

যে ভুলগুলো সবাই করে

“Patriot missile ব্যর্থতা ছিল একটা float == তুলনার bug।”

এটা এই ঘটনার সবচেয়ে সাধারণ ভুল বর্ণনা। প্রকৃত mechanism ছিল সম্পূর্ণ ভিন্ন: system clock-এর 0.1 সেকেন্ড টিককে সেকেন্ডে রূপান্তর করতে একটা ২৪-বিট fixed-point approximation ব্যবহার হতো (IEEE 754 float নয়), যেটা 0.1-এর অসীম বাইনারি expansion-কে truncate (chop) করত। এই ছোট্ট error প্রায় ১০০ ঘণ্টা একটানা চালু থাকা অবস্থায় জমতে জমতে 0.34 সেকেন্ড হয়ে যায় — কোনো তুলনা operation জড়িত ছিল না, শুধু একটা systematic, একদিকের truncation error দীর্ঘ সময় ধরে accumulate হয়েছিল। এটা এই লেসনের “naive summation”-এর ধারণার একটা বাস্তব, ভয়াবহ উদাহরণ, তুলনা-bug-এর নয়।

“Float তুলনার জন্য abs(a - b) \< একটা fixed ছোট epsilon (যেমন 1e-9) ব্যবহার করাই যথেষ্ট।”

এই pattern খুব জনপ্রিয়, কিন্তু magnitude-এর দুই প্রান্তেই ভেঙে পড়ে।

খুব বড় সংখ্যায় খুব কড়া: a, b \approx 10^{20}-এর কাছাকাছি হলে, সেই magnitude-এ ULP নিজেই 10^{20} \times 2^{-52} \approx 2.2 \times 10^4 — মানে দুই “বাস্তবিকভাবে সমান” (এক ULP দূরত্বের) সংখ্যার পার্থক্যও 1e-9-এর চেয়ে কয়েক অর্ডার অফ ম্যাগনিটিউড বড়। Fixed epsilon 1e-9 এখানে ভুলভাবে “সমান নয়” বলবে, সবসময়।

খুব ছোট সংখ্যায় খুব ঢিলা: a, b \approx 10^{-15}-এর কাছাকাছি হলে, 1e-9 epsilon সেই সংখ্যাগুলোর চেয়ে বহুগুণ বড় — সম্পূর্ণ ভিন্ন দুইটা ছোট সংখ্যাও abs(a-b) \lt 1e-9 শর্ত সহজেই মিটিয়ে “সমান” বলে ভুল রায় দেবে।

সঠিক pattern — relative + absolute tolerance একসাথে (Python-এর math.isclose-এর ঠিক এই সূত্র):

def is_close(a, b, rel_tol=1e-9, abs_tol=0.0):
    return abs(a - b) \<= max(rel_tol * max(abs(a), abs(b)), abs_tol)

rel_tol বড় সংখ্যায় scale করে (magnitude-এর সমানুপাতিক সীমা), আর abs_tol শূন্যের কাছাকাছি তুলনার জন্য একটা floor দেয় (যেখানে relative tolerance অর্থহীন হয়ে যায়, কারণ 0-এর কাছে সবকিছুর relative error বিশাল দেখায়)।

“0.1 + 0.2 ≠ 0.3 — এটা নিশ্চয়ই সেই নির্দিষ্ট ভাষা বা compiler-এর একটা bug।”

প্রতিটা ভাষা যা IEEE 754 double precision অনুসরণ করে (প্রায় সব আধুনিক ভাষা) — Python, JavaScript, Java, C, C++, Rust, Go, Swift — সবাই ঠিক একই ফলাফল দেবে, কারণ সমস্যাটা representation-এর, কোনো নির্দিষ্ট ভাষার implementation-এর নয়। 0.1-এর বাইনারি expansion অসীম (লেসন ৭), তাই কোনো সসীম বিট বাজেট — ৩২ হোক বা ৬৪ হোক, Python হোক বা C — সেটা exactly ধরতে পারবে না।

“Money-র হিসাবে float ব্যবহার করলেও সাবধানে round করলেই যথেষ্ট।”

এই approach fragile — এটা কাজ করে বলে মনে হয় ছোট পরীক্ষায়, কিন্তু scale-এ ভেঙে পড়ে। সমস্যাটা শুধু “round করতে ভুলে যাওয়া” নয় — এমনকি ঠিকঠাক round করলেও, প্রতিটা round নিজেই একটা ছোট error, আর লক্ষ লক্ষ লেনদেনে (এই লেসনের naive summation তত্ত্ব অনুযায়ী) সেই error জমতে পারে। আরও খারাপ — কোন সময়ে round করবেন (প্রতিটা intermediate ধাপে, নাকি শুধু শেষে) তার উপর নির্ভর করে ফলাফল বদলে যেতে পারে, যা audit/reconciliation-এ অসঙ্গতি তৈরি করে।

সঠিক নিয়ম, ব্যতিক্রমহীন: টাকার হিসাবে কখনো binary floating point ব্যবহার করবেন না। integer smallest-unit (cents) অথবা decimal.Decimal/সমতুল্য base-10 arithmetic ব্যবহার করুন — “সাবধানতা” কোনো স্থায়ী সমাধান নয়, সঠিক representation-ই একমাত্র স্থায়ী সমাধান।

“Kahan summation ব্যবহার করলেই summation error-এর সমস্যা চিরতরে শেষ।”

Kahan summation O(n\varepsilon)-কে O(\varepsilon)-এ নামায় — বিশাল উন্নতি, কিন্তু “চিরতরে শেষ” নয়, তিনটা কারণে।

১. এটা n-নির্ভরতা কমায়, শূন্য করে না। O(\varepsilon) মানে এখনো একটা non-zero error bound আছে, শুধু n বাড়লে সেটা বাড়ে না। খুব উচ্চ-নির্ভুলতার প্রয়োজনে (যেমন কিছু বৈজ্ঞানিক simulation) এটাও যথেষ্ট নাও হতে পারে।

২. Compiler optimization ভেঙে দিতে পারে। এই লেসনের “build” অংশে দেখানো হয়েছে — -ffast-math-এর মতো flag Kahan-এর মূল কৌশলকেই algebraically “সরল করে” মুছে দিতে পারে।

৩. চরম magnitude-ব্যবধানে এর নিজস্ব দুর্বলতা আছে — এই লেসনের Question ৬-এ দেখানো হয়েছে, যেখানে Neumaier summation প্রয়োজন হতে পারে।

সঠিক মানসিকতা: Kahan/Neumaier summation একটা উল্লেখযোগ্য উন্নতি, একটা গ্যারান্টি নয়। Critical numerical code-এ সবসময় একটা independent reference (যেমন decimal module, বা higher-precision computation) দিয়ে ফলাফল যাচাই করা ভালো অভ্যাস — এই লেসনের প্রতিটা Experiment ঠিক এই কাজটাই করেছে।

বুঝেছেন কি না দেখুন

লেসন ৮-এর bit-level জ্ঞান ব্যবহার করে ব্যাখ্যা করুন — 0.1 আর 0.2-এর exact stored মান যোগ করলে ফলাফল 0.3-এর exact stored মানের সাথে কেন মেলে না, শুধু “উভয়ই আনুমানিক” বলা ছাড়া নির্দিষ্টভাবে।

যুক্তি

তিনটা আলাদা সংখ্যাই স্বাধীনভাবে round হয়েছে, প্রতিটা তার নিজের নিকটতম representable double-এ:

0.1 round হয়ে store হয় 0.1000000000000000055511151231257827... (প্রকৃত 0.1-এর চেয়ে সামান্য বেশি)
0.2 (= 2 \times 0.1, exact দ্বিগুণকরণ) store হয় 0.2000000000000000111022302462515654... (আবার সামান্য বেশি)
এই দুইটা exact stored মান যোগ করলে গাণিতিকভাবে পাওয়া যায় 0.3000000000000000166533...-এর কাছাকাছি একটা মান, যেটা আবার নিকটতম double-এ round হয়ে দাঁড়ায় 0.3000000000000000444089209850062616... — প্রকৃত 0.3-এর চেয়ে সামান্য বেশি
কিন্তু 0.3 সরাসরি লিখলে (x = 0.3), সেই literal-টা স্বাধীনভাবে round হয় সবচেয়ে কাছের double-এ, যা হলো 0.2999999999999999888977697537484345... — প্রকৃত 0.3-এর চেয়ে সামান্য কম

সমস্যাটা হলো তিনটা ভিন্ন rounding ঘটনা, তিনটা ভিন্ন direction-এ (দুইটা উপরে, একটা নিচে) মিলে একটা যোগফল আর একটা সরাসরি-লেখা মান — এই দুইয়ের মধ্যে এক ULP পার্থক্য তৈরি করেছে (2^{-54} \approx 5.55 \times 10^{-17})। এটা কোনো একক বড় error নয় — এটা তিনটা independent, প্রতিটাই অত্যন্ত ছোট rounding সিদ্ধান্তের সমষ্টি যা কাকতালীয়ভাবে বিভিন্ন দিকে গেছে।

যদি round-off সব সংখ্যায় একই দিকে যেত (ধরুন সবসময় নিচে), তাহলে হয়তো যোগফল আর 0.3-এর সরাসরি representation একই মানে মিলে যেত (দুটোই সমান পরিমাণ নিচের দিকে সরে থাকত)। এখানে এমনটা হয়নি — exactly কেন সেটা bit-level হাতে-হিসাব ছাড়া আগে থেকে বলা কঠিন, যা-ই এই লেসনের মূল শিক্ষা: floating-point error-এর দিক আগে থেকে অনুমান করা যায় না, শুধু হিসাব করে বা experiment করে জানা যায়।

একটা প্রোগ্রামে আপনি একটা function-এর derivative আনুমানিক করছেন f'(x) \approx (f(x+h) - f(x))/h সূত্র দিয়ে। ছোট h (যেমন h = 10^{-15}, double precision-এ) ব্যবহার করলে ফলাফলের নির্ভুলতার কী হবে, আর কেন?

প্রয়োগ

স্বজ্ঞা বলে “যত ছোট h, তত ভালো approximation” — গণিতের limit-এর সংজ্ঞা অনুযায়ী এটা সত্য infinite precision-এ। কিন্তু floating-point-এ h খুব ছোট হলে catastrophic cancellation ঘটে, ঠিক এই লেসনের কেন্দ্রীয় প্যাটার্নে।

যখন h অত্যন্ত ছোট (10^{-15}-এর কাছাকাছি, double-এর machine epsilon \approx 2.2 \times 10^{-16}-এর কাছাকাছি), f(x+h) আর f(x) প্রায় সমান — এই দুইটা প্রায়-সমান সংখ্যা বিয়োগ করলে এই লেসনের “hood” অংশের নীতি অনুযায়ী: absolute error প্রায় অপরিবর্তিত (machine epsilon-এর মাত্রার), কিন্তু ফলাফল (f(x+h)-f(x)) নিজে খুব ছোট হয়ে যায় — relative error বিস্ফোরিত হয়। তারপর সেই ইতিমধ্যে-নষ্ট সংখ্যাকে আরও ছোট h দিয়ে ভাগ করলে error magnify হয় আরও বেশি।

ব্যবহারিক ফল: h কমানোর সাথে সাথে error প্রথমে কমে (truncation error কমছে — গাণিতিক approximation ভালো হচ্ছে), কিন্তু একটা বিন্দুর পর error আবার বাড়তে শুরু করে (cancellation error প্রাধান্য পাচ্ছে) — একটা U-আকৃতির error curve। সর্বোত্তম h সাধারণত \sqrt{\varepsilon} এর কাছাকাছি (double-এ প্রায় 10^{-8}), শূন্যের যত কাছে যাওয়া যায় তত ছোট নয়।

এটাই সেই একই নীতি যা quadratic formula উদাহরণে দেখা গেছে — “গাণিতিকভাবে ছোট h মানে বেশি নির্ভুল” স্বজ্ঞাটা floating-point-এ একটা নির্দিষ্ট সীমার পর ভুল হয়ে যায়।

Patriot missile সিস্টেমের chopping error প্রতি 0.1 সেকেন্ড tick-এ 9.5 \times 10^{-8} সেকেন্ড ছিল। যদি একটা battery ৭২ ঘণ্টা (৩ দিন) একটানা চালু থাকত (রিবুট ছাড়া), মোট accumulated timing error কত হতো, আর একটা ২,০০০ m/s গতির target-এ এটার position error কত দাঁড়াত?

প্রয়োগ

ধাপ ১ — মোট tick সংখ্যা:

$72 \text{ ঘণ্টা} = 72 \times 3600 = 259,200 \text{ সেকেন্ড}$

$259,200 \text{ সেকেন্ড} \div 0.1 \text{ সেকেন্ড/tick} = 2,592,000 \text{ টা tick}$

ধাপ ২ — মোট accumulated timing error:

$2,592,000 \times 9.5 \times 10^{-8} \approx 0.246 \text{ সেকেন্ড}$

ধাপ ৩ — position error:

$0.246 \text{ s} \times 2,000 \text{ m/s} = 492 \text{ মিটার}$

প্রায় ৪৯২ মিটার — মূল ঘটনার (~100 ঘণ্টা, ~573 মিটার) তুলনায় কম, কারণ কম সময় চালু ছিল, কিন্তু তবুও radar-এর range gate-এর সাধারণ প্রস্থের (কয়েক মিটার থেকে কয়েক দশ মিটার) তুলনায় বহুগুণ বড় — অর্থাৎ এমনকি ৭২ ঘণ্টার নিরবচ্ছিন্ন operation-ও সম্ভবত একই ধরনের ব্যর্থতার দিকে নিয়ে যেত।

এই হিসাবটাই দেখায় কেন সমাধানটা এত সহজ ছিল (আর কেন সেটা এতটা বিয়োগান্তক যে বাস্তবায়িত হয়নি সময়মতো) — শুধু system-কে নিয়মিত পুনরায় চালু করলেই (error counter শূন্যে রিসেট হয়ে) accumulated error কখনো বিপজ্জনক মাত্রায় পৌঁছাতে পারত না। যুদ্ধকালীন পরিস্থিতিতে এই নিয়মিত reboot ব্যবহারিকভাবে অবহেলিত হয়েছিল।

একটা e-commerce checkout system ডিজাইন করছেন, যেটা প্রতিটা item-এর দামের উপর tax হিসাব করে, তারপর সব item যোগ করে total বের করে। কেন এখানে binary float/double ব্যবহার করা বিপজ্জনক, আর সঠিক implementation কেমন হবে?

ডিজাইন

সমস্যাটা তিনটা স্তরে দেখা দিতে পারে:

১. একক মান representation-এ ভুল। $19.99 কে float-এ store করলে সেটা exactly 19.99 নয় (লেসন ৭-৮-এর মূল সত্য — 0.99-এর মতো ভগ্নাংশও বাইনারিতে সসীম নয়)। প্রতিটা মূল্যই একটা সামান্য আনুমানিক সংখ্যা।

২. Tax হিসাবে rounding। price \times tax\_rate একটা নতুন floating-point operation, নতুন rounding সহ। যদি প্রতিটা item-এ আলাদা করে tax round করা হয় (round(price \times rate, 2)), আর তারপর সেই rounded মানগুলো যোগ করা হয়, চূড়ান্ত total প্রকৃত “সঠিক” পদ্ধতি (সব যোগ করে একবারে round করা)-র চেয়ে কয়েক পয়সা আলাদা হতে পারে — regulatory audit-এ এই অসঙ্গতি সমস্যা তৈরি করতে পারে।

৩. Accumulated error — এই লেসনের naive summation তত্ত্ব। একটা bulk order-এ (হাজার হাজার item) সব দাম যোগ করলে, প্রতিটা addition-এর ছোট rounding জমে total-এ একটা লক্ষণীয় (যদিও সাধারণত ছোট) বিচ্যুতি তৈরি করতে পারে — বিশেষত যদি summation order প্রতিবার ভিন্ন হয় (যেমন distributed/parallel processing-এ), তখন একই cart, একই item, ভিন্ন run-এ ভিন্ন total আসতে পারে — একটা payment system-এর জন্য অগ্রহণযোগ্য।

সঠিক ডিজাইন:

সব monetary মান integer smallest-unit-এ রাখুন — cents (বা যে currency-র যা প্রযোজ্য)। $19.99 → integer 1999।
Tax rate-ও careful ভাবে হ্যান্ডেল করুন — সাধারণত একটা নির্দিষ্ট, দলিলভুক্ত rounding rule (round-half-up, বা round-half-even — যেটাই আইনত নির্ধারিত) ব্যবহার করে integer cents-এ ফলাফল বের করুন, প্রতিটা ধাপে।
সব addition integer arithmetic-এ — কোনো floating-point rounding-ই জড়িত নয়, তাই কোনো accumulated error নেই, কোনো order-নির্ভরতা নেই। Integer addition সবসময় exact আর deterministic।
শুধু display-এর সময় (UI-তে দেখানোর জন্য) integer cents-কে $19.99-স্টাইলে ফরম্যাট করুন — কোনো internal হিসাবে float জড়ান না।
বিকল্প (যদি fractional-cent precision দরকার হয়, যেমন কিছু crypto/forex system-এ): decimal.Decimal বা সমতুল্য arbitrary-precision base-10 type ব্যবহার করুন, কিন্তু কর্মক্ষমতার খরচ মাথায় রেখে।

এই পুরো ডিজাইনটাই এই লেসনের একটা বাক্যে সংক্ষিপ্ত: টাকার হিসাবে কখনো binary floating point-কে “সত্যের উৎস” (source of truth) বানাবেন না — শুধু display-এর জন্য ব্যবহার করুন, হিসাবের জন্য নয়।

একটা scientific simulation a = 1,000,000,000.0 (১০⁹) আর b = 1,000,000,000.0000001 (প্রায় একই, 10^{-7} পার্থক্য) — double precision-এ এই দুইটাকে abs(a-b) \lt 1e-9 দিয়ে তুলনা করলে কী ফলাফল আসবে, আর এটা কি “সঠিক” আচরণ?

যুক্তি

প্রথমে magnitude আর ULP হিসাব করি। 10^9 আছে 2^{29} \approx 5.37 \times 10^8 আর 2^{30} \approx 1.07 \times 10^9-এর মাঝে, তাই binade exponent e = 29। Double-এ ULP:

$\text{ULP} = 2^{(29-52)} = 2^{-23} \approx 1.19 \times 10^{-7}$

a আর b-এর প্রকৃত পার্থক্য 10^{-7} — এটা এই magnitude-এ ULP-এর (\approx 1.19 \times 10^{-7}) কাছাকাছি, তার মানে a আর b সম্ভবত একই representable double-এ round হয়ে গেছে অথবা মাত্র এক-দুই ULP দূরে আছে — তারা “একই রকম নির্ভুল” এই magnitude-এ যতটা সম্ভব।

এখন 1e-9 epsilon দিয়ে তুলনা করলে:

1e-9 এই magnitude-এর ULP-এর (1.19 \times 10^{-7}) চেয়ে প্রায় ১২০ গুণ ছোট। তার মানে abs(a-b) \lt 1e-9 প্রায় নিশ্চিতভাবে False দেবে — যদিও a আর b floating-point-এ “যতটা কাছে থাকা সম্ভব” ততটাই কাছে (হয়তো ০ বা ১ ULP দূরত্বে)।

এটা কি সঠিক আচরণ? এটা নির্ভর করে প্রশ্নটা কী — যদি “এরা কি bit-identical” জানতে চান, তাহলে এই তুলনা কাকতালীয়ভাবে প্রায় সঠিক উত্তর দিতে পারে। কিন্তু বেশিরভাগ ব্যবহারিক ক্ষেত্রে (convergence check, “essentially equal after floating-point noise” পরীক্ষা) উদ্দেশ্য হলো “এরা কি numerically সমতুল্য, floating-point-এর অনিবার্য imprecision বাদ দিয়ে” — আর সেই উদ্দেশ্যে 1e-9 ভুল উত্তর দেবে, কারণ এটা এই magnitude-এর জন্য অবাস্তবভাবে কড়া একটা মানদণ্ড প্রয়োগ করছে।

সঠিক পদ্ধতি — relative tolerance:

math.isclose(a, b, rel_tol=1e-9)

rel_tol=1e-9 মানে “একে অপরের 10^{-9} আপেক্ষিক ভগ্নাংশের মধ্যে” — 10^9 magnitude-এ এটা প্রায় 1.0 absolute tolerance দেয় (10^9 \times 10^{-9} = 1.0), যা 10^{-7} পার্থক্যকে সহজেই “close” হিসেবে সঠিকভাবে চিনবে।

সাধারণ নীতি: যখনই তুলনা করা সংখ্যার magnitude 1.0 থেকে বহুদূরে (বড় বা ছোট) হতে পারে, একটা fixed absolute epsilon কখনোই নির্ভরযোগ্য নয় — magnitude-স্কেলড (relative) tolerance, প্রয়োজনে একটা ছোট absolute floor সহ, ব্যবহার করতে হবে।

Kahan summation-এর y = x - c লাইনটা implicitly ধরে নেয় যে sum-এর magnitude সবসময় x-এর চেয়ে বড় বা কাছাকাছি। এই assumption কখন বাস্তবে ভেঙে পড়তে পারে, আর Neumaier summation কীভাবে সেটা ঠিক করে?

যুক্তি

Kahan-এর মূল সূত্র c = (t - sum) - y তখনই সবচেয়ে নির্ভুলভাবে কাজ করে যখন t - sum (গাণিতিকভাবে) y-এর কাছাকাছি নির্ভুলতায় representable হয় — যেটা সাধারণত ঘটে যখন |sum| \ge |x|, কারণ তখন t মূলত sum-এর magnitude-এর কাছাকাছি থাকে, আর t-sum সেই magnitude-এ যথেষ্ট precision নিয়ে y-এর হারানো অংশ ধরতে পারে।

সমস্যাটা কখন দেখা দেয়: ধরুন loop-এর শুরুতে sum এখনো ছোট (যেমন sum = 0.001), কিন্তু নতুন x অনেক বড় (x = 1000)। তখন |x| \gg |sum| — এখানে t = sum + x মূলত x-এর কাছাকাছি একটা মান, আর t - sum গাণিতিকভাবে প্রায় x-ই, কিন্তু এই subtraction-এ sum-এর ছোট মানটাই rounding-এ হারিয়ে যাওয়ার ঝুঁকিতে পড়ে (ঠিক এই লেসনের catastrophic cancellation নীতির একটা সূক্ষ্ম রূপ) — Kahan-এর সূত্র তখন sum-এর হারানো অংশ সঠিকভাবে ধরতে পারে না।

Neumaier-এর সমাধান — প্রতি ধাপে |sum| বনাম |x| পরীক্ষা করে, যেটা ছোট তার হারানো অংশ আলাদাভাবে হিসাব করে:

যদি |sum| ≥ |x|:  c += (sum - t) + x    # x-এর হারানো অংশ উদ্ধার
নাহলে:             c += (x - t) + sum    # sum-এর হারানো অংশ উদ্ধার

এই শর্তসাপেক্ষ যুক্তিটাই এই লেসনের “hood” অংশে দেখানো Neumaier কোডের মূল যোগ — এটা নিশ্চিত করে যেই ছোট অপারেন্ডের তথ্যই হারিয়ে যাক না কেন, সেটা সঠিকভাবে c-তে ধরা পড়বে, শুধু “sum সবসময় বড়” এই একমুখী অনুমানের উপর নির্ভর না করে।

ব্যবহারিক তাৎপর্য: যদি আপনার input sequence-এ magnitude বারবার ওঠানামা করে (কখনো ছোট সংখ্যার পর বড়, কখনো উল্টো — যা বাস্তব sensor data বা মিশ্র-scale financial data-তে সাধারণ), Neumaier summation Kahan-এর চেয়ে systematically বেশি নির্ভরযোগ্য। যদি input সবসময় একই দিকে বাড়তে থাকে (monotonically), দুটোর মধ্যে পার্থক্য practically নগণ্য।

এরপর কী

এরপর কী — সংখ্যা থেকে অক্ষরে

লেসন ৭-এ যে প্রশ্ন দিয়ে শুরু হয়েছিল — “ভগ্নাংশ কীভাবে bit-এ লেখা যায়?” — তার একটা সম্পূর্ণ, বাস্তবসম্মত উত্তর এখন আপনার কাছে আছে, শুধু ধারণাগত নয়, ব্যবহারিকও।

তিনটা লেসনে আমরা একটা সম্পূর্ণ পথ পাড়ি দিয়েছি — কেন integer যথেষ্ট নয় (fixed-point), scientific notation কীভাবে বাইনারিতে রূপান্তরিত হয় (IEEE 754-এর bit layout), আর সেই representation-এর অনিবার্য imprecision কীভাবে arithmetic-এর মধ্য দিয়ে বাতিল হয়, জমে, বা বিপর্যয়কর হয়ে ওঠে (এই লেসন)।

মূল সুতোটা মনে করুন — এই পুরো module-এর থিসিস: একটা bit pattern নিজে কিছু বলে না, interpretation বলে। আমরা দেখেছি সংখ্যার জন্য এই থিসিসটা কতটা গভীর — একই bit pattern integer, fixed-point, বা IEEE 754 float হিসেবে সম্পূর্ণ ভিন্ন মান বহন করতে পারে, আর floating-point-এর ক্ষেত্রে সেই “মান”-টাও সবসময় একটা compromise, একটা approximation।

এই module-এর পরের ধাপ — টেক্সট। 65 সংখ্যাটা 'A' অক্ষর হয় কীভাবে? পৃথিবীর হাজারো ভাষার (বাংলা-সহ) হাজারো অক্ষর মাত্র কয়েকটা byte-এ কীভাবে ধরা হয়? ASCII কেন যথেষ্ট ছিল না, আর Unicode/UTF-8 কীভাবে সেই সমস্যা সমাধান করে — এটাই সেই একই “bit pattern-এর অর্থ কে ঠিক করে” প্রশ্নের আরেকটা অধ্যায়, এবার সংখ্যার বদলে ভাষার জগতে।

একটা শেষ পর্যবেক্ষণ, যা এই তিনটা লেসনের সবচেয়ে গুরুত্বপূর্ণ ব্যবহারিক পাঠ হয়ে থাকুক: floating-point কোনো “খারাপ” প্রযুক্তি না — এটা একটা অসাধারণ সমাধান একটা মৌলিকভাবে কঠিন সমস্যার (সসীম বিটে অসীম সংখ্যা ধরা)। সমস্যাটা হয় তখনই যখন আমরা ভুলে যাই এটা একটা approximation, আর তার সাথে সেই ভুলে-যাওয়া অনুযায়ী কোড লিখি — == দিয়ে তুলনা করি, বিয়োগের ক্রম নিয়ে না ভেবেই সূত্র লিখি, বা টাকার হিসাবে সরাসরি float ব্যবহার করি। এই তিনটা লেসনের প্রতিটা misconception, প্রতিটা historical বিপর্যয় — সবগুলোরই মূলে এই একটাই ভুলে-যাওয়া।

Floating-Point Precision — যখন ছোট্ট Error জমে বিপর্যয় হয়

এই লেসন শেষে আপনি পারবেন

আগে যা বোঝা থাকা দরকার

আগে এটা বুঝি

মূল ধারণা

`0.1 + 0.2 ≠ 0.3` — সম্পূর্ণ প্রমাণ

Error-এর ভাষা — absolute বনাম relative

Machine epsilon

ভেতরে কী ঘটছে

যেখানে error জমে, বাতিল হয়, বা বিস্ফোরিত হয়

১. Catastrophic cancellation — তত্ত্ব

২. Naive summation — জমতে থাকা error

৩. Kahan summation — compensation দিয়ে error পুনরুদ্ধার

৪. Patriot missile (১৯৯১) — জমতে থাকা error-এর একটা প্রকৃত বিপর্যয়

উদাহরণ

সম্পূর্ণ hand-worked উদাহরণ — Catastrophic Cancellation

নিজে চালিয়ে দেখুন

0.1 + 0.2 বিট-বাই-বিট verify করা

Catastrophic cancellation নিজে দেখুন — quadratic formula

নিজে বানান

Kahan Summation — accumulated error পুনরুদ্ধার করুন

বাস্তব সিস্টেমে

যেখানে এই precision সমস্যা সত্যিই ঘটেছে

যে ভুলগুলো সবাই করে

বুঝেছেন কি না দেখুন

এরপর কী

এরপর কী — সংখ্যা থেকে অক্ষরে

আরও পড়ুন

এই লেসন শেষে আপনি পারবেন

আগে যা বোঝা থাকা দরকার

◔আগে এটা বুঝিWhy this matters

◆মূল ধারণাCore concept

0.1 + 0.2 ≠ 0.3 — সম্পূর্ণ প্রমাণ

Error-এর ভাষা — absolute বনাম relative

Machine epsilon

⚙ভেতরে কী ঘটছেUnder the hood

যেখানে error জমে, বাতিল হয়, বা বিস্ফোরিত হয়

১. Catastrophic cancellation — তত্ত্ব

২. Naive summation — জমতে থাকা error

৩. Kahan summation — compensation দিয়ে error পুনরুদ্ধার

৪. Patriot missile (১৯৯১) — জমতে থাকা error-এর একটা প্রকৃত বিপর্যয়

▶উদাহরণExample

সম্পূর্ণ hand-worked উদাহরণ — Catastrophic Cancellation

⌘নিজে চালিয়ে দেখুনExperiment

0.1 + 0.2 বিট-বাই-বিট verify করা

Catastrophic cancellation নিজে দেখুন — quadratic formula

⚒নিজে বানানBuild it yourself

Kahan Summation — accumulated error পুনরুদ্ধার করুন

◉বাস্তব সিস্টেমেReal-world usage

যেখানে এই precision সমস্যা সত্যিই ঘটেছে

⚠যে ভুলগুলো সবাই করেMisconceptions

?বুঝেছেন কি না দেখুনCheck understanding

→এরপর কীNext step

এরপর কী — সংখ্যা থেকে অক্ষরে

আরও পড়ুন

আগে এটা বুঝি

মূল ধারণা

`0.1 + 0.2 ≠ 0.3` — সম্পূর্ণ প্রমাণ

ভেতরে কী ঘটছে

উদাহরণ

নিজে চালিয়ে দেখুন

নিজে বানান

বাস্তব সিস্টেমে

যে ভুলগুলো সবাই করে

বুঝেছেন কি না দেখুন

এরপর কী