จำนวนที่น้อยที่สุดที่เมื่อบวกกับ 1 แล้วได้ผลลัพธ์ไม่ใช่ 1 (การทำวิทยานิพนธ์ภาคปฏิบัติตอนที่ ๙๔) MO Memoir : Tuesday 27 November 2561

เวลาที่คอมพิวเตอร์เก็บค่าตัวเลขในรูปแบบเลขจำนวนจริง (Real) มันก็จะเก็บในทำนอง ±0.ddddd.. x 10^±yyy (อันที่จริงมันเก็บในรูปแบบเลขฐาน 2 นะ แต่ขอเขียนในรูปแบบเลขฐาน 10 เพื่อให้มองเห็นภาพได้ง่าย) เมื่อ ddddd... คือตัวเลขที่อยู่หลังจุดทศนิยม ส่วน yyy ก็เป็นส่วนของเลขยกกำลัง และก็มีบิทที่เป็นตัวกำหนดเครื่องหมายของส่วน ddddd... ว่าจะให้เป็นบวกหรือลบ และส่วนของ yyy ว่าจะให้เป็นยกกำลังที่เป็นบวกหรือลบ ส่วนตัวเลขที่อยู่หลังจุดทศนิยมจะมีได้กี่ตำแหน่ง และ yyy จะมีค่าได้มากที่สุดเท่าใดนั้นก็ขึ้นอยู่กับจำนวนบิท (bit) ที่ใช้ในการเก็บข้อมูล ถ้าใช้บิทมากก็จะเก็บได้ละเอียดมาก แต่นั่นก็หมายถึงการสิ้นเปลืองหน่วยความจำมากขึ้น และการคำนวณที่กินระยะเวลายาวนานขึ้นด้วย
 

เวลาที่เราให้มันเก็บตัวเลขเช่น 12.34 มันก็จะเก็บเป็น 0.12340000 x 10² และถ้าให้เก็บตัวเลข 0.001234 มันก็จะเก็บค่าเป็น 0.12340000 x 10^-2 กล่าวคือมันจะทำการปรับในส่วนของตัวเลขยกกำลัง เพื่อให้ตัวเลขตัวแรกที่อยู่หลังจุดทศนิยมนั้นเป็นตัวเลขที่ไม่ใช่ศูนย์ ส่วนตัวเลขที่อยู่หลังเลข 4 ก็จะมีค่าเป็นศูนย์ไปทั้งหมด ส่วนจะมีกี่ตัวนั้นก็ขึ้นอยู่กับว่ามันเก็บทิศนิยมได้กี่ตำแหน่ง อย่างเช่นในที่นี้สมมุติว่ามันเก็บได้แค่ 8 ตำแหน่ง ตัวเลขก็จะเป็นดังที่แสดง (คือมี 0 ต่อท้ายอีก 4 ตัว)
 

เวลาที่ทำการบวกเลขสองจำนวนเข้าด้วยกัน สิ่งที่คอมพิวเตอร์ทำก็คือการปรับให้ส่วนที่เป็นเลขยกกำลังของตัวเลขทั้งสองจำนวนนั้นเท่ากันก่อน ด้วยการเขยิบตัวเลขส่วนทศนิยมนั้นมาข้างหลัง (หรือเพิ่ม 0 ไปข้างหน้า) และเพิ่มตัวเลขในส่วนที่เป็นเลขยกกำลังให้สูงขึ้น (1 หน่วยทุก ๆ การเขยิบถอยหลัง 1 หลัก) อย่างเช่นถ้าเราต้องการบวก 12.34 เข้ากับ 0.001234 มันก็จะปรับ 12.34 เป็น 0.12340000 x 10² และปรับ 0.001234 เป็น 0.00001234 x 10² จากนั้นจึงบวกส่วนที่เป็นเลขทศนิยมเข้าด้วยกัน ผลที่ได้คือ 0.12341234 x 10² 
 

ทีนี้ถ้าเราต้องการบวกตัวเลข 0.00001234 เข้ากับ 12.34 บ้าง เลข 0.00001234 ถูกเก็บเป็น 0.12340000 x 10^-4 จากนั้นก็ทำการปรับให้ส่วนที่เป็นเลขยกกำลังเท่ากัน ในคราวนี้ 0.00001234 จะกลายเป็น 0.00000012 x 10^-4 คือเลข 3 และ 4 สองตัวท้ายจะหายไปอันเป็นผลจากการที่มันเก็บทศนิยมไว้ได้เพียงแค่ 8 ตำแหน่ง ดังนั้นผลรวมที่ได้จึงกลายเป็น 0.12340012 x 10² ซึ่งมีความคลาดเคลื่อนไปจากคำตอบที่ควรเป็น
 

(หมายเหตุ : ในที่นี้เวลากล่าวถึงการ "บวก" ก็จะครอบคลุมถึงการ "ลบ" ด้วย และในทำนองเดียวกันเวลากล่าวถึงการ "คูณ" ก็จะครอบคลุมถึงการ "หาร" ด้วย)
 

แล้วปัญหาเรื่องนี้มันส่งผลต่อการคำนวณมากน้อยแค่ไหน ตรงนี้ก็ขึ้นอยู่กับว่าเราเขียนลำดับคำสั่งการคำนวณอย่างไร ตัวเช่นการคำนวณค่า 1 + e - 1, 1 - 1 + e และ 1 - (1 - e) เมื่อ e มีค่าต่าง ๆ กัน ซึ่งในทางทฤษฎีแล้วการเขียนทั้งสามรูปแบบควรให้คำตอบเดียวกัน แต่พอลองคำนวณด้วยโปรแกรม spread sheet ของ OpenOffice 4.1.5 กลับได้ผลดังแสดงในตารางที่ ๑ ซึ่งจะเห็นได้ว่าตั้งแต่ e มีค่าจาก 0.001 ไปจนถึง 10^-14 นั้น ทั้งสามรูปแบบให้คำตอบที่แตกต่างกัน และเมื่อ e มีค่าตั้งแต่ 10^-15 หรือต่ำลงไป การเขียนในรูปแบบ 1 + e - 1 และ 1 - (1 - e) ให้คำตอบที่เป็น 0 ในขณะที่การเชียนแบบ 1 - 1 + e นั้นยังให้คำตอบที่ถูกต้องอยู่
 

(หมายเหตุ : ผลการคำนวณที่ได้นั้นยังขึ้นอยู่กับคอมพิวเตอร์ที่ใช้ว่าเป็นระบบกี่บิท ดังนั้นถ้าลองทำกับเครื่องที่ตนเองใช้อยู่แล้วได้ผลที่ไม่เหมือนที่นำมาแสดงนี้ก็ไม่ใช่เรื่องแปลก แต่ถ้าใช้บิทเยอะ (เช่น 64 บิท) จะเห็นปัญหาเช่นนี้เกิดขึ้นช้ากว่ากรณีของเครื่อง 16 บิทหรือ 32 บิท)
 

ตารางที่ ๑ ผลการคำนวณค่า 1 + e - 1, 1 - 1 + e และ 1 - (1 - e) เมื่อ e มีค่าต่าง ๆ กัน การคำนวณนี้กระทำบนโปรแกรม spread sheet ของ OpenOffice 4.1.5 ในทางทฤษฎีแล้วมันควรจะให้ค่าที่เท่ากัน แต่ในทางปฏิบัติเมื่อทำการคำนวณด้วยเครื่องดิจิตอลคอมพิวเตอร์จะเห็นว่าค่าที่ได้นั้นแตกต่างกันอยู่

e	1 + e - 1	1 - 1 + e	1 - (1 - e)
1.00E-01	0.10000000000000000000	0.10000000000000000000	0.10000000000000000000
1.00E-02	0.01000000000000000000	0.01000000000000000000	0.01000000000000000000
1.00E-03	0.00099999999999989000	0.00100000000000000000	0.00100000000000000000
1.00E-04	0.00009999999999998900	0.00010000000000000000	0.00009999999999998900
1.00E-05	0.00001000000000006550	0.00001000000000000000	0.00000999999999995449
1.00E-06	0.00000099999999991773	0.00000100000000000000	0.00000100000000002876
1.00E-07	0.00000010000000005839	0.00000010000000000000	0.00000009999999994736
1.00E-08	0.00000000999999993923	0.00000001000000000000	0.00000001000000005025
1.00E-09	0.00000000100000008274	0.00000000100000000000	0.00000000099999997172
1.00E-10	0.00000000010000000827	0.00000000010000000000	0.00000000010000000827
1.00E-11	0.00000000001000000083	0.00000000001000000000	0.00000000001000000083
1.00E-12	0.00000000000100008890	0.00000000000100000000	0.00000000000099997788
1.00E-13	0.00000000000009992007	0.00000000000010000000	0.00000000000010003109
1.00E-14	0.00000000000000999201	0.00000000000001000000	0.00000000000000999201
1.00E-15	0.00000000000000000000	0.00000000000000100000	0.00000000000000000000
1.00E-16	0.00000000000000000000	0.00000000000000010000	0.00000000000000000000
1.00E-17	0.00000000000000000000	0.00000000000000001000	0.00000000000000000000
1.00E-18	0.00000000000000000000	0.00000000000000000100	0.00000000000000000000
1.00E-19	0.00000000000000000000	0.00000000000000000010	0.00000000000000000000
1.00E-20	0.00000000000000000000	0.00000000000000000001	0.00000000000000000000

(๑) พึงสังเกตว่า ณ ตำแหน่งค่า e = 0.001 การเขียนแบบ 1 + e - 1 จะให้ผลลัพธ์ที่มีความคลาดเคลื่อนเกิดขึ้น ในขณะที่การเขียนแบบ 1 - (1 - e) จะเกิดปัญหาเดียวกันที่ค่า e = 0.0001

(๒) ณ ตำแหน่งค่า e = 1E-15 การเขียนแบบ 1 + e - 1 และการเขียนแบบ 1 - (1 - e) จะให้ผลลัพธ์ที่เป็น 0 ในขณะที่การเขียนแบบ 1 - 1 + e ยังคงให้คำตอบที่ถูกต้องอยู่

ตัวอย่างที่ยกมานี้ทำการคำนวณโดยใช้ spread sheet สำหรับผู้ที่เขียนโปรแกรมแล้วให้ compiler แปลเป็นไฟล์ EXE แล้วนำไปคำนวณนั้นอาจได้ผลที่แตกต่างกันออกไปได้ เพราะแต่ก่อนก็เคยเจอเหมือนกันสำหรับ compiler ที่ฉลาดบางตัว เวลาที่เขียนคำสั่งทำนอง A = 1 + e - 1 หรือ A = 1 - 1 + e มันจะแปลงเป็น A = e ให้โดยอัตโนมัติ ดังนั้นไม่ว่าจะใช้ค่า e เป็นเท่าใดก็ได้ค่า A เท่ากับ e เสมอ 
  

(compiler ในที่นี้หมายถึงตัวแปลภาษา จากภาษาที่เราเขียนที่เป็นภาษาที่เราพอจะอ่านเข้าใจได้ ให้กลายเป็นภาษาที่คอมพิวเตอร์เข้าใจได้เพื่อมันจะได้ทำงานได้และสามารถสั่ง run ได้ เช่นเราเขียนโปรแกรมด้วยไฟล์ work แต่ตัว compiler จะแปลงเป็นไฟล์ .exe ให้)
 

แต่ก่อนนั้น ก่อนที่จะเริ่มทำ simulation สิ่งหนึ่งที่ต้องทำกันก็คือต้องหาว่าเครื่องคอมพิวเตอร์ที่ใช้นั้นมันเก็บเลขในส่วนที่เป็นเลขทศนิยมได้กี่ตำแหน่ง และส่วนที่เป็นเลขยกกำลังนั้นติดลบได้มากที่สุดเท่าใด จำนวนหลักหลังจุดทศนิยมที่มันเก็บได้เป็นตัวบ่งบอกถึงค่า "Machine precision" (ที่แปลเป็นไทยว่าค่าความเที่ยงของเครื่อง) หรือจำนวนที่น้อยที่สุดที่เมื่อบวกเข้ากับ 1 แล้วได้ผลลัพธ์ที่ไม่ใช่ 1 ส่วนเลขยกกำลังที่ติดลบได้มากที่สุดนั้นเป็นตัวบ่งบอกให้ถึงจำนวนที่น้อยที่สุดที่มากกว่าศูนย์ที่เครื่องนั้นสามารถเก็บค่า ค่าหลังนี้เรียกว่า "Machine accuracy" (ที่แปลเป็นไทยว่าค่าความแม่นของเครื่อง) อย่างเช่นในตัวอย่างที่ยกมาแสดงในตารางที่ ๑ นั้นค่า Machine precision ของเครื่องที่ผมใช้คำนวณจะอยู่ที่ประมาณ 10^-16 หรือแปลได้ว่า แต่จำนวนที่น้อยที่สุดที่ไม่ใช่ศูนย์ที่เครื่องมันรับค่าได้จะอยู่ที่ประมาณ 10^-307 (คือถ้าป้อนค่าน้อยกว่านี้มันจะบันทึกค่าเป็น 0) 
  

จากประสบการณ์ส่วนตัวที่เคยต้องเขียนโปรแกรม simulation ในยุคที่ต้องเขียน source code กันเอง การที่ต้องรู้ค่า Machine percision ของเครื่องคอมพิวเตอร์ที่ใช้คำนวณนั้นเป็นเรื่องสำคัญ (ส่วน Machine accuracy ไม่เคยใช้) เพราะมันส่งผลต่อการคำนวณค่าฟังก์ชันบางฟังก์ชันที่ไม่สามารถหาได้ด้วยเทคนิคการ analyse โดยตรง ตัวอย่างที่เห็นได้ชัดก็คือค่าอนุพันธ์ของฟังก์ชันที่ซับซ้อน การหาค่าอนุพันธ์ของฟังก์ชันเหล่านี้เราต้องกลับไปคำนวณโดยใช้สมการพื้นฐานคือ

ในทางทฤษฎีนั้น ค่าอนุพันธ์ของฟังก์ชัน ณ จุดใดก็เสมือนกับความชันของเส้นสัมผัสโค้ง ณ ตำแหน่งนั้น ความชันของเส้นสัมผัสโค้งประมาณได้จากผลต่างของค่า y และค่า x ที่ 2 ตำแหน่ง และเมื่อระยะห่างระหว่างจุด x 2 จุดนี้แคบลง (∆x วิ่งเข้าหา 0) ค่าความชันที่ได้ก็จะเป็นค่าความชันที่ถูกต้อง ณ ตำแหน่งนั้น 
  

ที่นี้เรามาลองคำนวณค่า df(x)/dx ของฟังก์ชัน f(x) = x² ด้วยการใช้สมการพื้นฐานดูบ้าง โดยเทคนิคการ analyse นั้นฟังก์ชันนี้มีค่า df(x)/dx = 2x ดังนั้นที่ x = 1 จะได้ค่า df(x)/dx = 2 ผลการคำนวณด้วยสมการพื้นฐานที่ค่า ∆x ต่าง ๆ กันแสดงไว้ในตารางที่ ๒ ข้างล่าง (เครื่องนี้มีค่า Machine precision อยู่ที่ระดับ 10^-16)

ตารางที่ ๒ ผลการคำนวณค่า df(x)/dx ของฟังก์ชัน y = x² ด้วยการใช้สมการพื้นฐานที่ค่า ∆x ต่าง ๆ กัน

dx	f(x+dx)	f(x)	df/dx (numerical)
1.00E-01	1.2100000000000000	1.0000000000000000	2.1000000000000000
1.00E-02	1.0201000000000000	1.0000000000000000	2.0100000000000000
1.00E-03	1.0020010000000000	1.0000000000000000	2.0009999999997000
1.00E-04	1.0002000100000000	1.0000000000000000	2.0000999999991700
1.00E-05	1.0000200001000000	1.0000000000000000	2.0000100000139300
1.00E-06	1.0000020000010000	1.0000000000000000	2.0000009999243700
1.00E-07	1.0000002000000100	1.0000000000000000	2.0000001010878100
1.00E-08	1.0000000200000000	1.0000000000000000	1.9999999878450600
1.00E-09	1.0000000020000000	1.0000000000000000	2.0000001654807400
1.00E-10	1.0000000002000000	1.0000000000000000	2.0000001654807400
1.00E-11	1.0000000000200000	1.0000000000000000	2.0000001654807400
1.00E-12	1.0000000000020000	1.0000000000000000	2.0001778011646800
1.00E-13	1.0000000000002000	1.0000000000000000	1.9984014443252800
1.00E-14	1.0000000000000200	1.0000000000000000	1.9984014443252800
1.00E-15	1.0000000000000000	1.0000000000000000	0.0000000000000000

จะเห็นว่าเมื่อลดค่า ∆x ลงเรื่อย ๆ ค่า df(x)/dx ที่คำนวณได้จะเข้าหาค่าที่ถูกต้องมากขึ้น แต่เมื่อ ∆x ลดลงถึงระดับหนึ่งกลับพบว่าค่าที่คำนวณได้นั้นจะผิดเพี้ยนมากขึ้น และมากขึ้นตามค่า ∆x ที่ลดต่ำลง และพอถึงระดับหนึ่งจะพบว่าได้ค่าอนุพันธ์เป็นศูนย์ ซึ่งเป็นค่าที่ผิด สาเหตุเป็นเพราะที่ระดับนี้ค่า x + ∆x จะมีค่าเท่ากับ x ทำให้ f(x + ∆x) = f(x)
 

จุดที่ให้คำตอบที่ถูกต้องมากที่สุดอยู่ที่ประมาณรากที่สองของค่า Machine precision กล่าวคือถ้าเครื่องคอมพิวเตอร์ที่ใช้นั้นคำนวณด้วยความละเอียดขนาด 16 ตำแหน่งหรือค่า Machine precision อยู่ที่ระดับ 10^-16 ถ้า x มีค่าอยู่ที่ระดับ 1 ค่า ∆x ที่ให้ผลการคำนวณที่ถูกต้องที่สุดก็อยู่ที่ประมาณ 10^-8 แต่ต้องย้ำนิดนึงว่าค่านี้เป็นค่าสัมพัทธ์ กล่าวคือถ้า x มีค่าอยู่ที่ระดับ 10³ ค่า ∆x ที่ให้ผลการคำนวณที่ถูกต้องที่สุดก็อยู่ที่ประมาณ 10^-5 และถ้า x มีค่าอยู่ที่ระดับ 10^-3 ค่า ∆x ที่ให้ผลการคำนวณที่ถูกต้องที่สุดก็อยู่ที่ประมาณ 10^-11 หรือพูดง่าย ๆ ก็คือค่า x มีค่าเท่าใด ก็เอาค่า 10^-8 (ซึ่งเป็นค่ารากที่สองของค่า Machine precision ของเครื่องที่ผมใช้คำนวณ) คูณเข้าไป