Float Calculator

কেন এই প্রজেক্ট

IEEE 754 লেসনে আমরা bit pattern থেকে দশমিক মান বের করেছি — decode। কিন্তু CPU-র ভেতরের FPU যা করে সেটা তার চেয়ে বেশি: এটা দুইটা bit pattern নিয়ে তৃতীয় একটা বৈধ bit pattern তৈরি করে, আর সেই তৈরির নিয়মটাই IEEE 754 standard-এর আসল শরীর। “IEEE 754 একটা storage format” — এই ধারণাটা এখানে ভেঙে যায়। IEEE 754 আসলে একটা hardware algorithm, আর addition তার সবচেয়ে basic অপারেশন।

এই প্রজেক্টে আপনি ঠিক সেই algorithm-টা হাতে লিখবেন — exponent align করা, mantissa যোগ করা, ফলাফল renormalize করা, round করা। প্রতিটা ধাপে আপনি এমন সব প্রশ্নের মুখোমুখি হবেন যেগুলো decode করার সময় আসে না: দুইটা float-এর exponent আলাদা হলে যোগ করব কীভাবে? যোগফল যদি ঠিক মাঝামাঝি পড়ে, কোন দিকে round করব? বিয়োগে দুইটা প্রায়-সমান সংখ্যা একে অপরকে প্রায় বাতিল করে দিলে কী ঘটে? — এগুলোই floating-point precision লেসনের abstract সতর্কতাগুলোকে (কেন 0.1 + 0.2 != 0.3) যান্ত্রিক বাস্তবতায় পরিণত করবে।

লক্ষ্য — একটা সহজ উদাহরণ হাতে করে

কোড লেখার আগে 1.5 + 0.25 হাতে করে দেখা যাক, algorithm-টা বোঝার জন্য:

  1.5  = 1.10000000000000000000000 × 2^0    (bits: sign=0 exp=127 mantissa→ 0x3FC00000)
  0.25 = 1.00000000000000000000000 × 2^-2   (bits: sign=0 exp=125 mantissa→ 0x3E800000)

  align   : exponent ফারাক 0-(-2)=2, তাই 0.25-এর mantissa ২ ঘর ডানে শিফট হবে
            0.25 এর mantissa, exponent-0 এর ভাষায়: 0.01000000000000000000000
  add     : 1.10000000000000000000000
          + 0.01000000000000000000000
          -------------------------
            1.11000000000000000000000   × 2^0
  normalize: ইতিমধ্যে 1.xxx রেঞ্জে আছে — শিফটের দরকার নেই
  round    : বাদ পড়া কোনো bit নেই (guard/round/sticky = 000) — round করার কিছু নেই

  ফলাফল = 1.75  →  bits = 0x3FE00000

0x3FC00000 (1.5), 0x3E800000 (0.25), 0x3FE00000 (1.75) — এই তিনটাই IEEE 754 single-precision-এর সুপরিচিত constant, যাচাই করা যায় হাতে-কলমে exponent bias (127) আর mantissa bit বসিয়ে। এই উদাহরণে কোনো rounding লাগেনি বলেই সহজ — নিচের কোডে আমরা এমন কেসও বানাব যেখানে rounding, renormalization, এমনকি overflow পর্যন্ত ঘটে।

ধাপে ধাপে

১. Decode — uniform representation-এ আনা

প্রতিটা operand-কে (sign, e, M, kind)-এ ভাঙব, যেখানে M সবসময় ২৪-bit significand (normal-এর জন্য implicit 1 bit-সহ, subnormal-এর জন্য ছাড়া), আর e এমন একটা exponent যাতে value = (-1)^sign × M × 2^(e-23) — normal আর subnormal দুইটার জন্যই এই একই সূত্র কাজ করে (subnormal-এর e সবসময় -126, normal-এর ন্যূনতম exponent-এর সমান — এই সমন্বয়টাই দুই case-কে একই কোডে ধরার চাবিকাঠি)।

def decode(bits_pattern: int):
    sign = (bits_pattern >> 31) & 1
    exp_field = (bits_pattern >> 23) & 0xFF
    mantissa = bits_pattern & 0x7FFFFF

    if exp_field == 0xFF:
        return sign, None, None, ('nan' if mantissa else 'inf')
    if exp_field == 0:
        if mantissa == 0:
            return sign, None, None, 'zero'
        return sign, -126, mantissa, 'subnormal'              # implicit bit = 0
    return sign, exp_field - 127, (1 \<\< 23) | mantissa, 'normal'   # implicit bit = 1


def encode(sign: int, exponent_field: int, mantissa_field: int) -> int:
    return (sign \<\< 31) | (exponent_field \<\< 23) | mantissa_field

২. Align — exponent মিলিয়ে mantissa শিফট করা

দুইটা ভিন্ন exponent-এর সংখ্যা সরাসরি যোগ করা যায় না — ছোট exponent-এর mantissa-কে বড়টার exponent-এ “নামিয়ে আনতে” হয়, ডানে শিফট করে। কিন্তু শুধু শিফট করলে নিচের bit-গুলো হারিয়ে যায়, আর সেটা rounding-এর সিদ্ধান্তকে ভুল করে দিতে পারে — তাই ৩টা অতিরিক্ত bit (guard, round, sticky) রেখে শিফট করি।

GRS_BITS = 3   # guard + round + sticky — শিফটে হারানো তথ্যের precision রাখার জন্য

def shift_right_sticky(value: int, shift: int) -> int:
    """value-কে shift বিট ডানে সরানো; ফেলে দেওয়া বিটের যেকোনোটা 1 হলে
       ফলাফলের সবচেয়ে নিচের বিট (sticky bit) 1 করে রাখা — rounding তথ্য না হারানোর জন্য।"""
    if shift \<= 0:
        return value
    dropped = value & ((1 \<\< shift) - 1)
    result = value >> shift
    if dropped:
        result |= 1
    return result


def align_and_add(sign_a, e_a, M_a, sign_b, e_b, M_b):
    Ma = M_a \<\< GRS_BITS
    Mb = M_b \<\< GRS_BITS

    if e_a >= e_b:
        e = e_a
        Ma_aligned = Ma
        Mb_aligned = shift_right_sticky(Mb, e_a - e_b)
    else:
        e = e_b
        Mb_aligned = Mb
        Ma_aligned = shift_right_sticky(Ma, e_b - e_a)

    val_a = Ma_aligned if sign_a == 0 else -Ma_aligned
    val_b = Mb_aligned if sign_b == 0 else -Mb_aligned
    total = val_a + val_b

    result_sign = 0 if total >= 0 else 1
    return result_sign, e, abs(total)

sticky bit — এই ধারণাটাই সবচেয়ে সূক্ষ্ম অংশ। শিফটে ফেলে দেওয়া bit-গুলোর ভেতর একটাও 1 থাকলে, বাকি সব 0 হলেও, ফলাফল mathematically “ঠিক মাঝামাঝি” থেকে সামান্য সরে যায় — সেই তথ্যটুকু একটা মাত্র bit-এ (OR করে) ধরে রাখাই যথেষ্ট সঠিক rounding-এর জন্য।

৩. Renormalize — carry-out ও cancellation সামলানো

যোগফল দুইভাবে “বেঠিক আকারে” আসতে পারে: (ক) দুইটা কাছাকাছি-বড় সংখ্যা যোগ হয়ে এক ঘর উপচে উঠতে পারে (carry-out), অথবা (খ) বিয়োগে দুইটা প্রায়-সমান সংখ্যা একে অপরকে প্রায় বাতিল করে দিতে পারে, ফলে অনেকগুলো leading zero তৈরি হয় (catastrophic cancellation)। দুইটাই renormalize করে ঠিক করতে হয়।

TARGET_BIT = 23 + GRS_BITS     # normalized mantissa-র top bit এখানে থাকা উচিত (=26)
MIN_EXP = -126                  # normal সংখ্যার সবচেয়ে ছোট exponent

def normalize(sign, e, M):
    if M == 0:
        return sign, 0, 0, 'zero'

    while M >> (TARGET_BIT + 1):        # carry-out: এক ঘর উপচে উঠেছে
        M = shift_right_sticky(M, 1)
        e += 1

    while M \< (1 \<\< TARGET_BIT) and e > MIN_EXP:   # leading zero: cancellation
        M \<\<= 1
        e -= 1

    kind = 'normal' if (e > MIN_EXP or M >= (1 \<\< TARGET_BIT)) else 'subnormal'
    return sign, e, M, kind

e > MIN_EXP শর্তটাই subnormal ফলাফলকে সঠিকভাবে ধরে — যদি leading zero সরাতে সরাতে exponent একদম -126-এ ঠেকে যায় কিন্তু mantissa তখনো “পুরো ২৪ bit” আকারে না পৌঁছায়, তাহলে সেটা আসলে একটা বৈধ subnormal ফলাফল, error না। এটাই gradual underflow — subnormal সংখ্যা কেন আছে তার আসল কারণ।

৪. Round-to-nearest-even ও re-encode

def round_and_pack(sign, e, M, kind):
    if kind == 'zero':
        return encode(sign, 0, 0)

    extra = M & ((1 \<\< GRS_BITS) - 1)     # guard/round/sticky bit তিনটা
    M = M >> GRS_BITS

    half = 1 \<\< (GRS_BITS - 1)             # ঠিক মাঝখান (=4)
    if extra > half or (extra == half and (M & 1)):
        M += 1                               # round up — tie হলে even mantissa-র দিকে
        if kind == 'normal' and (M >> 24):
            M >>= 1                          # round করার ফলে আবার carry-out
            e += 1
        elif kind == 'subnormal' and (M >> 23):
            kind = 'normal'                  # round করে subnormal থেকে normal-এ উঠে গেল

    if e > 127:
        return encode(sign, 0xFF, 0)         # overflow → ±∞

    if kind == 'subnormal':
        return encode(sign, 0, M & 0x7FFFFF)

    return encode(sign, e + 127, M & 0x7FFFFF)

extra == half and (M & 1) — এটাই round-to-even নিয়ম। ফলাফল যদি ঠিক দুইটা representable মানের মাঝামাঝি পড়ে (extra ঠিক অর্ধেক), তাহলে যেদিকে round করলে শেষ bit জোড় (even) হয় সেদিকেই যাওয়া হয় — এলোমেলোভাবে সবসময় উপরে বা নিচে round করলে বহু হিসাবে systematic bias তৈরি হতো, IEEE 754 এই bias এড়াতেই এই নিয়ম বেছে নিয়েছে।

৫. Special value handle করে সব জোড়া লাগানো

def make_nan() -> int:
    return 0x7FC00000     # quiet NaN-এর প্রচলিত bit pattern


def add_f32(bits_a: int, bits_b: int) -> int:
    sign_a, e_a, M_a, kind_a = decode(bits_a)
    sign_b, e_b, M_b, kind_b = decode(bits_b)

    if kind_a == 'nan' or kind_b == 'nan':
        return make_nan()
    if kind_a == 'inf' and kind_b == 'inf':
        return make_nan() if sign_a != sign_b else encode(sign_a, 0xFF, 0)
    if kind_a == 'inf':
        return bits_a
    if kind_b == 'inf':
        return bits_b
    if kind_a == 'zero' and kind_b == 'zero':
        return encode(sign_a & sign_b, 0, 0)   # -0 + -0 = -0, বাকি সব +0
    if kind_a == 'zero':
        return bits_b
    if kind_b == 'zero':
        return bits_a

    result_sign, e, M = align_and_add(sign_a, e_a, M_a, sign_b, e_b, M_b)
    result_sign, e, M, kind = normalize(result_sign, e, M)
    return round_and_pack(result_sign, e, M, kind)

৬. হার্ডওয়্যারের সাথে মেলানো

import struct

def f32_bits(value: float) -> int:
    return int.from_bytes(struct.pack('\<f', value), 'little')

def bits_to_f32(bits_pattern: int) -> float:
    return struct.unpack('\<f', bits_pattern.to_bytes(4, 'little'))[0]

def hardware_add_f32(a: float, b: float) -> int:
    """a, b কে আগে float32-এ round করে, তারপর সেই মান দুইটা যোগ করা।
       double-এর mantissa (52 bit) float32-এর (23 bit) চেয়ে অনেক বড়, তাই দুইটা
       float32-representable মান যোগ করলে double-এ কখনো নতুন rounding error ঢোকে না —
       ফলে struct.pack দিয়ে শেষে round করলেই সেটা সঠিকভাবে-rounded আসল float32 যোগফল।"""
    rounded_sum = bits_to_f32(f32_bits(a)) + bits_to_f32(f32_bits(b))
    return f32_bits(rounded_sum)

hardware_add_f32-এর docstring-এর যুক্তিটা গুরুত্বপূর্ণ — এটা কোনো shortcut না, বরং প্রমাণযোগ্য কারণ কেন Python-এর double দিয়েই আসল float32 hardware addition-এর সমতুল্য ফলাফল পাওয়া সম্ভব, কোনো CPU intrinsic ছাড়াই।

৭. Test harness — নিজেকে নিজেই যাচাই

def run_tests():
    FLT_MAX = (2 - 2 ** -23) * 2 ** 127
    FLT_MIN_NORMAL = 2 ** -126
    FLT_MIN_SUBNORMAL = 2 ** -149

    cases = [
        (1.5, 2.25, 'সাধারণ যোগ — exact, কোনো rounding লাগে না'),
        (1.0, -(1.0 - 2 ** -24), 'catastrophic cancellation — বহু bit renormalize'),
        (FLT_MIN_NORMAL, -FLT_MIN_SUBNORMAL, 'normal থেকে subnormal সীমানায় নামা'),
        (FLT_MAX, FLT_MAX, 'overflow — ফলাফল সীমার বাইরে, +∞ হওয়া উচিত'),
        (0.1, 0.2, 'দুইটা মানই ইতিমধ্যে float32-এ rounded, যোগফলও rounded'),
    ]

    for a, b, desc in cases:
        bits_a, bits_b = f32_bits(a), f32_bits(b)
        mine = add_f32(bits_a, bits_b)
        hw = hardware_add_f32(a, b)
        status = 'match ✓' if mine == hw else 'MISMATCH ✗'

        print(desc)
        print(f'  a={a!r}  b={b!r}')
        print(f'  mine: {mine:08X}  ({bits_to_f32(mine)!r})')
        print(f'  hw  : {hw:08X}  ({bits_to_f32(hw)!r})')
        print(f'  {status}\n')


if __name__ == '__main__':
    run_tests()

প্রতিটা টেস্ট কেস কী পরীক্ষা করছে, খেয়াল করুন:

1.0 + -(1.0 - 2⁻²⁴): 1.0 - 2⁻²⁴ হলো float32-তে 1.0-এর ঠিক নিচের representable মান (exponent-এর সীমানায় ULP বদলে যায় বলে এটা exactly representable)। তাই গাণিতিকভাবে ফলাফল হওয়া উচিত ঠিক 2⁻²⁴ — কিন্তু হিসাবটা exponent -1 থেকে শুরু হয়ে renormalize করে -24-এ নামতে হয়, প্রায় পুরো mantissa-র জায়গা জুড়ে শিফট। এখানে normalize()-এর leading-zero loop সবচেয়ে বেশি পরীক্ষিত হয়।
FLT_MIN_NORMAL - FLT_MIN_SUBNORMAL: normal-এর সবচেয়ে ছোট মান থেকে সবচেয়ে ছোট subnormal বিয়োগ করলে ফলাফল হয় ঠিক সবচেয়ে বড় subnormal — normal/subnormal সীমানা পার হওয়ার test।
FLT_MAX + FLT_MAX: FLT_MAX-এর দ্বিগুণ কোনো normal float32-তেই ধরে না (exponent 127-এর বেশি লাগবে, যেটা inf-এর জন্য সংরক্ষিত) — তাই এটা নিশ্চিতভাবে +∞-এ overflow করবে।

নিজেকে চ্যালেঞ্জ করুন

বিয়োগ, গুণ যোগ করুন — গুণ আসলে সহজ: exponent-গুলো যোগ হয়, mantissa-গুলো গুণ হয় (কিন্তু rounding logic প্রায় একই)
Double precision (64-bit) ভার্সন বানান — শুধু bit-width আর bias বদলাবে, algorithm এক
Round-toward-zero আর round-toward-+∞ যোগ করুন, চারটা rounding mode-এর ফলাফল পাশাপাশি দেখান
Fuzz testing — হাজার হাজার random float32 pair generate করে mine == hw কতবার মেলে না দেখুন, mismatch পেলে সেই exact input রিপোর্ট করুন
Double rounding ইচ্ছাকৃতভাবে দেখান — x87-এর মতো 80-bit extended precision-এ প্রথমে round করে তারপর 32-bit-এ round করলে কখনো কখনো সরাসরি 32-bit-এ round করার চেয়ে ভিন্ন ফল আসে, এমন একটা কেস খুঁজে বের করুন

এটা যেখানে গিয়ে মিশবে

এখানে যা শিখলেন	পরে কোথায় লাগবে
Exponent align, mantissa add, renormalize	Level 2 — FPU/floating-point adder সার্কিট ডিজাইন
Round-to-nearest-even bit-level-এ	Level 3 — SIMD floating-point instruction (SSE/AVX), FPU pipeline
Guard/round/sticky bit technique	Level 11 — Numerical performance টিউনিং, error accumulation বোঝা
Catastrophic cancellation দেখা ও measure করা	Level 6 — Numerically stable algorithm ডিজাইন
Hardware বনাম manual implementation cross-check	Level 13 — Formal verification, hardware correctness প্রমাণ

মাইলস্টোন

আগে যা পড়া দরকার