Floating point-এর হিসাব — যোগ, গুণ, rounding আর ভুলের ফাঁদ

১ · একটু ঝালিয়ে নেওয়া — float কীভাবে রাখা থাকে, আর হিসাব কেন কঠিন

পর্ব ২৩-এ (Floating point — IEEE 754) আপনি দেখেছেন floating point সংখ্যা কীভাবে বিটে রাখা হয়। এখানে শুধু সেটুকু মনে করিয়ে দিই, যা হিসাবের জন্য লাগবে। IEEE 754 মানদণ্ডে একটি normal সংখ্যার মান

এখানে s হলো চিহ্নের বিট, e হলো biased exponent (excess-K রূপ, পর্ব ০৫), আর 1.f হলো significand — সামনের ১-টি লুকানো (hidden bit), তাই লেখা থাকে শুধু ভগ্নাংশ f. ৩২ বিটের single precision (float)-এ ১ + ৮ + ২৩ বিট, bias 127, নির্ভুলতা p = ২৪ বিট। ৬৪ বিটের double-এ ১ + ১১ + ৫২ বিট, bias 1023, p = ৫৩ বিট। Subnormal, ∞, NaN, machine epsilon আর ULP-এর সংজ্ঞাও পর্ব ২৩-এ আছে; এখানে আমরা সেগুলো ব্যবহার করব।

Integer যোগে (পর্ব ২০) একটি adder-ই যথেষ্ট। Floating point যোগে কিন্তু সংখ্যা দুটি আলাদা “স্কেলে” থাকে: একটি হয়তো 2¹⁰-এর ঘরে, অন্যটি 2⁻³-এর ঘরে। তাই আগে স্কেল মেলানো, তারপর যোগ, তারপর ফলকে আবার নিয়মের রূপে আনা (normalization), তারপর বাড়তি বিট ফেলে দেওয়ার নিয়ম (rounding). এই পর্বের বেশিরভাগই এই ক্রমের ভেতরের খুঁটিনাটি।

২ · যোগ-বিয়োগের algorithm — exponent মেলানো থেকে সীমা যাচাই

Floating point যোগ (আর বিয়োগ, যা আসলে চিহ্ন উল্টে যোগ) পাঁচটি ধাপে হয়। প্রথমে দশমিকে, ৪ অঙ্কের নির্ভুলতায় — Patterson-Hennessy-র বইয়ের পরিচিত উদাহরণ: 9.999 × 10¹ + 1.610 × 10⁻¹.

1. Exponent মেলানো (alignment): ছোট exponent-এর সংখ্যাটিকে বড়টির exponent-এ আনা হয়। 1.610 × 10⁻¹ = 0.01610 × 10¹. দুই ঘর ডানে সরল। পর্দায় ৪ অঙ্ক রাখা গেলে থাকে 0.016 × 10¹ — শেষের “10” পড়ে গেল।
2. Significand যোগ: 9.999 + 0.016 = 10.015, অর্থাৎ 10.015 × 10¹.
3. Normalise: দশমিক বিন্দুর আগে একটিমাত্র অশূন্য অঙ্ক রাখতে হবে: 1.0015 × 10². Significand এক ঘর ডানে সরল, তাই exponent এক বাড়ল। বিয়োগে উল্টোটাও হয় — ফল ছোট হয়ে গেলে বাঁয়ে সরিয়ে exponent কমাতে হয়।
4. Round: ৪ অঙ্কে আনতে হবে। 1.0015-এর বাড়তি অঙ্ক “5”, ঠিক মাঝখান; nearest-even নিয়মে 1.001 বিজোড়, তাই ওপরে: 1.002 × 10². (আসল যোগফল 100.151; ৪ অঙ্কে সবচেয়ে কাছের মানও 1.002 × 10² — এবার ভাগ্য ভালো ছিল, কেন সবসময় থাকে না তা অংশ ৩-এ।)
5. সীমা যাচাই: exponent বেশি বড় হলে overflow (ফল ±∞), বেশি ছোট হলে underflow (subnormal বা ০)। Round-এর পরে আবার normalise লাগতে পারে (যেমন 9.9999 → 10.000), তাই যাচাই শেষে।

একই কাজ binary-তে: 0.5 + (−0.4375)

এবার বইয়ের binary উদাহরণ। 0.5 = 1.000₂ × 2⁻¹, আর 0.4375 = 7/16 = 0.0111₂ = 1.110₂ × 2⁻². ৪ বিট নির্ভুলতায় (1.xxx):

• Align: ছোট exponent −2, বড় −1; পার্থক্য ১. তাই −1.110 × 2⁻² = −0.111 × 2⁻¹.
• যোগ: 1.000 − 0.111 = 0.001, অর্থাৎ 0.001 × 2⁻¹.
• Normalise: সামনে তিনটি শূন্য, তাই significand তিন ঘর বাঁয়ে, exponent তিন কম: 1.000 × 2⁻⁴.
• সীমা আর round: single precision-এ normal exponent −126 থেকে 127; −4 ভেতরে। Significand ঠিক ৪ বিটে ধরে, round-এর দরকার নেই। ফল 1.000 × 2⁻⁴ = 1/16 = 0.0625 — আর 0.5 − 0.4375 সত্যিই 0.0625.

এখানে একটি নতুন জিনিস লক্ষ করুন: বিয়োগে সামনের অনেকগুলো বিট বাতিল হয়ে যেতে পারে (এখানে তিনটি), তখন বড় করে বাঁয়ে shift লাগে। কতগুলো শূন্য আছে গোনার জন্য hardware-এ একটি leading-zero counter থাকে। এই “বাতিল হয়ে যাওয়া” বিটের কথা অংশ ৭-এ আবার আসবে — সেটিই catastrophic cancellation-এর মূল।

আসল float-এ: 0.1f + 0.2f, বিট ধরে ধরে

পর্ব ২৩-এ দেখেছেন double-এ 0.1 + 0.2 ≠ 0.3. Float-এ কী হয়? নিচের stepper-এ একটি আসল single-precision যোগ বিট ধরে দেখানো হলো — exponent মেলানো, যোগ, normalise আর round. সবগুলো বিট একটি script দিয়ে যাচাই করা। “পরে” চাপুন:

0.2f = 0x3E4CCCCD  s=0  e=124 (2^-3)
       1.1001 1001 1001 1001 1001 101
0.1f = 0x3DCCCCCD  s=0  e=123 (2^-4)
       1.1001 1001 1001 1001 1001 101

d = 124 - 123 = 1
smaller number: 0.1f  -> shift it right by 1
result exponent (for now): 2^-3

                                     G R S
0.2f : 1.1001 1001 1001 1001 1001 101 0 0 0
0.1f : 0.1100 1100 1100 1100 1100 110 1 0 0
       (the last 1 of 0.1f fell into G)

   1.1001 1001 1001 1001 1001 101 000
 + 0.1100 1100 1100 1100 1100 110 100
 ------------------------------------
  10.0110 0110 0110 0110 0110 011 100   x 2^-3

   1.0011 0011 0011 0011 0011 001 1 1 0   x 2^-2
                                  G R S

GRS = 110  (more than half)  -> add 1 ulp
   1.0011 0011 0011 0011 0011 001
 +                              1
   1.0011 0011 0011 0011 0011 010

e = -2 + 127 = 125   (1..254: no overflow/underflow)
0 01111101 0011 0011 0011 0011 0011 010
= 0x3E99999A = 0.300000011920929
same bits as 0.3f  -> (0.1f + 0.2f == 0.3f) is true

৩ · Guard, round আর sticky — তিনটি বাড়তি বিট কেন লাগে

Align করার সময় ছোট সংখ্যার যে বিটগুলো ডানে সরে significand-এর বাইরে চলে যায়, সেগুলো কি সোজা ফেলে দেওয়া যায়? না। তাতে correct rounding-এর প্রতিশ্রুতি ভাঙে। আবার সবগুলো বিট রাখাও অসম্ভব — exponent-এর পার্থক্য ২০০ হলে ২০০টি বাড়তি বিট লাগত। মজার ব্যাপার হলো, মাত্র তিনটি বাড়তি বিট দিয়েই সব ক্ষেত্রে সঠিক rounding পাওয়া যায়:

কেন guard লাগে — একটি দশমিক উদাহরণ

৩ অঙ্কের নির্ভুলতায় 2.34 × 10² + 2.56 × 10⁰. Align করলে 2.56 হয় 0.0256 × 10². বাড়তি অঙ্ক না রাখলে থাকে 0.02, যোগফল 2.36 × 10². দুটি বাড়তি অঙ্ক (guard 5, round 6) রাখলে 2.3400 + 0.0256 = 2.3656, round করে 2.37 × 10². আসল যোগফল 236.56; ৩ অঙ্কে সঠিক উত্তর 237, অর্থাৎ বাড়তি অঙ্ক ছাড়া উত্তর ভুল।

কেন sticky লাগে — যেখানে শুধু G আর R ভুল করে

ধরুন একটি খেলনা format: significand 1.ffff (৪ বিট ভগ্নাংশ)। যোগ করুন 1.0000 × 2⁰ + 1.0001 × 2⁻⁵. দ্বিতীয়টিকে ৫ ঘর ডানে সরালে পাই 0.0000 1 0 001 — significand-এর পরে G = 1, R = 0, আর তার পরে 001, যার মধ্যে একটি ১ আছে, তাই S = 1.

• G, R, S তিনটিই থাকলে: যোগফল 1.0000 | 101. বাদ পড়া অংশ 101 > 100, অর্থাৎ আধা ULP-এর বেশি, তাই ওপরে: 1.0001 = 1.0625.
• শুধু G আর R থাকলে: যোগফল 1.0000 | 10 — দেখতে ঠিক মাঝখান। Ties-to-even-এ 1.0000 জোড়, তাই নিচে: 1.0.

আসল যোগফল 1 + 0.033203125 = 1.033203125. দুই প্রতিবেশী 1.0 আর 1.0625-এর মাঝখান 1.03125; আসল মান মাঝখানের সামান্য ওপরে, তাই সঠিক উত্তর 1.0625. Sticky ছাড়া hardware ভেবেছে tie, অথচ tie ছিল না — ফলে ভুল দিকে round হয়েছে।

তিনটি বিটই যথেষ্ট কেন — যুক্তি

পরীক্ষায় এর যুক্তি চাওয়া হয়। দুটি ক্ষেত্র ভাবুন। (ক) exponent-এর পার্থক্য d ≤ ১: ছোট সংখ্যা সর্বোচ্চ ১ ঘর সরে, তাই G-তেই সব বাড়তি বিট ধরে যায় — হিসাব সম্পূর্ণ নির্ভুল। বিয়োগে অনেক বিট বাতিল হয়ে বড় বাঁ-shift লাগতে পারে, কিন্তু তখন হারানোর মতো কোনো বিটই ছিল না, ফল exact. (খ) d ≥ ২: বড় সংখ্যাটির exponent E হলে সেটি অন্তত 2^E, আর ছোটটি 2^E−1-এর কম — তাই বিয়োগ করলেও ফল 2^E−1-এর বেশি থাকে, অর্থাৎ normalise-এ বাঁয়ে সর্বোচ্চ ১ ঘর shift. সেই এক ঘরে G ফলের ভেতরে ঢোকে, R হয় নতুন rounding বিট, আর বাকি সব বিটের খবর লাগে শুধু “শূন্য কি না” — সেটি S দেয়। যোগে ফল ডানে সর্বোচ্চ ১ ঘর সরে, তখন G হয় rounding বিট আর R ও S মিলে নতুন sticky. দুই ক্ষেত্রেই “আধা ULP-এর চেয়ে কম, সমান, না বেশি” প্রশ্নের উত্তর তিন বিটে পুরো জানা যায়।

নিচের simulator এই খেলনা format-এর (৪ বিট ভগ্নাংশ, exponent −6 … 7) একটি adder — align, যোগ-বিয়োগ, normalise, round আর সীমা যাচাই করে, আর প্রতিটি যোগ দুবার চালায়: G, R, S দিয়ে একবার, sticky ছাড়া একবার। শেষে Math.fround দিয়ে একটি আসল float32 যোগ করে field খুলে দেখায়। Run চাপুন, তারপর cases-এ নিজের সংখ্যা বসিয়ে দেখুন:

// Toy float: sign s, exponent E, significand m = 1.ffff (16..31) plus G R S bits.
const b = (v, n) => v.toString(2).padStart(n, "0");
const show = r => b(r >> 3, 5).replace(/^(.)/, "$1.") + " " + b(r & 7, 3);
const val = t => (t.s ? -1 : 1) * t.m / 16 * 2 ** t.E;
function add(x, y, sticky) {
  if (y.E > x.E || (y.E === x.E && y.m > x.m)) [x, y] = [y, x];   // |x| >= |y|
  const d = x.E - y.E, kept = (y.m << 2) >> d;            // align: keep G, R
  const S = sticky && ((y.m << 2) & ((1 << d) - 1)) ? 1 : 0;   // OR of the rest
  const yy = (kept << 1) | S;
  let E = x.E, r = x.s === y.s ? (x.m << 3) + yy : (x.m << 3) - yy;
  if (r >= 256) { r = (r >> 1) | (r & 1); E++; }          // 1x.xxxx -> 1.xxxxx
  while (r && r < 128) { r <<= 1; E--; }                  // cancellation
  let m = r >> 3, grs = r & 7;
  if (grs > 4 || (grs === 4 && (m & 1))) m++;             // nearest, ties to even
  if (m === 32) { m = 16; E++; }
  const flag = E > 7 ? " OVERFLOW" : E < -6 ? " UNDERFLOW" : "";   // range check
  return { d, r, E, out: val({ s: x.s, m, E }) + flag };
}
const cases = [[{ s: 0, E: -1, m: 0b10000 }, { s: 1, E: -2, m: 0b11100 }],
               [{ s: 0, E: 0, m: 0b10000 }, { s: 0, E: -5, m: 0b10001 }],
               [{ s: 0, E: 0, m: 0b10000 }, { s: 1, E: -3, m: 0b10011 }]];
for (const [x, y] of cases) {
  console.log(`${val(x)} + (${val(y)})   exact = ${val(x) + val(y)}`);
  for (const st of [true, false]) {
    const t = add(x, y, st);
    console.log(`  ${st ? "G,R,S" : "G,R  "}  shift ${t.d}  normalised ${show(t.r)}` +
                ` x2^${t.E}  ->  ${t.out}`);
  }
}
const f32 = v => new Uint32Array(new Float32Array([v]).buffer)[0];   // raw bits
const a = Math.fround(0.1), c = Math.fround(0.2);
for (const v of [a, c, Math.fround(a + c), Math.fround(0.3)]) {   // float32 a + c
  const u = f32(v); console.log(`float32 ${v.toPrecision(9)}: s=${u >>> 31} e=${(u >>> 23) & 255}` +
              ` f=${b(u & 0x7FFFFF, 23)} (0x${u.toString(16)})`);
}

0.5 + (-0.4375)   exact = 0.0625
  G,R,S  shift 1  normalised 1.0000 000 x2^-4  ->  0.0625
  G,R    shift 1  normalised 1.0000 000 x2^-4  ->  0.0625
1 + (0.033203125)   exact = 1.033203125
  G,R,S  shift 5  normalised 1.0000 101 x2^0  ->  1.0625
  G,R    shift 5  normalised 1.0000 100 x2^0  ->  1
1 + (-0.1484375)   exact = 0.8515625
  G,R,S  shift 3  normalised 1.1011 010 x2^-1  ->  0.84375
  G,R    shift 3  normalised 1.1011 100 x2^-1  ->  0.875
float32 0.100000001: s=0 e=123 f=10011001100110011001101 (0x3dcccccd)
float32 0.200000003: s=0 e=124 f=10011001100110011001101 (0x3e4ccccd)
float32 0.300000012: s=0 e=125 f=00110011001100110011010 (0x3e99999a)
float32 0.300000012: s=0 e=125 f=00110011001100110011010 (0x3e99999a)

৪ · IEEE-এর চারটি rounding mode — আর ties-to-even কেন

বাদ পড়া বিট দেখে ফল ওপরে যাবে না নিচে — এই নিয়মকে বলে rounding mode. IEEE 754 binary format-এর জন্য চারটি mode দেয় (২০০৮ সালের সংস্করণে পঞ্চম একটি, ties-to-away, যোগ হয়েছে — মূলত decimal format-এর জন্য):

শেষ তিনটিকে একসাথে বলে directed rounding. এগুলোর বড় ব্যবহার interval arithmetic-এ: একই হিসাব −∞ আর +∞ দুই দিকে round করে চালালে আসল উত্তর নিশ্চিতভাবে দুটির মাঝখানে থাকে। Float থেকে int-এ রূপান্তরে C ভাষা নিজে truncate করে ((int)2.7 = 2), যা rounding mode-এর ওপর নির্ভর করে না।

Bit-এর ভাষায় চার mode

Hardware-এ কাজটি সরল: বাদ পড়া অংশকে আধা ULP-এর সাথে তুলনা করা, GRS দিয়ে। Nearest-even-এর নিয়ম: GRS < 100 হলে নিচে; > 100 হলে ওপরে; ঠিক 100 হলে শেষ বিট ০ থাকলে নিচে, ১ থাকলে ওপরে। Toward zero-তে সবসময় কেটে ফেলা। Toward +∞-তে ধনাত্মক সংখ্যায় বাদ পড়া অংশ শূন্য না হলেই ওপরে, ঋণাত্মক সংখ্যায় কেটে ফেলা (কারণ sign-magnitude রূপে পরম মান কমালে মান বাড়ে); toward −∞ এর ঠিক উল্টো। নিচের simulator কয়েকটি binary significand-কে ৪ বিট ভগ্নাংশে চার mode-এ round করে, তারপর ties-to-even-এর পক্ষপাতহীনতা হাতে-কলমে দেখায়:

// Round binary significands to K fraction bits under the four IEEE 754 modes.
const K = 4;
const modes = {
  "nearest-even": (q, rem, half, neg) =>
      rem > half || (rem === half && q % 2n === 1n) ? q + 1n : q,
  "toward 0":     (q) => q,                                   // chop
  "toward +inf":  (q, rem, half, neg) => rem > 0n && !neg ? q + 1n : q,
  "toward -inf":  (q, rem, half, neg) => rem > 0n && neg ? q + 1n : q,
};
function roundBin(str, mode) {                 // str like "-1.0101 1"
  const neg = str.startsWith("-"), [ip, fp] = str.replace(/[-+ ]/g, "").split(".");
  const extra = fp.length - K, all = BigInt("0b" + ip + fp);
  const q = all >> BigInt(extra), rem = all - (q << BigInt(extra));
  const half = 1n << BigInt(extra - 1);
  const r = modes[mode](q, rem, half, neg).toString(2).padStart(K + 1, "0");
  return (neg ? "-" : "+") + r.slice(0, -K) + "." + r.slice(-K);
}
const inputs = ["+1.0100 01", "+1.0100 1", "+1.0101 1", "+1.0100 11",
                "-1.0100 11", "-1.0101 1", "+1.1111 1"];
console.log("exact        " + Object.keys(modes).map(m => m.padEnd(13)).join(""));
for (const s of inputs)
  console.log(s.padEnd(13) + Object.keys(modes)
    .map(m => roundBin(s, m).padEnd(13)).join(""));

// Bias check: round 0.5, 1.5, ..., 7.5 to integers and add them up.
const ties = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5];
const even = x => { const f = Math.floor(x); return x - f === 0.5 ? (f % 2 ? f + 1 : f) : Math.round(x); };
const sum = a => a.reduce((s, x) => s + x, 0);
console.log("\nexact sum        :", sum(ties));
console.log("half-up rounding :", sum(ties.map(x => Math.floor(x + 0.5))));
console.log("ties-to-even     :", sum(ties.map(even)), ties.map(even).join(" "));

exact        nearest-even toward 0     toward +inf  toward -inf  
+1.0100 01   +1.0100      +1.0100      +1.0101      +1.0100      
+1.0100 1    +1.0100      +1.0100      +1.0101      +1.0100      
+1.0101 1    +1.0110      +1.0101      +1.0110      +1.0101      
+1.0100 11   +1.0101      +1.0100      +1.0101      +1.0100      
-1.0100 11   -1.0101      -1.0100      -1.0100      -1.0101      
-1.0101 1    -1.0110      -1.0101      -1.0101      -1.0110      
+1.1111 1    +10.0000     +1.1111      +10.0000     +1.1111      

exact sum        : 32
half-up rounding : 36
ties-to-even     : 32 0 2 2 4 4 6 6 8

কেন এটি গুরুত্বপূর্ণ: Binary-তে হিসাবের ফল প্রায়ই ঠিক মাঝখানে পড়ে — বিশেষ করে exponent-এর পার্থক্য ছোট হলে, বা ২ দিয়ে ভাগে। লাখো operation-এ প্রতিটি tie ওপরে গেলে গড় ভুল আর শূন্য থাকে না, একদিকে জমতে থাকে। Ties-to-even দীর্ঘমেয়াদে গড় ভুল শূন্যের কাছে রাখে। এজন্য একে “banker's rounding”-ও বলে।

৫ · গুণ, ভাগ আর FP unit-এর hardware

অবাক করার মতো হলেও floating point গুণ যোগের চেয়ে সহজ — exponent মেলাতে হয় না। নিয়ম পাঁচটি: (১) biased exponent দুটি যোগ করে একবার bias বাদ (কারণ দুটিতেই bias আছে, যোগফলে দুবার); (২) significand দুটি গুণ (unsigned integer গুণ, পর্ব ২১); (৩) normalise — দুটি 1.x গুণ করলে ফল ১ থেকে ৪-এর কম, তাই সর্বোচ্চ ১ ঘর ডানে; (৪) round আর সীমা যাচাই; (৫) চিহ্ন = s₁ ⊕ s₂ (একই চিহ্ন হলে ধনাত্মক)।

দুটি worked উদাহরণ

দশমিকে (৪ অঙ্ক): 1.110 × 10¹⁰ × 9.200 × 10⁻⁵. Exponent: 10 + (−5) = 5. Significand: 1.110 × 9.200 = 10.212. Normalise: 1.0212 × 10⁶. Round: 1.021 × 10⁶. চিহ্ন: দুটিই ধনাত্মক, ফল ধনাত্মক।

Binary-তে (বইয়ের উদাহরণ): 0.5 × (−0.4375) = 1.000₂ × 2⁻¹ × (−1.110₂ × 2⁻²).

ফল −1.110₂ × 2⁻³ = −0.21875, আর 0.5 × (−0.4375) সত্যিই −0.21875. আসল single precision-এ একটি উদাহরণ: 1.5f (0x3FC00000, e = 127) × 2.75f (0x40300000, e = 128). Exponent 127 + 128 − 127 = 128; significand 1.1₂ × 1.011₂ = 10.0001₂ (1.5 × 1.375 = 2.0625); normalise করে 1.00001 × 2¹, তাই biased exponent 129. ফল 4.125 = 0x40840000.

ভাগ — সংক্ষেপে

ভাগে উল্টো নিয়ম: exponent বিয়োগ করে bias আবার যোগ (e₁ − e₂ + bias), significand ভাগ, normalise (ফল ০.৫ থেকে ২-এর কম, তাই সর্বোচ্চ ১ ঘর বাঁয়ে), round, আর চিহ্ন আবার XOR. Significand ভাগের পদ্ধতিগুলো — restoring, non-restoring, SRT, আর Newton–Raphson দিয়ে reciprocal বের করা — আর SRT টেবিলের ভুলে ১৯৯৪ সালের বিখ্যাত Pentium FDIV বাগ পর্ব ২২-এ (ভাগ) পাবেন। ভাগ আর বর্গমূল গুণের চেয়ে অনেক ধীর, আর প্রায়ই পুরোপুরি pipelined নয়।

FP adder-এর hardware

অংশ ২-এর পাঁচ ধাপ hardware-এ একটি পাইপের মতো সাজানো থাকে। নিচের ছবিটি Patterson-Hennessy-র FP adder-এর সরল রূপ:

Multiplier-এর গঠন আরও সরল: exponent-এর জন্য একটি adder (bias বাদসহ), significand-এর জন্য একটি বড় integer multiplier (single-এ ২৪ × ২৪, double-এ ৫৩ × ৫৩ বিট — partial product-এর গাছ, পর্ব ২১), তারপর ১ ঘরের normalise, round আর চিহ্নের একটি XOR গেট।

কত দ্রুত — latency আর pipelining

প্রতিটি ধাপ একটি clock cycle-এ শেষ হয় না বলে আধুনিক CPU এগুলোকে কয়েকটি pipeline ধাপে ভাগ করে। মোটামুটি হিসেবে, সাম্প্রতিক x86-64 আর বড় ARM core-এ FP যোগ, গুণ বা FMA-র latency প্রায় ৩ থেকে ৫ cycle, কিন্তু pipelined বলে প্রতি cycle-এ নতুন একটি (প্রায়ই দুটি unit-এ দুটি) operation শুরু করা যায় — অর্থাৎ throughput প্রতি cycle-এ ১–২টি। ভাগ আর বর্গমূলে latency সাধারণত ১০ থেকে কয়েক দশক cycle. সঠিক সংখ্যা chip-ভেদে বদলায় — নির্মাতার optimization manual দেখুন। Pipelining-এর পুরো গল্প পর্ব ৩৫-এ, আর একটি instruction-এ অনেকগুলো float একসাথে হিসাব (SIMD) পর্ব ৪২-এ।

# উদাহরণ — RISC-V F extension (single precision) instructions
fadd.s   ft2, ft0, ft1          # ft2 = ft0 + ft1, rounding mode from fcsr.frm
fadd.s   ft2, ft0, ft1, rtz     # same add, but round toward zero for this one
fmul.s   ft3, ft0, ft1          # ft3 = ft0 * ft1
fmadd.s  ft4, ft0, ft1, ft2     # ft4 = ft0 * ft1 + ft2 with ONE rounding
frflags  t0                     # read the 5 sticky flags NV DZ OF UF NX

; x86-64 with FMA3 (Intel syntax)
vfmadd231sd xmm0, xmm1, xmm2    ; xmm0 = xmm1 * xmm2 + xmm0, one rounding

Fused multiply-add (FMA) — এক rounding-এ a × b + c

উপরের fmadd.s একটি বিশেষ instruction: fused multiply-add (FMA). সাধারণভাবে a × b + c হিসাব করলে দুবার round হয় — গুণফলে একবার, যোগফলে একবার। FMA পুরো গুণফল (single-এ ৪৮ বিট) না কেটে রেখে c-র সাথে যোগ করে, তারপর একবারই round করে। লাভ দুটি: নির্ভুলতা (একটি rounding কম, আর গুণফলের নিচের বিট হারায় না) আর গতি (দুটি instruction-এর কাজ একটিতে, প্রায় একটির latency-তে)। Dot product, matrix গুণ, বহুপদীর Horner পদ্ধতি — সবই “গুণ করে যোগ” এর লম্বা শিকল, তাই বৈজ্ঞানিক হিসাব আর AI-এর matrix কাজে FMA-ই মূল ইট। IBM-এর POWER1 (১৯৯০) FMA hardware-এ আনে বলে ধরা হয়; IEEE 754-2008 একে মানদণ্ডের অংশ করে; x86-এ FMA3 এসেছে মোটামুটি ২০১২–২০১৩ সালে (AMD Piledriver, Intel Haswell). অংশ ৬-এর C প্রোগ্রামে FMA আর সাধারণ a × b + c-এর পার্থক্য চোখে দেখবেন।

৬ · পাঁচটি exception আর flag — প্রোগ্রাম থামে না, শুধু খাতায় লেখা থাকে

Integer-এ শূন্য দিয়ে ভাগ করলে x86-এ প্রোগ্রাম সাধারণত crash করে। Floating point-এ কিন্তু 1.0 / 0.0 করলে প্রোগ্রাম চলতেই থাকে, ফল আসে +∞. কারণ IEEE 754-এর নকশা: পাঁচ রকম floating point exception-এর প্রতিটির একটি ডিফল্ট ফল আছে, আর ঘটনাটি একটি status flag-এ লিখে রাখা হয়।

Flag-গুলো sticky: একবার উঠলে নিজে থেকে নামে না, প্রোগ্রাম মুছে না দেওয়া পর্যন্ত থাকে। ফলে লাখো হিসাবের শেষে একবার দেখে নেওয়া যায় — “এর মধ্যে কোথাও কি NaN জন্মেছিল?” চাইলে trap চালু করা যায়, তখন exception হলেই CPU interrupt দিয়ে প্রোগ্রাম থামায় (debug-এর জন্য কাজের); কিন্তু ডিফল্টে সব trap বন্ধ (masked). C ভাষায় এই flag দেখার মানক পথ <fenv.h>: feclearexcept(FE_ALL_EXCEPT) দিয়ে মোছা, fetestexcept(FE_INVALID) দিয়ে দেখা। Linux-এর glibc-তে এই function-গুলো math library-তে থাকে, তাই নিজের কম্পিউটারে gcc prog.c -lm লাগে।

নিচের প্রোগ্রাম কোনো library ছাড়াই সরাসরি hardware-এর flag register পড়ে — x86-64-এর MXCSR, যার নিচের ছয়টি বিট এই flag (বিট ১, denormal-input, IEEE-র অংশ নয়, তাই বাদ)। তাই এটি শুধু x86-64-এ চলবে; বর্গমূলের জন্য সরাসরি SSE2-র sqrtsd instruction, আর শেষে FMA বনাম a × b + c:

#include <stdio.h>
#include <float.h>
/* x86-64 only: MXCSR status bits IE=1 ZE=4 OE=8 UE=16 PE=32 */
static void show(const char *what, double r) {
    unsigned f = __builtin_ia32_stmxcsr();
    printf("%-11s = %-12g flags:%s%s%s%s%s\n", what, r,
           (f & 1) ? " invalid" : "", (f & 4) ? " div-by-zero" : "",
           (f & 8) ? " overflow" : "", (f & 16) ? " underflow" : "",
           (f & 32) ? " inexact" : "");
    __builtin_ia32_ldmxcsr(f & ~63u);            /* clear the sticky flags */
}
static double hw_sqrt(double x) {               /* the SSE2 sqrtsd instruction */
    double r; __asm__("sqrtsd %1, %0" : "=x"(r) : "x"(x)); return r;
}
__attribute__((target("fma")))
static double hw_fma(double a, double b, double c) { return __builtin_fma(a, b, c); }
int main(void) {
    volatile double zero = 0.0, one = 1.0, big = DBL_MAX, tiny = DBL_MIN;
    volatile double p1 = 0.1, p2 = 0.2, a = 1.0 + 0x1p-27;
    show("1.0 / 0.0", one / zero);
    show("0.0 / 0.0", zero / zero);
    show("sqrt(-1.0)", hw_sqrt(-one));
    show("DBL_MAX * 2", big * 2.0);
    show("DBL_MIN / 3", tiny / 3.0);
    show("0.1 + 0.2", p1 + p2);
    show("0.5 + 0.25", one * 0.5 + 0.25);
    /* exact a*a = 1 + 2^-26 + 2^-54; a*a - 1 rounds the product first */
    printf("a*a - 1       = %.17g\n", a * a - 1.0);
    if (__builtin_cpu_supports("fma"))
        printf("fma(a, a, -1) = %.17g\n", hw_fma(a, a, -1.0));
    else
        puts("no FMA instruction on this CPU");
    printf("2^-26 + 2^-54 = %.17g\n", 0x1p-26 + 0x1p-54);
    return 0;
}

1.0 / 0.0   = inf          flags: div-by-zero
0.0 / 0.0   = -nan         flags: invalid
sqrt(-1.0)  = -nan         flags: invalid
DBL_MAX * 2 = inf          flags: overflow inexact
DBL_MIN / 3 = 7.41691e-309 flags: underflow inexact
0.1 + 0.2   = 0.3          flags: inexact
0.5 + 0.25  = 0.75         flags:
a*a - 1       = 1.4901161193847656e-08
fma(a, a, -1) = 1.4901161249358807e-08
2^-26 + 2^-54 = 1.4901161249358807e-08

৭ · যখন গণিতের নিয়ম ভাঙে — associativity, cancellation আর quadratic সূত্র

বাস্তব সংখ্যায় (a + b) + c = a + (b + c) — এটি associativity. Floating point-এ প্রতিটি যোগের পরে round হয় বলে ক্রম বদলালে ফল বদলাতে পারে; একে বলে non-associativity. Double-এ 10¹⁶-এর আশেপাশে পাশাপাশি দুটি মানের দূরত্ব (ULP) ২. তাই:

• (10¹⁶ + 1) − 10¹⁶: প্রথম যোগফল 10¹⁶ + 1 ঠিক দুই প্রতিবেশী 10¹⁶ আর 10¹⁶ + 2-এর মাঝখানে; ties-to-even 10¹⁶-এ ফেরায়। তারপর বিয়োগে 0.
• (10¹⁶ − 10¹⁶) + 1 = 0 + 1 = 1.

একই কারণে বড় x-এর জন্য x + 1 == x সত্য: float-এ 2²⁴ = 16777216-এর পরে ULP ২, তাই 16777216 + 1 আবার 16777216 (16777217 float-এ রাখাই যায় না); double-এ একই সীমা 2⁵³. লুপের counter float রাখলে এক সময় সেটি বাড়া বন্ধ করে দেয় — চিরকাল চলা লুপ। Distributivity-ও ভাঙে: double-এ 10 × (0.1 + 0.2) = 3.0000000000000004, কিন্তু 10 × 0.1 + 10 × 0.2 = 3. আর (0.1 + 0.2) + 0.3 = 0.6000000000000001, অথচ 0.1 + (0.2 + 0.3) = 0.6.

Catastrophic cancellation — প্রায় সমান দুটি সংখ্যার বিয়োগ

বিয়োগ নিজে খুব নির্ভুল: দুটি সংখ্যা পরস্পরের দুই গুণের মধ্যে থাকলে তাদের পার্থক্য floating point-এ নির্ভুলভাবে রাখা যায় (Sterbenz lemma). বিপদ অন্যখানে। ধরুন a আর b দুটিই আগের হিসাব থেকে এসেছে, প্রতিটিতে শেষের এক-দুই বিট ভুল। a − b করলে সামনের মিলে যাওয়া বিটগুলো সব বাতিল হয়ে যায় — থাকে শুধু শেষের সেই ভুল-মেশানো বিটগুলো, আর normalise সেগুলোকে টেনে সামনে আনে। আপেক্ষিক ভুল হঠাৎ বিশাল। একে বলে catastrophic cancellation.

পাঠ্যবইয়ের সবচেয়ে বিখ্যাত উদাহরণ দ্বিঘাত সমীকরণের সূত্র x = (−b ± √(b² − 4ac)) / 2a, যখন b² ≫ 4ac. তখন √(b² − 4ac) প্রায় |b|-এর সমান, তাই −b + √(…)-এ প্রায় সমান দুটি সংখ্যা বিয়োগ হয়। প্রতিকার একটি পুরোনো বীজগণিত: দুটি মূলের গুণফল c/a (Vieta-র সূত্র)। তাই প্রথমে বাতিলহীন মূলটি বের করুন — q = −(b + sign(b)·√(b² − 4ac)) / 2, x₁ = q / a — তারপর অন্যটি x₂ = c / q. নিচের stepper-এ x² + 10⁸x + 1 = 0, double-এ:

true roots:  x1 = -99999999.99999999  (about -1e8)
             x2 = -1.0000000000000000e-08
small root is what we care about

b*b       = 1e16              (exact)
4*a*c     = 4
b*b - 4ac = 9999999999999996  (exact: ulp here is 2)

true sqrt = 99999999.99999998000...
double    = 99999999.99999998510  = 1e8 - 2^-26
            (ulp near 1e8 is 2^-26 = 1.49e-8)

-1e8 + (1e8 - 2^-26) = -2^-26   (cancellation!)
x2 = -2^-26 / 2 = -7.4505805969238281e-09
relative error = 25.5 %

q  = -(b + d) / 2 = -1e8      (same signs: no cancellation)
x1 = q / a        = -1e8
x2 = c / q        = 1 / -1e8 = -1e-08

naive  x2 = -7.45e-09   error 2.5e-01
stable x2 = -1.00e-08   error ~1e-16
check: x1 * x2 = c / a = 1

৮ · যোগফলে ভুল জমা — Kahan আর pairwise summation

লম্বা তালিকা যোগ করা — গড়, মোট বিক্রি, dot product — কম্পিউটারের সবচেয়ে সাধারণ কাজগুলোর একটি। সাধারণ লুপ s += x[i]-এ প্রতিটি যোগে নতুন rounding ভুল ঢোকে, আর যোগফল যত বড় হয়, প্রতিটি নতুন ছোট সংখ্যার তত বেশি বিট align-এর সময় ডানে পড়ে যায়। একে বলে accumulation error. n টি সংখ্যার সাধারণ যোগে ভুলের সীমা মোটামুটি (n − 1) × u × Σ|x_i|, যেখানে u আধা ULP-এর আপেক্ষিক মান (float-এ প্রায় 6 × 10⁻⁸). n = এক কোটি হলে এই সীমা ভয়ংকর।

বাস্তবে দেখুন: float-এ 0.1f এক কোটি বার যোগ করলে কত আসে? 0.1f নিজেই 0.100000001490116…, তাই আসল উত্তর প্রায় 1000000.0149. কিন্তু যোগফল 2¹⁹ = 524288 পেরোলে float-এর ঘর 0.0625 দূরে দূরে, 2²⁰ পেরোলে 0.125 দূরে — প্রতিটি যোগের rounding ভুল বড়, আর তা একদিকে জমে। নিচের C প্রোগ্রামে দেখবেন উত্তর 1087937 — প্রায় ৮.৮% বেশি।

Kahan summation — হারানো বিটগুলো আলাদা খাতায়

William Kahan ১৯৬৫ সালে এই পদ্ধতি প্রকাশ করেন; একে বলে Kahan (compensated) summation. যোগফল sum-এর পাশে একটি ছোট সংশোধনী c রাখা হয়:

/* উদাহরণ — Kahan compensated summation */
float sum = 0.0f, c = 0.0f;
for (i = 0; i < n; i++) {
    float y = x[i] - c;      /* first fix the next number with the old error  */
    float t = sum + y;       /* big + small: low bits of y are lost here      */
    c = (t - sum) - y;       /* (t - sum) is what was really added; minus y
                                gives the lost part (with a minus sign)       */
    sum = t;
}

কৌশলটি হলো (t - sum) - y: t আর sum কাছাকাছি, তাই t − sum নির্ভুলভাবে বলে দেয় y-এর কতটুকু আসলে যোগ হয়েছে; সেখান থেকে y বাদ দিলে পাওয়া যায় হারানো অংশ। গণিতে c সবসময় শূন্য — তাই -ffast-math দিলে compiler সত্যিই একে মুছে দিতে পারে! নিচের stepper-এ ৪ অঙ্কের দশমিক যন্ত্রে 1000-এর সাথে 3.142 চারবার যোগ:

sum = 1000    c = 0      naive = 1000
add 3.142 four times
exact answer = 1012.568

y   = 3.142 - 0          = 3.142
t   = 1000 + 3.142       = 1003.142 -> 1003
c   = (1003 - 1000) - 3.142 = -0.142
sum = 1003               naive = 1003

y   = 3.142 - (-0.142)   = 3.284
t   = 1003 + 3.284       = 1006.284 -> 1006
c   = (1006 - 1003) - 3.284 = -0.284
sum = 1006               naive = 1006

y   = 3.142 + 0.284      = 3.426
t   = 1006 + 3.426       = 1009.426 -> 1009
c   = (1009 - 1006) - 3.426 = -0.426
sum = 1009               naive = 1009

y   = 3.142 + 0.426      = 3.568
t   = 1009 + 3.568       = 1012.568 -> 1013
c   = (1013 - 1009) - 3.568 = +0.432
sum = 1013               naive = 1009 + 3.142 -> 1012

exact          = 1012.568  (4 digits: 1013)
naive          = 1012      (lost 0.568)
Kahan sum      = 1013      (correctly rounded)
sum - c        = 1013 - 0.432 = 1012.568  (c keeps the rest)

Pairwise summation — যোগের ক্রম বদলে দেওয়া

আরেকটি সরল কৌশল: তালিকাকে দুই ভাগ করে প্রতিটি অর্ধেক আলাদা যোগ, তারপর দুই যোগফল যোগ — আর একইভাবে ভেতরে ভেতরে (একটি binary tree). প্রতিটি সংখ্যা মাত্র log₂ n বার rounding-এর মধ্য দিয়ে যায়, তাই ভুলের সীমা n-এর বদলে log n-এর অনুপাতে বাড়ে। একে বলে pairwise summation. একটি আলাদা script-এ float32-এ এক কোটি 0.1f-কে এভাবে tree ধরে যোগ করে পাওয়া গেছে ঠিক 1000000.0. বাড়তি খরচ প্রায় নেই, আর এটি parallel (দুই অর্ধেক দুই core-এ) করাও সহজ — তাই NumPy-র মতো library ভেতরে এর একটি রূপ ব্যবহার করে বলে জানা যায়। (Python-এর math.fsum আরও এগিয়ে — correctly rounded যোগফল দেয়।)

নিচের প্রোগ্রাম এই অংশ আর আগের অংশের তিনটি দাবি একসাথে যাচাই করে: ক্রম বদলালে ফল বদলায়, float-এ বড় x-এ x + 1 == x, এক কোটি 0.1f-এর তিন রকম যোগ, আর quadratic সূত্রের দুই রূপ (বর্গমূলের জন্য আবার সরাসরি sqrtsd, তাই x86-64):

#include <stdio.h>
static double hw_sqrt(double x) {               /* x86-64 sqrtsd, no libm needed */
    double r;
    __asm__("sqrtsd %1, %0" : "=x"(r) : "x"(x));
    return r;
}
int main(void) {
    volatile double big = 1e16, one = 1.0;
    printf("(1e16 + 1) - 1e16 = %g\n", (big + one) - big);
    printf("(1e16 - 1e16) + 1 = %g\n", (big - big) + one);
    volatile float f24 = 16777216.0f;           /* 2^24 */
    printf("2^24f + 1 == 2^24f : %d\n", f24 + 1.0f == f24);
    const int N = 10000000;
    float naive = 0.0f, sum = 0.0f, c = 0.0f;   /* c = Kahan compensation */
    double dsum = 0.0;
    for (int i = 0; i < N; i++) {
        naive += 0.1f;
        dsum  += 0.1f;
        float y = 0.1f - c, t = sum + y;
        c = (t - sum) - y;                      /* the low bits that were lost */
        sum = t;
    }
    printf("float naive = %.1f\ndouble      = %.6f\nfloat Kahan = %.1f\n",
           naive, dsum, sum);
    volatile double a = 1.0, b = 1e8, cc = 1.0;   /* roots near -1e-8 and -1e8 */
    double d = hw_sqrt(b * b - 4 * a * cc);
    double naive_x = (-b + d) / (2 * a);        /* subtracts two nearly equal numbers */
    double q = -(b + d) / 2, stable_x = cc / q; /* no cancellation */
    printf("naive  small root = %.17g\nstable small root = %.17g\n", naive_x, stable_x);
    return 0;
}

(1e16 + 1) - 1e16 = 0
(1e16 - 1e16) + 1 = 1
2^24f + 1 == 2^24f : 1
float naive = 1087937.0
double      = 1000000.014901
float Kahan = 1000000.0
naive  small root = -7.4505805969238281e-09
stable small root = -1e-08

৯ · ভিন্ন যন্ত্রে ভিন্ন উত্তর, float তুলনা, টাকা আর AI-এর mixed precision

একই কোড, ভিন্ন উত্তর — x87 আর FMA

IEEE 754 প্রতিটি operation-এর ফল নির্দিষ্ট করে, তবু একই C প্রোগ্রাম দুই কম্পিউটারে শেষ অঙ্কে আলাদা উত্তর দিতে পারে। দুটি প্রধান কারণ ইতিহাসে আর বর্তমানে:

• x87 আর ৮০ বিট: Intel-এর 8087 coprocessor (১৯৮০) আর তার পরের x87 FPU ভেতরে সব হিসাব করত ৮০ বিটের extended precision-এ (৬৪ বিট significand). ৩২ বিট x86-এর পুরোনো compiler-এ double-এর মধ্যবর্তী ফল register-এ ৮০ বিটে থাকত, আর মেমরিতে লেখার (spill) সময় ৬৪ বিটে round হতো। কোন মান কখন মেমরিতে যাবে তা optimization-এর ওপর নির্ভর করে, তাই -O0 আর -O2-এ ফল বদলাত; আর দুবার round হওয়া (double rounding) কখনো ১ ULP সরিয়ে দিত। x86-64-এর ABI-তে float আর double SSE2-তে হয়, তাই এই সমস্যা আজ অনেক কম (C-তে FLT_EVAL_METHOD 0).
• FMA contraction: a*b + c লিখলে compiler চাইলে একে একটি FMA instruction বানাতে পারে (C মানদণ্ড এই “contraction” অনুমতি দেয়)। FMA-ওয়ালা CPU-তে, বা -march=native দিলে, ফল এক rounding-এর; FMA ছাড়া দুই rounding-এর — অংশ ৬-এর C প্রোগ্রামে যেমন দেখেছেন, পার্থক্য শেষ বিটে। GCC-র নথি অনুযায়ী GNU mode-এ ডিফল্ট -ffp-contract=fast, কিন্তু -std=c11-এর মতো ISO mode-এ বন্ধ। নিশ্চিত মিল চাইলে -ffp-contract=off.

Library function-এর শেষ-বিটের পার্থক্য আর parallel যোগের ক্রম — আরও দুটি কারণ — ভাবনার প্রশ্ন ৩-এ।

Float তুলনা — কখনো সরাসরি == নয়

0.1 + 0.2 আর 0.3 double-এ ঠিক এক ULP দূরে (5.55 × 10⁻¹⁷), তাই == মিথ্যা বলে। তিনটি সঠিক উপায়:

1. Absolute tolerance: |a − b| ≤ ε_abs. শূন্যের কাছে কাজের, কিন্তু বড় সংখ্যায় অর্থহীন — 10⁹-এর কাছে double-এর ঘরই প্রায় 1.2 × 10⁻⁷ দূরে, তাই ε_abs = 10⁻⁹ সেখানে ==-এর মতোই কঠোর।
2. Relative tolerance: |a − b| ≤ ε_rel × max(|a|, |b|). মানের আকারের সাথে মানিয়ে নেয়; কিন্তু একটি মান ঠিক ০ হলে কাজ করে না। তাই বাস্তবে দুটি একসাথে: absolute দিয়ে শূন্যের কাছে, relative দিয়ে বাকি জায়গায়। Python-এর math.isclose ঠিক এভাবে কাজ করে (ডিফল্ট rel_tol = 10⁻⁹).
3. ULP দূরত্ব: একই চিহ্নের দুটি float-এর বিটকে integer হিসেবে পড়লে তাদের পার্থক্যই বলে মাঝখানে কয়টি float আছে — কারণ IEEE-র বিন্যাসে বিট-প্যাটার্নের ক্রম আর মানের ক্রম এক (excess-K exponent-এর সুফল, পর্ব ০৫)। 0.1 + 0.2 আর 0.3-এর দূরত্ব ঠিক ১ ULP. “৪ ULP-এর মধ্যে হলে সমান” ধরনের পরীক্ষা গাণিতিক library-র test-এ প্রচলিত।

Catastrophic cancellation পার হয়ে আসা মানে লাখো ULP ভুল থাকতে পারে, তাই tolerance ঠিক করুন হিসাবের ইতিহাস দেখে। NaN-এর কথাও মনে রাখুন: NaN == NaN মিথ্যা, তাই x != x দিয়ে NaN ধরা যায় (তবে -ffast-math-এ এই পরীক্ষাও মুছে যেতে পারে; isnan() নিরাপদ)।

টাকার হিসাব — float নয়, পূর্ণসংখ্যা পয়সা

১০ পয়সা মানে 0.10 টাকা, আর 0.1 binary-তে অসীম ভগ্নাংশ — লাখো লেনদেনে এক-দুই পয়সার অমিল জমতে পারে। নিয়ম: টাকা রাখুন পূর্ণসংখ্যা পয়সায় (১,২৩৪.৫০ টাকা = 123450 পয়সা, ৬৪ বিট integer-এ) বা decimal type-এ, আর rounding-এর নিয়ম স্পষ্ট লিখে রাখুন। বিস্তারিত পাবেন পর্ব ২৩-এর Floating point অংশে।

বাস্তবের ঘটনা

আরও দুটি ঘটনা আগের পর্বগুলোতে আছে, এখানে শুধু সংযোগটুকু: Patriot ক্ষেপণাস্ত্র-প্রতিরোধ ব্যবস্থা (১৯৯১) — ০.১ সেকেন্ডকে কেটে রাখার ছোট ভুল ঘণ্টার পর ঘণ্টা জমে বড় হয়েছিল; বিবরণ পর্ব ০৪-এ (সংখ্যা পদ্ধতি)। Ariane 5 (১৯৯৬) — ৬৪ বিট floating point থেকে ১৬ বিট integer-এ রূপান্তরের overflow; বিবরণ পর্ব ০৫-এ। তিনটির শিক্ষা একই: ছোট, একমুখী ভুল সময়ের সাথে জমে; আর সীমা ছাড়ানোর পরিকল্পনা না থাকলে বিপদ।

AI-তে mixed precision — নিচু precision-এ গুণ, উঁচুতে যোগ

বড় neural network-এর প্রশিক্ষণ আর চালানোর মূল কাজ বিশাল matrix গুণ — কোটি কোটি FMA. মেমরি আর শক্তি বাঁচাতে আধুনিক AI hardware গুণ করে ১৬ বিটে — FP16 বা bfloat16-এ (বিট-বিন্যাস পর্ব ২৩-এ) — কিন্তু যোগফল জমায় ৩২ বিট FP32-এ। কেন? অংশ ৮-ই উত্তর: dot product মানে হাজারো গুণফলের যোগ। bfloat16-এ 256 + 1 = 256 (৮ বিটে 257 রাখাই যায় না), FP16-এ 2048 + 1 = 2048. ছোট precision-এ জমাতে থাকলে যোগফল একটু বড় হলেই নতুন গুণফলগুলো আর যোগই হয় না (“stagnation”). তাই NVIDIA-র Tensor Core আর Google-এর TPU-র মতো hardware-এর মূল নকশা: নিচু precision-এ গুণ, FP32-এ জমা — একে বলে mixed precision (বিস্তারিত ভাবনার প্রশ্ন ৪-এ; GPU-র ভেতরের গঠন পর্ব ৫৫-এ)।

১০ · AI-এর যুগে floating point-এর হিসাব

AI টুল floating point-এর সংজ্ঞা ঝরঝরে বলতে পারে, কিন্তু বিট-স্তরের হিসাবে প্রায়ই ভুল করে — আর সেই ভুল দেখতে একদম আত্মবিশ্বাসী। কয়েকটি প্রচলিত ভুল:

AI-এর উত্তর হাতে পেলে — যাচাইয়ের চেকলিস্ট

• Format (float/double/bfloat16), bias আর rounding mode স্পষ্ট বলা আছে
• প্রতিটি operand-এর hex নিজে script দিয়ে (Float32Array, struct.pack) মিলিয়েছি
• Align-এর shift সংখ্যা = exponent-এর পার্থক্য, আর G, R, S হাতে মিলিয়েছি
• চূড়ান্ত ফল আসল মানের সবচেয়ে কাছের float কি না, প্রতিবেশী দুটি float বের করে দেখেছি
• AI-এর লেখা সংখ্যাগত কোডে == তুলনা, বাতিল-প্রবণ সূত্র বা লম্বা naive যোগ আছে কি না খুঁজেছি
• Chip-এর latency বা instruction-এর দাবি নির্মাতার manual বা ISA manual-এ মিলিয়েছি