8.5 math.h และฟังก์ชันคณิตศาสตร์ที่ใช้บ่อย

ภาษา C มีฟังก์ชันคณิตศาสตร์สำเร็จรูปให้ใน math.h เช่น รากที่สองและยกกำลัง ใช้คำนวณระยะทาง พื้นที่ ดอกเบี้ย ซึ่งออกในโจทย์บ่อย แต่มีกับดักใหญ่เรื่องจำนวนเต็มที่ทำให้คนเสียคะแนนกันมาก ครึ่งหลังของหน้าเป็นฟังก์ชันที่ไม่มีในไลบรารี แต่น้องจะได้เขียนเองซ้ำแล้วซ้ำอีกในโจทย์แข่งขัน คือตรวจจำนวนเฉพาะและหา ห.ร.ม.

sqrt รากที่สอง

บ้านอยู่ที่พิกัด (0, 0) โรงเรียนอยู่ที่ (3, 5) ระยะทางเป็นเส้นตรงหาได้จากทฤษฎีบทพีทาโกรัส คือรากที่สองของ dx² + dy² ฟังก์ชัน sqrt (มาจาก square root) รับ double และส่ง double กลับ

#include <stdio.h>
#include <math.h> // (1)!

double distance(double x1, double y1, double x2, double y2)
{
    double dx = x2 - x1;
    double dy = y2 - y1;
    return sqrt(dx * dx + dy * dy);
}

int main()
{
    double home_x, home_y, school_x, school_y, tutor_x, tutor_y;
    scanf("%lf %lf", &home_x, &home_y);
    scanf("%lf %lf", &school_x, &school_y);
    scanf("%lf %lf", &tutor_x, &tutor_y);
    double first = distance(home_x, home_y, school_x, school_y);
    double second = distance(school_x, school_y, tutor_x, tutor_y);
    printf("home to school %.2f\n", first);
    printf("school to tutor %.2f\n", second);
    printf("total %.2f\n", first + second);
    return 0;
}

ฟังก์ชันใน math.h ต้อง include ไฟล์นี้ก่อน เหมือน string.h ในบทที่ 6

0 0
3 5
7 2

home to school 5.83
school to tutor 5.00
total 10.83

ห่อสูตรไว้ในฟังก์ชัน distance แล้ว main ก็ไม่ต้องเขียนสูตรยาวๆ ซ้ำ โจทย์ที่ต้องหาระยะระหว่างจุดหลายคู่จะสั้นลงมาก

เทียบระยะไม่ต้องใช้ sqrt

ถ้าโจทย์แค่ถามว่าจุดไหนใกล้กว่า ให้เทียบ dx² + dy² ตรงๆ เพราะระยะที่มากกว่าก็มีกำลังสองมากกว่าเสมอ ถ้าพิกัดเป็นจำนวนเต็ม ค่านี้เป็นจำนวนเต็มด้วย เทียบได้แม่นยำ ไม่มีปัญหาทศนิยมคลาดแบบหน้า 1.4

pow ยกกำลัง

pow(x, y) คือ x ยกกำลัง y รับ double สองตัว ส่ง double กลับ เหมาะกับงานที่เป็นทศนิยมอยู่แล้ว เช่น ฝากเงินได้ดอกเบี้ยทบต้นปีละ r เปอร์เซ็นต์ ผ่านไป y ปีจะมีเงิน เงินต้น × (1 + r/100)^y

#include <stdio.h>
#include <math.h>

int main()
{
    double money, rate;
    int years;
    scanf("%lf %lf %d", &money, &rate, &years);
    double result = money * pow(1 + rate / 100, years);
    printf("%.2f\n", result);
    printf("%f\n", pow(2, 10));
    return 0;
}

1000 2 10

1218.99
1024.000000

บรรทัดสุดท้ายสังเกตว่า pow(2, 10) ได้ 1024.000000 ไม่ใช่ 1024 แม้จะส่งจำนวนเต็มเข้าไป เพราะ pow ส่ง double กลับเสมอ ต้องพิมพ์ด้วย %f ถ้าพิมพ์ด้วย %d จะได้คำเตือนเรื่อง format แบบหน้า 1.4

กับดัก pow กับจำนวนเต็มใหญ่

เพราะ pow คำนวณเป็น double ซึ่งแม่นได้ราว 15 ถึง 16 หลักตามหน้า 1.4 ถ้าคำตอบเป็นจำนวนเต็มที่ยาวกว่านั้น หลักท้ายๆ จะเพี้ยน ลองเทียบกับฟังก์ชันยกกำลังที่เขียนเองด้วยการคูณซ้ำแบบ long long

#include <stdio.h>
#include <math.h>

long long power(long long base, int exponent)
{
    long long result = 1;
    for (int i = 0; i < exponent; i++) {
        result *= base;
    }
    return result;
}

int main()
{
    int base, exponent;
    scanf("%d %d", &base, &exponent);
    long long by_pow = (long long)pow(base, exponent);
    long long by_loop = power(base, exponent);
    printf("pow  = %lld\n", by_pow);
    printf("loop = %lld\n", by_loop);
    printf("%.1f\n", pow(base, exponent));
    return 0;
}

3 39

pow  = 4052555153018976256
loop = 4052555153018976267
4052555153018976256.0

3³⁹ เป็นเลข 19 หลัก long long เก็บได้พอดี แต่ double เก็บได้แค่ค่าที่ใกล้ที่สุด บรรทัดสุดท้ายแสดงค่าใน double ตรงๆ จะเห็นว่าหลักท้ายเพี้ยนไปตั้งแต่ก่อนแปลงเป็น long long แล้ว ผลต่างแค่ 11 จากเลขหลักล้านล้านล้านดูเล็กมาก แต่ grader เทียบทุกหลัก ผิดหลักเดียวก็ผิดทั้งข้อ

กฎจึงเหมือนเรื่อง double ในหน้า 1.4 คำตอบเป็นจำนวนเต็ม ห้ามคำนวณผ่าน pow ให้เขียน power เองแบบข้างบน และระวังว่าผลต้องไม่เกิน long long ด้วย ถ้าโจทย์ให้เลขชี้กำลังใหญ่ระดับ 10⁹ ลูปนี้จะช้าเกินไป วิธีที่เร็วกว่าจะเรียนในบทที่ 9 และบทที่ 12

abs กับ fabs ค่าสัมบูรณ์

ค่าสัมบูรณ์คือระยะห่างจาก 0 โดยไม่สนเครื่องหมาย ใช้หาว่าสองค่าห่างกันเท่าไร โดยไม่ต้องรู้ว่าตัวไหนมากกว่า มีสองฟังก์ชันที่ชื่อคล้ายกันแต่ใช้คนละชนิด

ฟังก์ชัน	ใช้กับ	อยู่ใน
`abs(x)`	`int`	`stdlib.h`
`fabs(x)`	`double` (f มาจาก floating point)	`math.h`

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main()
{
    int floor_from, floor_to;
    double temp_morning, temp_noon;
    scanf("%d %d", &floor_from, &floor_to);
    scanf("%lf %lf", &temp_morning, &temp_noon);
    printf("floors = %d\n", abs(floor_to - floor_from));
    printf("change = %.1f\n", fabs(temp_noon - temp_morning));
    return 0;
}

7 2
31.5 26.8

floors = 5
change = 4.7

ลิฟต์ลงจากชั้น 7 ไปชั้น 2 ผลลบคือ -5 แต่จำนวนชั้นที่ผ่านคือ 5 อุณหภูมิลดลง 4.7 องศาก็เช่นกัน ใช้ผิดตัวแล้วจะเกิดอะไรขึ้น ดูแบบที่ 2 ในส่วนลองทำให้พัง

is_prime ตรวจจำนวนเฉพาะ

ในหน้า 4.5 เราตรวจจำนวนเฉพาะด้วยลูปที่หยุดที่ √n โจทย์จริงมักถามต่อจากนั้น เช่น นับจำนวนเฉพาะในช่วง หาจำนวนเฉพาะคู่แฝด หรือตรวจว่าเลขกลับหลังยังเป็นจำนวนเฉพาะไหม ทุกแบบต้องตรวจซ้ำหลายครั้ง จึงควรเป็นฟังก์ชัน

#include <stdio.h>

int is_prime(int n)
{
    if (n < 2) { // (1)!
        return 0;
    }
    for (long long i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            return 0; // (2)!
        }
    }
    return 1;
}

int main()
{
    int limit;
    scanf("%d", &limit);
    for (int p = 2; p + 2 <= limit; p++) {
        if (is_prime(p) && is_prime(p + 2)) { // (3)!
            printf("%d %d\n", p, p + 2);
        }
    }
    return 0;
}

0 กับ 1 ไม่ใช่จำนวนเฉพาะ และลูปข้างล่างจะไม่ทำงานเลยกับสองค่านี้ ต้องแยกไว้ตั้งแต่ต้น
เจอตัวหารตัวแรกก็ตอบได้เลย return จบฟังก์ชันและออกจากลูปไปพร้อมกัน ไม่ต้องใช้ break กับตัวแปร flag แบบหน้า 4.5
ฟังก์ชันคืน 1 เมื่อจริง และ 0 เมื่อเท็จ ตามหน้า 3.1 จึงใช้ในเงื่อนไขได้ตรงๆ อ่านออกเสียงได้ว่า "ถ้า p เป็นจำนวนเฉพาะ และ p + 2 เป็นจำนวนเฉพาะ"

ฟังก์ชันที่ตอบใช่หรือไม่ใช่แบบนี้ นิยมตั้งชื่อขึ้นต้นด้วย is_ เช่น is_prime is_leap is_vowel อ่านในเงื่อนไขแล้วเป็นประโยคเลย เทียบกับการเขียนลูปตรวจจำนวนเฉพาะสองชุดซ้อนอยู่ใน main แบบนี้ทั้งสั้นและผิดยากกว่ามาก

gcd หา ห.ร.ม.

ห.ร.ม. (ตัวหารร่วมมาก) หรือ gcd (greatest common divisor) ใช้ทำเศษส่วนให้เป็นอย่างต่ำ หา ค.ร.น. และโจทย์แบ่งของเป็นกองเท่าๆ กัน วิธีลองหารทุกตัวก็ได้ แต่มีวิธีที่เร็วกว่ามากเรียกว่า วิธีของยุคลิด (Euclid) ซึ่งใช้แค่ %

ความคิดคือ ตัวเลขใดที่หารทั้ง a และ b ลงตัว จะหาร a % b ลงตัวด้วย เพราะ a % b คือ a ลบด้วย b หลายๆ เท่า และกลับกันก็จริง ห.ร.ม. ของ a กับ b จึงเท่ากับ ห.ร.ม. ของ b กับ a % b ซึ่งเป็นเลขที่เล็กลง ทำซ้ำไปจนตัวหลังเป็น 0 ตัวหน้าคือคำตอบ ลองกับ 84 และ 36

a	b	a % b
84	36	12
36	12	0
12	0	หยุด ได้ 12

#include <stdio.h>

int gcd(int a, int b)
{
    while (b != 0) {
        int remainder = a % b;
        a = b; // (1)!
        b = remainder;
    }
    return a;
}

int main()
{
    int top, bottom;
    scanf("%d %d", &top, &bottom);
    int g = gcd(top, bottom);
    printf("gcd = %d\n", g);
    printf("fraction = %d/%d\n", top / g, bottom / g);
    printf("lcm = %d\n", top / g * bottom); // (2)!
    return 0;
}

แก้ค่า a กับ b ได้ตามสบาย เพราะเป็นสำเนาตามหน้า 8.2 top กับ bottom ใน main ยังอยู่ครบ
ค.ร.น. คือ a × b ÷ ห.ร.ม. แต่หารก่อนแล้วค่อยคูณ ผลระหว่างทางจะได้ไม่ใหญ่เกิน ดูแบบที่ 3 ในส่วนลองทำให้พัง

84 36

gcd = 12
fraction = 7/3
lcm = 252

84/36 ทำเป็นอย่างต่ำได้ 7/3 ด้วยการหารทั้งเศษและส่วนด้วย ห.ร.ม. วิธีของยุคลิดเร็วมาก คู่ 1836311903 กับ 1134903170 ซึ่งเป็นคู่ที่ต้องวนนานที่สุดในช่วงของ int ยังวนแค่ 44 รอบ เหตุผลว่าทำไมเร็ว และเรื่องจำนวนเฉพาะกับ ห.ร.ม. ที่ลึกกว่านี้ อยู่ในบทที่ 12

ลองทำให้พัง

แบบที่ 1 ลืม include math.h

ในโปรแกรม route ลบบรรทัด #include <math.h> ออก

route.c:7:12: error: implicit declaration of function 'sqrt' [-Wimplicit-function-declaration]
    7 |     return sqrt(dx * dx + dy * dy);
      |            ^~~~
route.c:2:1: note: include '<math.h>' or provide a declaration of 'sqrt'
    1 | #include <stdio.h>
  +++ |+#include <math.h>
    2 | 
route.c:7:12: warning: incompatible implicit declaration of built-in function 'sqrt' [-Wbuiltin-declaration-mismatch]
    7 |     return sqrt(dx * dx + dy * dy);
      |            ^~~~
route.c:7:12: note: include '<math.h>' or provide a declaration of 'sqrt'

เป็น error แบบเดียวกับแบบที่ 1 ในหน้า 8.1 เพราะไม่มี prototype ของ sqrt ต่างกันตรงที่ compiler รู้จักฟังก์ชันมาตรฐานตัวนี้ จึงบอกเลยว่าต้อง include ไฟล์ไหน และชี้ตำแหน่งที่ควรเติม บรรทัดที่มี + คือบรรทัดที่มันแนะนำให้เพิ่ม

แบบที่ 2 ใช้ abs กับ double

ในโปรแกรม abs เปลี่ยน fabs เป็น abs และพิมพ์ด้วย %d เพราะคิดว่า abs ใช้ได้กับทุกอย่าง

#include <stdio.h>
#include <stdlib.h>

int main()
{
    double temp_morning, temp_noon;
    scanf("%lf %lf", &temp_morning, &temp_noon);
    printf("change = %d\n", abs(temp_noon - temp_morning));
    return 0;
}

31.5 26.8

change = 4

abs รับ int ค่า -4.7 จึงถูกตัดทศนิยมเหลือ -4 ตั้งแต่ตอนส่งเข้าตามหน้า 8.2 ได้ผล 4 แทน 4.7 และ -Wall ไม่เตือนเลย ถ้าเปิด -Wextra แบบที่เห็นในหน้า 3.4 จะได้

abs.c:8:29: warning: using integer absolute value function 'abs' when argument is of floating-point type 'double' [-Wabsolute-value]
    8 |     printf("change = %d\n", abs(temp_noon - temp_morning));
      |                             ^~~

แบบที่ 3 คูณก่อนหารจนล้น

ในโปรแกรม gcd เขียน ค.ร.น. ตามสูตรตรงๆ เป็น top * bottom / g

#include <stdio.h>

int gcd(int a, int b)
{
    while (b != 0) {
        int remainder = a % b;
        a = b;
        b = remainder;
    }
    return a;
}

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

60000 90000

multiply first = 36834
divide first = 180000

ค.ร.น. ของ 60000 กับ 90000 คือ 180000 ซึ่ง int เก็บได้สบาย แต่ top * bottom เป็นห้าพันสี่ร้อยล้าน เกิน int ไปแล้วตั้งแต่ก่อนหาร ล้นแบบหน้า 1.3 คราวนี้ค่าที่วนกลับบังเอิญเป็นบวก ผลจึงดูเหมือนตัวเลขปกติ ซึ่งอันตรายกว่าค่าติดลบที่เห็นแล้วรู้ทันทีว่าผิด การหารด้วย g ก่อนทำให้ผลระหว่างทางไม่เคยใหญ่กว่าคำตอบ และหารลงตัวเสมอเพราะ g หาร top ลงตัวแน่นอน

แบบที่ 4 ลืมเช็ก n < 2

ในโปรแกรม twin ลบสามบรรทัดแรกของ is_prime ออก แล้วลองตรวจเลข 0 ถึง 5

#include <stdio.h>

int is_prime(int n)
{
    for (long long i = 2; i * i <= n; i++) {
        if (n % i == 0) {
            return 0;
        }
    }
    return 1;
}

int main()
{
    for (int n = 0; n <= 5; n++) {
        printf("%d %d\n", n, is_prime(n));
    }
    return 0;
}

0 กับ 1 ได้ 1 ซึ่งแปลว่าเป็นจำนวนเฉพาะ ผิดทั้งคู่ เพราะ 2 * 2 <= 0 และ 2 * 2 <= 1 เป็นเท็จ ลูปไม่ทำงานเลย แล้วตกลงมาที่ return 1; ชุดทดสอบที่มีเลข 1 อยู่ในช่วงจะจับบั๊กนี้ได้ทันที ทุกครั้งที่เขียนฟังก์ชันเสร็จ ให้ลองเรียกกับค่าน้อยสุดที่โจทย์อนุญาตก่อน

สรุป

math.h มี sqrt pow fabs ทุกตัวรับและส่ง double ส่วน abs สำหรับ int อยู่ใน stdlib.h
pow ใช้กับงานทศนิยม คำตอบเป็นจำนวนเต็มใหญ่ให้เขียน power เองด้วย long long
เทียบระยะทางใช้กำลังสองของระยะ ไม่ต้อง sqrt
is_prime ต้องจัดการ n < 2 แยก และลูปถึงแค่ i * i <= n
gcd แบบยุคลิดคือแทน (a, b) ด้วย (b, a % b) จน b เป็น 0 ค.ร.น. คือ a / gcd(a, b) * b หารก่อนคูณ
ฟังก์ชันที่ตอบใช่หรือไม่ ตั้งชื่อขึ้นต้นด้วย is_ แล้วใช้ในเงื่อนไขได้ตรงๆ

กลับไปเช็กรายการ ผ่านบทนี้เมื่อ แล้วไปทำแบบฝึกหัด