Worked examples — Floating point gotchas — catastrophic cancellation, associativity failure

5.4.22 · D3 · Coding › Scientific Computing (Python) › Floating point gotchas — catastrophic cancellation, associat

Kuch bhi shuru karne se pehle: teen words jo hum baar baar use karenge.

---

The scenario matrix

Is topic par har trouble in cells mein se ek mein aati hai. Har row neeche ek worked example hai.

#	Case class	Trigger	What breaks	Fix family
A	Near-equal positives subtract karo	$a\approx b$, dono $>0$	catastrophic cancellation	algebraic reformulation
B	Quadratic, same-sign root	$b^2\gg 4ac$, safe sign lo	theek hai (koi fix nahi chahiye)	—
C	Quadratic, opposite-sign root	doosra root	cancellation in $-b+\sqrt{\cdot }$	$x=c/(a{x}_{\text{safe}})$
D	Degenerate: term rounds to zero	$1-\cos x$, $x\to 0$	total loss (answer = 0)	trig identity
E	Absorption — small + huge	$1+10^{16}$	chota wala gayab ho jaata hai	reorder / Kahan
F	Associativity flip	$(a+b)+c\ne a+(b+c)$	grouping se answer badal jaata hai	pehle large terms cancel karo
G	Summation order over a list	many tiny + one huge	error $n$ ke saath badhti hai	sort ascending / Kahan
H	Equality comparison	0.1+0.2 == 0.3	unexpectedly False	tolerance compare
I	Real-world word problem	GPS / distance difference	cancellation disguise mein chhupta hai	physics reformulate karo
J	Exam twist	"higher precision fixes it?"	conditioning ≠ precision	reformulate karo, upgrade mat karo

---

Example A — near-equal positives subtract karo

Steps.

1. True value. $\sqrt{x+1}-\sqrt{x}\approx \frac{1}{2\sqrt{x}}=\frac{1}{2\times 10^8}=5\times 10^{-9}$. Ye step kyun? Hume honest target chahiye taaki hum uske against error measure kar sakein. Approximation $\sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}}\approx \frac{1}{2\sqrt{x}}$ se aati hai.

2. Naive compute. Doubles mein, $\sqrt{10^{16}}=10^8$ exactly hai, aur $\sqrt{10^{16}+1}$ — lekin $10^{16}+1$ rounds to $10^{16}$ (gap near $10^{16}$ $2$ hai, isliye $1$ add karna edge se gir jaata hai, cell E in disguise). Toh dono square roots $10^8$ ke barabar hain, aur subtraction $\hat{f}=0$ deti hai. Ye step kyun? Double failure dikhata hai: pehle absorption, phir do identical numbers subtract karo toh exactly $0$ milta hai.

3. Damage measure karo. Relative error $=\left|\frac{0-5\times 10^{-9}}{5\times 10^{-9}}\right|=1$ — 100% galat. Ye step kyun? Confirm karta hai "catastrophic": answer completely lost ho gaya.

4. Reformulate. Top aur bottom ko conjugate se multiply karo: $\sqrt{x+1}-\sqrt{x}=\frac{(\sqrt{x+1}-\sqrt{x})(\sqrt{x+1}+\sqrt{x})}{\sqrt{x+1}+\sqrt{x}}=\frac{1}{\sqrt{x+1}+\sqrt{x}}.$ Ab $\frac{1}{\sqrt{10^{16}+1}+\sqrt{10^{16}}}=\frac{1}{2\times 10^8}=5\times 10^{-9}$ compute karo. Near-equals ki koi subtraction nahi — hum sirf same-sign roots add karte hain aur divide karte hain. Ye step kyun? Division aur same-sign addition kabhi relative error amplify nahi karte (parent note, §2), toh hum villain ke around route le lete hain.

Verify: reformulated value $5\times 10^{-9}$ true $\frac{1}{2\sqrt{x}}$ se full precision tak match karta hai. Sanity: large $x$ ke liye consecutive roots ke beech gap $1/\sqrt{x}$ ki tarah shrink hona chahiye — karta hai. ✓

---

Examples B & C — do quadratic roots

Hum ek quadratic lete hain aur dono roots nikalte hain, kyunki dono roots alag matrix cells mein rehte hain: ek safe hai (B), ek cancel hota hai (C).

Steps.

1. Discriminant. $b^2-4ac=4\times 10^{10}-4=39999999996$, toh $\sqrt{\cdot }\approx 199999.99999$. Ye step kyun? Dono roots ye square root share karte hain; ek baar compute karo.

2. Same-sign root (Cell B — safe). Kyunki $b<0$ hai, term $-b=+200000$ positive hai; safe root do positives add karta hai: $x_1=\frac{-b+\sqrt{b^2-4ac}}{2a}=\frac{200000+199999.99999}{2}\approx \mathrm{199999.999995.}$ Ye step kyun? Do same-sign near-equal numbers add karna cancel nahi karta — koi amplification nahi. Ye root trustworthy hai.

3. Opposite-sign root, naive (Cell C — catastrophic). $x_2^{\text{naive}}=\frac{200000-199999.99999}{2}=\frac{0.00001}{2}\approx 5\times 10^{-6}.$ Yahan cancellation do operands $p=200000$ aur $q=199999.99999$ ke beech hoti hai. Parent note ka amplifier, generic subtraction $p-q$ ke liye likha, $\frac{|p|+|q|}{|p-q|}$ hai (parent ke §2 mein letters $a,b$ sirf placeholder operands hain — ye quadratic coefficients nahi hain). Plug in karo: $\frac{|p|+|q|}{|p-q|}\approx \frac{400000}{10^{-5}}=4\times 10^{10}$ — precision annihilated. Ye step kyun? Cancellation aur uski size expose karta hai, operands $p,q$ ko coefficients $a,b,c$ se confuse nahi karte.

4. Stable second root. $x_1{x}_2=c/a$ (product of roots) use karo, toh $x_2=\frac{c}{a{x}_1}=\frac{1}{199999.999995}\approx 5.00000000013\times 10^{-6}.$ Ye step kyun? Hamare paas pehle se safe root $x_1$ hai; product identity $x_2$ sirf division use karke deta hai — near-equals ki koi subtraction nahi. Dekho Quadratic formula numerical issues.

Verify: $x_2$ back plug karo: $x_2^2-200000x_2+1\approx 0$. Product check: $x_1{x}_2=c/a=1$. Sum check: $x_1+x_2=-b/a=200000$. ✓

---

Example D — ek term zero round ho jaata hai (degenerate)

Steps.

1. True value. Taylor: $\cos x=1-\frac{x^2}{2}+\dots$, toh $1-\cos x=\frac{x^2}{2}$ aur $f\to \frac{1}{2}$. Ye step kyun? Target $0.5$ establish karta hai.

2. Machine kya store karta hai. $\cos (10^{-8})=1-5\times 10^{-17}+\dots$. Lekin $5\times 10^{-17}<u\approx 1.11\times 10^{-16}$, toh $\cos (10^{-8})$ exactly $1.0$ round ho jaata hai. Ye step kyun? Woh information jo $f$ ko nonzero banati thi machine epsilon ke neeche rehti hai — subtraction se pehle hi gone ho jaati hai.

3. Naive result. $1-1=0$, toh $\hat{f}=0/x^2=0$. Relative error $=1$ (total loss). Ye step kyun? Degeneracy confirm karta hai: answer zero collapse ho jaata hai.

4. Trig identity se reformulate karo. $1-\cos x=2{\sin }^2(x/2)$, hence $f(x)=\frac{2{\sin }^2(x/2)}{x^2}=\frac{1}{2}{\left(\frac{\sin (x/2)}{x/2}\right)}^2.$ Chhote argument ka $\sin$ accurately compute hota hai (woh $\approx$ argument hai, koi cancellation nahi), aur $\frac{\sin t}{t}\to 1$, toh $f\to \frac{1}{2}$. Ye step kyun? Humne "subtract near-equals" ko "$\sin$ of a small number" se replace kiya, jise library faithfully evaluate karti hai.

Verify: $x=10^{-8}$ par, $\sin (5\times 10^{-9})/5\times 10^{-9}\approx 1$, toh $f\approx \frac{1}{2}(1{)}^2=0.5$. Limit se match karta hai. ✓

Figure neeche dono naive aur reformulated versions ko kaafi $x$ values par plot karta hai. Aise padho: dashed gray line true limit $0.5$ hai; green curve (reformulated) seedha usse hug karta rehta hai neeche tak; red curve (naive) $0.5$ pe ride karta hai jab tak $x\approx 10^{-8}$, phir cliff se gir ke $0$ par aa jaata hai jis waqt $\cos x$ $1$ round ho jaata hai. Woh cliff hi catastrophic cancellation hai — picture exactly woh $x$ dikhati hai jahan naive formula marr jaata hai.

---

Example E — absorption (small + huge)

Steps.

1. $10^{16}$ par gap. $2^{53}\approx 9\times 10^{15}$ magnitude ke paas, consecutive doubles $2$ se differ karte hain (mantissa ka last bit ab $1$ nahi, $2$ weight karta hai). Toh $10^{16}$ ke neighbours $2$ apart hain. Ye step kyun? Us magnitude par resolution jaanna zaroori hai — precision relative hoti hai, isliye bade numbers coarse spacing rakhte hain (dekho IEEE-754 floating point representation).

2. $1$ add karo. $10^{16}+1=10000000000000001$ (odd) — lekin woh representable nahi; nearest doubles $10^{16}$ aur $10^{16}+2$ hain. Rounding-to-nearest $10^{16}$ rakhta hai. Ye step kyun? $+1$ gap ke aadhe se chhhota hai, isliye absorbed ho jaata hai — completely vanish ho jaata hai.

3. Consequence. $\widehat{(1+10^{16})}=10^{16}$. $1$ ne koi trace nahi chhoda. Ye step kyun? Ye Example F mein associativity break ka mechanism hai.

Verify: Python mein, 1e16 + 1 == 1e16 is True; 1e16 + 2 == 1e16 is False. Gap exactly $2$ hai. ✓

---

Example F — associativity flip

Example se pehle, ek chhhota sa lemma jis par hum lean karenge.

Steps.

1. Left grouping $(a+b)+c$. $a+b=1+10^{16}\to 10^{16}$ (absorption, Example E). Phir $+c=10^{16}-10^{16}=0$. Ye step kyun? Dikhata hai ki pehle add mein $1$ mar jaata hai, isliye result $0$ hai.

2. Right grouping $a+(b+c)$. $b+c=10^{16}-10^{16}=0$ exactly — ye Sterbenz lemma upar hai: $b$ aur $-c$ magnitude mein equal hain, isliye subtraction koi rounding introduce nahi karta. Phir $+a=0+1=1$. Ye step kyun? Yahan huge numbers pehle cancel ho jaate hain (exactly), tiny $1$ ko protect karte hain.

3. Compare. Left $=0$, right $=1$. $100\%$ differ karte hain. Ye step kyun? Demonstrate karta hai ki floating addition associative nahi hai — grouping decide karta hai kaun se bits survive karte hain.

4. Rule of thumb. Large opposite-sign terms ko pehle cancel karo, chote ones add karne se pehle. Ye step kyun? Chote term ko absorb hone se bachata hai.

Verify: (1.0 + 1e16) + (-1e16) == 0.0 aur 1.0 + (1e16 + (-1e16)) == 1.0. ✓

Bar chart neeche flip ko vivid banata hai. Aise padho: left red bar $(a+b)+c$ grouping hai — woh $0$ tak collapse ho jaata hai kyunki $1$ ko $10^{16}$ ne absorb kar liya $c$ ke cancel karne se pehle; right green bar $a+(b+c)$ hai — woh $1$ tak pahunchta hai kyunki do huge terms pehle annihilate ho jaate hain (Sterbenz-exact) aur $1$ survive karta hai. Numbers ke do identical sets, do alag heights: woh height gap hi associativity failure hai.

---

Example G — list par summation order

Steps.

1. Naive left-to-right. Start $s=10^{16}$. Pehla $1$ add karo: $10^{16}+1\to 10^{16}$ (absorbed). Har subsequent $+1$ similarly absorbed ho jaata hai. Final $\hat{s}=10^{16}$ — saare $10^6$ ones lost. Ye step kyun? Running partial sum pehle step se hi huge hai, isliye har tiny addend uske gap ke neeche gir jaata hai.

2. Ascending order. Sort karo taaki huge number last aaye: pehle million ones sum karo, $10^6$ milega (har add comparable numbers ke beech hai — safe), phir $10^{16}$ add karo. Ab $10^{16}+10^6$: gap at $10^{16}$ $2$ hai, aur $10^6\gg 2$, isliye survive karta hai: $\hat{s}=10000000001000000$. Ye step kyun? Pehle smallest add karne se woh itne bade number mein accumulate ho jaate hain ki absorb na ho sakein. Dekho Round-off error propagation.

3. Kahan option. Kahan compensated summation lost low bits ko ek compensator c mein carry karta hai, jo error $\sim u$ deta hai order aur $n$ se independent. Ye step kyun? Jab sort nahi kar sakte (streaming data), Kahan waise bhi accuracy recover karta hai.

Verify: naive $10^{16}$ deta hai; ascending $10^{16}+10^6$ deta hai; dono true $10000000001000000$ se compare karo. ✓

---

Example H — equality trap

Steps.

1. Kyun woh exact nahi hain. $0.1,0.2,0.3$ finite binary fractions nahi hain (jaise decimal mein $1/3$). Har ek nearest double ke roop mein store hota hai, pehle se $\sim u$ off. Ye step kyun? Premise "math mein equal ⇒ floats mein equal" galat hai kyunki inputs khud rounded hain.

2. Sum. $\text{fl}(0.1)+\text{fl}(0.2)=0.30000000000000004\dots$, jabki $\text{fl}(0.3)=0.29999999999999998\dots$. Dono $\sim 5.6\times 10^{-17}$ se differ karte hain. Ye step kyun? Dikhata hai ki dono sides alag doubles par land karte hain, isliye == False return karta hai.

3. Fix: tolerance compare. math.isclose(x, y) use karo, yaani $|x-y|\le \text{atol}+\text{rtol}\cdot \max (|x|,|y|)$ default rtol=1e-9 ke saath. Ye step kyun? Hum poochhte hain "kya woh round-off ke andar equal hain?" bajaaye "bit-identical?".

Verify: 0.1 + 0.2 == 0.3 is False; difference lagbhag $5.55\times 10^{-17}$ hai; abs((0.1+0.2)-0.3) < 1e-9 is True. ✓

---

Example I — real-world word problem

Steps.

1. True difference. $\Delta =6378137.250-6378137.000=0.250\text{ m}$. Ye step kyun? Target set karta hai: ek 25 cm difference.

2. Digit accounting. Har altitude $\sim 16$ significant digits (double) ke liye good hai. Har ek ka absolute error lagbhag $|R|\times u\approx 6.378\times 10^6\times 1.11\times 10^{-16}\approx 7\times 10^{-10}\text{ m}$ hai. Subtraction exact hai (Sterbenz — do values factor of 2 ke andar hain), toh woh do absolute errors seedha $\Delta$ mein carry ho jaate hain: worst case $\sim 1.4\times 10^{-9}\text{ m}$. $\Delta$ ka relative error tab $\frac{1.4\times 10^{-9}}{0.25}\approx 5.6\times 10^{-9}$ hai. Ye step kyun? Relative-error amplifier $\frac{|R_1|+|R_2|}{|\Delta |}\approx \frac{1.2756\times 10^7}{0.25}\approx 5.1\times 10^7$ batata hai ki hum lagbhag ${\log }_{10}(5.1\times 10^7)\approx 7.7$ digits lose karte hain. $\sim 16$ good digits se start karke, roughly $8$ bachte hain — toh $0.25\text{ m}$ yahan bhi theek hai. Lekin agar do altitudes 25 nanometres apart hote, amplifier $10^{15}$ hota aur kuch bhi survive na karta. Ye kyun matter karta hai: danger scale karta hai ki do bade numbers kitne close hain.

3. Engineering fix. Difference ko ek common reference se store karo: har reading se ek base altitude (maano $6378137\text{ m}$) rounding se pehle subtract karo, ya ek local coordinate frame mein kaam karo jahan chhhoti quantity seedha represent ki jaati hai na ki do bade numbers ke difference ke roop mein. Tab tum 25 cm quantity par full precision rakhte ho, kyunki woh kabhi big-minus-big se paida nahi hoti. Ye step kyun? Interesting chhhoti number ko kabhi do nearly-equal large numbers ke difference ke roop mein paida mat hone do — seedha represent karo. Dekho Numerical stability and conditioning.

Verify: $\Delta =0.25$; amplifier $=\frac{R_1+R_2}{\Delta }\approx 5.1\times 10^7$; digits lost $={\log }_{10}(\text{amplifier})\approx 7.7$. ✓

---

Example J — exam twist

Steps.

1. float128 kya deta hai. `np.float128` (x86 80-bit extended) mein $\sim 18$–$19$ significant digits hain, $u\approx 5.4\times 10^{-20}$. Ye step kyun? Claim karne se pehle naya resolution jaano ki woh help karta hai ya nahi.

2. Kya cancellation vanish ho jaati hai? $\cos (10^{-8})=1-5\times 10^{-17}$; ab $5\times 10^{-17}>u$, toh term is particular $x$ par survive karti hai, aur tumhe $\approx 0.5$ milta hai. Lagta hai fixed ho gaya. Ye step kyun? Dikhata hai ki higher precision sirf cliff ko postpone karta hai.

3. $x$ chhhota karo. $x=10^{-10}$ par, $1-\cos x\approx 5\times 10^{-21}<u_{128}$, toh float128 bhi $0$ round kar deta hai. Catastrophe wapas aata hai — tumne sirf threshold move ki. Ye step kyun? Prove karta hai ki problem formula ki conditioning hai, machine ki precision nahi.

4. Correct fix. $\frac{1}{2}{\left(\frac{\sin (x/2)}{x/2}\right)}^2$ ke roop mein reformulate karo (Example D). Ye plain float64 mein sabhi small $x$ ke liye kaam karta hai. Ye step kyun? Expression fix karna hardware upgrade karne se better hai — parent note ka mistake "more precision always fixes accuracy", cell J.

Verify: $x=10^{-8}$ par, float64 naive $=0$ lekin reformulated $\approx 0.5$; reformulation ko koi extra precision nahi chahiye. ✓

---

Recall

Recall Answers cover karo — kaun sa fix kaun se cell ke liye?

Kaun sa cell hai "do bade same-sign numbers subtract karo"? ::: A / C — catastrophic cancellation, fix by reformulation. Safe quadratic root kaun sa sign of $b$ use karta hai? ::: Woh sign jo $-b$ ko add karta hai (same sign as $-b\pm \sqrt{\cdot }$), toh koi cancellation nahi. Dangerous quadratic root safely kaise nikaalte hain? ::: $x_2=c/(a{x}_1)$ product of roots se. Kyun $1-\cos (10^{-8})$ float64 mein exactly 0 deta hai? ::: $5\times 10^{-17}$ correction $u\approx 1.1\times 10^{-16}$ ke neeche hai, toh $\cos$ $1$ round ho jaata hai. One-huge-many-tiny ke liye best summation order? ::: Ascending (smallest first) — ya Kahan. Kya float128 ek cancelling formula fix karta hai? ::: Nahi, ye sirf threshold move karta hai; bajaaye reformulate karo. Sterbenz lemma kya guarantee karta hai? ::: Agar do same-sign floats factor of 2 ke andar hain, unki subtraction exactly compute hoti hai (koi rounding nahi).