Showing posts with label Python. Show all posts
Showing posts with label Python. Show all posts

Mar 23, 2016

ใบเฟิร์น Barnsley

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


การสร้าง fractal รูปใบเฟิร์นอย่างง่าย สามารถอธิบายคร่าวๆ ได้ด้วยขั้นตอนวิธีนี้


  1. ลากเส้นตรงความยาว $l$ ไม่ต้องเยอะมากนัก แล้วเลือกปลายข้างหนึ่งไว้เป็นปลายที่จะต่อยอดออกไป
  2. ที่ปลายยอดของเส้นตั้งต้น ทำสองขั้นตอนต่อไปนี้เพื่อลากเส้นเพิ่ม 3 เส้น
    • ลากเส้นตรงความยาว $0.2$ เท่าของความยาวเส้นตั้งต้น ในทิศทางตั้งฉากกับเส้นตั้งต้นทั้งสองข้าง
    • ลากเส้นตรงความยาว $0.8$ เท่าของความยาวเส้นตั้งต้น ในทิศทางต่อขึ้นไปตรงๆ จากเส้นตั้งต้น
  3. สำหรับแต่ละเส้นที่ลากต่อออกมาในข้อ 2 นั้น เลือกปลายด้านที่ว่างๆ ไว้เป็นปลายต่อยอด แล้ววนกลับไปทำขั้นตอนที่ 2 โดยเปลี่ยนมาใช้เส้นตั้งต้นเป็นเส้นที่ได้จากข้อ 2 แทน (คิดซะว่าเส้นใหม่นี้ยาว $l$) ทำไปเรื่อยๆ จนกว่าจะพอใจผลลัพธ์

วิธีการข้างต้นจะให้ผลลัพธ์เป็นรูปใบเฟิร์นเรียบๆ ตรงๆ ไม่หวือหวานัก แต่ก็มีนักคณิตศาสตร์นามว่า Michael Barnsley ได้เขียนสมการกำหนดมุมและจัดตำแหน่งของการต่อยอดให้สวยงามมีลูกเล่น เนื่องจากเขาเป็นคนแรกที่เขียนเรื่องนี้ไว้ตั้งแต่ปี 1993 ทุกคนเลยพากันเรียกใบเฟิร์นที่สวยงามเป็นเอกลักษณ์นี้ว่า ใบเฟิร์น Barnsley



พูดให้รัดกุมเป็นภาษาคณิตศาสตร์ก็คือ การสร้างใบเฟิร์น Barnsley เริ่มจากลากเส้นแรกยาว $1.6$ หน่วยจากจุด $(0,0)$ ขึ้นไปตามแกน $y$ แล้วเอาเส้นนั้นมาวาดซ้ำๆ อีกครั้ง จากเลื่อนและปรับขนาดด้วยสมการดังนี้


f_1(x, y) = [ 0.85, 0.04; -0.04, 0.85 ] [ x; y ] + [ 0; 1.6 ]
f_2(x, y) = [ 0.20 & -0.26 \\ 0.23 & 0.22 ] [ x; y ] + [ 0; 1.6 ]
f_3(x, y) = [ -0.15, 0.28; 0.26 & 0.24 ] [ x; y ] + [ 0; 0.44 ]

โดย $f_1, f_2, f_3$ คือสมการสำหรับคำนวณรูปร่างก้านที่ต่อยอดขึ้นข้างบน ก้านที่เอียงไปด้านซ้าย และก้านที่พลิกตัวแล้วเอียงไปด้านขวา ตามลำดับ


ผมแนบโค้ดในภาษา Python มาให้ด้วยสำหรับใครที่ต้องการนำไปเล่นต่อ ซึ่งวิธีการที่แนบนี้เป็นการสร้างแบบ recursive ต่อยอดออกไปเรื่อยๆ ตามแนวทางขั้นตอนวิธีที่ได้กล่าวไว้ข้างต้น หากไปค้นดูวิธีอื่นเพิ่มเติมจะพบว่าวิธีการสุ่มเลือกจุดจะให้ผลลัพธ์ที่รวดเร็วกว่า


from PIL import Image, ImageDraw

def transform(matrix, xy):
    return [sum(e * p for e, p in zip(row, xy)) for row in matrix]

def translate(vector, xy):
    return [e + p for e, p in zip(vector, xy)]

class BarnsleyFern(object):
    affines = [([[0.85, 0.04], [-0.04, 0.85]], [0, 1.6]),
               ([[0.20, -0.26], [0.23, 0.22]], [0, 1.6]),
               ([[-0.15, 0.28], [0.26, 0.24]], [0, 0.44])]

    def __init__(self, depth=10, size=(800,800)):
        self.image = Image.new('RGB', size)
        self.draw = ImageDraw.Draw(self.image)
        self.iterate([0.0, 0.0], [0.0, 1.6], depth)

    def canvas_coordinate(self, xy):
        width, height = self.image.size
        return [width/2 + xy[0]*width/10, height - xy[1]*height/10]

    def linespec(self, xy0, xy1):
        return self.canvas_coordinate(xy0) + self.canvas_coordinate(xy1)

    def draw_too_small(self, xy0, xy1, pos=2):
        return all(round(a, pos) == round(b, pos) for a, b in zip(xy0, xy1))

    def iterate(self, xy0, xy1, depth, affine=None):
        if depth == 0:
            return
        if affine is not None:
            xy0 = translate(affine[1], transform(affine[0], xy0))
            xy1 = translate(affine[1], transform(affine[0], xy1))
        if self.draw_too_small(xy0, xy1):
            return
        self.draw.line(self.linespec(xy0, xy1), fill='#0C3')
        for affine in self.affines:
            self.iterate(xy0, xy1, depth-1, affine)



fractal ยังมีรูปร่างอื่นๆ อีกมากมาย ทั้งที่พบได้ในธรรมชาติ เช่น ดอกทานตะวัน บรอคโคลี เกล็ดหิมะ หรือมากจากการสังเคราะห์ขึ้นอย่าง สามเหลี่ยม Sierpinski หรือ เซ็ต Mandelbrot


แต่ความมหัศจรรย์ของ fractal ที่แท้จริง คือเราไม่สามารถรู้ได้เลยว่าตอนนี้เรากำลังอยู่ตรงไหนของมันกันแน่ ลองคิดภาพว่าตัวเองเป็น Ant-Man ที่จะปรับขนาดตัวให้ไม่ใหญ่ไปกว่าก้านใบเฟิร์นที่ยืนอยู่ได้ แล้วก็ลองไปยืนอยู่บนใบเฟิร์นข้างต้นดู ณ ขณะหนึ่งเราอาจจะคิดว่าตัวเองอยู่บนก้านที่ใหญ่ที่สุดแล้ว แต่นอกจากขนาดที่แตกต่างกัน ก้านที่เราคิดว่าใหญ่ที่สุดนั้น ก็ไม่ได้มีอะไรที่แตกต่างจากก้านย่อยก้านอื่นๆ เลย หากเราเดินย้อนกลับไปเรื่อยๆ อาจพบว่าก้านที่เราคิดว่าใหญ่ที่สุดนั้น เป็นเพียงแค่ก้านย่อยที่แตกแขนงแยกออกมาก็เป็นได้


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


ก็เปรียบเหมือนความรู้ เราจะรู้ได้ไงว่าอะไรคือพื้นฐานที่แท้จริงกันแน พันปีก่อนเราอาจคิดว่าถ้าได้รู้จักอะตอมก็คือรู้ถึงพื้นฐานแห่งทุกสรรพสิ่งแล้ว แต่ไม่กี่ร้อยปีก่อนเราเพิ่งพบว่าอะตอมยังแบ่งแยกย้อยลงไปได้อีกเป็นโปรตอน นิวตรอน อิเล็กตรอน และปฏิรูปพื้นฐานความรู้ด้านเคมีขึ้นใหม่ ส่วนปัจจุบันเรายังแบ่งย่อยสิ่งต่างๆ ลงไปได้อีกเป็นควาร์ก หรือแม้กระทั่งเสนอทฤษฎีสตริงซึ่งเป็นพื้นฐานของทุกอนุภาค


อีกนัยหนึ่ง พื้นฐานความรู้ที่ลึกเกินไปอาจเป็นสิ่งไม่จำเป็น หรือยิ่งไปกว่านั้น มันอาจเป็นตัวถ่วงให้เราไม่กล้าที่จะเดินหน้าสำรวจโลกใหม่ๆ ก็เป็นได้


Dec 29, 2015

เจองานจับฉลากแลกของขวัญที่มีคนได้ของตัวเอง ไม่ต้องตกใจไปโอกาสมีสูงเกินครึ่ง

วันก่อนเห็นเพื่อนแลกของขวัญปีใหม่แล้วดันจับได้ของตัวเองก็ฮาดี มองเผินๆ เหมือนว่าเรื่องแบบนี้จะเกิดขึ้นได้ยาก แต่เท่าที่เคยเล่นแลกของขวัญมา เรื่องประมาณนี้เกิดขึ้นบ่อยกว่าที่คิดนะ เลยลองมานั่งโซ้วสมการความน่าจะเป็นดู ว่าโอกาสที่ในงานเลี้ยงหนึ่งๆ จะมีคนจับของขวัญได้ของตัวเองเป็นเท่าไหร่

แต่ก่อนอื่นต้องเข้าใจกฎการแลกของขวัญ ซึ่งมีขั้นตอนประมาณนี้
  1. ทุกคนซื้อของขวัญมาร่วมสนุกในงานคนละอย่าง โดยเอาของขวัญใส่กล่องไว้ไม่ให้ใครรู้ว่าด้านในเป็นอะไร
  2. พอแต่ละคนไปถึงงาน คนจัดงานจะแปะป้ายชื่อเจ้าของของขวัญไว้กับกล่อง พร้อมหย่อนฉลากรายชื่อของคนนั้นลงไห
  3. เมื่อถึงเวลาแลกของขวัญ ให้ประธานในพิธีหรือใครซักคนที่ใหญ่ที่สุดในงาน เป็นคนเริ่มจับฉลากรายชื่อจากไหเป็นคนแรก
  4. ประธานจับได้ชื่อใครก็เอากล่องของขวัญของคนนั้นไป แล้วก็ให้คนที่โดนจับชื่อเมื่อกี้ เป็นคนเริ่มจับฉลากในไหวนต่อไปเรื่อยๆ
  5. ถ้าเกิดมีคนจับได้ชื่อประธาน จะถือว่าเป็นการแลกของขวัญที่ครบรอบวง ก็ให้คนนั้นเอาของขวัญของประธานไป แล้วผู้จัดงานเลือกใครซักคนจากคนที่ยังไม่ได้จับของขวัญ มาเป็นหัวจับของขวัญรอบต่อไป
  6. การแลกของขวัญจะสิ้นสุดลงเมื่อไม่มีฉลากรายชื่อเหลือในไห
จากกฎนี้เห็นได้ชัดว่า เมื่อถึงตอนใกล้จบการแลกของขวัญที่มีฉลากรายชื่อในไหเหลือแค่ 2 ใบ ใบนึงจะเป็นของคนแรกและอีกใบเป็นของคนที่ยังไม่ได้จับฉลาก ดังนั้นโอกาสที่คนจับก่อนจะจับฉลากได้ชื่อคนแรกและบังคับให้คนสุดท้ายจับชื่อตัวเอง โอกาสนี้จะสูงถึง 50% เลย

ยกตัวอย่างให้ชัดขึ้นไปอีก สมมตินายประธานเป็นคนจับฉลากคนแรก แล้วก็มีคนจับตามๆ กันมาอีกหลายคน จนตอนสุดท้ายเหลือคนที่ยังไม่ได้จับฉลาก 2 คนคือสมศักดิ์กับมานี ถ้าคนก่อนหน้าจับฉลากได้ชื่อสมศักดิ์ แปลว่าสมศักดิ์จะเป็นคนจับฉลากคนต่อไป แต่รายชื่อในไหตอนนี้จะเหลือแค่มานีกับประธานเท่านั้น ที่โอกาสครึ่งๆ ถ้าสมศักดิ์จับได้ชื่อประธาน ก็จะเหลือคนไม่ได้จับฉลากคือมานี และฉลากที่เหลือในไหก็เป็นชื่อมานี นั่นก็คือมานีจะต้องจับฉลากได้ของขวัญของตัวเองแน่นอน

โอกาสที่ในงานแลกของขวัญงานหนึ่ง จะมีคนได้ของขวัญเป็นของตัวเองสูงถึง 50% ซึ่งจริงๆ แล้วโอกาสน่าจะมากกว่านี้อีกนิดหน่อยด้วยซ้ำ เพราะยังไม่ได้คิดโอกาสที่คนกลางๆ จะจับได้ของตัวเองทันทีหลังจากเกิดการจับฉลากครบรอบวงเลย



คิดได้แค่นี้ก็พอใจแล้ว แต่ก็โดน @NungNing ถามต่อว่า ถ้าไม่จับฉลากกันต่อไปเรื่อยๆ แบบนี้หละ แต่เปลี่ยนเป็นให้ทุกคนหยิบฉลากจากไหพร้อมกันเลย แล้วค่อยเปิดดูทีเดียวว่าใครได้ของขวัญจากใคร ยังจะมีโอกาสสูงอยู่มั้ยที่จะมีคนจับฉลากได้ของขวัญของตัวเอง?

ตอนนั้นหมดพลังแล้ว ให้คิดต่ออย่างเดียวคงไม่ออกแน่ๆ เลยเขียนโปรแกรมง่ายๆ มาดูค่าเอาเลยนั่นแหละ
#!/usr/bin/env python3
from random import shuffle
from collections import Counter

shuffled = lambda ls: shuffle(ls) or ls
own_gift = lambda n: [i for i, j in enumerate(shuffled(list(range(n)))) if i == j]
สมมติให้มีคนแลกของขวัญกัน 10 คน ทดลองไปล้านครั้ง พบว่ามีประมาณสามแสนหกหมื่นครั้งเท่านั้น (36%) ที่การแลกของขวัญนั้นจะไม่มีคนได้ของตัวเอง

ตอนนั้นก็ตกใจว่าเลข 36% มาจากไหน นึกว่าโอกาสไม่ได้ของตัวเองมันจะมีค่าสูงกว่า 50% ที่คิดได้จากการแลกของขวัญวนไปเรื่อยๆ ซะอีก นี่กลายเป็นว่าโอกาสไม่ได้ของตัวเองมันยิ่งเหลือน้อย เลยทดลองเพิ่มโดยนับแยกว่า ในแต่ละการทดลองแลกของขวัญนั้น มีคนกี่คนที่ได้ของขวัญของตัวเอง

ทดลองแลกของขวัญกัน 10 คนล้านครั้งเหมือนเดิม แล้วก็ตกใจเพิ่มกับหน้าตาผลลัพธ์ที่ไม่โกหกดังนี้

จำนวนคนที่จับได้ของขวัญของตัวเอง จำนวนครั้งของเหตุการณ์ที่เกิด ข้อคาดการณ์โอกาสการเกิด
0 368,243 1/0! x 1/e
1 367,640 1/1! x 1/e
2 183,800 1/2! x 1/e
3 61,406 1/3! x 1/e
4 15,269 1/4! x 1/e
5 3,014 1/5! x 1/e
6 545 1/6! x 1/e
7 67 1/7! x 1/e
8 16 1/8! x 1/e

ถึงตอนนี้ก็พอจะเข้าใจแล้วว่ามันเกิดอะไรขึ้น



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

เริ่มจากกรณีที่ไม่มีใครจับของขวัญได้ของตัวเองก่อน จะได้ว่า
  • คนแรกมีโอกาสจะจับไม่ได้ของตัวเองที่ (n-1)/n
  • คนที่สองมีโอกาสที่จะจับไม่ได้ของตัวเอง แบ่งได้ 2 แบบ คือ
    • คนแรกจับได้ของคนที่สองไปแล้ว (โอกาส 1/n) ดังนั้นคนที่สองจับไม่ได้ของตัวเองแน่ๆ (โอกาส 1) รวมโอกาสได้เป็น 1/n
    • คนแรกจับไม่ได้ของคนที่สอง (โอกาส (n-1)/n) คนที่สองต้องจับหลบให้ไม่ได้ของตัวเอง (โอกาส (n-2)/(n-1)) รวมโอกาสได้เป็น ((n-1)/n)((n-2)/(n-1))
    • สรุปว่าโอกาสทั้งหมดที่คนที่สองจะจับไม่ได้ของตัวเองคือ 1/n + ((n-1)/n)((n-2)/(n-1))
  • โอกาสคนที่สามจะคล้ายๆ กับของคนที่สอง คือ 2/n + ((n-2)/n)((n-3)/(n-2))
  • คนที่ i มีโอกาส (i-1)/n + ((n-i-1)/n)((n-i)/(n-i-1)) ที่จะจับไม่ได้ของตัวเอง

เนื่องจากเหตุการณ์ทั้งหมด ต้องเกิดขึ้นพร้อมกัน ดังนั้นสรุปได้ว่า
P(n คนไม่มีใครได้ของตัวเอง) = (n-1)/n x (1/n + ((n-1)/n)((n-2)/(n-1))
                                  x ... x ((i-1)/n + ((n-i-1)/n)((n-i)/(n-i-1))
                                  x ... x ((n-n)/n + 0)
                        = (n-1)/n x (1/n + (n-2)/n) x .... x ((i-1)/n + (n-i)/n) x ... x 1
                        = Product  (n-1)/n  for integer i from 1 to n-1
                        = Product  1 - 1/n  for integer i from 1 to n-1
                        = (1-1/n)^(n-1)
จากสมการข้างต้น ถ้าลองแทนตัวเลขน้อยๆ เช่น n=1,2,3,4 เข้าไป จะได้คำตอบเป็น 1, 0.5, 0.445, 0.422 ตามลำดับ แต่เมื่อลองแทนเลขมากๆ แล้ว ผลลัพธ์จะลู่เข้าหาค่าเดียวกันคือ 0.367879... นั่นก็เพราะ
limit  (1-1/n)^(n-1)  when n -> infinity = limit  (1-1/n)^n  when n -> infinity
                                         = limit  ((n-1)/n)^n  when n -> infinity
                                         = limit  (n/(n+1))^n  when n -> infinity
                                         = limit  1/((n+1)/n)^n  when n -> infinity
                                         = 1/(limit  ((n+1)/n)^n  when n -> infinity)
                                         = 1/(limit  (1+1/n)^n  when n -> infinity)
                                         = 1/e
สวยมั้ย อยู่ดีๆ ก็ได้ค่า e ออกมาเฉยเลย

ส่วนการพิสูจน์กรณีมีคนได้ของขวัญของตัวเองแค่คนเดียว จะมีฟอร์มที่ต่างไปเล็กน้อยเนื่องจากเหตุการณ์จะไม่เกิดขึ้นต่อเนื่องกันแล้ว เขียนสมการออกมาได้เป็น
P(n คนมี 1 คนได้ของตัวเอง) = (1/n)P(n-1 คนไม่มีใครได้ของตัวเอง)
                              + ((n-1)/n)(1/(n-1))P(n-2 คนไม่มีใครได้ของตัวเอง)
                              + ...
                              + ((n-(n-1))/n)(1/(n-(n-1)))P(n-n คนไม่มีใครได้ของตัวเอง)
                        = (1/n)P(n-1 คนไม่มีใครได้ของตัวเอง)
                              + (1/n)P(n-2 คนไม่มีใครได้ของตัวเอง)
                              + ...
                              + (1/n)P(n-n คนไม่มีใครได้ของตัวเอง)
                        = (1/n)(Sum  P(n-i คนไม่มีใครได้ของตัวเอง)  for integer i from 1 to n)
สมมติว่ามีคนมาร่วมงานแลกของขวัญเยอะ อาจจะมองว่า P(n-i คนที่ไม่มีใครได้ของตัวเอง) มีค่าไม่แตกต่างกันสำหรับแต่ละค่า i เลยก็ได้ ทำให้เราสามารถยุบผลรวมยุ่งๆ ด้านหลัง ให้เหลือแค่เอาความน่าจะเป็นไปคูณด้วย n ก็พอ ซึ่งมันก็จะโดนหักล้างทิ้งไปกับ (1/n) ด้านหน้าทันที

ดังนั้นจะได้ว่า ค่าความน่าจะเป็นของงานเลี้ยงแลกของขวัญที่จะมีคนหนึ่งคนได้จับได้ของตัวเอง เป็น 1/e เท่ากันกับกรณีไม่มีใครจับได้ของตัวเองเลย

สำหรับกรณีที่มีคนสองคนจับได้ของตัวเอง จะเริ่มพิสูจน์ยุ่งยากมากขึ้นไปอีก แต่ใจความจะคล้ายๆ กับการพิสูจน์คนเดียวได้ของตัวเอง ซึ่งได้แก่การมองแยกกรณี โดยแบบที่หนึ่งจะให้คนแรกจับได้ของตัวเองไปแล้ว หลังจากนั้นจึงหาโอกาสของคนที่เหลือที่จะมีหนึ่งคนจับได้ของตัวเองเป็นเท่าไหร่ กับแบบที่สองที่ให้คนแรกจับไม่ได้ของตัวเอง แล้วหาโอกาสของคนที่เหลือที่จะมีสองคนจับได้ของตัวเองเป็นเท่าไหร่ (ถ้าคิดด้วยวิธีนี้ไม่ผิด จะได้หน้าตาสมการเป็นการหาผลรวมของเทอม ((n-k)/n)((n-k+1)/(n-1))((n-k+2)/(n-2))...(1/(n-k+1))P(n-k คนมี 1 คนได้ของตัวเอง))

กรณีอื่นๆ ที่เหลือเนื่องด้วยพื้นที่กระดาษไม่พอเขี.... เฮ้ย ไม่ใช่แฟร์มาต์! จริงๆ แล้วต้องบอกว่าการพิสูจน์ด้วยท่านี้มันเริ่มยุ่งยากซับซ้อนจนไปต่อไม่ไหวแล้ว แต่ถ้ายังจำกันได้ นี่มันการแจกแจงปัวซองที่มีค่า λ=1 ชัดๆ ก็เอาเป็นว่าใครอยากพิสูจน์ต่อก็ลองใช้ท่านี้ดูนะ อาจจะได้บทพิสูจน์ที่สั้นและสวยงามกว่าด้วย



เห็นแล้วว่าการสุ่มจับฉลากทั้งหมดทีเดียวพร้อมกันมันไม่เวิร์ค กลับไปดูวิธีดั้งเดิมที่ค่อยๆ ให้จับฉลากทีละคน จับได้ชื่อใครก็ให้คนนั้นไปจับฉลากต่อ วิธีนี้เขียนเป็นโค้ดออกมาได้ว่า
from random import randint

def exchange(ls):
    owned = [None for _ in ls]
    current = 0
    while ls:
        if owned[current] is not None:
            current = owned.index(None)
        gone = randint(0, len(ls)-1)
        owned[current] = ls[gone]
        current = ls[gone]
        del ls[gone]
    return owned
เอาฟังก์ชัน exchange ไปแทนที่ฟังก์ชัน shuffled ข้างบนแล้วทดลองใหม่ ก็ได้ผลลัพธ์ที่ไม่ต่างไปจากตารางแรกซักเท่าไหร่เลย



ดังนั้นไม่ว่าจะให้จับฉลากพร้อมกันหมดในทีเดียว หรือจะค่อยๆ จับทีละคน พอจับได้ชื่อใครก็ให้คนนั้นจับต่อ ทั้งสองวิธีนี้ยังไงก็จะมีคนซวยโชคดีอย่างน้อยหนึ่งคน ที่จับได้ของขวัญของตัวเองด้วยโอกาสสูงเท่ากันที่ 1-1/e หรือประมาณ 63.21% ครับ โดยที่(แทบ)ไม่ต้องสนใจด้วยว่ามีคนมาแลกของขวัญกันกี่คน

แต่วิธีแก้ก็ไม่ได้ยากเกินไป ใช้วิธีค่อยๆ ไล่จับฉลากไปทีละคนนั่นแหละ (จะได้เล่นมุกคั่นเวลา ดูรีแอคชันของคนได้ของขวัญด้วย) แล้วเพิ่มกฎเข้าไปสองข้อดังนี้
  1. ถ้าจับได้ชื่อตัวเอง ให้จับฉลากเพิ่มอีกใบเพื่อจะได้เอาของขวัญจากคนนั้น แต่อย่าลืมหย่อนฉลากชื่อตัวเองที่จับขึ้นมาเมื่อกี้ใส่กลับลงไหก่อนเอาไปให้จับฉลากต่อ
  2. ตอนเหลือคนยังไม่ได้จับฉลาก 2 คน ไม่ต้องให้จับฉลากแล้ว ให้คนที่กำลังจะจับฉลากไปเอาของขวัญจากคนสุดท้ายได้เลย ส่วนคนสุดท้ายก็ไปเอาของขวัญจากคนแรก
อย่างไรก็ตาม หากรู้สึกว่ากฎพวกนี้มันวุ่นวายจนทำให้งานกร่อยก็ไม่ต้องไปทำตามหรอก เพราะการมีคนจับได้ของขวัญของตัวเองก็ไม่ใช่เรื่องเลวร้ายแต่อย่างใด หนำซ้ำมันยังเป็นประสบการณ์สุดประหลาดที่หาได้ยาก และเป็นเรื่องขำขันชั้นดีที่สามารถเก็บเอาไว้เล่นได้เรื่อยๆ ยันลูกบวชเลยหละ

สวัสดีปีใหม่ 2016 ล่วงหน้าครับ

May 22, 2015

Voronoi จากภาพ

วันก่อนไปเจอคำถามใน Code Golf มา เค้าให้สร้าง Voronoi เพื่อเลียนแบบภาพต้นฉบับครับ

อ่านแล้วสนุกดี ผลงานของคนอื่นที่ได้ๆ ก็ดูเหมือนงานศิลป์พวก Divisionism/Pointillism ด้วย เลยลองกลับมาเขียนเองแบบมั่วๆ บ้าง 555

import random
from PIL import Image
from PIL.ImageFilter import FIND_EDGES, GaussianBlur, SHARPEN

def coordinate(width, height):
    yield from ((x, y) for y in range(height) for x in range(width))

def normalized(choices):
    lower = min(weight for _, weight in choices)
    upper = max(weight for _, weight in choices)
    norm = lambda weight: (weight-lower) + (upper-lower)//4
    return [((x, y), norm(weight)) for (x, y), weight in choices]

def random_by_weight(choices):
    rand_val = random.uniform(0, sum(weight for _, weight in choices))
    index = 0
    count = 0
    while count < rand_val:
        count += choices[index][1]
        index += 1
    return choices.pop(index-1)[0]

def init_centroids(image, cells):
    width, height = image.size
    edge_img = image.filter(FIND_EDGES).filter(GaussianBlur).filter(SHARPEN)
    weight = lambda x, y: 256 - max(edge_img.getpixel((x, y)))
    choices = [((x, y), weight(x, y)) for x, y in coordinate(width, height)]
    return [random_by_weight(normalized(choices)) for _ in range(cells)]

def init_rgbs(image, centroids):
    rgb_im = image.convert('RGB')
    return [rgb_im.getpixel((x, y)) for x, y in centroids]

def simulate_voronoi(image_path, cells=25, scale=None):
    image = Image.open(image_path)
    if scale is not None:
        image.thumbnail((scale, scale), Image.ANTIALIAS)
    centroids = init_centroids(image, cells)
    rgbs = init_rgbs(image, centroids)
    return image.size, list(zip(centroids, rgbs))
ผลลัพธ์ที่ออกมาก็ประมาณนี้ครับ (ที่ 500 เซลล์ Voronoi)

เอาไปสู้เขาไม่ได้หรอก เพราะเขียนไปมั่วแบบคนไม่มีพื้นฐาน image processing อะไรเลย แต่ก็สนุกดีได้ลองใช้ PIL เป็นครั้งแรกด้วย

Apr 16, 2015

Code Jam 2015 รอบคัดเลือก


ปีนี้ประมาทไปเยอะครับ คือตอนแรกตั้งใจว่าจะถ่างตามาเริ่มเล่นตั้งแต่ปล่อยโจทย์ตอน 6 โมงเช้าเลย แต่สุดท้ายก็เผลอหลับไปตอนตี 4 ตื่นมารู้ตัวอีกทีก็บ่ายแล้ว เลยพักสมองออกเล่นบอร์ดเกมกะ @tamwb และ @aibig ตามนัดปาร์ตี้สงกรานต์ แถมกลับมาก็ยังไม่มีสมาธิแก้โจทย์มัวแต่ดู MV เกาหลีอีก กว่าได้เริ่มทำแข่งจริงจังก็เหลือเวลาไม่กี่ชั่วโมงแล้ว

ตอนอ่านโจทย์ผ่านๆ ทุกข้อนี่คิดว่าง่าย ปีนี้ได้แก้มือเก็บคะแนนเต็มแน่ๆ ... ที่ไหนได้ bug ซ่อนเร้น 1 ข้อ แล้วก็ตีความโจทย์อีก 1 ข้อ สุดท้ายแก้ไม่ทัน ได้ส่งแค่ 2 ข้อ ดีแค่ไหนแล้วที่ยังไม่ตกรอบ 555



ข้อแรกง่ายฮะ แต่ถ้ายังอ่านโจทย์แล้วไม่อินขอให้ดูคลิปนี้ตอนนาทีที่ 7:05 และ 7:23 (เรียบร้อยแล้วก็กลับไปดูตั้งแต่ต้นด้วยจะดีมาก) :P



ส่วน code ก็ประมาณ
import Text.Printf (printf)
import Data.List.Split (splitOn)

invites audience =
    let aux []     _   _   inv = inv
        aux (p:ps) shy acc inv = aux ps (shy+1) (acc+p+more) (inv+more)
            where more = max 0 (shy-acc)
    in  aux audience 0 0 0

test t = do
    it <- getLine
    let [_, ps] = splitOn " " it
    printf "Case #%d: %d\n" t $ invites [read [p] :: Integer | p <- ps]

main = do
    loop <- getLine
    sequence_ [test t | t <- [1..read loop :: Integer]]
(อย่างไรก็ตาม เท่าที่เคยไปคอนเสิร์ตหลายครั้ง รู้สึกว่าคนไทยจะมีค่าความอายเป็น infinity นะ)


ข้อสองข้ามไปนะฮะ ทำไม่ได้ (ข้อสี่ด้วย)

ส่วนข้อสามนี่สนุกดี ถ้ายังจำเนื้อหาม.ปลายได้ เราจะมีเลขมหัศจรรย์ตัวนึงเรียกว่า i ซึ่งมีนิยามที่สวยงามเรียบง่าย คือ
i^2 = -1
สังเกตเพิ่มอีกหน่อยก็จะเห็นว่า
i^4 = 1    (multiplicative identity)
โจทย์เพิ่มความสนุกโดยให้ตัวเลขมหัศจรรย์ j และ k มาด้วย ทั้งสองตัวที่เพิ่มเข้ามานี้ มีนิยามแบบเดียวกันกะ i เด๊ะๆ เลย เพียงแต่มันไม่ใช่ i เพราะเรานิยามการคูณระหว่างตัวเลขเหล่านี้ว่า
ij = k
jk = i
ki = j
ความซับซ้อนของตัวเลขชุดนี้ยังเพิ่มขึ้นไปอีก เมื่อเรานิยามให้มันไม่มีสมบัติสลับที่ด้วย ซึ่งในที่นี้คือ
ji = -k
kj = -i
ik = -j
เพื่อความสะดวก สมมติชื่อให้เลขชุดนี้ว่า ℍu ละกันครับ โชคดีที่อย่างน้อยในเซ็ต ℍu = {1, i, j, k} มันยังมี
  • สมบัติปิดบนการคูณ คือ ถ้าใช้การคูณอย่างเดียวบนสมาชิกของเซ็ต ℍu ผลลัพธ์ที่ได้ก็จะเป็นสมาชิกของเซ็ต ℍu ด้วย
  • สมบัติจัดกลุ่มการคูณ คือ a(bc) = (ab)c สำหรับทุก a, b, c ∈ ℍu
สรุปง่ายๆ ถ้าใครถนัดเรื่องการคูณไขว้เวกเตอร์ ก็นึกถึงกฎมือขวาได้เลย

จาก X, LS ที่โจทย์ให้ การแก้ปัญหาแบบตรงๆ ก็คือวิ่งผ่านชุดตัวเลข LS เป็นจำนวน X ครั้ง โดยเลขแต่ละตัวเลข s ∈ LS ที่วิ่งผ่าน ก็จัดการคูณเก็บไว้ในใจไปเรื่อยๆ ถ้าเจอตัวเลขครบทั้ง i, j, k และหางที่เหลือคูณกันออกมาเป็น 1 ก็ตอบว่าถูกได้เลย

ถ้าอ่านมาจนถึงตรงนี้แล้วยังงงไม่หาย ลองดูการทำงานตามขั้นตอนวิธีด้วยตัวอย่างนี้
X  = 10
LS = ji

sequence: ji ji ji ji ji ji ji ji ji ji
running:  --->--->-------->----------->
found:       i   j        k           1
โจทย์ข้อนี้ฉลาดตรงที่เอาตัวเลข Xขนาดใหญ่มากๆ มาขู่ ทั้งที่เราสามารถ modulo มันทิ้งให้เหลือแค่หลักสิบได้ (ซึ่งก็คือ คำถามชุดใหญ่จะมีความยากแทบไม่ต่างกับคำถามชุดเล็กเลย)

การลดขนาด X ทำได้โดยใช้ข้อสังเกตที่ว่า i^4 = j^4 = k^4 = 1 ซึ่งบอกเป็นนัยอยู่ 2 อย่างคือ
  • ถ้าวิ่งคูณหา i, j, k แต่ละตัวเกิน 4 รอบแล้วยังไม่เจอ ไม่ต้องวิ่งต่อ เพราะแต่ละรอบที่วิ่งเพิ่ม จะไม่ได้ผลลัพธ์ใหม่แล้ว
  • หลังจากหา i, j, k ได้ครบแล้ว หางที่เหลือสามารถ mod 4 ได้
ดังนั้น เราจะวิ่งผ่านชุดตัวเลขนี้ 3 ครั้ง (เพื่อหา i, j, k) โดยแต่ละครั้งวิ่งไม่เกิน 4 รอบ ซึ่งก็เท่ากับวิ่ง 12 รอบ แล้วรอบที่เหลือก็จับ mod 4 ได้เลย

หรือพูดจาภาษา code ก็คือ
def correct_misspelling(ijks, x):
    x = min(x, 12+(x%4))
    accumulate = Quaternion('+1')
    searching = [Quaternion(it) for it in 'kji']
    for _ in range(x):
        for it in ijks:
            accumulate *= Quaternion(it)
            if searching and searching[-1] == accumulate:
                accumulate = Quaternion('+1')
                searching.pop()
    return not searching and accumulate == Quaternion('+1')

for case in range(int(input())):
    _, x = [int(n) for n in input().split()]
    ijks = input().strip()
    ans = 'YES' if correct_misspelling(ijks, x) else 'NO'
    print('Case #{}: {}'.format(case+1, ans))
ข้อนี้ถ้าจะ optimize เพิ่มก็ยังทำได้อีก (code จริงที่เอาไปส่งก็เขียนสวยน้อยกว่านี้เพราะมัวแต่ optimize มากเกินไป) โดยเรามีข้อสังเกตดังนี้
  • aba = b สำหรับทุก a, b ∈ ℍu
  • เราสามารถทดเก็บค่าสุดท้ายที่ได้จากการวิ่งผ่านเพื่อหาผลคูณของทุกตัวใน LS ไว้ได้ แล้วพอจะวิ่งผ่าน LS อีกรอบก็ใช้ค่าที่ทดไว้ได้เลย
  • ในการทดลองจริง จำนวน X ที่จำเป็นต่อการคำนวณไม่เคยมีค่าเกิน 8 เลย ซึ่งก็คือ x = min(x, 4+(x%4))
ถ้ายังอยากรู้ว่าข้อสังเกตเหล่านี้จริงหรือไม่ ถ้าจริงแล้วจะช่วยการคำนวณอย่างไร proof เพิ่มดูเองได้เลยครับ

Jun 1, 2014

การสอนด้วยตัวอย่างที่ดีนั้นยาก

ถ้าช่วงนี้ใครเข้าเว็บ Python.org คงจะเห็นตัวอย่างโค้ดนี้
list = [2, 4, 6, 8]
sum = 0
for num in list:
    sum = sum + num
print("The sum is:", sum)
เรียบง่ายและงดงาม...

แต่ถ้าไปถามคนที่เขียน Python ประจำนี่ด่าตรึมเลย เพราะโค้ดด้านบนนี้มีจุดอ่อนอยู่ 2 ประการ
  1. มีการประกาศชื่อตัวแปร list กับ sum ซึ่งทั้ง 2 ชื่อนี้ เป็นชื่อของฟังก์ชันมาตฐาน ถึงแม้ว่าจะไม่มีข้อห้ามให้ประกาศชื่อซ้ำได้เหมือนกับพวก reserved keyword แต่การทำเช่นนี้ก็บาปมากเพราะ 2 ฟังก์ชันนี้จะไม่สารมารถถูกเรียกใช้ได้โดยง่ายอย่างที่ตัวภาษาออกแบบไว้
  2. Python มีฟังก์ชัน sum สำหรับหาผลรวมของชุดตัวเลขให้ใช้งาน นอกจากจะเขียนได้ง่ายและทำงานเร็วกว่าแล้ว มันยังสร้างสุขลักษณะที่ดีต่อการเขียนโปรแกรมสำหรับใช้งานในโลกจริงด้วย
  3. (ข้อนี้เพิ่มให้ในฐานะนักทฤษฎีคณิตศาสตร์และคอมพิวเตอร์) ขั้นตอนวิธีข้างต้นสามารถพิสูจน์ลดรูปให้เป็น O(1) ได้ ถ้าลำดับเลขนำเข้าเป็นเลขจำนวนเต็มเรียงกันทั้งหมด หรือเป็นเลขคู่เรียงกัน หรือเป็นเลขคี่เรียงกัน
ตัวอย่างนี้จึงเป็นตัวอย่างที่ไม่ค่อยดีในการนำมาสอนเขียนโปรแกรมที่จริงจังหน่อย เพราะว่าโค้ดแบบนี้ไม่ควรเกิดขึ้นในโลกจริงเลย (แน่นอนว่าคนที่ยกตัวอย่างนี้มา อาจแก้ต่างว่าโค้ดนี้มีเจตนาที่จะยกตัวอย่างลูปให้มือใหม่หัดเขียนโปรแกรม ซึ่งโดยส่วนตัวผมไม่ค่อยเห็นด้วยเพราะนี่คือหน้าเว็บที่มีทุกระดับฝีมือมาเยี่ยมเยียน ไม่ใช่บทที่ 1 ในหนังสือ Programming 101)

ถ้าเช่นนั้น แล้วจะมีตัวอย่างแบบไหนที่ดีพอหละ?

ผมคิดว่าตัวอย่างที่ดี เข้าเกณฑ์แค่ข้อเดียวก็พอ คือต้องไม่ทำหน้าที่ซ้ำซ้อนกับสิ่งที่มีอยู่แล้ว (do not reinvent the wheel)

ปัญหาสำคัญสำหรับแนวคิดนี้ คือของที่ง่ายพอจะเอามาเป็นตัวอย่างมักมีอยู่เกือบหมดแล้ว จะให้สอนโดยไม่ยกตัวอย่างซ้ำกับของที่มีเลยคงยาก

โชคดี (?) ที่ Python ยังไม่มีฟังก์ชันผลคูณรวม (product) ทำให้โค้ดในโลกความจริง ที่สามารถนำมาสอนเขียนได้ดังนี้
numbers = [2, 4, 6, 8]
product = 1
for number in numbers:
    product = product * number
print("The product is:", product)
เพราะการสร้างตัวอย่างที่ดีไม่ใช่เรื่องง่ายๆ

May 28, 2014

อย่าเสียเวลาให้กับลูปที่เรียก System บ่อยๆ

ช่วงนี้มีวิเคราะห์ string บนไฟล์จำนวนมหาศาล (หลายแสนไฟล์) ตอนแรกก็เขียน Python ง่ายๆ ตามความเคยชิน
counter = 0
for file in glob('*.cpp'):
    if '//' in file.read():
        counter += 1
print(counter)
โค้ดนี้ทำงานได้ดีเลยแหละ ติดตรงที่ไฟล์ที่นำมาวิเคราะห์บางไฟล์ก็ใช้ encoding แปลกๆ แล้วเปิดอ่านไม่ได้ (ยังงงอยู่ว่า Python ไม่เช็คว่าเป็น ASCII หรือ UTF-8 ให้อัติโนมัติ?) ก็แก้ๆ โค้ดด้านบนนี้โดยใช้โมดูล codecs เข้าช่วย ไม่มีปัญหาอะไร

แต่แล้วก็คิดขึ้นมาได้ว่า งานประเภทนี้ถ้าทำใน Bash ไปเลยก็น่าจะแจ่มเหมือนกัน เพราะปรัชญาของ Unix คือให้โปรแกรมจัดการกับ I/O Stream แบบข้อความอยู่แล้ว เลยได้โค้ดออกมาเป็น
counter=0
for file in *.cpp
do
    grep -q '//' "$file" && : $((counter++))
done
echo $counter
ผลลัพท์ออกมาไม่ต่างกับ Python ด้านบน แต่ว่า Bash ทำงานช้ามากกกกก (ก.ไก่ 140 ตัว) ตอนแรกก็คิดว่า Bash มันน่าจะทำลูปไม่เก่งมั้ง เพราะเท่าที่อ่านโค้ดของโครงการ Linux ก็ไม่ค่อยเห็นเขียนลูปกันเท่าไหร่ แต่พอลองเขียนลูปเปล่าแสนรอบใน Bash มันก็ไม่เห็นจะช้า เลยคิดว่าน่าจะเป็นที่ grep มั้ง? แต่ก็ไม่อยากจะเชื่ออยู่ดีเพราะเคยอ่านงานของ Russ Cox ก็ไม่เห็นว่า grep มันจะช้าซักเท่าไหร่

อาจารย์ @juggapong ผ่านมาเห็นเลยยิงคำถามมาว่า มันช้าเพราะมัวแต่มุดเข้ามุดออก system หรือเปล่า?

ได้ยินดังนั้นเลยใช้ time จับเวลาดู โอ้ว user กับ sys ใช้เวลาพอๆ กันเลยแฮะ

ถึงตอนนี้ก็รู้แล้วว่าสาเหตุที่มันช้าเพราะ grep แต่ละครั้งเป็น process ใหม่เลย ต้องรอให้ system จองทรัพยากรเสร็จก่อนจึงจะใช้งานได้ แถมพอจบงานก็ทำลาย process นี้ทิ้งซะอีก ไม่ได้จัดสรรทรัพยากรอะไรกันระหว่างลูปเลย ... เทียบง่ายๆ คงเหมือนกับการที่มีเตารีดแสนอัน เสียบปลั๊กเตารีดตัวแรกรอร้อนแล้วจึงรีดผ้าหนึ่งชิ้น พอผ้าชิ้นนี้เรียบเนี๊ยบจนพอใจก็ถอดปลั๊กเตารีดซะงั้น แล้วไปเสียบปลั๊กเตารีดอีกตัวมารีดผ้าอีกชิ้นนั่นเอง

คิดได้ดังนี้ก็เปลี่ยนโค้ดไปเป็น (ใช้เตารีดตัวเดียว เสียบปลั๊กรอร้อนครั้งเดียว)
grep -l '//' *.cpp | wc -l
คราวนี้เร็วส์เสียยิ่งกว่า Python อีก :P

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

Apr 13, 2014

Code Jam 2014 รอบคัดเลือก


ปีนี้ปัญหาเยอะนิดหน่อยเพราะวันแข่งตรงกับวันเดินทางออกต่างจังหวัดพอดี และเนื่องจาก #แล้วรถไฟความเร็วสูงกูหละ ทำให้เสียเวลาเดินทางสั้นๆ แค่ 300 กว่ากิโลเมตรไป 8 ชั่วโมงครึ่ง แถมพอถึงที่หมายแล้วก็หาเวลาส่วนตัวไม่ค่อยได้เลย กว่าจะได้ลงมือจริงๆ ก็ตีหนึ่งแล้ว (มีเวลา 8 ชั่วโมงจาก 27 ชั่วโมง) ผลงานที่ได้เทียบกับ effort ที่สามารถลงให้มันไปถือว่าน่าผิดหวังเล็กน้อย :\

อย่างไรก็ตาม ปีนี้ก็บรรลุเป้าข้อนึงคือใช้ภาษาอย่างน้อย 3 ภาษาเขียนแก้โจทย์ ก็ได้แก่ Python, Haskell และ CoffeeScript (อันที่จริง CoffeeScript นี้แปลจาก Python มาตรงๆ นาทีเดียวเสร็จเลยนะ :P) ปีต่อไปพยายามตั้งเป้าไว้ที่ Racket/Scheme กับ matrix processing language ซักตัวครับ

มาดูข้อแรกข้อ A ที่ตามธรรมเนียมจะง่ายที่สุดกัน ยิ่งใช้ Python ยิ่งง่ายใหญ่เลยเพราะมี set intersection ให้ใช้ เราก็แค่หาว่า intersect ของตัวเลขที่นักมายากลถามมานั้นคืออะไร ถ้ามีตัวเดียวก็คือตัวนั้นแหละ ถ้ามีมากกว่าหนึ่งแปลว่านักมายากลแม่งกาก แต่ถ้าไม่มีเลยนี่ผู้ชมโกงแล้วหละ

def foo(answer_1, matrix_1, answer_2, matrix_2):
    candidate_1 = matrix_1[answer_1-1]
    candidate_2 = matrix_2[answer_2-1]
    final = set(candidate_1) & set(candidate_2)
    if len(final) == 0:
        return 'Volunteer cheated!'
    elif len(final) > 1:
        return 'Bad magician!'
    else:
        return final.pop()

for case in range(int(input())):
    answer_1 = int(input())
    matrix_1 = [[int(n) for n in input().split()] for _ in range(4)]
    answer_2 = int(input())
    matrix_2 = [[int(n) for n in input().split()] for _ in range(4)]
    answer = foo(answer_1, matrix_1, answer_2, matrix_2)
    print('Case #{}: {}'.format(case+1, answer))



ข้อ B ถัดมาเริ่มจะต้องสังเกตละ แต่เนื่องจากผมตั้งใจเขียนข้อนี้ในสไตล์แบบ clean code พอสมควร (ไม่มี comment, ตั้งชื่อตัวแปรให้รู้เรื่อง) ดังนี้ข้อนี้ขออธิบายด้วย code อย่างเดียวนะครับ

def optimal_cookie(new_farm_cost, prodct_cookie, finale_cookie):
    income_cookie = 2.0
    optimal_spent = 0
    while True:
        new_farm_time = new_farm_cost / income_cookie
        expect_cookie = prodct_cookie + income_cookie
        wait_win_time = finale_cookie / income_cookie
        upgd_win_time = finale_cookie / expect_cookie + new_farm_time
        if wait_win_time > upgd_win_time:
            income_cookie = expect_cookie
            optimal_spent += new_farm_time
        else:
            optimal_spent += wait_win_time
            break
    return optimal_spent

for case in range(int(input())):
    answer = optimal_cookie(*[float(n) for n in input().split()])
    print('Case #{}: {:.7f}'.format(case+1, answer))

หมายเหตุว่าข้อนี้ถ้าไปดูคำตอบผมในเว็บ จะเห็นคำตอบสำหรับข้อเล็กเขียนด้วย Haskell ส่วนข้อใหญ่เขียนด้วย CoffeeScript ครับ แต่ในการทดลองนี้จะใช้อันไหนทดสอบที่ระดับใดก็ได้



กระโดดข้ามมาข้อ D ก่อน ข้อนี้ยากตรงที่เราต้อง proof จนได้ strategy ที่ดีที่สุด เราอาจสังเกตว่า ในเกม war ธรรมดามีวิธีเล่นที่ดีที่สุดดังนี้
  • ถ้า Naomi เล่นกล่องที่หนักกว่าทุกกล่องของ Ken ทางเลือกที่ดีที่สุดสำหรับ Ken คือเล่นกล่องที่เบาที่สุด เพื่อให้กล่องที่เหลือๆ มีประสิทธิภาพในการสู้สูงที่สุด
  • ถ้า Naomi เล่นกล่องที่เบากว่ากล่องบางกล่องของ Ken แล้ว Ken จะตอบโต้ด้วยกล่องที่หนักกว่ากล่องของ Naomi แต่เป็นกล่องที่เบาที่สุดที่ยังจะหนักกว่ากล่องของ Naomi เพราะการเล่นแบบนี้ทำให้ได้คะแนน และยังเก็บกล่องหนักๆ ไว้ใช้สู้ในอนาคตได้อีก
พอเกมเปลี่ยนมาเป็น deceitful war ทาง Naomi ก็ใช้ข้อได้เปรียบจากการล่วงรู้น้ำหนักกล่องของ Ken ทั้งหมดมาเล่น โดยยึดกลยุทธ์ดังนี้
  • Naomi สามารถเล่นกล่องที่เบาที่สุดของตน แต่หลอกว่ากล่องนี้มีน้ำหนักน้อยกว่ากล่องที่หนักที่สุดของ Ken เนื่องจาก Ken ยังยึดกฎเดิม เขาจะตอบโต้ด้วยกล่องที่หนักที่สุด (กินคะแนน)
  • Naomi สามารถเล่นกล่องที่เบาที่สุดของตน แต่หลอกว่ากล่องนี้มีน้ำหนักมากกว่ากล่องที่หนักที่สุดของ Ken เช่นกันเพราะ Ken ยึดกฎเดิม เขาจะตอบโต้ด้วยกล่องที่เบาที่สุด (ทิ้งของกากสุด)
เทคนิคสำคัญคือรอจังหวะที่แต่ละกล่องของ Naomi เมื่อเรียงน้ำหนักและนำไปจับคู่กับ Ken แล้ว ทุกกล่องมีน้ำหนักมากกว่า เมื่อนั้น Naomi จะใช้กลยุทธ์ล่างทันที แต่ถ้ายังไม่มีโอกาสก็จะยึดกลยุทธ์แรกไปก่อน เพื่อกำจัดกล่องหนักๆ ของ Ken ทิ้งไป

from collections import deque

def normal_war(ns, ks):
    c = 0
    while ns:
        n = ns.pop()
        if ks[-1] < n:
            ks.pop(0)
            c += 1
        else:
            i = 0
            while n > ks[i]:
                i += 1
            ks.pop(i)
    return c

def deceitful_war(ns, ks):
    while any(n < k for n, k in zip(ns, ks)):
        ns.popleft() and ks.pop()
    return len(ns)

for case in range(int(input())):
    _ = input()
    ns = sorted(float(n) for n in input().split())
    ks = sorted(float(n) for n in input().split())
    dw_answer = deceitful_war(deque(ns), deque(ks))
    nw_answer = normal_war(ns[:], ks[:])
    print('Case #{}: {} {}'.format(case+1, dw_answer, nw_answer))

ข้อ C ถามว่าในอุดมคติแล้ว เราสามารถกดปุ่มเดียวผ่าน minesweeper ได้มั้ย ซึ่งข้อนี้มีจุดน่าสนใจอยู่ที่การแตกเคสต่างๆ ออกมาให้ครอบคลุม ดังนี้
  • ถ้ากระดานมีความกว้างแค่ 1 หน่วย ยังไงก็กดครั้งเดียวผ่านแน่ เพราะเราก็ให้ระเบิดทั้งหมดไปกองด้านนึง อีกด้านก็ไม่มีอะไร
  • ถ้ากระดานมีความกว้าง 2 หน่วย เราจะได้ว่าต้องมีพื้นที่ว่าง 1 ช่อง (สำหรับกดครั้งเดียวผ่านเลย) หรือที่ว่างต้องเป็นเลขคู่ที่มากกว่า 2 (ถ้ามี 2 ช่องพอดีแม้ช่องนั้นจะอยู่ติดกันก็ต้องกด 2 ครั้ง)
  • กระดานแบบอื่นๆ ทั้วไป พยายามกั้นช่องว่างให้ตัวเองเป็นจัตุรัสก่อน แล้วค่อยๆ เวนคืนที่ว่างเพิ่มเติม ดังนี้
    • ที่ว่าง 1 ช่องเป็นจัตุรัสขนาด 1x1 ง่ายสุด คือกดครั้งเดียวผ่านเลย อย่างไรก็ตาม ที่ว่างแบบนี้ทำให้เราไม่สามารถมีที่ว่าง 2, 3 ได้
    • ที่ว่าง 4 ช่องเป็นจัตุรัสขนาด 2x2 สามารถมีติ่งเพิ่มออกมาอีก 2, 4 ช่องได้ คือจะได้ ที่ว่างรวมเป็น 6, 8 ช่อง ส่วนติ่งที่เพิ่มแค่ 1, 3 ช่องเพิ่มไม่ได้
    • ที่ว่าง 9 ช่องเป็นจัตุรัสขนาด 3x3 และใหญ่กว่า คราวนี้จะมีติ่งเป็นเศษแค่ 1 ช่องก็ได้ ทำโดยไปยืมพื้นที่ว่างหัวมุมมาจับคู่เข้ากับติ่งนี้ซะ
    • สรุปว่าท่านี้ไม่รองรับการเหลือที่ว่างแค่ 2, 3, 5, 7 ครับ
    • ท่านี้ต้องพลิกแพลงนิดหน่อยด้วย ตรงที่ถ้าเกิดขนาดจัตุรัสใหญ่เกินด้านใดด้านหนึ่ง ต้องเปลี่ยนไปใช้โหมดเติมเต็มทีละบรรทัดแทน
อ่านข้อความอย่างเดียวอาจคิดไม่ออก ดูตัวอย่างการ fill แบบต่างๆ บนตาราง 4x4 นี้น่าจะง่ายกว่า
Case #0:      Case #1:      Case #2:      Case #3:
c...          c...          c...          c...
....          ....          ....          ....
....          ....          ....          ...*
....          ...*          ..**          ..**

Case #4:      Case #5:      Case #6:      Case #7:
c..*          c..*          c..*          c..*
...*          ...*          ...*          ...*
...*          ...*          ..**          ...*
...*          ..**          ..**          ****

Case #8:      Case #9:      Case #10:     Case #11:
c.**          Impossible    c.**          Impossible
..**                        ..**
..**                        ..**
..**                        ****

Case #12:     Case #13:     Case #14:     Case #15:
c.**          Impossible    Impossible    c***
..**                                      ****
****                                      ****
****                                      ****
ส่วน code ก็หน้าตาประมาณ
def transpose(matrix):
    return [list(line) for line in zip(*matrix)]

def fill_line(r, c, left):
    width, length = sorted([r, c])
    matrix = [['*' for _ in range(length)] for _ in range(width)]
    for i in range(width):
        matrix[i][:left//width] = ['.' for _ in range(left//width)]
    matrix[0][0] = 'c'
    return matrix if r == len(matrix) else transpose(matrix)

def fill_rect(r, c, left):
    width, length = sorted([r, c])
    sq_length = int(left**0.5)
    matrix = [['*' for _ in range(width)] for _ in range(length)]
    if sq_length > width:
        i = 0
        while left - width > 0:
            matrix[i] = ['.' for _ in range(width)]
            left -= width
            i += 1
        if left == 1:
            matrix[i-1][-1:] = ['*']
            matrix[i][:2] = ['.' for _ in range(2)]
        else:
            matrix[i][:left] = ['.' for _ in range(left)]
    else:
        i = 0
        while i < sq_length:
            matrix[i][:sq_length] = ['.' for _ in range(sq_length)]
            left -= sq_length
            i += 1
        if left > sq_length:
            if left == sq_length + 1:
                matrix[i][:sq_length-1] = ['.' for _ in range(sq_length-1)]
                left = 2
            else:
                matrix[i][:sq_length] = ['.' for _ in range(sq_length)]
                left -= sq_length
            if i == len(matrix)-1:
                matrix = transpose(matrix)
            else:
                i += 1
            matrix[i][:left] = ['.' for _ in range(left)]
        elif left > 0:
            if left == 1:
                matrix[i-1][sq_length-1:sq_length] = ['*']
                left = 2
            matrix[i][:left] = ['.' for _ in range(left)]
    matrix[0][0] = 'c'
    return matrix if r == len(matrix) else transpose(matrix)

def mine_fill(r, c, m):
    width = min(r, c)
    left = r * c - m
    if width == 1:
        if left == 0:
            return None
        return fill_line(r, c, left)
    if width == 2:
        if left in [0, 2] or left % 2 == 1 and left != 1:
            return None
        return fill_line(r, c, left)
    else:
        if left in [0, 2, 3, 5, 7]:
            return None
        return fill_rect(r, c, left)

for case in range(int(input())):
    r, c, m = [int(n) for n in input().split()]
    answer = mine_fill(r, c, m)
    print('Case #{}:'.format(case+1))
    if answer is None:
        print('Impossible')
    else:
        for line in answer:
            print(''.join(line))
ข้อนี้เสียดายมากที่หา bug ไม่เจอ มาเจออีกทีตอนแข่งจบไปแล้ว ไม่งั้นรอบนี้คงได้ทำคะแนนเต็มให้เป็นเกียรติเป็นศรี \(TwT)/

Nov 10, 2013

Division ใน Relational Algebra

การหาร (DIVISION) ใน relational algebra อธิบายง่ายๆ คือการเอาตารางมา join กันแล้วนับว่าแต่ละ pk ของตารางตัวตั้ง มี row เท่ากับตารางตัวหารหรือไม่ ถ้าครบก็ตัด column ที่ปรากฏในตารางตัวหารทิ้งแล้วยุบรวม row ซ้ำด้วย pk เช่นเรามีข้อมูล
table: students
------+------+-------
 name | subj | grade 
------+------+-------
 bar  | code | A
 foo  | code | C+
 foo  | trek | A
 qux  | code | B+
 qux  | math | A
 qux  | trek | A
------+------+-------

table: colleges
-----------+------
  campus   | subj 
-----------+------
 leetademy | code
 leetademy | trek
 hack-cool | code
 hack-cool | math
-----------+------
คำสั่ง SQL ของการหารก็จะมีหน้าตาประมาณ
SELECT *
FROM   students
DIVISION (
    SELECT subj
    FROM   colleges
    WHERE  campus = 'leetademy'
) AS leetademy_courses
ซึ่งให้ผลลัพท์ออกมาคือ
------
 name 
------
 foo
 qux
------



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

อย่างไรก็ตามปัญหาจริงๆ คือมันเป็นท่าที่อาจารย์ database ส่วนใหญ่ชอบเอาไปออกสอบ ทำให้เราต้องหาท่าอ้อมๆ ที่ให้ผลลัพท์แบบเดียวกันมาใช้แทน

การจะเข้าใจมันได้ ถ้าคิดแบบ programming ก็แค่ดูว่า set ที่เป็นตัวตั้งหารนั้นเป็น superset ของ set ที่ถูกหารหรือเปล่า
students_set = defaultdict(set)
for name, subj, grade in students:
    students_set[name].add(subj)

colleges_set = defaultdict(set)
for campus, subj in colleges:
    colleges_set[campus].add(subj)

output = set()
for name, subj_set in students_set.items():
    if subj_set >= colleges_set['leetademy']:
        output.add(name)
print(output)
แต่การทำงานบน database เรามักเขียนพวกนี้ทั้งหมดให้เป็น one-liner เพื่อที่จะรับประกันว่าระหว่างการ query จะไม่เกิด race condition กับ query อื่นๆ สิ่งที่ต้องทำก่อนอื่นใดคือพยายามกำจัดการสร้าง set ของข้อมูลล่วงหน้าซะ

ในที่นี้ ขั้นแรกเราจะกำจัดการสร้าง colleges_set โดยวน loop ใน colleges ทุกตัวดูว่าแต่ละ subj มีใน students_set หรือเปล่า
def learnt_all(name):
    for campus, subj in colleges:
        if campus == 'leetademy':
            if subj not in students_set[name]:
                return False
    return True

output = set()
for name, _, _ in students:
    if learnt_all(name):
        output.add(name)
print(output)
ขั้นต่อมาก็กำจัด students_set โดย loop เข้า students ดูว่านักเรียนที่ชื่อ name ได้เรียนวิชา subj แล้วหรือยัง
def learnt(name, subj):
    for student_name, student_subj, _ in students:
        if student_name == name and student_subj == subj:
            return False
    return True

def learnt_all(name):
    for campus, subj in colleges:
        if campus == 'leetademy':
            if learnt(name, subj):
                return False
    return True
    
output = set()
for name, _, _ in students:
    if learnt_all(name):
        output.add(name)
print(output)
เนื่องจากในฟังก์ชัน learnt กับ learnt_all จะคืนค่า False ทันทีที่เข้าเงื่อนไขครบ แต่จะคืนค่า True หลังจากจบ loop เท่านั้น เราสามารถเขียนตรงนี้ใหม่ด้วย generator แล้ว wrap ด้านนอกด้วย not any(...) เช่นนี้
def inverse_learnt(name, subj):
    for student_name, student_subj, _ in students:
        if student_name == name and student_subj == subj:
            yield True

def inverse_learnt_all(name):
    for campus, subj in colleges:
        if campus == 'leetademy':
            if not any(inverse_learnt(name, subj)):
                yield True

output = set()
for name, _, _ in students:
    if not any(inverse_learnt_all(name)):
        output.add(name)
print(output)
พอทุกอย่างเป็น single-statement และ generator คราวนี้ก็ง่ายที่จะทำ one-liner แล้ว ... จับมาเขียนสวยๆ ได้แบบนี้
print({
    name
    for name, _, _ in students
    if not any(
        True
        for campus, subj in colleges
        if campus == 'leetademy'
        and not any(
            True
            for student_name, student_subj, _ in students
            if student_name == name
            and student_subj == subj
        )
    )
})
เห็น code นี้แล้วคงเดาได้ไม่ยากว่า พอเอาเขียนเป็น SQL คงหนีไม่พ้น
SELECT DISTINCT name
FROM   students AS s
WHERE  NOT EXISTS (
    SELECT *
    FROM   colleges AS c
    WHERE  campus = 'leetademy'
    AND    NOT EXISTS (
        SELECT *
        FROM   students AS ss
        WHERE  ss.name = s.name
        AND    ss.subj = c.subj
    )
)

อ้างอิง


Apr 14, 2013

Code Jam 2013 รอบคัดเลือก


เนื่องจากไม่อยากให้เกิดเหตุการณ์แบบรอบ 1 ของปีที่แล้ว ปีนี้เลยนอนเอาแรงไม่บ้าพลังถ่างตาทำตั้งแต่เช้ามืด ทำให้กว่าจะได้เริ่มอ่านโจทย์ก็บ่าย 2 เข้าไปแล้ว

ปีนี้ความตั้งใจแรกคือจะพยายามเขียนส่งให้ได้หลายๆ ภาษา คือฝึก ML, Lisp, Bash อะไรพวกนี้จนคิดว่าน่าจะคล่องพอเอามาใช้แก้โจทย์ปัญหาจริงๆ ได้แล้ว น่าเสียดายที่รับ input ไม่เป็น เลยกลับมาตายรังด้วย Python เหมือนเดิม (แต่ก็มี Haskell โผล่มาครึ่งข้อ) ไม่แน่ใจว่ารอบต่อๆ ไปจะเขียน C++ ทันหรือเปล่า เพราะโดนจำกัดเวลาเหลือแค่ 2.5 ชั่วโมงแล้ว

ข้อแรกโจทย์ extended เกม OX ให้ใหญ่ขึ้นเป็นตาราง 4x4 แถมยังอนุญาตให้มี wildcard (เป็นได้ทั้ง X และ O) อย่างมาก 1 ที่บนกระดาน ก็ถามว่าเกมตอนนี้ใครชนะ หรือว่าเสมอ หรือว่ายังเล่นไม่จบ

เคยเจอคำถามเกี่ยวกับ grid แบบนี้มาพอสมควร อันที่จริงก็ต้องไล่เช็คทั้งแนวตั้ง แนวนอน และก็แนวทะแยงซ้ายขวาเลย แต่เพิ่งมาสำนึกได้ว่าเขียนแบบนั้นมี code ซ้ำกันเยอะๆ จะเริ่ม maintenance ยากละ เลยสังเกตว่าถ้าหมุน grid ให้ได้เป็นแบบนี้

1 2 3      7 4 1
4 5 6  =>  8 5 2
7 8 9      9 6 3

จะสามารถ reuse code ส่วนที่เช็คแนวนอนมาเช็คแนวตั้งได้ด้วย เช่นเดียวกับ code เช็คแนวทะแยง การ implement แค่ส่วนนี้คงมีวิธีเจ๋งๆ ให้เลือกใช้เยอะอยู่ ส่วนผมที่มาจากสาย functional เลือกใช้ transpose . reversed ครับ

def transpose(grid):
    return [''.join(line) for line in zip(*grid)]

def chk_hori(grid):
    for line in grid:
        if all(c in 'XT' for c in line):
            return 'X won'
        if all(c in 'OT' for c in line):
            return 'O won'

def chk_diag(grid):
    if all(grid[i][i] in 'XT' for i in range(len(grid))):
        return 'X won'
    if all(grid[i][i] in 'OT' for i in range(len(grid))):
        return 'O won'

def chk_won(grid):
    for _ in range(2):
        for chk in [chk_hori, chk_diag]:
            answer = chk(grid)
            if answer:
                return answer
        grid = transpose(reversed(grid))

def chk_full(grid):
    return not any('.' in line for line in grid)

for case in range(int(input())):
    grid = [input() for _ in range(4)]
    input()
    answer = chk_won(grid)
    if not answer:
        answer = 'Draw' if chk_full(grid) else 'Game has not completed'
    print('Case #{}: {}'.format(case+1, answer))



ข้อถัดมาถามว่าสามารถใช้เครื่องตัดหญ้าอัตโนมัติ (ที่ดันวิ่งตรงไปข้างหน้าได้อย่างเดียว) ตัดหญ้าให้แต่ละช่องมีความสูงตาม pattern ที่วางไว้ได้หรือเปล่า

พบว่า algorithm ง่ายกว่าที่คิด คือดูว่าแต่ละจุดบนสนามหญ้านั้น ต้องมีความสูงเท่ากับความสูงหญ้าที่สูงที่สุดในแนวนอนหรือแนวตั้งเท่านั้น เพราะถ้าเกิดมีจุดที่สูงกว่าทั้งแนวนอนและแนวตั้งพร้อมกัน จะทำให้ไม่สามารถตัดหญ้าในช่องนั้นให้ต่ำกว่าค่าสูงสุดได้

def transpose(lawn):
    return [list(line) for line in zip(*lawn)]
    
def beautiful_garden(lawn, n, m):
    tlawn = transpose(lawn)
    def chk(y, x):
        return lawn[y][x] in [max(tlawn[x]), max(lawn[y])]
    return all(all(chk(y, x) for x in range(m)) for y in range(n))
            
for case in range(int(input())):
    n, m = [int(i) for i in input().split()]
    lawn = [[int(i) for i in input().split()] for _ in range(n)]
    answer = 'YES' if beautiful_garden(lawn, n, m) else 'NO'
    print('Case #{}: {}'.format(case+1, answer))



ข้อที่ 3 ให้ว่ามีเลข palindrome ที่รากที่สองของมันก็ยังเป็น palindrome ด้วยทั้งหมดกี่ตัวในช่วงที่กำหนด ซึ่งมีความรู้สึกว่าถ้าเล่นกับ palindrome แล้ว Haskell จะสวยงามมาก แต่ก็ไปตายที่ test กลางเลยกลับไปใช้ Python แทน

import Data.List.Split (splitOn)

boolAsNum b = if b then 1 else 0

isSquare x = root^2 == x
    where root = sqrt x

isPalindrome x = y == reverse y
    where y = show x

countFair x y
    | x > y     = 0
    | otherwise = (boolAsNum $ all isPalindrome [x,x^2]) + countFair (x+1) y

eachLoop nosLoop = do
    raw <- getLine
    let rawSqrtNum = [sqrt $ read x | x <- splitOn " " raw]
        [start, stop] = [f x | (f,x) <- zip [ceiling, floor] rawSqrtNum]
    putStrLn $ "Case #" ++ (show nosLoop) ++ ": " ++ (show $ countFair start stop)

main = do
    allLoop <- getLine
    sequence_ [eachLoop n | n <- [1..read allLoop]]

ส่วนนี่คือ Python ที่ optimize ไปนิดหน่อยสำหรับความยากปานกลาง โดยส่วนที่วนสร้างเลข palindrome ตัวถัดไปยังไม่ได้ optimize สำหรับโจทย์นี้โดยเฉพาะเลย (ว่าแล้วก็ต้องเขียนเก็บเข้า lib ตัวเองซะแล้ว)

odd = lambda x: x % 2 == 1
square = lambda n: int(n ** 0.5) ** 2 == n
palindrome = lambda n: str(n) == str(n)[::-1]

def int_sqrt(n):
    if n == 0: 
        return 0
    a, b = divmod(n.bit_length(), 2)
    x = 2 ** (a + b)
    while True:
        y = x + n // x 
        y //= 2
        if y >= x: 
            return x
        x = y

def next_palindrome(n):
    s = str(n)
    size = len(s)

    fst = s[:size//2+1] if odd(size) else s[:size//2]
    lst = s[-size//2:]

    size_old_fst = len(fst)
    if fst <= lst[::-1]:
        fst = str(int(fst) + 1)

    lst = fst[:-1] if len(fst) > size_old_fst else fst
    return int(fst[:-1] + lst[::-1]) if odd(size) else int(fst + lst[::-1])

def palindrome_range(start, stop):
    while start < stop:
        if palindrome(start):
            yield start
        start = next_palindrome(start)

for case in range(int(input())):
    start, stop = [int(n) for n in input().split()]

    start = int_sqrt(start) + (0 if square(start) else 1)
    stop = int_sqrt(stop)

    answer = sum(palindrome(pal**2) for pal in palindrome_range(start, stop+1))
    print('Case #{}: {}'.format(case+1, answer))



ส่วนข้อ 4 ไม่ได้ทำ เพราะสังเกตเห็น pattern ของข้อ 3 เลยกะจะแก้ความยากระดับโหดสุด (input ใหญ่ได้ถึง 10^100) แต่ลอง optimize ดูแล้วทำทันแค่ 10^80 เท่านั้น ก็เลยต้องตัดใจไปตามระเบียบ แต่แต้มแค่นี้ก็เพียงพอสำหรับรอบ 1 แล้วครับ ;)

Apr 7, 2013

Bit Flag

ในภาษาสมัยใหม่ เราคงคุ้นเคยกับการเขียนฟังก์ชั่น (และเรียกใช้) แบบนี้

def foo(n, m, a=False, b=False, c=False):
    if a:
        n += 42
    if b:
        n *= 7
    if c:
        n **= 2
    return n % m

...

foo(5, 99)
foo(5, 99, a=True)
foo(5, 99, c=True)
foo(5, 99, c=True, b=True, a=True)

โชคร้ายที่ระบบ keyword argument ไม่ได้เป็นแบบนี้ทุกภาษา อย่างเช่นใน PHP ถ้าจะเปิดตัว flag c เพียงตัวเดียว ก็ยังคงต้องบอกว่า a, b ถูกปิดอยู่ด้วย (ต้องบอกสถานะของ flag ทุกตัวที่อยู่ก่อนหน้า c)

ทางออกง่ายๆ แต่ทำให้ syntax ดูรุงรังหน่อยคือการให้ argument ส่วนที่เป็น flag ใช้ datatype แบบ hash table เช่นนี้

function foo($n, $m, $options=array()) {
    if (array_key_exists('a', $options) && $options['a'])
        $n += 42;
    if (array_key_exists('b', $options) && $options['b'])
        $n *= 7;
    if (array_key_exists('c', $options) && $options['c'])
        $n *= $n;
    return $n % $m;
}   

...

foo(5, 99);
foo(5, 99, array('a' => true));
foo(5, 99, array('c' => true));
foo(5, 99, array('c' => true, 'b' => true, 'a' => true));

ทางออกที่คลาสสิก (แต่ปวดหัว) กว่าคือการใช้ bit flag แบบนี้

const int A = 0x01;
const int B = 0x02;
const int C = 0x04;

int foo(int n, int m, int options) {
    if (options & A)
        n += 42;
    if (options & B)
        n *= 7;
    if (options & C)
        n *= n;
    return n % m;
}

...

foo(5, 99, 0);
foo(5, 99, A);
foo(5, 99, C);
foo(5, 99, C|B|A);

หลายคนอาจบอกว่าการใช้ bit flag นั้นมีข้อดีเหนือกว่า keyword argument ตรงที่สามารถ store สถานะของ flag เก็บไว้ใน data ตัวหนึ่ง (เพื่อที่จะนำไปใช้ซ้ำในที่อื่น) ซึ่งข้อได้เปรียบนี้ไม่เป็นความจริงเลยสำหรับภาษาที่แตก array, hash table ลงไปเป็น argument ได้

args = [5, 99]
flags = {'b': True, 'c': True}

foo(*args, **flags) == foo(5, 99, b=True, c=True)

ก็เลือกใช้ความสามารถเหล่านี้ตามความสะดวกของแต่ละภาษานะครับ อย่าไปฝืนมันมากถ้าตัวภาษามันไม่ได้ออกแบบให้ทำเรื่องยากๆ บางเรื่องไว้แต่แรก

Mar 2, 2013

Swap ตัวแปร

โปรแกรมที่เล็กและง่ายรองจาก hello world ที่โปรแกรมเมอร์หลายคนต้องผ่านตามาบ้าง คือโปรแกรม (ฟังก์ชัน) สำหรับสลับค่าตัวแปรนั่นเอง

แต่ก่อนจะลงลึกที่รายละเอียด มาดู simple yet best practice ใน Pyhton กันก่อน :P

a, b = b, a



อธิบายก่อนว่าการเข้าถึง data ในคอมพิวเตอร์ จะคล้ายๆ กับคนที่มีแขนข้างเดียว คือหยิบของได้ทีละ 1 อย่าง ไม่สามารถหยิบของพร้อมกัน 2 อย่างแล้ววางสลับที่กันทันทีได้เหมือนในชีวิตจริง

วิธีแก้ปัญหาที่ตรงไปตรงมาภายใต้ข้อจำกัดนี้ คือหยิบของชิ้นแรกไปวางไว้ในที่ว่างก่อน หยิบของชิ้นที่สองไปวางแทนที่ๆ ของชิ้นแรกเคยอยู่ แล้วค่อยกลับไปหยิบของชิ้นแรกวางแทนที่ๆ ของชิ้นสองเคยอยู่

จะอธิบายให้คอมพิวเตอร์เข้าใจได้ ก็ต้องบอกด้วยว่าที่ว่างนั้นคือตรงไหน (พูดลอยๆ ไม่ได้ เดี๋ยวกลับไปหาไม่เจอ) เช่นนี้

int t;

t = a;
a = b;
b = t;

ปัญหาของการย้ายแบบนี้ คือมันเป็นการมองในเชิง object โลกมนุษย์ ซึ่งไม่ได้ใกล้เคียงกับกิจกรรมที่เกิดขึ้นใน memory เลย มันไม่ใช่การย้ายของด้วยซ้ำ แต่เป็นการทำสำเนาสิ่งของแล้วย้ายที่สลับไปมาต่างหาก

อนึ่ง วิธีนี้ยังมีข้อเสียตรงที่มันใช้พื้นที่ว่างเพิ่มเติมสำหรับเก็บข้อมูลระหว่างทางอีก แล้วมันจะมีทางมั้ยที่จะสลับที่ตัวแปรโดยไม่ต้องอาศัยตัวแปรอื่นเข้ามาทำหน้าที่เป็นที่พักข้อมูล?

ลองนั่งคิดเล่นๆ ซักพักคงพบว่ามันก็ไม่ได้ยากอะไร เริ่มจากเอาค่าของตัวแปรสองไปเก็บรวมกับตัวแปรแรก แล้วเอาตัวแปรแรก (ที่ตอนนี้มีค่าของตัวแปรแรกและตัวแปรสองรวมกัน) กลับไปหักลบกับตัวแปรสอง เท่านี้ตัวแปรที่สองก็จะมีค่าเท่ากับตัวแปรแรกแบบดั้งเดิมแล้ว (ขั้นที่เหลือคงไม่ต้องบอกว่าทำยังไงต่อ)

a = a + b;
b = a - b;
a = a - b;

แบบนี้ก็เกือบดีแล้ว แต่ยังมีข้อควรระวังคือ overflow เนื่องจากมันเป็นวิธีทางคณิตศาสตร์ แถมโปรแกรมแบบนี้ก็ยังไม่ค่อยสวยงามเท่าไหร่ด้วย

ถ้าจะมองให้เป็นไปในทางคอมพิวเตอร์จริงๆ ต้องเปลี่ยนการดำเนินการ +, - ไปใช้ xor แทน เช่นนี้

a ^= b;
b ^= a;
a ^= b;

ซึ่งสามารถลดรูปได้อีกขั้นจนเหลือเพียงบรรทัดเดียวว่า

a ^= b ^= a ^= b;

อย่างไรก็ตาม การมองว่าคอมพิวเตอร์เปรียบเหมือนคนที่มีแขนข้างเดียวนั้น ก็คงไม่ถูกต้องอีกต่อไปแล้ว (เป็นการมองแบบ Turing machine ซึ่งตกยุคไปชาติเศษๆ) เพราะในปัจจุบัน คอมพิวเตอร์มีแขนเหล่านี้เพิ่มมากขึ้น และสามารถเข้าถึงข้อมูลได้หลายตำแหน่งในเวลาเดียวกัน แต่การเขียนโปรแกรมให้ทำงานแบบ parallel ให้ได้ดีนั้น ก็ไม่ใช่เรื่องง่ายซักเท่าไหร่นัก หน้าตาของโปรแกรมอาจจะออกมาเช่นนี้

process [1, 2]:
  goto address [a, b]
  read into process memory
  sync with process [2, 1]
  goto address [b, a]
  write from process memory



ส่วนภาษาเชิง functional ที่ไม่ยอมให้มี mutable อย่างการสลับเปลี่ยนค่าตัวแปร ก็ยังสามารถทำท่านี้ได้โดยใช้เทคนิคตั้งชื่อตัวแปรสลับที่กันใน environment ที่ต่างกัน (ภาษาเค้าเรียกว่า Monad) ต่อไปนี้คือตัวอย่างใน Haskell

do (a,b) <- return (b,a)
   someMagic a b
   ...

หรือแบบนี้ใน Lisp (สังเกตการใช้คำสั่ง let เพื่อสลับตัวแปร)

(let ([a b]
      [b a])
  (magic-happens a b ...))



สุดท้ายนี้กลับไปดู Python ที่เรียบง่ายอีกรอบ แท้จริงแล้วมันคือการสั่ง

container = b, a
a, b = container

แปลได้ว่า แรกสุดมันจะเก็บค่าตัวแปร b และ a ตามลำดับลง container ชนิดหนึ่ง (ใน Python มันคือ tuple) แล้วหลังจากนั้นจึงดึงตัวแปรออกจาก container นี้มาให้ตัวแปร a และ b ตามลำดับ แนวคิดแบบนี้ค่อนข้าง general มาก และสามารถนำไปใช้กับกรณีที่มีตัวแปรให้สลับเยอะกว่านี้ได้ด้วย เช่น

a, b, c, d = d, -b, -c, a

Feb 25, 2013

ลำดับ Farey และ วงกลม Ford

ปีที่แล้วเล่น Project Euler ไปได้หลายข้อ แล้วก็สะดุดตากับข้อที่น่าสนใจ เลยเอามาขยายผลต่อจนได้เป็นเรื่องนี้ขึ้นมา

อธิบายสั้นๆ คือลำดับ Farey เป็นลำดับของเศษส่วนอย่างต่ำตั้งแต่ 0/1 ไปจนถึง 1/1 เพราะฉะนั้นตัวเลขอย่าง 2/4 จึงไม่นับ (หรือถ้าจะนับ ก็ทำให้มันเป็น 1/2 ก่อน)

อย่างไรก็ดี ถ้านิยามอย่างนี้เราจะมีปัญหาว่าลำดับ Farey มันดันเป็นลำดับอนันต์ซะหนิ ดังนั้นเลยเพิ่มข้อจำกัดไปว่า ลำดับ Farey อันดับ n จะบรรจุเศษส่วนอย่างต่ำ ที่ตัวส่วนมีค่าน้อยกว่าเท่ากับ n เท่านั้น

พอมันเป็นลำดับจำกัดแล้ว ก็สามารถเพิ่มข้อกำหนดสุดท้ายได้ว่า สมาชิกแต่ละตัวในลำดับนี้ ต้องเรียงค่าจากน้อยไปมากด้วย

ส่วนวงกลม Ford ก็เป็นวงกลมที่นำเอาลำดับ Farey ไปวาดแสดงผลได้อย่างสวยงาม ;)

รายละเอียดทั้งหมดอยู่ในเอกสารฉบับนี้ ไม่ต้องกลัวว่าเป็นภาษาคณิตศาสตร์แล้วจะเข้าใจยาก เพราะมีโปรแกรมที่เขียนด้วย Python ให้ไปแกะโค้ดเล่นกันด้วย

ปล. blog ตอนนี้ เขียนเป็นพิเศษให้ @flurrywong @FordAntiTrust ฮะ <3

Jan 2, 2013

ทำ Optional Argument แบบ LaTeX

command ใน LaTeX นั้นมีวิธีเรียกใช้ที่ต่างจากภาษาโดยทั่วไปอยู่บ้าง ลองดู
\sqrt{x}
\sqrt[3]{x}
อันด้านบนจะวาด √x (รากที่สองของ x -- โดยปรกติ รากที่สองไม่จำเป็นต้องมีเลข 2 กำกับ) ส่วนอันล่างจะวาด ∛x

รูปร่างของ syntax ที่น่าสนใจคือ ส่วนที่เป็น optional argument จะถูกเรียกขึ้นมาก่อน main argument แถมยังใช้วงเล็บแยกกันอีก นี่ทำให้การ currying เป็นไปได้อย่างง่ายดาย
\newcommand{\cbrt}{\sqrt[3]}
...
\cbrt{x}
จริงๆ แล้วทาง math ก็มีฟังก์ชันในแนวคิดนี้อยู่เยอะ อย่างเช่น σx(n) ซึ่งเขียนแบบนี้แล้วเข้าใจง่ายกว่า σ(n, x) และยังไม่สับสนอีกว่า argument ตัวไหนที่ควรเป็น n หรือ x กันแน่

แล้วถ้าอยากได้แบบนี้ใน Python บ้าง? ง่ายนิดเดียวเพราะใช้เทคนิคเดียวกับ Infinite List นั่นเอง
@singleton
class sigma:
    def __call__(self, n, x=1):
        if n == 1:
            return 1
        return product( sum((k**x)**i for i in range(v+1))
                            for k, v in Counter(factor(n)).items() )
    def __getitem__(self, x):
        def partial_sigma(n):
            return self(n, x)
        return partial_sigma
คราวนี้จะเขียน
6 == sigma[1](6) - 6
หรือกระทั่ง
tau = sigma[0]
...
[tau(i) for i in range(return 1, 10)]
ก็ย่อมได้

Dec 14, 2012

9/3(2+1) = ?

เวลามีคนถามว่า

9/3(2+1) = ?

จะตอบโดยไม่ลังเลเลยว่า 1 แถมถ้าใครมาบอกว่าต้องคิดหารทางด้านซ้ายมือ แล้วค่อยมาคูณแบบคอมพิวเตอร์นะ จะตอบกลับไปเลยว่ามัน syntax error โว้ย ภาษาคอมพิวเตอร์ที่ไหนเค้ายอมรับการเขียนเลขติดกันว่าเป็นการคูณเลขบ้าง?

อนึ่ง ถ้ามองว่าการเขียนเลขติดกัน (คั่นด้วยวงเล็บ) เป็นการคูณจริงๆ ตัวที่ติดกับวงเล็บก็ต้องมี precedence สูงกว่า (ทำงานก่อน) เครื่องหมายหารอยู่แล้ว นั่นหมายความว่าต้องคิดส่วน 3(2+1) ให้เสร็จก่อนเอาไปคำนวณต่ออยู่ดี

เจอบ่อยๆ ก็เริ่มขี้เกียจอธิบาย เลยเขียนโปรแกรมที่มันแปลสมการด้านบนนี้ได้ถูกต้องเลยละกัน -> จิ้มโลด

ตัวอย่างผลลัพท์จากโปรแกรม

>>> 9/3(2+1)
1.0
>>> a, b, c, d = 48, 2, 9, 3
>>> a/b(c+d)
2.0
>>> 10(9(8(7(6(5(4(3(2(1)))))))))
3628800
>>> (1+2j)(3-4j)
(11+2j)
>>> 1 + 1/phi
1.618033988749895
>>> phi(43)
42

อันที่จริงก็คิดว่าจะทำมาได้ซักพักละ แต่ด้วยความที่ขี้เกียจวางกฎ grammar ด้วย flex bison ทั้งที่หลายๆ กฎนั้นภาษาส่วนใหญ่ก็ทำไว้ให้แล้ว เลยดองมาเรื่อยๆ จนทำ Infinite List เสร็จ และไปเจอ Python's Magic Methods เลยได้ฤกษ์เริ่มจากตรงนั้น

แนวคิดก็ง่ายนิดเดียว แค่เพิ่มเมธอด __call__ เข้าไปให้ class ของตัวเลข (เป็น Monkey Patching เหมือนใน Ruby) โดย define มันกว่าถ้าเรียก a.__call__(b) แล้วจะมีค่าเท่ากับเอามันมาคูณกันซะ

ปัญหาที่ต้องไล่เก็บก็คือ
  1. Python ไม่อนุญาติให้แก้ไขพวก builtins -- ที่จริงปัญหานี้ก็เหมือนกับคราวที่ทำ Infinite List นั่นแหละ ซึ่งแก้ไม่ยาก แค่ inherit class นั้นๆ ออกมาแล้วไล่ define magic method ซะ
  2. ปัญหาจริงๆ คือ Python parser มันจะแปลตัวเลขไปเป็น class ตัวเลขใน builtins เท่านั้น ดังนั้นก็ต้องแปลง token เองโดยหาว่าตรงไหนคือเลข แล้วก็เอา class ตัวเลขที่ inherit มาแล้วครอบไว้
  3. งานนี้ได้โมดูล tokenize มาช่วยแปลตัวเลขให้ ไม่ต้องเขียน regex เอง เพราะมันจะมีท่าประกาศตัวเลขประหลาดๆ ทำให้ tokenize เองพลาดได้
  4. ที่เหลือก็เอา code ที่แก้ token เรียบร้อยไปทำงาน ตรงนี้ดึงโมดูล code ที่ทำทุกอย่างเตรียมไว้แล้วมาใช้ฟังก์ชันเดียวจบ
  5. แต่ตอนที่เอา tokenize มาใช้ร่วมกับ code จะมีปัญหาที่พิมพ์ได้ทีละบรรทัด เนื่องจาก tokenize มันปรับแต่งอะไรไม่ได้เลย ก็ต้องหลอกมันเอานิดหน่อย
  6. ครบถ้วนกระบวนความแล้วก็ clean up โดยจุดที่เอา code ออกไปได้เยอะมากๆ คือส่วน define magic method อันนี้ใช้ meta-programming บอกแค่ว่าจะเพิ่ม method ไหนบ้างก็พอ แล้วปล่อยให้มัน generate ตัวเองไปซะ
  7. อยู่ๆ ก็เขียน function decorator เป็นเฉยเลย งงตัวเอง 555+
อย่างไรก็ตาม แนวคิดของ callable number นี้ไม่ควรเอาไปใช้เขียนโปรแกรมในโลกจริง เพราะนักคณิตศาสตร์ชอบตังชื่อเดียว แต่ใช้ในคนละบริบท เช่น φ (phi) ที่อาจหมายถึง golden radio หรือ Euler's totient ก็ได้ ซึ่งแม้ว่าทั้ง 2 รูปของ φ นี้จะไม่ overlap กัน (เช่น 1+1/φ จะถูกมองว่าเป็นตัวเลข ในขณะที่ φ(43) คือการใช้ฟังก์ชัน) แต่ด้วยความที่ Python เป็น 1st class function ที่ยอมให้เราส่งผ่านฟังก์ชันเป็นตัวแปรของฟังก์ชันอื่นๆ ได้ ดังนั้นถ้าเราเจออะไรอย่าง δ(φ,ε) เราจะไม่สามารถแน่ใจได้เลยว่า φ ในบริบทนี้คืออะไรครับ

Oct 17, 2012

โปรแกรมที่เขียนมันช้า...

เวลาเขียน Python แล้วต้อง implement อะไรประมาณนี้
numbers = [random() for r in range(1000)]

for i in range(100):
    print(max(numbers)**i)
จะระลึกตัวเสมอว่า implement ตรงๆ แบบนี้แล้วมันช้านะ แต่ถ้า input ไม่ได้ใหญ่จนน่าเป็นห่วง ก็ไม่มีเหตุผลที่จะเขียนให้ยากขึ้น เพราะการแก้โดยเอา max(numbers) ไปไว้เป็นค่าคงที่นอก loop แม้จะทำให้โปรแกรมเร็วขึ้นก็ตาม แต่ก็จะเสีย readability ไปบางส่วน เวลากลับมาอ่านอาจจะสับสนใน logic ว่าส่วนนี้มีไว้ทำอะไรกันแน่ แถมถ้าภายหลังโปรแกรมถูกพัฒนาจนซับซ้อนมากๆ แล้ว ค่า max(numbers) อาจไม่เท่ากันทุกรอบที่วน loop ก็ได้ ทำให้เกิดเป็น bug ซ่อนเร้นอีก

ส่วนถ้าต้องการประสิทธิภาพที่ดีกว่าเดิมมากๆ การแก้ Python จนอ่านยากอาจไม่ใช่ความคิดที่เข้าท่าเท่าไหร่ ลองถอดโปรแกรมเดิมมาเขียนบน C++ แบบเป๊ะๆ ก็คงได้ประมาณนี้
int main(void) {
    float numbers[1000];
    fill_with_random(numbers);

    for (int i=0; i<100; i++) {
        printf("%.32f\n", pow(max(numbers), i));
    }
    return 0;
}
(ตัวเลขขนาด array / loop ควรเปลี่ยนไปยังค่าสูงกว่านี้ เพื่อให้เห็นได้ชัดว่าโปรแกรมทำงานได้ช้าลง)

ถึงตอนนี้จะเห็นว่า max(numbers) มีผลอย่างมากต่อประสิทธิภาพโปรแกรม การย้ายมันออกมาไว้นอก loop จะเป็นการกระทำที่มีเหตุผลขึ้นทันที

แต่ถ้าโปรแกรมเรามีอะไรแบบนี้เป็นสิบเป็นร้อยที่หละ การย้ายมือเองไม่สนุกแน่

ทางออกก็ง่ายๆ บอกให้ gcc มันทำแทนเราไปซะ
$ gcc -O2 exp-max.cpp
เพียงเท่านี้ก็จะได้ code ที่ยังคง logic + readability เดิมเอาไว้อยู่ พร้อมทั้งโปรแกรมที่เร็วกว่าเดิมเยอะ

แต่ถ้าโปรแกรมใหญ่มากๆ ก็ต้องรอ compile กันนานหน่อยนะ :P



Sep 11, 2012

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)

Sep 3, 2012

แก้โจทย์ป้ายสมัครงาน 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 นะครับ ^^

Aug 30, 2012

ทำ Infinity List ใน Python

พอดีไปขุด Project Euler มาเล่น แล้วมันต้องได้ยุ่งกับพวก infinity sequence อย่างจำนวนเฉพาะ, ฟีโบนักชีบ่อยๆ ตอนแรกก็ใช้แค่ list ธรรมดาเนี่ยแหละ แล้วถ้า index ที่อยากได้ไม่มีก็ค่อยไปเรียก function แยกมาคำนวณแล้วเก็บข้อมูลกลับเข้า list เอา

ทีนี้เขียนไปเขียนมารู้สึกว่า code มันไม่สวยเอาซะเลย เพราะไอ้จำนวนเหล่านี้เนี่ย เวลาหามันจะกระโดดข้ามไป index ที่ต้องการไม่ได้อยู่แล้ว บวกกับจำนวนพวกนี้มันหาเพิ่มได้เรื่อยๆ ไม่มีที่สิ้นสุด ดังนั้นในสายตาของนักคณิตศาสตร์อย่างผม ถ้าเขียน prime[100] แล้วเจอฟ้องว่า index out of range นี่คงเซ็งน่าดู เลยคิดว่ามันน่าจะมีวิธี hack ให้ใช้ syntax นี้หาค่าใน index ที่ต้องการแม้ต้อนนี้จะยังไม่มีได้

ค้นไปค้นมาก็เจอ list.__getitem__ เลยจัดการลงมือลุย
class InfinityList(list):
    def __iter__(self):
        n = 0
        while True:
            yield self[n]
            n += 1

    def __repr__(self):
        return super(InfinityList, self).__repr__()[:-1] + ', ...]'


class fibonacci(InfinityList):
    def __getitem__(self, n):
        while True:
            try:
                return super(InfinityList, self).__getitem__(n)
            except IndexError:
                self.append( super(InfinityList, self).__getitem__(-1) +
                             super(InfinityList, self).__getitem__(-2) )

    def __init__(self):
        super(InfinityList, self).__init__([1, 1])


class prime(InfinityList):
    def __getitem__(self, n):
        while True:
            try:
                return super(InfinityList, self).__getitem__(n)
            except IndexError:
                c = super(InfinityList, self).__getitem__(-1)
                while True:
                    c += 2
                    for p in self[:]:
                        if not c % p:
                            break
                    else:
                        self.append(c)
                        break

    def __init__(self):
        super(InfinityList, self).__init__([2, 3])

fibonacci = fibonacci()
prime = prime()
กรอเวลาไปข้างหน้า 8 ชั่วโมง ก็มีของเล่นใหม่ตามที่แปะเป็นตัวอย่างไว้ด้านบน
  1. ใน __getitem__ นี่ลืมไปว่าถ้าแก้ตรงนี้แล้วมันจะทำ recursion กับตัวเอง ทำให้ไม่สามารถใช้ syntax อย่าง fibonacci[-1] เพื่อดึงเอาตัวเลขออกมาได้ ก็เลยต้องเลี่ยงไปเรียก method จาก super class เอา
  2. จะให้ syntax แบบ list[n] ทำงานโดยเรียก __getitem__ ต้อง inherit class มาแล้ว define ไว้ตอนสร้าง class เท่านั้นด้วย มาสั่ง class.__getitem__ = outer_function ไม่ได้ (ป้องกัน injection เข้าระบบ)
  3. ตัว __getitem__ มันรับ parameter ได้อันเดียวก็จริง แต่ไม่ใช่แค่ตัวเลขเท่านั้น เพราะยังมี slice object อีกด้วย (นึกถึงเวลาเขียน some_list[4:-4] ไอ้ [4:-4] นั่นแหละ slice object) ดังนั้นงานนี้ play safe ใส่ while True ไปดีกว่า ช้าหน่อยยอมรับได้
  4. สุดท้ายก็คืออยากให้ class พวกนี้เป็น static ไปซะ (ตอนแรกออกแบบว่าจะให้เขียนเป็น prime.numbers[n] ด้วยซ้ำ) แต่เหนื่อยมากขี้เกียจทำแล้ว เลยจับประกาศชื่อ fibonacci, prime ทับกับชื่อ class ตัวมันเองไปซะเลย ยังไงก็มีได้แค่ instance เดียวอยู่แล้วหนิ 555+
เพียงเท่านี้ ถ้าอยากได้ prime ตำแหน่งที่ 101 ก็แค่สั่ง
prime[100]
จบเลย ง่ายมั้ย? อยากลองเล่นเองแล้ว? ไปจิ้ม code ได้จากที่นี่เลย XD

Aug 22, 2012

อะไรคือ Y Combinator

ถึงตอนนี้เราคงทำ recursion บน lambda เป็นกันแล้ว (ถ้ายังไม่เป็น แนะนำให้อ่านตอนก่อน)

กลับไปดูอีกรอบ จะเห็นว่าตอนเรียก recursion ต้องใส่ชื่อฟังก์ชันตัวเองลงไปเป็นตัวแปรหนึ่งเสมอ
lambda f, x: 1 if x == 0 else x * f(f, x-1)
แล้วถ้าเราไม่อยากทำอย่างนี้หละ อยากจะเรียกแค่ f(x-1) เหมือนการเรียกใช้งานฟังก์ชันตามปรกติทั่วไป จะทำได้มั้ย?

ถอยหนึ่งก้าว กลับไปเขียน factorial แบบธรรมดาๆ ก่อน
f = lambda x: 1 if x == 0 else x * f(x-1)
ฟังก์ชันนี้ยังคงสามารถเรียกใช้ได้ เนื่องจากว่าเราจองชื่อ f ไว้ให้มันแล้ว (ซึ่งไม่ใช่สิ่งที่เราต้องการ) นั่นหมายความว่าชื่อ f ต้องมีปรากฏอยู่ใน scope นี้จึงจะทำ recursion ได้ และมันมีอีกวิธีนึงที่จะทำให้ f ปรากฏอยู่ภายใน scope นี้โดยไม่ต้องจองชื่อหรือส่งผ่านเป็นตัวแปรใน scope เดียวกัน คือ
lambda f: lambda x: 1 if x == 0 else x * f(x-1)
แต่การประกาศเช่นนี้จะทำให้เราไม่สามารถใช้ฟังก์ชันได้ทันที เพราะเราต้องเอาฟังก์ชั่นนี้ส่งไปเป็นตัวแปรให้ตัวมันเองเรื่อยๆ เช่น ถ้าต้องการหา 5 factorial ต้องเขียนออกมาอย่างนี้
(lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
  (lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
    (lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
      (lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
        (lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
          (lambda f: lambda x: 1 if x == 0 else x * f(x-1))(
            (lambda whatever: 42)
          )
        )
      )
    )
  )
)(5)
นั่นหมายความว่าเราต้องการอะไรซักอย่าง ที่จะทำหน้าที่ generate ฟังก์ชันของเรานี้ออกมาเรื่อยๆ ไม่ว่าจะมีอะไรเกิดขึ้นก็ตาม นี่คือแนวคิดของ fixed-point combinator โดยมีตัวหนึ่งที่เป็นมีชื่อเสียงเป็นที่รู้จักกันมากคือ Y combinator ซึ่งเขียนได้ดังนี้
Y = lambda f: (lambda x: f(lambda v: x(x)(v)))(lambda x: f(lambda v: x(x)(v)))
(นิยามเต็มๆ คือ Y = λf.(λx.f (x x)) (λx.f (x x)) ซึ่งไม่จำเป็นต้องติดตัวแปร v มาด้วย)

ลองมาสำรวจการทำงานของมันดีกว่า จะเห็นได้ว่า
Y(g)(n)
-> (lambda f: (lambda x: f(lambda v: x(x)(v)))(lambda x: f(lambda v: x(x)(v))))(g)(n)
   # reduce: g -> f
-> (lambda x: g(lambda v: x(x)(v)))(lambda x: g(lambda v: x(x)(v)))(n)
   # reduce: (lambda x: g(lambda v: x(x)(v))) -> x
-> g(lambda v: (lambda x: g(lambda v: x(x)(v)))(lambda x: g(lambda v: x(x)(v)))(v))(n)
   # reduce: n -> v
-> g((lambda x: g(lambda v: x(x)(v)))(lambda x: g(lambda v: x(x)(v)))(n))
   # for inner g function, expand back to lambda: f <- g
-> g(lambda f: (lambda x: f(lambda v: x(x)(v)))(lambda x: f(lambda v: x(x)(v)))(g)(n))
   # by definition
-> g(Y(g)(n))
   # thus
-> g(g(Y(g)(n)))
-> g(g(g(Y(g)(n))))
-> g(g(g(g(...))))
วิธีเอามาใช้งานก็ไม่ยุ่งยาก เพียง
(lambda f: (lambda x: f(lambda v: x(x)(v)))(lambda x: f(lambda v: x(x)(v))))(  # Y-comb
  lambda f:                                  # recursive function name
    lambda x: 1 if x == 0 else x * f(x-1)    # recursive function definition
  )(5)                                       # recursive function argument
ลองเปลี่ยนมาหา fibonacci บ้าง ก็ง่ายนิดเดียว
(lambda f: (lambda x: f(lambda v: x(x)(v)))(lambda x: f(lambda v: x(x)(v))))(
  lambda f:
    lambda x: x if x <= 1 else f(x-1) + f(x-2)
  )(5)
ง่ายจริงป่ะ?

Aug 21, 2012

Recursive Lambda ความงามบนความงาม...

พอถล่ำลึกลงไปใน functional programming ซักพักหนึ่ง จะเริ่มเกิดคำถามขึ้นมาว่า
"แล้วเราจะทำ recursion บน lambda ได้หรือเปล่า?"
ฟังดูเผินๆ ไม่น่าเป็นไปได้ เพราะ lambda ทำให้ฟังก์ชันไม่มีชื่อ แต่การ recursion ต้องอาศัยชื่อของฟังก์ชันนั้นๆ เพื่อเรียกมันขึ้นมาทำงานเสมอ สมมติเช่นฟังก์ชันง่ายๆ ที่จะทำการยกกำลังสองเลขที่ได้รับมาทันที
(lambda x: x**2)(5)
นี่คือฟอร์มที่เรียบง่ายสุด และถ้าเราทดลองเพิ่มตัวแปรเข้าไปให้ฟังก์ชันนี้ (โดยไม่จำเป็นต้องใช้มัน) อาจจะเขียนฟังก์ชันนี้ออกมาได้ว่า
(lambda w, x: x**2)(999, 5)
แน่นอนว่ากรณีนี้เราไม่สนใจเลข 999 ที่ถูกนำไปเก็บไว้ยังตัวแปร w เลย ...แล้วเราก็คงจะถึงทางตัน ถ้าไม่สังเกตว่าตัวแปร w นั้นหนะ รับค่าอะไรเข้ามาก็ได้ (เพราะไม่ได้ใช้) นี่รวมไปถึงเราอาจจะยัด lambda เข้าไปให้มันก็ได้
(lambda w, x: x**2)(lambda u: 42, 5)
นั่นหมายความว่า lambda ไม่จำเป็นต้องไร้ชื่อเสียทีเดียว แน่หละมันอาจจะไร้ชื่อใน scope หนึ่งๆ แต่กับ scope ที่ใหญ่ขึ้นกว่าตัวมันแล้ว เราสามารถตั้งชื่อให้มันได้ ถึงตอนนี้จะเปลี่ยนชื่อตัวแปรกันซักหน่อย พร้อมกับเอาฟังก์ชันสำหรับคำนวณจากด้านบนมาใส่แทน lambda u: 42 ไปซะ
(lambda g, y: y**2)(lambda x: x**2, 5)
จะเห็นว่าเลข 5 จะถูกเก็บไว้ในตัวแปร y และคำนวณ y**2 ตรงๆ โดยที่ lambda x: x**2 ไม่ได้แสดงบทบาทอะไร แต่เนื่องจากฟังก์ชัน lambda x: x**2 ถูก scope ด้านหน้าแปะป้ายชื่อให้ว่า g เป็นที่เรียบร้อยแล้ว ดังนั้น เราสามารถส่งผ่านเลข 5 ให้ไปทำงานกับ lambda x: x**2 ได้ดังนี้
(lambda g, y: g(y))(lambda x: x**2, 5)
แม้ว่าตอนนี้เราจะสามารถเรียกชื่อของ lambda x: x**2 ใน scope ที่ใหญ่กว่าตัวมันได้ แต่มันก็ไม่มีประโยชน์เลยเพราะการจะทำ recursion นั้นต้องสามารถเรียกชื่อตัวมันเองใน scope ที่มันอยู่ได้ ดังนั้น เราจะใช้ท่าเดิมคือส่งผ่านเจ้า lambda x: x**2 นี้ต่อไปเรื่อยๆ ไม่ให้มันหายไป โดย
(lambda g, y: g(g, y))(lambda f, x: x**2, 5)
เท่านี้ ใน scope ของ lambda f, x: x**2 ก็สามารถเรียกตัวเองได้แล้ว (ชื่อว่า f) ถึงตอนนี้ก็ได้เวลาเขียน recursion ให้มัน อาจเริ่มจากฟังก์ชันง่ายๆ อย่าง factorial ก็ได้
(lambda g, y: g(g, y))(lambda f, x: 1 if x == 0 else x * f(f, x-1), 5)
สังเกตว่าการเรียก recursion ภายใน lambda นี้ ต้องส่งผ่านชื่อฟังก์ชันตัวเองเป็นตัวแปรพ่วงลงไปด้วยเสมอนะ

จบแล้ว แค่นี้แหละ recursion บน lambda ไม่ยากอย่างที่คิดเลยใช่มั้ยครับ ;)