12.2 ห.ร.ม. ค.ร.น. และจำนวนเฉพาะ

ห.ร.ม. กับ ค.ร.น. ใช้กับโจทย์แบ่งของให้ลงตัว และโจทย์หาเวลาที่หลายเหตุการณ์มาตรงกัน ส่วนจำนวนเฉพาะเป็นฐานของโจทย์ตัวเลขแทบทุกแบบ หน้านี้ใช้ ห.ร.ม. แบบยุคลิดที่น้องเขียนเป็นแล้วในบทที่ 9 มาหา ค.ร.น. ที่ไม่ล้น แล้วต่อจากการตรวจจำนวนเฉพาะในหน้า 4.5 ไปสู่ตะแกรงที่หาจำนวนเฉพาะทั้งหมดถึงสิบล้านได้ในเวลาไม่ถึงหนึ่งในสามสิบวินาที ปิดด้วยการแยกตัวประกอบ

ห.ร.ม. แบบยุคลิด

ทบทวนสั้นๆ ห.ร.ม. แบบยุคลิดเขียนเป็นลูปในหน้า 8.5 และเป็น recursion ในหน้า 9.2 หน้านี้ใช้แบบ recursion ที่รับ long long ตามโปรแกรม lights ข้างล่าง เลขขนาด 10¹⁸ ก็ยังใช้ไม่ถึงร้อยรอบ เพราะทุกสองรอบตัวเลขเหลือไม่ถึงครึ่ง

ถ้ามีหลายตัว ให้หา ห.ร.ม. ทีละคู่ไปเรื่อยๆ g = gcd(g, x) เพราะ ห.ร.ม. ของสามตัวเท่ากับ ห.ร.ม. ของ (ห.ร.ม. ของสองตัวแรก) กับตัวที่สาม

ค.ร.น. ต้องหารก่อนคูณ

ไฟสองดวงเปิดพร้อมกันที่วินาทีที่ 0 ดวงแรกกะพริบทุก a วินาที ดวงที่สองทุก b วินาที จะกะพริบพร้อมกันอีกครั้งตอนไหน คำตอบคือ ค.ร.น. (ตัวคูณร่วมน้อย least common multiple หรือ lcm) ของ a กับ b

มีความจริงข้อหนึ่งว่า ห.ร.ม. × ค.ร.น. = a × b เช่น 12 กับ 18 มี ห.ร.ม. 6 และ ค.ร.น. 36 ซึ่ง 6 × 36 = 216 = 12 × 18 ค.ร.น. จึงเท่ากับ a × b / ห.ร.ม. แต่ในโปรแกรมให้เขียนเป็น a / g * b คือหารก่อนแล้วค่อยคูณ เหตุผลอยู่ในส่วนลองทำให้พัง

#include <stdio.h>

long long gcd(long long a, long long b)
{
    if (b == 0) {
        return a;
    }
    return gcd(b, a % b);
}

int main()
{
    long long a, b;
    scanf("%lld %lld", &a, &b);
    long long g = gcd(a, b);
    printf("gcd = %lld\n", g);
    printf("lcm = %lld\n", a / g * b);
    return 0;
}

1000000000 999999999

gcd = 1
lcm = 999999999000000000

a / g หารลงตัวเสมอ เพราะ g เป็นตัวหารของ a ไม่มีเศษหายไป และ ค.ร.น. ของเลขไม่เกิน 10⁹ สองตัวใหญ่ได้ถึงราว 10¹⁸ อย่างตัวอย่างนี้ จึงต้องใช้ long long

ค.ร.น. ของหลายตัวแบบสะสม

ค.ร.น. ของหลายตัวทำแบบสะสมเหมือน ห.ร.ม. เริ่มจาก l = 1 แล้วรวมเลขเข้าไปทีละตัวด้วย l = l / gcd(l, x) * x ซึ่งก็คือ ค.ร.น. ของ l กับ x นั่นเอง ลองกับ 6 10 15 4

x	gcd(l, x)	l ใหม่
6	1	1 / 1 × 6 = 6
10	2	6 / 2 × 10 = 30
15	15	30 / 15 × 15 = 30
4	2	30 / 2 × 4 = 60

15 ไม่ทำให้ l เปลี่ยน เพราะ 15 หาร 30 ลงตัวอยู่แล้ว ต่างจาก ห.ร.ม. ที่มีแต่เล็กลง ค.ร.น. ที่สะสมมีแต่โตขึ้นหรือเท่าเดิม จึงใหญ่ได้เร็วมาก การหารก่อนคูณจึงสำคัญกว่าตอนมีแค่สองตัว เพราะ l อาจใหญ่ระดับ 10¹⁷ แล้ว ถ้าคูณ l * x ก่อนก็ล้น long long ได้ แม้คำตอบจะไม่เปลี่ยนเลย ดูแบบที่ 1 ในส่วนลองทำให้พัง

ถ้าต้องรู้ว่าทุกเลขเป็นจำนวนเฉพาะไหม

หน้า 4.5 ตรวจเลขตัวเดียวด้วยการลองหาร 2 ถึง √n ตัวละไม่เกิน √n รอบ ถ้าโจทย์ถามว่าเลขตั้งแต่ 1 ถึง 10⁷ มีจำนวนเฉพาะกี่ตัว หรือถามหลายแสนคำถามว่าเลขนั้นเลขนี้เป็นจำนวนเฉพาะไหม การตรวจทีละตัวต้องทำซ้ำเป็นสิบล้านครั้ง ข้างล่างจะวัดให้เห็นว่าช้าแค่ไหน ทางออกคือตะแกรงเอราทอสเทนีส (sieve of Eratosthenes) ซึ่งหาคำตอบของทุกเลขพร้อมกันในรอบเดียว

แนวคิดของตะแกรง

แทนที่จะถามเลขแต่ละตัวว่ามีอะไรหารลงตัวไหม ให้กลับด้าน เอาจำนวนเฉพาะแต่ละตัวไปขีดฆ่าพหุคูณของมันทิ้ง เลขที่รอดจากการขีดฆ่าคือจำนวนเฉพาะ ลองกับเลข 2 ถึง 30

i	ยังไม่ถูกขีด	ขีดฆ่า i × i, i × i + i, ... จนถึง 30
2	ใช่ 2 เป็นจำนวนเฉพาะ	4 6 8 10 12 14 16 18 20 22 24 26 28 30
3	ใช่ 3 เป็นจำนวนเฉพาะ	9 12 15 18 21 24 27 30
4	ไม่ ถูก 2 ขีดไปแล้ว	ข้าม
5	ใช่ 5 เป็นจำนวนเฉพาะ	25 30
6	6 × 6 = 36 เกิน 30 หยุดได้

เลขที่เหลือคือ 2 3 5 7 11 13 17 19 23 29 รวม 10 ตัว มีสองจุดที่ทำให้ตะแกรงทำงานน้อยลง

เริ่มขีดที่ i × i พหุคูณที่เล็กกว่านั้น เช่น 3 × 2 = 6 ถูกจำนวนเฉพาะที่เล็กกว่าขีดไปแล้ว
หยุดเมื่อ i × i เกิน n เพราะจุดเริ่มขีดเกิน n แล้ว ไม่มีอะไรให้ขีด เหตุผลเดียวกับการลองหารถึง √n

#include <stdio.h>

char is_composite[10000001]; // (1)!

int main()
{
    int n;
    scanf("%d", &n);
    for (int i = 2; i * i <= n; i++) {
        if (!is_composite[i]) { // (2)!
            for (int j = i * i; j <= n; j += i) {
                is_composite[j] = 1;
            }
        }
    }
    int count = 0;
    for (int i = 2; i <= n; i++) {
        if (!is_composite[i]) {
            count++;
        }
    }
    printf("%d\n", count);
    return 0;
}

ช่องที่ i เป็น 1 แปลว่า i ถูกขีดฆ่าแล้ว คือเป็นจำนวนประกอบ (composite) อาเรย์อยู่นอก main ตามหน้า 5.5 จึงเริ่มต้นเป็น 0 ทุกช่อง ใช้ char เพราะเก็บแค่ 0 กับ 1 ช่องละ 1 byte สิบล้านช่องใช้ราว 10 MB ถ้าเป็น int จะใช้ 40 MB
ขีดเฉพาะพหุคูณของเลขที่ยังรอดอยู่ ถ้า i ถูกขีดไปแล้ว พหุคูณของ i ก็ถูกตัวประกอบของ i ขีดไปหมดแล้วเหมือนกัน

10000000

พอสร้างตะแกรงเสร็จ จะถามว่าเลข x ใดๆ ที่ไม่เกิน n เป็นจำนวนเฉพาะไหม ก็ดูแค่ is_composite[x] ช่องเดียว ตอบได้ทันทีกี่คำถามก็ได้ อย่าลืมว่า 0 กับ 1 ก็เป็น 0 ในอาเรย์ แต่ไม่ใช่จำนวนเฉพาะ ต้องเช็ก x >= 2 เอง

วัดเวลาเทียบกับไล่หาร

นับจำนวนเฉพาะไม่เกิน n สองวิธี วิธีไล่หารคือเอาโปรแกรมตรวจจำนวนเฉพาะของหน้า 4.5 ที่ลองหารถึง i * i <= x มาครอบด้วยลูป x ตั้งแต่ 2 ถึง n วัดด้วย clock() บนเครื่องเดียวกับหน้า 12.1 รันแบบละ 5 ครั้ง ใช้ค่ากลาง

วิธี	n	แบบ Code::Blocks	แบบ -O2
ไล่หารถึง √x ทุกตัว	10⁶	0.096 วินาที	0.092 วินาที
ไล่หารถึง √x ทุกตัว	10⁷	2.450 วินาที	2.329 วินาที
ตะแกรง	10⁶	0.002 วินาที	0.001 วินาที
ตะแกรง	10⁷	0.027 วินาที	0.012 วินาที

ทั้งสองวิธีได้ 664579 ตรงกัน แต่ที่ n = 10⁷ ตะแกรงเร็วกว่าราว 90 เท่าแบบ Code::Blocks ส่วนวิธีไล่หารเกินเวลา 1 วินาทีไปไกล สังเกตด้วยว่า n เพิ่ม 10 เท่า วิธีไล่หารช้าลงราว 25 เท่า แต่ตะแกรงช้าลงราว 13 เท่า ยิ่ง n ใหญ่ ยิ่งห่างกัน เหตุผลเชิงคณิตศาสตร์ว่าทำไมเป็นแบบนี้จะเรียนในบทที่ 18

เลือกวิธีไหนดี

ตรวจเลขไม่กี่ตัว หรือเลขใหญ่เกินกว่าจะสร้างอาเรย์ เช่น 10¹² ใช้ไล่หารถึง √n ตรวจเลขจำนวนมากที่ไม่เกินราว 10⁷ สร้างตะแกรงครั้งเดียวแล้วตอบทุกคำถามจากอาเรย์

แยกตัวประกอบ

ทุกจำนวนเต็มที่มากกว่า 1 เขียนเป็นผลคูณของจำนวนเฉพาะได้แบบเดียว เช่น 360 = 2³ × 3² × 5 วิธีแยกคือลองหารด้วย p = 2, 3, 4, ... ตัวไหนหารลงตัวก็หารซ้ำจนหารไม่ลง แล้วนับว่าหารไปกี่ครั้ง

#include <stdio.h>

int main()
{
    long long n;
    scanf("%lld", &n);
    for (long long p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            int power = 0;
            while (n % p == 0) {
                n /= p;
                power++;
            }
            printf("%lld^%d\n", p, power);
        }
    }
    if (n > 1) {
        printf("%lld^1\n", n);
    }
    return 0;
}

2^3
3^2
5^1

คำถามที่น้องน่าจะสงสัยมีสามข้อ

p = 4 หรือ 6 ไม่ใช่จำนวนเฉพาะ จะพิมพ์ออกมาไหม ไม่ เพราะตอนถึง p = 4 ตัวประกอบ 2 ถูกหารออกจาก n หมดแล้ว 4 จึงหาร n ไม่ลงตัว ตัวที่หารลงตัวได้จึงเป็นจำนวนเฉพาะเสมอ
เงื่อนไข p * p <= n ใช้ n ที่เล็กลงเรื่อยๆ ได้หรือ ได้ และยิ่งดี n ที่เหลือไม่มีตัวประกอบที่เล็กกว่า p แล้ว ถ้ามันยังมีตัวประกอบสองตัวขึ้นไป ตัวที่เล็กสุดต้องไม่เกิน √n ลูปจึงหยุดเร็วขึ้นเมื่อ n ลดลง
บรรทัดที่ 17 มีไว้ทำไม พอลูปจบ ถ้า n ยังเหลือมากกว่า 1 ที่เหลือคือจำนวนเฉพาะตัวใหญ่ที่เกิน √ ของ n ตั้งต้น มีได้อย่างมากตัวเดียว เช่น 2000006 = 2 × 1000003 ลูปหาร 2 ออกแล้วเหลือ 1000003 ซึ่ง 1001 × 1001 เกินแล้ว ลูปจึงจบโดยไม่เจอมัน

เลข 10¹² วนไม่เกิน 10⁶ รอบ เสร็จทันที ลองป้อน 999999999989 ซึ่งเป็นจำนวนเฉพาะ ลูปจะวนจนสุดทาง ได้ 999999999989^1 และป้อน 1000000000000 ได้ 2^12 กับ 5^12

ลองทำให้พัง

แบบที่ 1 คูณก่อนหารใน ค.ร.น.

เขียนโปรแกรมไฟกะพริบด้วย int เพราะโจทย์บอกว่า a กับ b ไม่เกิน 10⁵ และคำตอบรับประกันว่าไม่เกิน int แล้วคำนวณ ค.ร.น. ทั้งสองแบบ

#include <stdio.h>

int gcd(int a, int b)
{
    if (b == 0) {
        return a;
    }
    return gcd(b, a % b);
}

int main()
{
    int a, b;
    scanf("%d %d", &a, &b);
    int g = gcd(a, b);
    printf("gcd = %d\n", g);
    printf("multiply first = %d\n", a * b / g);
    printf("divide first   = %d\n", a / g * b);
    return 0;
}

100000 100000

gcd = 100000
multiply first = 14100
divide first   = 100000

ไฟสองดวงที่กะพริบทุก 100000 วินาทีเท่ากัน ย่อมกะพริบพร้อมกันทุก 100000 วินาที คำตอบเล็กนิดเดียว แต่ a * b คือ 10¹⁰ ล้น int ไปก่อนจะได้หาร แบบหารก่อนค่าระหว่างทางไม่เคยเกินคำตอบ จึงไม่ล้น เรื่องเดียวกันเกิดกับ long long เมื่อหา ค.ร.น. ของหลายตัว ค.ร.น. ที่สะสมมาใหญ่ระดับ 10¹⁷ คูณกับเลขถัดไปก็ล้นได้ แม้ ค.ร.น. ใหม่จะไม่เปลี่ยนเลย

แบบที่ 2 ปล่อยลูปนอกของตะแกรงวิ่งถึง n

ในโปรแกรม sieve เปลี่ยนบรรทัดที่ 9 เป็น for (int i = 2; i <= n; i++) { คิดว่าวนเกินไปก็แค่เสียเวลานิดหน่อย คอมไพล์ผ่านโดยไม่มีคำเตือน แต่ป้อน 100000 แล้วโปรแกรมปิดตัวไปโดยไม่พิมพ์อะไรเลย บนเครื่องที่เขียนเอกสารนี้ Code::Blocks ขึ้นว่า

Process returned -1073741819 (0xC0000005)

พอ i เกิน 46340 ค่า i * i เกิน 2147483647 และล้น int กลายเป็นเลขติดลบ j จึงเริ่มที่ช่องติดลบซึ่งไม่มีอยู่จริง เงื่อนไข j <= n ก็ยังเป็นจริง โปรแกรมจึงเขียนลงหน่วยความจำที่ไม่ใช่ของตัวเอง 0xC0000005 คือรหัสของ Windows ที่แปลว่าเข้าถึงหน่วยความจำที่ไม่มีสิทธิ์ เหมือนอาการที่ 3 ในหน้า 5.5 เงื่อนไข i * i <= n จึงไม่ได้มีไว้ให้เร็วอย่างเดียว แต่กันไม่ให้ i * i ล้นด้วย

แบบที่ 3 ลืมตัวประกอบที่เหลือ

ตัดสามบรรทัด if (n > 1) ท้ายโปรแกรม factorize ออก แล้วป้อน 2000006

#include <stdio.h>

int main()
{
    long long n;
    scanf("%lld", &n);
    for (long long p = 2; p * p <= n; p++) {
        if (n % p == 0) {
            int power = 0;
            while (n % p == 0) {
                n /= p;
                power++;
            }
            printf("%lld^%d\n", p, power);
        }
    }
    return 0;
}

2^1

1000003 หายไปเงียบๆ ทดสอบกับ 360 หรือ 1000000000000 จะไม่เห็นปัญหานี้เลย เพราะตัวประกอบทุกตัวเล็ก ชุดทดสอบที่ดีต้องมีเลขที่เป็นจำนวนเฉพาะตัวใหญ่เอง และเลขที่เป็นจำนวนเฉพาะเล็กคูณจำนวนเฉพาะใหญ่

สรุป

ห.ร.ม. ใช้ gcd(b, a % b) เร็วมากแม้เลขใหญ่ระดับ 10¹⁸ หลายตัวก็ทำทีละคู่
ค.ร.น. คือ a / gcd(a, b) * b หารก่อนคูณเสมอ ค.ร.น. ของเลข 10⁹ สองตัวต้องใช้ long long
ตะแกรงเอราทอสเทนีสขีดพหุคูณตั้งแต่ i * i และหยุดเมื่อ i * i > n นับจำนวนเฉพาะถึง 10⁷ ได้ใน 0.027 วินาทีบนเครื่องที่เขียนเอกสารนี้ ส่วนไล่หารใช้ 2.45 วินาที
แยกตัวประกอบด้วยการหารซ้ำตั้งแต่ 2 จนถึง √n และอย่าลืมจำนวนเฉพาะตัวใหญ่ที่เหลืออยู่ตอนจบ