Category Archives: Math Geek

วาดภาพขาวดำให้เป็นจำนวนเฉพาะ

ผมไปเห็นโพสต์ตลกๆใน Reddit มา การ์ตูนบอกว่าสามารถเอาตัวเลข 0 และ 1 มาเรียงกันบนจอโทรศัพท์ให้เป็นรูปยีราฟและเลข 0 และ 1 ที่เรียงกันนั้นถ้าตีความเป็นเเลขฐานสองแล้ว มันจะเป็นเลขจำนวนเฉพาะอีกด้วย (เลขจำนวนเฉพาะคือจำนวนเต็มบวกที่มีตัวหารที่เป็นบวกแค่ 1 และตัวมันเอง)

ภาพการ์ตูนจาก https://www.reddit.com/r/math/comments/7qpfls/does_there_exist_a_prime_number_whose/ ครับ
ภาพการ์ตูนจาก https://www.reddit.com/r/math/comments/7qpfls/does_there_exist_a_prime_number_whose/ ครับ

 

ใต้การ์ตูนมีคนทำรูปยีราฟที่วาดด้วย 0  และ  1 เรียงกัน 4096 ตัว ซึ่งถ้ามองเป็นเลขฐานสองจะเป็นจำนวนเฉพาะด้วยครับ:

ภาพยีราฟวาดด้วย 0 และ 1
ภาพยีราฟวาดด้วย 0 และ 1

ถ้าเราสังเกตที่ด้านล่างของภาพยีราฟจะเห็นว่ามันไม่เป็น 0 ทั้งหมด ทำให้น่าสงสัยว่าคนสร้างภาพเปลี่ยนค่า 0 บางตัวเป็น 1 จนกระทั่งทั้งภาพกลายเป็นจำนวนเฉพาะ ผมก็เลยคิดสงสัยว่าเราสามารถสร้างรูปขาวดำอะไรก็ได้ให้เป็นจำนวนเฉพาะแบบนี้ได้ไหม 

ก่อนอื่นเราต้องหาตรวจสอบว่าเลขฐานสองยาวๆเป็นจำนวนเฉพาะหรือไม่ก่อน ถ้าเราใช้ภาษาไพธอน แล้วเรามีสตริงชื่อ x ที่มีเลข 0/1 เรียงกันเป็นเลขฐานสองเราใช้คำสั่ง int(x,2) เปลี่ยนเป็นค่าตัวเลขได้ แล้วใช้คำสั่ง gmpy2.is_prime(num) มาตรวจว่าเป็นจำนวนเฉพาะหรือเปล่า 


x = "0100011"
num = int(x,2) #num = 35
gmpy2.is_prime(num) #False เพราะ 35 = 1x5x7

y = '0100101'
num = int(x,2) #num = 37
gmpy2.is_prime(num) #True เพราะ 37 = 1x37

เวลาตัวเลขของเรามีขนาดไม่ใหญ่นัก มีไม่กี่หลัก เราสามารถทดลองหาตัวหารไปเรื่อยๆได้ แต่ถ้าตัวเลขใหญ่มีเป็นร้อยเป็นพันหลัก เราจะลองอย่างนั้นไม่ได้เพราะจะใช้เวลามากเกินไป ความวิเศษของ gmpy2.is_prime( ) ก็คือมันใช้เทคนิค Miller-Rabin ที่สามารถตรวจสอบเลขที่มีหลายๆพันหลักได้ด้วยวิธีสุ่ม โดยที่เราเลือกลดความน่าจะเป็นที่มันจะตอบผิดว่าตัวเลขเป็นจำนวนเฉพาะทั้งๆที่ไม่ใช่ให้มีค่าน้อยๆใกล้ศูนย์ได้ตามใจ 

พอเรารู้วิธีตรวจสอบจำนวนเฉพาะขนาดใหญ่ๆได้แล้ว เราก็แค่เขียนโปรแกรมให้อ่านภาพเข้ามาแล้วแปลงเป็นภาพขาวดำ เปลี่ยนสีดำเป็น 0 เปลี่ยนสีขาวเป็น 1 

โค้ดส่วนนี้อ่านภาพเข้ามาแล้วทำให้ขนาดไม่ใหญ่เกินไป:


def scale_image(file, size = 64):
    """
    Read an image file, change it to grayscale (0-255),
    and scale it so that its dimension is at most size x size,
    keeping the image's aspect ratio.
    """

    im = PIL.Image.open(file)
    im = im.convert('L') #change image to grayscale (0 to 255)
   
    width, height = im.size
    if width > height:
        max_size = width
    else:
        max_size = height
       
    scale_factor = size/max_size
    new_width = int(scale_factor*width)
    new_height = int(scale_factor*height)
    im = im.resize((new_width, new_height))
   
    return im

ส่วนนี้เปลี่ยนภาพให้เป็นสตริง 0/1


def image_to_binary(image):
    """
    Convert a grayscale (0-255) image to its binary string representation.
    Change the last pixel to 1 to make it a prime number more easily.
    """

    width, height = image.size
   
    pixels = np.asarray(image) #store pixel value in numpy array
    binary_pixels = (pixels > 128).astype(np.int) #low values go to 0, high values go to 1
    binary_pixels[(-1,-1)] = 1 #make sure last pixel is 1 to make it a prime more easily
    binary_string = "".join(binary_pixels.flatten().astype(np.str).tolist())
   
    return width, height, binary_string #return width and height along with 1-dimension binary representation

แล้วเราก็ตรวจสอบว่าสตริง 0/1 เป็นจำนวนเฉพาะด้วยฟังก์ชั่นนี้:


def is_binary_string_prime(x, n_tests = 50):
    """
    Check if the binary string x is prime using gmpy2's Miller-Rabin
    algorithm with n_tests steps.
    """

    num = int(x,2) #convert binary string to decimal number
    return gmpy2.is_prime(num, n_tests)

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


def mutate(x):
    """
    Given a binary string x, flip one random bit and return the result.
    """

    i = random.randint(0, len(x)-1) #location to flip a bit

    if x[i] == '0':
        flip = '1'
    else:
        flip = '0'
       
    result = x[:i] + flip + x[i+1:] #copy x to result, except one flipped bit
   
    return result

ปัญหาต่อไปก็คือถ้าภาพยังไม่เป็นจำนวนเฉพาะ เราจะทำการเปลี่ยนจุด 1 จุดกี่ครั้งดี เราอาศัยความจริงที่ว่าจำนวนเฉพาะที่น้อยกว่า N จะมีประมาณ N/Log(N) ตัว ดังนั้นความน่าจะเป็นที่เราสุ่มเลือกเลขที่อยู่ระหว่าง 1 ถึง N มาแล้วมันเป็นจำนวนเฉพาะก็จะประมาณ 1/Log(N)

เปรียบเทียบจำนวนของจำนวนเฉพาะที่ไม่เกิน N และ N/Log(N)
เปรียบเทียบจำนวนของจำนวนเฉพาะที่ไม่เกิน N และ N/Log(N)

ถ้าภาพของเรามีจุด 0/1 ทั้งหมด K จุด ความน่าจะเป็นที่มันเป็นจำนวนเฉพาะก็จะประมาณ 1/Log(2K) = 1/(K Log(2)) = 1.44/K  เช่นถ้าภาพมีขนาด 64×64 แสดงว่า K = 64×64 = 4096 และความน่าจะเป็นที่ภาพจะเป็นจำนวนเฉพาะก็ประมาณ 1.44/4096 เท่ากับประมาณ 1/2800 แปลว่าถ้าเราเลือกเลขสุ่มๆมาสักสองสามเท่าของ 2800 ครั้ง เราน่าจะได้เลขจำนวนเฉพาะมาสองสามตัวถ้าเราไม่โชคร้ายเกินไป

จากนั้นเราก็แค่ลองเปลี่ยนเพียง 1 จุดในภาพไปเรื่อยๆจนภาพกลายเป็นจำนวนเฉพาะด้วยฟังก์ชั่นหน้าตาแบบนี้ครับ:


def mutate_to_prime(x, max_tries = 12000):
    """
    Given a binary string x, flip just one bit to make it a prime number.
    max_tries is the maximum number of tries before giving up.
    If x is N bits long, probability that x is a prime is about 1/(N log(2)),
    so max_tries about a few times N log(2) should be OK.
   
    The first prime found will be returned.
    If no prime is found, None will be returned.
    """


    if is_binary_string_prime(x):
        return x  #if x is already prime, just return it.
   
    found_prime = False
    for k in range(max_tries):
        m = mutate(x) #try flipping one bit
        if is_binary_string_prime(m): #if it's a prime, return it
            return m
   
    return None

ผมประกอบฟังก์ชั่นต่างๆเป็นโปรแกรมแบบนี้ครับ โหลดได้ที่นี่:


import PIL.Image
import gmpy2
import random
import numpy as np
import math

def scale_image(file, size = 64):
    """
    Read an image file, change it to grayscale (0-255),
    and scale it so that its dimension is at most size x size,
    keeping the image's aspect ratio.
    """

    im = PIL.Image.open(file)
    im = im.convert('L') #change image to grayscale (0 to 255)
   
    width, height = im.size
    if width > height:
        max_size = width
    else:
        max_size = height
       
    scale_factor = size/max_size
    new_width = int(scale_factor*width)
    new_height = int(scale_factor*height)
    im = im.resize((new_width, new_height))
   
    return im

def image_to_binary(image):
    """
    Convert a grayscale (0-255) image to its binary string representation.
    Change the last pixel to 1 to make it a prime number more easily.
    """

    width, height = image.size
   
    pixels = np.asarray(image) #store pixel value in numpy array
    binary_pixels = (pixels > 128).astype(np.int) #low values go to 0, high values go to 1
    binary_pixels[(-1,-1)] = 1 #make sure last pixel is 1 to make it a prime more easily
    binary_string = "".join(binary_pixels.flatten().astype(np.str).tolist())
   
    return width, height, binary_string #return width and height along with 1-dimension binary representation

def mutate(x):
    """
    Given a binary string x, flip one random bit and return the result.
    """

    i = random.randint(0, len(x)-1) #location to flip a bit

    if x[i] == '0':
        flip = '1'
    else:
        flip = '0'
       
    result = x[:i] + flip + x[i+1:] #copy x to result, except one flipped bit
   
    return result

def is_binary_string_prime(x, n_tests = 50):
    """
    Check if the binary string x is prime using gmpy2's Miller-Rabin
    algorithm with n_tests steps.
    """

    num = int(x,2) #convert binary string to decimal number
    return gmpy2.is_prime(num, n_tests)

def mutate_to_prime(x, max_tries = 12000):
    """
    Given a binary string x, flip just one bit to make it a prime number.
    max_tries is the maximum number of tries before giving up.
    If x is N bits long, probability that x is a prime is about 1/(N log(2)),
    so max_tries about a few times N log(2) should be OK.
   
    The first prime found will be returned.
    If no prime is found, None will be returned.
    """


    if is_binary_string_prime(x):
        return x  #if x is already prime, just return it.
   
    found_prime = False
    for k in range(max_tries):
        m = mutate(x) #try flipping one bit
        if is_binary_string_prime(m): #if it's a prime, return it
            return m
   
    return None

def print_pic(x, width):
    """
    Given a binary string x, print it
    so that each line contains width characters.
    """

    height = int(len(x)/width)
    for row in range(height):
        print(x[row*width : row*width+width])

def pic_to_prime(file, size = 64, max_tries = 12000):
    """
    Open an image file and attempt to
    print it as a binary string picture
    that is a prime number.
    """

    im = scale_image(file, size)
   
    width, height, binary = image_to_binary(im)
   
    x = mutate_to_prime(binary, max_tries)
    if x:
        print_pic(x,width)
        num = int(x,2)
        print("\nThe number is {} in base 10.".format(num))
        print("It has {} digits".format(len(str(num))))
        print("It's likely a prime with probability = {}".format(1-math.pow(0.5,50)))

def main():
    import sys
    usage = """
    python pic_to_prime.py image_file [size] [max_tries]
    Will attempt to convert image_file into a black and white picture whose ASCII representation is
    (black = 0, white = 1) and the string of 0's and 1's forms a binary number which is likely prime with
    an extremely high probablity.

    size is an optional parameter.  The converted image's dimension will be at most size x size
    Default value for size is 64

    max_tries is an optional parameter. It's the maximum number of tries to flip one bit in the converted image
    until it's a prime number.  If max_tries is reached without find a prime, the program does not return any Ascii image.
    Typically, setting max_tries to be a few times size x size should be sufficient.
    Default value for max_tries is 12,000
    """


    if len(sys.argv) < 2:
        print(usage)
        exit(1)

    try:
        file_name = sys.argv[1]
        print("Converting {}...".format(file_name))
        if len(sys.argv) == 2:
            pic_to_prime(file_name)
        if len(sys.argv) == 3:
            size = int(sys.argv[2])
            pic_to_prime(file_name, size)
        if len(sys.argv) &gt;= 4:
            size = int(sys.argv[2])
            max_tries = int(sys.argv[3])
            pic_to_prime(file_name, size, max_tries)      
       
       
    except:
        print("An error has occured:", sys.exc_info())
   
if __name__ == '__main__':
    main()

วิธีใช้ก็เรียก python pic_to_prime.py filename โดยที่ filename คือชื่อไฟล์รูปภาพที่เราต้องการเปลี่ยนเป็นจำนวนเฉพาะครับ 

ตัวอย่างภาพที่ได้จากการเปลี่ยนภาพ Darth Vader:

ภาพ Darth Vader ข้างบนคือเลขฐานสิบที่มี 1195 หลัก หน้าตาแบบนี้ครับ: 

3069183072727172278860337378214364811175034036025244450162785708918361127976622462244955616465542073310944718954147895761643252918537495912453992270956943332625163100410236971285347282168904965253931503727481796300351734500671920231110149740644603732320920567437139574018949692454814303373410337477406380375065946238083359795848929370044692287614885238710013203784493630576976222535714074987122348018208566276582108014408460395720661886870781674435952140263967388925330430468219482326166783846446434209037095617083753020217147213702448717413427206823791407299246422960343709648378621123745736208899112612922595944987425466506332619408247083609772049609869403316890387132071258208841211446813145804008085966690045193746244730410209439050722697005342344505437504795775660537102148374631191587907498516228228789209160814532098859180914649594125752152645342655118749520955226752989697313945409163186662411266021539069921213862430095726119501306030813846080124361789827996639620644316501126994186024697789684863191124018231390812596919287605853413555802043957703338340785628246023479885637324716903027104728159535177737667950181169447546144659362616465956401475098937299932469625956492313055325061119

 

สร้าง Bifurcation Diagram แบบง่ายๆ


พอดีมีนักเรียนตั้งกระทู้ถามเรื่องการวาดรูปเรื่อง Bifurcation Diagram เมื่อคืน ผมเลยทดลองวาดบน Mathematica เห็นว่าวิธีวาดง่ายดี เลยมาบันทึกไว้เผื่อมีใครค้นหาอีกในอนาคต ถ้าจะลองก็คัดลอกเอาไปลองได้เลย แล้วเปลี่ยนนู่นเปลี่ยนนี่เล่นดูเอง


(* กำหนด mapping ที่เราสนใจ อันนี้เรียกว่า logistic map *)
f[r_, x_] := r x (1 – x)

(* ทำการ iterate ด้วย ‘r’ ไป ‘iterations’ ครั้ง เริ่มด้วย ‘x0’ แล้วตัดมาดู ‘count’ ตัว *)
longtermValues[r_, count_, iterations_, x0_] :=
Map[{r, #} &, Take[NestList[f[r, #] &, x0, iterations], -count]]

(* เราเปลี่ยน r ตั้งแต่ 2.6 ไปจนถึง 4 โดยขยับทีละ 0.001 สำหรับแต่ละ r เรา iterate 500 ครั้งแล้วเอา 200 ตัวสุดท้่ายมาใช้ เราต้องใช้ Flatten[…,1] เพื่อให้คู่ลำดับทั้งหมดอยู่ในลิสท์ระดับเดียวกัน *)
allValues = Flatten[Table[longtermValues[r, 200, 500, 0.2], {r, 2.6, 4, 0.001}], 1];

(* วาดรูป *)
ListPlot[allValues, PlotStyle -> PointSize[0.001]]

A Barnsley’s Fern In 7 Lines of Mathematica

I used to draw a Barnsley’s fern with a program written in Pascal when I was 19 years old. Yesterday someone asked about it in a forum I visited, so I drew another one using Mathematica. The code is much shorter this time. (I’m sure that many people can shorten it even more.)
Here’s the code to draw the fern with 10,000 points. You can copy and paste and run it in Mathematica:

          

ifsFern[p_] := Module[{i},
i = Random[Integer, 99];
If[i < 1, Return[{{0., 0.}, {0., .16}}.p ]];
If[i >= 1 && i < 86, Return[{{0.85, 0.04}, {-0.04, 0.85}}.p + {0., 1.6}]];
If[i >= 86 && i < 93, Return[{{0.20, -0.26}, {0.23, 0.22}}.p + {0., 1.6}]];
If[i >= 93, Return[{{-0.15, 0.28}, {0.26, 0.24}}.p + {0., 0.44}]]]


Graphics[{RGBColor[0, 0.5, 0], Point[NestList[ifsFern, {0, 0}, 10000]]}]


The result looks like this: