Where Do I Begin?

where do i begin? .... ตอบยากมาก เริ่มไม่ถูกเลยแฮะ พอรู้ตัวอีกทีก็หลงรักเสียงของปู่ไปแล้ว

หลับให้สบายฮะปู่ :'(

P.S. Love means never having to say you're sorry

Einleitung

ปรกติแล้ว เพลงโหมโรงในโลกของ classical music จะค่อนข้างยาว บางเพลงก็เป็น 10 นาทีเลย ซึ่งมันก็ไม่ได้ peak ตลอดเวลาทั้งเพลง มีส่วนที่ฟังสนุกๆ แค่ไม่เท่าไหร่

แต่ Einleitung ใน Also sprach Zarathustra ของ R. Strauss กลับไม่ใช่อย่างนั้น มันเป็นเพลงโหมโรงสั้นๆ แค่ 2 นาทีที่ทรงพลังอย่างมาก

ยิ่งเอาไปใช้คู่กับบทเปิดตัวกับอะไรก็ตามนี่ สุดยอดมาก >w<)b

1000^1000 กับ 1001^999 อะไรมากกว่ากัน?

สืบเนื่องจาก blog ตอนที่แล้ว มีหลายท่านท้วงเข้ามาว่าโจทย์ผิดหรือเปล่า ซึ่งผมได้ไปเช็คกับคำถามอีกรอบแล้วก็พบว่า โจทย์ที่ได้รับมานั้นไม่ผิดจริงๆ แต่ก็น่าสงสัยอย่างมากว่าคนตั้งโจทย์นั่นแหละที่เขียนเลขตกหล่นไปเอง

จริงๆ ตอนก่อนก็แอบเกรียนไปหน่อย คือรู้แหละว่าตั้งใจจะให้พิสูจน์ทางคณิตศาสตร์ งั้นมาไถ่บาปกับคำถามที่ว่า 1000^1000 กับ 1001^999 (เวอร์ชั่นแก้คำผิดแล้ว) อะไรมากกว่ากันอีกรอบดีกว่า ;)

---

ก่อนอื่น สังเกตว่า

1000^1000 = 1000 * 1000^999

ที่เราต้องการคือ จัดรูปของ 1001^999 ให้กลายเป็น 1000^999 ให้ได้ ซึ่งก็ง่ายๆ โดนทำการกระจายทวินาม

1001^999 = (1000 + 1)^999
         = summation   C(999, i) * 1000^(999-i) * 1^i   for i from 0 to 999
         = summation   C(999, i) * 1000^(999-i)         for i from 0 to 999
         = 1000^999 + C(999, 1) * 1000^998 + C(999, 2) * 1000^997 + ... + 1

จะเห็นว่า พอกระจาย 1000^999 ออกมาแล้ว จะมี 1000 พจน์พอดี (index ไล่จาก 0 จนถึง 999) ซึ่งมันสามารถเอาไปจับคู่กับ 1000*1000^999 ได้ ดังนี้

pair(1000^1000, 1001^999) = [ (1000^999, 1000^999),
                              (1000^1 * 1000^998, C(999, 1) * 1000^998),
                              (1000^2 * 1000^997, C(999, 2) * 1000^997),
                              ...
                              (1000^i * 1000^(999-i), C(999, i) * 1000^(999-i)),
                              ...
                              (1000^999 * 1000^0, C(999, 999) * 1000^0) ]

จะเห็นว่า สิ่งที่ต้องการตรวจสอบคือ 1000^i กับ C(999, i) อะไรมากกว่ากัน

จาก

P(n, r) = n!/r!
C(n, r) = n!/(r!(n-r)!) = P(n, r)/(n-r)!

จะเห็นว่า P(n, r) มีค่ามากกว่าหรือเท่ากับ C(n, r) เสมอ

ถึงตรงนี้ก็เห็นชัดแล้วว่า

1000^0 = P(999, 0) = 999!/999! = 1
1000^1 > P(999, 1) = 999!/998! = 999
1000^2 > P(999, 2) = 999!/997! = 999 * 998
1000^3 > P(999, 3) = 999!/997! = 999 * 998 * 997
...
1000^n >= P(999, n) = 999!/(999-n)! = 999 * 998 * ... * (999-n)

ดังนั้น

1000^i >= P(999, i) >= C(999, i)   for all i in [0..999]

เพราะฉะนั้น

for i in [1..999]:
    summation 1000^i * 1000^(999-i) >= summation C(999, i) * 1000^(999-i)
    summation 1000^(i+999-i) >= summation C(999, i) * 1000^(999-i) * 1^i
    1000 * 1000^999 >= (1000 + 1)^999
    1000^1000 >= 1001^999

แต่เนื่องจาก มีบางพจน์ 1000^i ที่มีค่ามากกว่า C(999, i) ดังนั้น

1000^1000 > 1001^999

จบการพิสูจน์ครับ

ยังไหว...

วันก่อนเจอคำถามงี่เง่าว่า 1000^1000 กับ 101^999 อันไหนมากกว่ากัน? แล้วก็บอกว่าห้ามใช้ sense ด้วยนะ

ห้ามใช้ sense ก็กดเครื่ิงคิดเลขไปสิโว้ย

คนตั้งคำถามคงคิดไม่ถึงเนาะ ว่าเครื่องคิดเลขบ้าอะไรจะหาคำตอบของเลขใหญ่ๆ ขนาดนั้นได้

เลยลองเขียนอะไรเล่นๆ เพื่อทดสอบประสิทธิภาพทางการคำนวณตัวเลขดู พบว่าตัวเลขระดับหมื่นหลักก็ยังไหว

Fibonacci ลำดับที่ 57075 เมื่อเขียนในฐานสิบจะต้องใช้ตัวอักษร 11928 ตัว

อันที่จริง ด้วย algorithm ด้านบนนี้ จะเอาไปใช้คำนวณหา Fibonacci ลำดับที่ห้าล้าน ซึ่งถ้าจะเขียนออกมาในเลขฐานสิบ ต้องใช้ตัวอักษรถึง 1044939 ตัวก็ยังไหวนะ (ภายในเวลาประมาณครึ่งนาทีเอง)

ก็อย่าประเมินประสิทธิภาพของคอมพิวเตอร์ต่ำไปละกัน ;)

ทำ Recursion บนการ Integral

ฟังก์ชัน factorial นั้นมีนิยามง่ายๆ ตรงไปตรงมาคือ ผลคูณของจำนวนเต็มบวกตั้งแต่ 1 ไปจนถึง n

factorial n = product [1..n]

ซึ่งถ้าสังเกตดูซักหน่อย จะพบว่ามันสามารถเขียนนิยามเป็น recursion ได้

factorial 1 = 1
factorial n = n * factorial (n - 1)

ตอนนี้เราอาจเพิ่มนิยามที่ว่า 0! = 1 เพื่อความสะดวกในการกระจายพจน์สำหรับคำนวณความน่าจะเป็น

แต่ทั้งหมดนี้จะมีปัญหาอยู่อย่างนึง ตรงที่ทุกอย่างเป็น discrete math หมดเลย นั่นคือเราไม่สามารถหาอะไรอย่าง 0.5! ได้

ถ้าดูจากที่มาและการใช้งานของมัน (ส่วนมากเป็นเรื่องความน่าจะเป็นนั่นแหละ) ก็อาจนับว่าไม่มีปัญหา แต่สำหรับ pure math แล้ว มันก็น่าจะมีอะไรซักอย่างมาตอบคำถามตรงนี้ได้นะ

---

โชคดี (?) ที่เรามีเลขมหัศจรรย์อย่าง e ซึ่งมีสมบัติประหลาดอย่าง

e^x = d/dx e^x

ดูๆ ไปแล้ว มันก็คล้ายกับ fixed-point combinator ที่เคยพูดไว้ในตอนก่อนๆ เพียงแต่เปลี่ยนจาก function call เป็นการ diff-integral แทน โดยมี terminate point อย่าง

integral 1/e^x dx from 0 to infinity = - 1/e^infinity + 1/e^0
                                     = 1

นี่ทำให้เราสามารถเขียน recursion ในรูปของ integral ได้ เช่น

Gamma(n) = integral 1/e^t t^(n-1) dt from 0 to infinity

ลองดูสมบัติของมันโดยทำการ integral เข้าไป จะได้ว่า

Gamma(n) = integral 1/e^t t^(n-1) dt from 0 to infinity
         = integral u dv from 0 to infinity            ; let u = 1/e^t, dv = t^(n-1) dt
         = (limit u*v for x from 0 to infinity) - integral v du from 0 to infinity
         = (1/n)(1/e^infinity infinity^n - 1/e^0 0^n) - integral v du from 0 to infinity
         = 0 - integral v du from 0 to infinity
         = (1/n) integral 1/e^t t^n dt from 0 to infinity
         = Gamma(n+1)/n

ดังนั้น จะเห็นว่า

Gamma(1) = 1
Gamma(n) = (n-1) * Gamma(n-1) = (n-1)(n-2) * Gamma(n-2) = ...

ซึ่งก็คือ

factorial(n) = Gamma(n+1)

---

ถึงตอนนี้ก็คงตอบคำถามได้แล้วว่า 0.5! นั้น สามารถหาค่าได้จาก

Gamma(3/2) = (1/2) * Gamma(1/2)

Gamma(1/2) = integral 1/e^t 1/sqrt(t) dt from 0 to infinity
           = 2 * integral 1/e^u^2 du from 0 to infinity      ; let du = 1/sqrt(t)
           = integral 1/e^u^2 du from -infinity to infinity
           = sqrt(pi)

Gamma(3/2) = sqrt(pi)/2

ตอนพยายามหาค่าของ Gamma(1/2) ต้องใช้ความรู้เรื่อง double integral และการแปลงระนาบเข้าช่วย ลองศึกษาเทคนิคที่นี่

ข้อดีของการเขียนในรูป integral นี่มีอีกอย่าง คือถ้ารู้สึกว่าพิสูจน์ห่าค่าตรงๆ แบบนี้มันยากไป จะเลี่ยงไปใช้ Riemann integral เพื่อประมาณค่าแทนก็ย่อมได้ครับ

Scope และชื่อตัวแปร

JavaScript เป็นตัวอย่างที่ดี (หรือแย่?) ในเรื่อง scope ลองดู

function starline(n) {
    stars = ''
    for (i=0; i<n; i++)
        stars += '*'
    return stars + '\n'
}

function square(n) {
    shape = ''
    for (i=0; i<n; i++)
        shape += starline(n)
    return shape
}

console.log(square(5))

ความคาดหวังของเราคือรูปสี่เหลี่ยมจตุรัสที่สร้างจากดาวขนาด 5x5 แต่ถ้าลองรันโปรแกรมนี้ดู จะเห็นผลลัพท์เป็นรูปดาวที่มีแค่แถวเดียว นั่นเป็นเพราะว่า i ของ square กับ i ของ starline มันคือตัวเดียวกัน ทำให้หลังจากทำงานที่ starline เสร็จแล้ว i ที่ square ก็จะกลายเป็น 5 ไปด้วย ส่งผลให้หลุดออกจาก loop ทันที

ทางแก้สำหรับ JavaScript นี้ก็ไม่ยาก แค่เพิ่ม var ไว้หน้า i เพื่อประกาศว่าเป็นตัวแปรใน scope นี้ก็พอ (อันที่จริงก็ควรจะประกาศตัวแปรทุกตัวที่จะใช้) เพียงเท่านี้ ตัวแปรหนึ่งๆ ก็จะไม่หลุดออกไปนอก scope ที่มันอยู่แล้ว

...แต่คำถามก็คือ มันมีความจำเป็นจริงๆ หรือ ที่เราต้องไล่ประกาศตัวแปรแต่ละตัวในแต่ละ scope เพียงเพื่อไม่ให้มันหลุดไปยัง scope ที่ใหญ่กว่า

ลองดูตัวอย่างแบบเดียวกันใน Python บ้าง

def starline(n):
    stars = ''
    for i in range(n):
        stars += '*'
    return stars + '\n'

def square(n):
    shape = ''
    for i in range(n):
        shape += starline(n)
    return shape

print(square(5))

ความสามารถในการละการประกาศตัวแปรของแต่ละ scope คงไม่ช่วยเรื่อง readability มากเท่าใดนัก แต่มันช่วยให้เขียนโปรแกรมด้วยความกังวลน้อยลง โฟกัสกับเรื่องของ logic มากขึ้น แถมยังไม่ต้องเจอกับ bug ประหลาดอย่างด้านบนนั้นอีก ... อันที่จริง เราสามารถทำอย่างนี้ได้ด้วยซ้ำ

def square(n):
    shape = ''
    for i in range(n):
        for i in range(n):
            shape += '*'
        shape += '\n'
    return shape

print(square(5))

พอลองย้อนกลับไปดูดีๆ บางทีเราก็ไม่จำเป็นต้องเรียกตัวแปรบางตัวมาใช้เสมอไป (แค่ประกาศมันไว้เป็น loop index) แล้วงั้นทำไมไม่ทำให้ตัวแปรนั้นไม่มีชื่อไปซะเลยหละ?

let square n = unlines [['*' | _ <- [1..n]] | _ <- [1..n]]

putStrLn $ square 5

maintenance ง่ายขึ้นเป็นกองเลย (หรือเปล่า? 555)

ผิดพลาด

บ่อยครั้งที่เรามักกลัวในสิ่งใหม่ๆ เพราะนั่นหมายถึงการเข้าไปอยู่ในสภาพแวดล้อมที่ไม่คุ้นเคย ส่งผลให้มีโอกาสทำผิดพลาดสูงขึ้นมากจนทนรับไม่ไหว

แต่บางครั้ง ความผิดพลาดก็นำมาซึ่งสิ่งใหม่ที่สวยงามนะ

อย่างเช่นรูปนี้ที่ตอนแรกตั้งใจว่าจะวาดออกมาเป็นท้องฟ้า (แบบ Google Sky Map) แต่ว่าแปลงสมการออกมาผิดนิดหน่อย เลยเละตุ้มเป๊ะอย่างนี้

ก็อย่าไปกลัวว่าจะทำพลาดนักเลย ถ้าผลลัพท์ที่เลวร้ายที่สุดไม่ได้ไปทำใครตายเข้า ก็หลับหูหลับตาแล้วลุยๆ ไปเถอะ อย่ามัวแต่ลีลานักเลยน่า :3

แก้โจทย์ป้ายสมัครงาน Google แบบบรรทัดเดียว

นอนไม่หลับ บวกกับช่วงนี้แก้โจทย์เล่นเยอะไปหน่อย เหมือนสมองไม่อยากหยุดแก้โจทย์ 555

แล้วก็นึกได้ว่า ไอ้โจทย์ป้ายสมัครงานของ Google อันนี้ ยังไม่เคยแก้เล่นเองเลยนี่หว่า

แต่ถ้าจะให้แก้ตามปรกติธรรมดา ก็ดูจะไม่ท้าทายซักเท่าไหร่ไปแล้ว (เพราะเล่นไปอ่านเฉลยมาเรียบร้อย) แบบนี้เลยต้องเพิ่ม standard ให้กับตัวเองด้วยการเขียนเป็น one-liner โดยไม่พึ่งพา stdlib ซะเลย

ตัวโจทย์คงไม่ต้องแปลแล้ว? (ถ้าตีโจทย์ไม่ออกจริงๆ แนะนำบล็อก @khajochi เลย) เอาเป็นว่าเริ่มกันเลยดีกว่านะ

---

โดย algorithm ที่จะใช้คำนวณหาค่า e คือ

e = 1 + 1/1! + 1/2! + 1/3! + 1/4! + ...

ต้องใช้ algorithm นี้เพราะว่ามันให้ค่าทศนิยมได้แม่นยำรวดเร็วพอสมควร และยัง implement ตามไม่ยากเท่าไหร่ ในที่นี้ใช้ถึงแค่พจน์ 1/200! ก็ให้ค่าออกมาแม่นยำเกือบ 400 ตำแหน่งแล้ว

มาเริ่มที่ส่วนที่เล็กที่สุด ฟังก์ชั่น factorial ที่ดูเผินๆ ก็ไม่มีอะไรประหลาด แต่ถ้าจะเขียนแบบ one-liner ก็เลี่ยงไม่ได้ที่จะต้องเขียนเป็น recursion และต้องไม่จองชื่อ function ไว้ก่อนด้วย ซึ่งก็ทำได้โดยใช้ combinator แบบนี้

factorial = lambda n: (lambda f, x: f(f, x))(lambda r, y: y*r(r, y-1) if y > 1 else 1, n)

(รายละเอียดอ่านได้จากตอนก่อนๆ)

ต่อมาคือการจับแต่ละพจน์มาบวกกัน ซึ่งจะ implement ตรงๆ แบบนี้ไม่ได้

sum(1/factorial(i) for i in range(200))

เพราะ float มันเก็บ precision ของทศนิยมได้เล็กนิดเดียว ดังนั้นเราต้องดัดแปลงสมการต้นแบบให้กลายเป็น

e = (1 + 1/1! + 1/2! + 1/3! + 1/4! + ...)(1!2!3!4!.../1!2!3!4!...)

เทคนิคหนึ่งที่ใช้ได้เสมอๆ เมื่อต้องการทศนิยม n ตำแหน่ง คือคูณด้วย 10 ยกกำลัง n เข้าไปนั่นเอง ดังนั้น

e*10**n = 10**n * (1 + 1/1! + 1/2! + 1/3! + 1/4! + ...)(1!2!3!4!.../1!2!3!4!...)
        = 10**n * (1 + 2!3!4!... + 1!3!4!... + 1!2!4!... + 1!2!3!... + ...)/1!2!3!4!...

จะเห็นว่าถ้าเรียงลำดับ operator ดีๆ เราจะไม่สูญเสีย precision เลย

ได้แนวคิดแล้วก็เริ่ม implement กัน เริ่มโดยหา factorial ของพจน์ต่างๆ ตามปรกติ

fact_each = [factorial(i) for i in range(200)]

ทีนี้เพื่อป้องกันการคำนวณซ้ำหลายรอบ ก็เอา lambda มาห่อหุ้มมันซะ ซึ่งสิ่งที่เราต้องการต่อจากนี้คือผลคูณรวมของค่าทั้งหมด (1!2!3!4!...) โดยเรายังคงต้องใช้ list เดิมนี้อยู่ ดังนั้น

prod_and_fact_each = (lambda l: [reduce(lambda x, y: x*y, l), l])(fact_each)

ถึงตอนนี้ก็ได้เวลาเอาผลคูณรวมของ factorial ทุกตัว มาหารแต่ละตัวใน list นั้นแล้ว (เรายังต้องใช้ผลคูณรวมอีกรอบ ดังนั้นส่งต่อมันด้วย)

prod_and_div_each = (lambda a, l: [a, [a / i for i in l]])(*prod_and_fact_each)

สุดท้ายก็จับทุกตัวมาบวกกัน คูณด้วย 10 ยกกำลัง n แล้วค่อยหารด้วยค่ารวมนั้น

e = str((lambda a, l: 10**400 * sum(l) / a)(*prod_and_div_each))

เท่านี้ก็ได้เลข e ที่มีทศนิยมตามต้องการแล้วนะ <3

---

สำหรับส่วนที่สองจะมาดูที่จำนวนเฉพาะ เนื่องจากเลข 10 หลักค่าที่ใหญ่ที่สุดคือ n = 9999999999 การที่จะบอกว่ามันเป็นจำนวนเฉพาะหรือไม่ ก็ทดสอบด้วยการเอาไปหารจำนวนเฉพาะทุกตัวที่น้อยกว่ามันไปจนถึง sqrt(n) ซึ่งหมายความว่าเราตรวจถึงแค่จำนวนเฉพาะที่น้อยกว่า sqrt(9999999999) = 99999 ก็พอแล้ว

เริ่มด้วยฟังก์ชั่นเพิ่มจำนวนเฉพาะเข้าไป

add_prime = lambda p, t: p + [t] if all(t % q for q in p) else p

ฟังก์ชั่นนี้จะรับตัวแปร 2 ตัวคือ list ของจำนวนเฉพาะก่อนหน้า และเลขที่จะตรวจว่าเป็นจำนวนเฉพาะหรือเปล่า ส่วนค่าที่คืนมานั้น ถ้าเป็นจำนวนเฉพาะจะคืน list ที่เพิ่มจำนวนที่ตรวจเข้าไปด้วย ดังนั้นสิ่งที่เราต้องทำคือใส่ list ว่างๆ เข้าไปก่อน แล้วส่งเลขที่น้อยที่สุดคือ 2 เข้าไปตรวจ หลังจากได้ list ที่มีเลข 2 ก็ส่งเลข 3,4,5,6,...,99999 เข้าไปตรวจทีละตัว โดย list ที่เป็นผลลัพท์แต่ละครั้งที่ได้มาก็จะเก็บไว้ใช้ต่อไปเรื่อยๆ นี่ก็คือแนวคิดของ reduce นั่นเอง ซึ่งก็จะเขียนได้ว่า

p = reduce(add_prime, [[], 2] + range(3, 100000, 2))

ถึงตอนนี้ก็มีจำนวนเฉพาะทั้งหมดที่ต้องการสำหรับงานนี้แล้ว :3

---

สุดท้าย ได้เวลาจับทุกอย่างมารวมกัน เริ่มจากการพยายาม loop ไปยังเลข 10 ตัวในทศนิยมของ e ก่อน ซึ่งก็คือ

[e[i:10+i] for e in [e] for i in range(len(e)-9)]

แต่เนื่องจากเราต้องการ filter ตัวที่ไม่ใช่จำนวนเฉพาะออกไป ซึ่งคราวนี้เราจะเปลี่ยนมาทดสอบด้วยฟังก์ชั่น

test_prime = lambda p, e: all(int(e) % q for q in p)

ดังนั้นก็จะได้

[e[i:10+i] for e, p in [[e, p]] for i in range(len(e)-9) if test_prime(p, e[i:10+i])]

แน่นอนว่างานนี้เราต้องการแค่จำนวนเฉพาะตัวแรกที่เจอ ดังนั้นสามารถเปลี่ยนไปเขียน

next(e[i:10+i] for e, p in [[e, p]] for i in range(len(e)-9) if test_prime(p, e[i:10+i]))

ก็จะได้คำตอบแล้วหละ \(^0^)/

---

สำหรับ code ฉบับเต็มๆ ที่เป็น one-liner (แต่ขึ้นบรรทัดใหม่เพื่อให้อ่านสะดวก) คือ

next(e[i:10+i] for e, p in [[
        str((lambda a, l: 10**400 * sum(l) / a)(*
            (lambda a, l: [a, [a / i for i in l]])(*
                (lambda l: [reduce(lambda x, y: x*y, l), l])(
                    [(lambda f, x: f(f, x))(
                           lambda r, y: y * r(r, y-1) if y > 1 else 1, i)
                        for i in range(200)])))),
        reduce(
            lambda ps, t: ps + [t] if all(t % p for p in ps) else ps,
            [[], 2] + range(3, 100000, 2))]]
    for i in range(len(e)-9)
        if all(int(e[i:10+i]) % q for q in p))

happy coding นะครับ ^^