ทุกสิ่งที่ควรรู้เกี่ยวกับอัลกอริทึม Fast Inverse Square Root

 (github.com/francisrstokes)
5 คะแนน โดย GN⁺ 2024-06-03 | 2 ความคิดเห็น | แชร์ทาง WhatsApp
  • fast inverse square root ซึ่งโด่งดังจาก Quake 3 คือวิธีแก้ปัญหาด้านประสิทธิภาพในยุคนั้น ที่ประมาณค่า 1 / sqrt(x) อย่างรวดเร็วด้วยการตีความบิตของ float ใหม่และการปรับแก้แบบ Newton-Raphson
  • แก่นสำคัญคือรูปแบบบิตแบบจำนวนเต็มของ IEEE-754 32-bit float สามารถถูกมองเหมือนเป็น ค่าประมาณ log2(x) ที่ถูกสเกลและเลื่อนตำแหน่ง ได้
  • 0x5f3759df - (i >> 1) คือรูปแบบที่ย้าย log2(x^-0.5) = -0.5 * log2(x) มาเป็นการเลื่อนบิตและการลบแบบจำนวนเต็ม ส่วนค่าคงที่มหัศจรรย์มาจาก 3/2 * 2^23 * (127 - σ)
  • จากนั้นใช้ การปรับแก้แบบ Newton-Raphson หนึ่งครั้งด้วย y = y * (1.5 - 0.5x * y * y) และในโค้ด Quake การวนซ้ำครั้งที่สองถูกคอมเมนต์ไว้
  • ในปี 1999 inverse square root จำเป็นต้องใช้หลายร้อยถึงหลายพันครั้งต่อวินาทีสำหรับการจัดแสงและการทำ normalization ของเวกเตอร์ 3D แต่บนฮาร์ดแวร์สมัยใหม่ ประโยชน์เชิงปฏิบัติของทริกเดียวกันนี้ลดลงเพราะมีหน่วยประมวลผล floating point เฉพาะทาง

โค้ด Quake ทำอะไร

float Q_rsqrt(float number) {
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = *(long*)&y;
  i  = 0x5f3759df - ( i >> 1 );
  y  = *(float*)&i;
  y  = y * ( threehalfs - ( x2 * y * y ) );

  return y;
}
  • ฟังก์ชันนี้คำนวณค่าประมาณของ inverse square root 1 / sqrt(number) สำหรับ number
  • ส่วนที่โด่งดังที่สุดคือการตีความค่า float เหมือนเป็น long แล้วทำการจัดการบิต 0x5f3759df - (i >> 1)
  • ในปี 1999 ที่ Quake 3 เปิดตัว inverse square root เป็นการคำนวณที่ช้าและมีต้นทุนสูง และจำเป็นต้องใช้หลายร้อยถึงหลายพันครั้งต่อวินาทีในการคำนวณเวกเตอร์ 3D ที่ต้องใช้สมการแสงและการทำ normalization
  • บนฮาร์ดแวร์สมัยใหม่ การคำนวณแบบนี้อาจไม่ได้ทำบน CPU หรือแม้จะรันบน CPU ก็เร็วขึ้นด้วยฮาร์ดแวร์ floating point เฉพาะทางที่พัฒนามาแล้ว

การแทนค่า IEEE-754 32-bit float

  • 32-bit float ประกอบด้วยสามส่วน
    • Sign: 1 บิต แสดงว่าเป็นค่าบวกหรือลบ
    • Exponent: 8 บิต กำหนดช่วงที่ค่าตกอยู่
    • Mantissa: 23 บิต แสดงตำแหน่งภายในช่วงนั้นแบบเชิงเส้น
  • ค่าทั่วไปถูกตีความในรูปแบบต่อไปนี้
N = (-1)^S * 2^(E - 127) * (1 + M / 2^23)
  • B = 127 คือค่า bias ที่ใช้กับ biased exponent และเลขชี้กำลังจริงคือ e = E - B
  • mantissa ไม่ได้ใช้โดยการคูณ m ตรง ๆ แต่ใช้ในรูป 1 + m
    • ถ้า m = 0 จะเป็น 2^e
    • เมื่อ m เข้าใกล้ 1 จะสามารถแทนค่าได้จนเกือบถึงช่วง exponent ถัดไปคือ 2^(e+1)
  • ถ้า exponent เป็น 0 ทั้งหมด จะเป็นตัวเลข sub-normal และสูตรจะเปลี่ยนไป
N = (-1)^S * 2^-126 * m
  • sub-normal จำเป็นสำหรับการแทนค่า 0 และค่าขนาดเล็กที่ใกล้ 0 มาก ๆ
  • ถ้า exponent เป็น 1 ทั้งหมด จะถูกจัดการเป็นค่าพิเศษ
    • ถ้า E = 255, M = 0 จะเป็น Infinity หรือ -Infinity
    • ถ้า M != 0 จะเป็น NaN

ความสัมพันธ์แบบลอการิทึมเมื่อมองบิตของ float เป็นจำนวนเต็ม

  • หากมองการแทนค่าภายในของ float เหมือนเป็นจำนวนเต็ม 32 บิต จะเขียนได้เป็นสูตรต่อไปนี้
I_x = 2^31 S + 2^23 E + M
  • inverse square root ใช้กับอินพุตที่เป็นบวก ดังนั้นถ้ากำหนด S = 0 สูตรจะง่ายขึ้น
L = 2^23
I_x = L E + M
  • ภายในช่วง exponent เดียวกัน mantissa แสดงตำแหน่งแบบเชิงเส้น แต่เมื่อ exponent เพิ่มขึ้น จำนวนขั้นของ mantissa ที่เท่ากันจะครอบคลุมช่วงบนเส้นจำนวนที่กว้างกว่า
    • E = 127 หรือ e = 0 อยู่ประมาณช่วง [1, 2)
    • E = 128 หรือ e = 1 อยู่ประมาณช่วง [2, 4)
    • ทั้งสองช่วงมีจำนวนขั้น mantissa เท่ากัน แต่ช่วงที่สองกว้างเป็นสองเท่า
  • เพราะโครงสร้างนี้ เมื่อมองรูปแบบบิตดิบของ float เป็นจำนวนเต็ม จึงเกิด ความสัมพันธ์แบบลอการิทึม ขึ้น

บิตดิบคือค่าประมาณของ log2(x)

  • หากตีความรูปแบบบิตของ float เป็นจำนวนเต็ม I_x จะมองได้เหมือนเป็นค่าประมาณเชิงเส้นแบบแบ่งช่วงของ log2(x)
  • ความสัมพันธ์นี้เขียนเป็นสูตรประมาณได้ดังนี้
log2(x) ≈ I_x / L - B
  • เมื่อนำจำนวนเต็มของบิตดิบไปหารด้วยขนาด mantissa L = 2^23 แล้วลบ exponent bias B = 127 จะได้ค่าที่ใกล้กับ log2(x)
  • ลอการิทึมภายในช่วง mantissa ถูกจัดการเป็นค่าประมาณเชิงเส้น
log2(1 + x) ≈ x + σ
  • σ คือ พารามิเตอร์ปรับจูน สำหรับปรับค่าประมาณ และ x แสดงตำแหน่งภายในช่วง exponent ในช่วง [0, 1]

เปลี่ยน inverse square root เป็นเอกลักษณ์ลอการิทึม

  • เป้าหมายคือการหาค่าต่อไปนี้
y = 1 / sqrt(x)
  • สามารถเปลี่ยนเป็นรูปเลขชี้กำลังได้ดังนี้
y = x^-0.5
  • เมื่อใช้เอกลักษณ์ลอการิทึม การคำนวณ inverse square root จะกลายเป็นความสัมพันธ์ต่อไปนี้
log2(1 / sqrt(x)) = log2(x^-0.5) = -0.5 * log2(x)
  • โดยใช้ข้อเท็จจริงที่ว่าบิตของ float ทำงานเหมือนค่าประมาณของ log2(x) จึงสามารถประมาณรูปแบบบิตแบบจำนวนเต็ม I_y ของ y ได้โดยตรงจากรูปแบบบิตแบบจำนวนเต็ม I_x ของ x
I_y ≈ -0.5 I_x + 1.5 L (B - σ)
  • สูตรนี้นำไปสู่บรรทัดหลักของโค้ด Quake
i = 0x5f3759df - ( i >> 1 );
  • i >> 1 ทำหน้าที่เลื่อนบิตจำนวนเต็มไปทางขวา 1 บิต ซึ่งเทียบได้กับการคูณด้วย 1/2
  • ค่าคงที่ด้านหน้า 0x5f3759df สอดคล้องกับ 1.5 * L * (B - σ)

ตัวตนของค่าคงที่ 0x5f3759df

  • ถ้ากำหนด σ = 0 ค่าคงที่จะคำนวณได้ดังนี้
1.5 * 2^23 * 127 = 1598029824
  • ค่าในรูปเลขฐานสิบหกคือ 0x5f400000
  • ค่านี้ต่างจากค่าคงที่จริงของ Quake คือ 0x5f3759df อยู่ 566817
  • จากความต่างนี้ สามารถคำนวณค่า σ ที่สอดคล้องกับโค้ด Quake ได้ดังนี้
σ = 377878 / 2^23
σ = 0.04504656
  • ใน C สามารถคำนวณค่าคงที่เดียวกันได้ดังนี้
int32_t compute_magic(void) {
  double sigma = 0.0450465;
  double expression = 1.5 * pow(2.0, 23.0) * (127.0 - sigma);
  int32_t i = expression;
  return i;
}

// -> 0x5f3759df
  • ในที่นี้ใช้ double และการแปลงเป็นจำนวนเต็มเป็นการ cast ตามปกติ ไม่ใช่การตีความบิตใหม่
  • ค่า σ นี้ถูกเลือกเพื่อปรับปรุงค่าประมาณ แต่ไม่ใช่ค่าที่เหมาะที่สุดจริง ๆ และยังไม่แน่ชัดว่าใครเป็นผู้สร้าง

เหตุผลที่ไม่ใช่แค่แฮ็กง่าย ๆ

  • 0x5f3759df - (i >> 1) คือสูตรที่ใช้ข้อเท็จจริงว่าบิตดิบของ float เป็นค่าประมาณลอการิทึม เพื่อสร้าง ค่าเริ่มต้นของ inverse square root
  • แม้จะอิงกับความสัมพันธ์ทางคณิตศาสตร์ที่ซับซ้อน แต่ในขั้นตอนการรันจริงใช้เพียงการดำเนินการที่เร็ว เช่น การเลื่อนบิตและการลบ
  • ในยุคนั้นต้องประมวลผลการคำนวณราคาแพงหลายพันครั้งต่อวินาที วิธีนี้จึงเป็นการออกแบบเชิงวิศวกรรมที่เหมาะกับข้อจำกัดของฮาร์ดแวร์
  • อย่างไรก็ตาม อัลกอริทึมนี้ทำงานเฉพาะกับ normal float เท่านั้น
    • ค่า sub-normal ไม่เป็นไปตามสมมติฐานของค่าประมาณ log2(1 + x) ≈ x + σ
    • ใน sub-normal จะมีรูปแบบที่ใกล้กับ 0 + x จริง ๆ ทำให้ค่าประมาณพังลง

ลดความคลาดเคลื่อนด้วยการปรับแก้แบบ Newton-Raphson

  • ค่าเริ่มต้นที่ได้จากการจัดการบิตค่อนข้างดี แต่ยังมีความคลาดเคลื่อนที่วัดได้อยู่
  • บรรทัดต่อไปนี้ช่วยปรับปรุงค่าประมาณได้อย่างมาก
y = y * ( threehalfs - ( x2 * y * y ) );
  • บรรทัดนี้เป็นรูปแบบที่ใช้ Newton-Raphson method
  • เพื่อปรับปัญหา inverse square root ให้เข้ากับ Newton method จึงเปลี่ยนเป็นปัญหาการหารากของฟังก์ชันต่อไปนี้
f(y) = 1 / y^2 - x = 0
  • Newton method จะสร้างค่าประมาณที่ดีกว่า y_(n+1) จากค่าประมาณปัจจุบัน y_n ดังนี้
y_(n+1) = y_n - f(y_n) / f'(y_n)
  • อนุพันธ์ของ f(y) = y^-2 - x คือ
f'(y) = -2y^-3 = -2 / y^3

สูตรปรับแก้ Newton ที่ไม่ใช้การหาร

  • หากใช้สูตร Newton ตรง ๆ จะมีการหาร floating point หลายครั้ง
  • หนึ่งในเหตุผลที่อัลกอริทึมนี้เร็วคือ การหลีกเลี่ยงการหาร floating point
  • เมื่อจัดรูปทางพีชคณิต จะได้รูปที่ใช้เฉพาะการคูณโดยไม่ต้องหาร
y_(n+1) = y_n * (1.5 - 0.5x * y_n^2)
  • ในโค้ด Quake มีการคำนวณ 0.5x ไว้ล่วงหน้าด้วย x2 = number * 0.5F และใช้ในบรรทัดต่อไปนี้
y = y * ( threehalfs - ( x2 * y * y ) );
  • หลังจากวนซ้ำหนึ่งครั้ง ค่าคลาดเคลื่อนสัมบูรณ์สูงสุดคือ 0.175% และในหลายกรณีค่าคลาดเคลื่อนต่ำกว่านั้น
  • ในโค้ดต้นฉบับมีการวนซ้ำ Newton ครั้งที่สอง แต่ถูกคอมเมนต์ไว้
// y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

ที่มาและอัลกอริทึมที่เกี่ยวข้อง

  • อัลกอริทึมนี้ไม่ได้ถูกคิดค้นโดย John Carmack และที่มาที่แน่ชัดก็ไม่มั่นใจได้ 100%
  • มีบทความของ Beyond3D ที่เกี่ยวข้อง: The truth is the exact origin is not 100% certain
  • Chris Lomont เขียนบทความวิชาการเพื่อหา ค่า sigma ที่เหมาะที่สุด ในขั้นตอนการประมาณลอการิทึม: InvSqrt.pdf
  • CORDIC เป็นอัลกอริทึมที่คำนวณ sine และ cosine ด้วยการบวกและการเลื่อนบิตเท่านั้นโดยไม่ใช้ floating point และมีรายละเอียดวิธีการต่างจาก fast inverse square root อย่างมาก
  • อัลกอริทึมทั้งสองมีจุดร่วมคือการนำข้อสังเกตทางคณิตศาสตร์มาประยุกต์ใช้อย่างมีประสิทธิภาพให้เข้ากับข้อจำกัดของฮาร์ดแวร์ในยุคนั้น

บทความที่เกี่ยวข้อง

2 ความคิดเห็น

 
joyfui 2024-06-03

โค้ดชวนทึ่งที่โผล่ขึ้นมาทุกครั้งพอเริ่มจะลืม..ฮ่า
 
GN⁺ 2024-06-03
ความคิดเห็นจาก Hacker News

  • คอมพิวเตอร์ที่ผลิตหลังปี 1999 โดยทั่วไปจะรองรับ ชุดคำสั่ง SSE และในนั้นมี _mm_rsqrt_ps ที่คำนวณส่วนกลับของรากที่สองได้ 4 ค่าในครั้งเดียวเร็วกว่า: https://www.intel.com/content/www/us/en/docs/intrinsics-guid...
    ถึงอย่างนั้น เทคนิคที่พูดถึงในที่นี้ก็ยังไม่ได้ไร้ความหมายไปเสียทีเดียว การแปลง float/int นั้นเร็ว แต่ยังมีฮาร์ดแวร์ที่ไม่มีคำสั่ง rsqrt, sqrt, pow, log อยู่ และการคำนวณเหล่านี้สามารถประมาณค่าได้ด้วยทริกนี้

    • คำสั่งส่วนกลับแบบ floating-point ของ SSE อาจให้ ผลลัพธ์ต่างกันเล็กน้อยระหว่าง Intel กับ AMD ดังนั้นถ้าคาดหวังผลลัพธ์แบบ deterministic ระหว่างพีซีหลายเครื่อง ก็อาจปวดหัวได้: https://robert.ocallahan.org/2021/09/rr-trace-portability-di...
    • น่าสนใจที่ SSE มีคำสั่งรากที่สองปกติด้วย แต่ ช้ากว่าส่วนกลับของรากที่สองมาก ดังนั้นถ้ายอมเสียความแม่นยำได้ การคำนวณ sqrt(x) เป็น x * 1/sqrt(x) จะเร็วกว่า
    • จริง ๆ แล้วคอมพิวเตอร์ส่วนใหญ่อย่างท่วมท้น ไม่ได้รองรับแม้แต่ SSE หรือชุดคำสั่ง i386/amd64 เอง และสัดส่วนที่รองรับแบบไม่ใช่การจำลองก็ลดลงเรื่อย ๆ
      ในชุดคำสั่ง GPU, ARM, RISC-V, AVR, PIC, 8051, FPGA ฯลฯ มักมีการดำเนินการประมาณค่าส่วนกลับของรากที่สองในตัว แต่ก็น่าจะเป็นไปได้สูงว่าใช้อัลกอริทึมทำนองนี้ในการ implement

  • ถ้าจะจับผิดบทความนิดหน่อย คำอธิบายที่สื่อว่าการคำนวณแบบนี้ไม่ได้เกิดขึ้นบน CPU ยุคนี้แล้วนั้นไม่ถูกต้อง เป็นความเข้าใจผิดที่พบบ่อยว่าเกมหรือแอปที่มีการคำนวณ floating-point หนัก ๆ อยากส่งการคำนวณ floating-point ทั้งหมดไปให้ GPU
    ในความเป็นจริง งานที่เหมาะจะส่งไป GPU คือเฉพาะงานขนาดใหญ่และสม่ำเสมอเท่านั้น ถ้าเป็น การ normalize เวกเตอร์แบบครั้งเดียว เช่น การสร้างเมทริกซ์หมุนเพื่อให้อ็อบเจกต์หนึ่งหันไปมองอีกอ็อบเจกต์หนึ่ง ปล่อยไว้บน CPU จะเร็วกว่า ต่อให้ไม่นับเวลาส่งข้อมูลไป GPU การคำนวณ floating-point เดี่ยว ๆ ก็ยังเร็วกว่าเมื่อทำบน CPU เพราะ GPU มักมี clock ต่ำกว่า และได้จำนวน FLOP สูงจาก parallelism

    • ตรงนี้น่าจะหมายถึง FPU ไม่ใช่ GPU เมื่อก่อน FPU คำนวณแบบ asynchronous และตอนนี้ถือเป็นส่วนหนึ่งที่รวมอยู่ใน CPU

  • ลองเขียน implementation สำหรับ MMIX แล้ว และตั้งสมมติฐานว่า input เดิมมีค่ามากกว่า 2^-1021

  • ถ้าสนใจ Wikipedia ก็มีคำอธิบายที่ดีเกี่ยวกับฟังก์ชันนี้และประวัติของมัน: https://en.wikipedia.org/wiki/Fast_inverse_square_root

  • ผมรวบรวมของพวกนี้ไว้หลายอย่าง: https://github.com/ncruces/fastmath/blob/main/fast.go
    มีโพสต์ StackOverflow ที่เกี่ยวข้องด้วย: https://stackoverflow.com/questions/32042673/optimized-low-a...

    • มีประโยชน์มาก เพราะผมกำลังคิดจะเริ่มรวบรวมเทคนิคแบบนี้เพื่อกลับไปเขียน 3D engine สไตล์ปลายยุค 80 ที่เคยทำไว้ใหม่
    • อยากเห็น benchmark ของแพ็กเกจ fastmath ด้วย

  • ได้เวลาจับผิดแล้ว มี typo ในสูตร float และควรเป็น (-1)^S ไม่ใช่ -1^S แบบหลังจะได้ -1 เสมอ
    คำอธิบายที่ว่าการตีความ raw bit pattern เป็นการประมาณเชิงเส้นเป็นช่วง ๆ ของลอการิทึมก็ไม่ถูกต้องเช่นกัน เส้นระหว่างจุดข้อมูลในกราฟสีน้ำเงินไม่ได้มีอยู่จริง และบิตไม่สามารถตั้งเป็น 1 ได้แค่ครึ่งเดียว มันใกล้เคียงกับ เวอร์ชันแบบไม่ต่อเนื่อง ของลอการิทึมมากกว่า และจุดข้อมูลที่มีอยู่จริง ซึ่งก็คือจุดที่เส้นสีแดงกับสีน้ำเงินตัดกันนั้น แท้จริงแล้วเท่ากับลอการิทึมที่ถูกสเกลและเลื่อนตำแหน่งพอดี นอกนั้นเป็นบทความที่ดี

    • ผมไม่ค่อยเข้าใจ ลองคิดถึง float ขนาดเล็กมากแบบ 6 บิต โดยมีบิตเครื่องหมาย 1 บิต, exponent 2 บิต, significand 3 บิต ช่วง [010000, 010111] จะบรรจุค่า 2, 2.25, 2.5, 2.75, 3, 3.25, 3.5, 3.75
      แต่ significand ที่ implied จากลอการิทึมฐาน 2 ของตัวเลขเหล่านี้คือ .0000000, .0010101, .0101001, .0111010, .1001010, .1011001, .1100111, .1110100 ตามลำดับ และยกเว้นตัวแรกแล้วไม่เหมือนกับ 001, 010 ฯลฯ ของ float เลย เพราะ float ในช่วง [2,4) มีระยะห่างแบบเชิงเส้น แต่ลอการิทึมที่สอดคล้องกันไม่เป็นแบบนั้น จึงมองได้ว่า float เป็น การประมาณเชิงเส้นเป็นช่วง ๆ ของลอการิทึมตามที่บทความว่าไว้
    • มันเป็น การประมาณเชิงเส้นเป็นช่วง ๆ แบบไม่ต่อเนื่อง ไม่ใช่การประมาณเชิงเส้นเป็นช่วง ๆ แบบต่อเนื่อง ประเด็นที่ว่าเส้นสีน้ำเงินไม่ต่อเนื่องนั้นถูกต้อง แต่การตีความนั้นผิด กราฟสีน้ำเงินประกอบด้วยจุดเดี่ยว ๆ 256 จุดที่วางเรียงอย่างสม่ำเสมอตามแกน x ไม่ใช่แค่จุดตัดไม่กี่จุด
      ถ้าเป็นกราฟเต็ม ๆ ก็จะมีตัวเลือก 2^32 ค่าอยู่ในแพตเทิร์นเชิงเส้นเป็นช่วง ๆ แต่สิ่งที่บทความต้นฉบับวาดไม่ใช่กราฟเต็มแบบนั้น ในเมื่อบทความพูดถึงจำนวนเต็ม 32 บิตและการคำนวณ IEEE-754 32-bit float ผมคิดว่าการละคำว่า “ไม่ต่อเนื่อง” ในคำอธิบายก็พอรับได้

  • เป็นบทความที่ดีและอธิบายแนวคิดน่าสนใจไว้มาก แต่ การแปลงพีชคณิต ในส่วนหนึ่งแย่อย่างน่าตกใจ
    หลังข้อความ “ขั้นตอนที่แน่นอนจากรูปแรกไปยังรูปนี้มีหลายขั้น แต่เพื่อความครบถ้วนจึงใส่ไว้ทั้งหมด” การแปลงมีขั้นตอนที่ไม่จำเป็นเยอะ และมีข้อผิดพลาดเรื่องเครื่องหมายที่หักล้างกันเองหลายจุด โดยเฉพาะตอนจากบรรทัดที่สองไปบรรทัดที่สาม เครื่องหมายลบไม่ได้ถูกแจกแจงอย่างถูกต้อง ถ้าเริ่มหลังบรรทัดที่สอง ก็สามารถเริ่มจาก y_n+1 = y_n + (1 - x * y_n^2) / y_n^2 * (y_n^3 / 2) แล้วไปถึง y_n+1 = y_n (1.5 * y_n - 0.5 * x * y_n * y_n) ได้สั้นกว่ามาก และขั้นตอนกลางก็ถูกต้องด้วย ผมคิดว่ามีแต่ขั้นตอนที่ชัดเจนอยู่แล้วสำหรับคนที่เข้าใจพีชคณิต

  • magic number ในโค้ดสั้น ๆ อันโด่งดังนั้นไม่ใช่ค่าคงที่ที่เหมาะที่สุด ถ้าใช้ค่าคงที่อื่น น่าจะลด relative error ได้อีกราว 0.5%
    ตอนนั้นการหาค่าที่เหมาะที่สุดแบบสัมบูรณ์อาจยาก แต่ตอนนี้ค่อนข้างง่าย ผมเองก็เคยหลงลง rabbit hole นี้อยู่พักหนึ่ง และมี Jupyter notebook สำหรับหา magic number ที่เหมาะที่สุดของ (1/x^2) กับ (1/x)

    • มีลิงก์บทความวิจัยที่สำรวจคำถามนั้นอยู่ด้านล่างของบทความ

  • สิ่งที่น่าสนใจที่สุดสำหรับผมในบทความนี้คือลิงก์ “How Java's Floating-Point Hurts Everyone Everywhere”: https://people.eecs.berkeley.edu/~wkahan/JAVAhurt.pdf
    ผู้เขียนคือ William Kahan ซึ่งเป็นที่รู้จักในชื่อ “Old Man of Floating-Point” ด้วย: https://news.ycombinator.com/item?id=29042853 - An Interview with the Old Man of Floating-Point (1998)

    • นอกเรื่องจากหัวข้อ แต่ผมเริ่มอ่าน PDF JAVAhurt แล้วพบว่า การจัดตัวอักษรแย่มาก เหมือนใช้แพ็กเกจ TeX ที่ถ่างช่องว่างระหว่างคำมากเกินไป แถมไม่สม่ำเสมอ หรือเหมือน OCR เอกสารอื่นแล้วมีช่องว่างเกินติดมา
      แม้แต่ส่วนที่เป็นฟอนต์ monospace ก็มีช่องว่างแปลก ๆ เพิ่มเข้ามา อ่านแล้วมีสมาธิยากมาก และถึงจะรู้ว่าไม่ใช่จริง ๆ แต่แทบให้ความรู้สึกเหมือนแถลงการณ์ของพวก geek สายวิทย์เลย

  • วิดีโอนี้ที่เคยดูเมื่อก่อนดีมากจริง ๆ: https://www.youtube.com/watch?v=p8u_k2LIZyo