Bỏ qua điều hướng, tới nội dung chính
Học C
Bài 9.717 phút đọc

Tổng và trung bình

Sau bài này bạn sẽ làm được

  • Chọn kiểu tích lũy đủ rộng để không tràn
  • Ép kiểu đúng chỗ khi tính trung bình
  • Hiểu vì sao cộng dồn số thực mất chính xác
  • Dùng thuật toán Kahan khi cần độ chính xác cao

Tính tổng và trung bình là bài toán đơn giản nhất trên mảng, và cũng là nơi ba loại lỗi số học khác nhau cùng gặp nhau: tràn số nguyên, chia nguyên mất phần lẻ, và sai số tích lũy của số thực.

#Hai hàm cơ bản

tong-tb.c
long long tong(const int *a, size_t n)
{
    long long s = 0;                    /* kiểu tích lũy rộng hơn kiểu phần tử */

    for (size_t i = 0; i < n; ++i)
        s += a[i];

    return s;
}

double trung_binh(const int *a, size_t n)
{
    if (n == 0) return 0.0;             /* tránh chia cho 0 */

    return (double)tong(a, n) / (double)n;
}

Bốn quyết định thiết kế trong mười dòng: kiểu tích lũy là long long chứ không phải int, kiểu trả về của trung bình là double, có kiểm tra mảng rỗng, và có ép kiểu rõ ràng trước phép chia. Ba mục dưới giải thích từng quyết định.

#Bẫy một: tràn số khi cộng

Tích lũy bằng int
int tong(const int *a, size_t n)
{
    int s = 0;

    for (size_t i = 0; i < n; ++i)
        s += a[i];          /* tràn khi tổng vượt 2 147 483 647 */

    return s;
}

/* Một triệu phần tử, mỗi phần tử 3000, tổng thật là 3 tỷ.
   Kết quả nhận được là một số âm. */
Tích lũy bằng long long
long long tong(const int *a, size_t n)
{
    long long s = 0;

    for (size_t i = 0; i < n; ++i)
        s += a[i];          /* long long chứa được tới hơn 9 tỷ tỷ */

    return s;
}
Kiểu phần tửKiểu tích lũy nên dùngGiới hạn an toàn
char, short, intlong longTổng tới hơn 9 tỷ tỷ
unsigned intunsigned long longKhông có hành vi không xác định, chỉ quay vòng
long longlong long, và phải tự kiểm traKhông có kiểu nào rộng hơn theo chuẩn
floatdoubleVừa rộng hơn vừa chính xác hơn
doubledouble, hoặc dùng thuật toán KahanVấn đề là độ chính xác chứ không phải tràn

Khi ngay cả long long cũng không đủ

#include <limits.h>

/* Cộng có kiểm tra tràn. Trả về 0 nếu an toàn, âm một nếu sẽ tràn. */
int cong_an_toan(long long a, long long b, long long *kq)
{
    if (b > 0 && a > LLONG_MAX - b) return -1;
    if (b < 0 && a < LLONG_MIN - b) return -1;

    *kq = a + b;

    return 0;
}

Điểm mấu chốt là kiểm tra trước khi cộng. Viết if (a + b < a) để phát hiện tràn là sai, vì phép cộng gây tràn đã xảy ra rồi và hành vi từ đó là không xác định.

#Bẫy hai: chia nguyên

Chia hai số nguyên
int a[3] = { 1, 2, 2 };

double tb = tong(a, 3) / 3;      /* 5 / 3 = 1 vì cả hai đều là số nguyên */

printf("%.2f\n", tb);            /* 1.00, sai */
Ép về double trước khi chia
int a[3] = { 1, 2, 2 };

double tb = (double)tong(a, 3) / 3;

printf("%.2f\n", tb);            /* 1.67, đúng */

Trong C, phép chia giữa hai số nguyên là phép chia nguyên, phần lẻ bị cắt bỏ chứ không làm tròn. Việc gán kết quả vào biến double không cứu được, vì phép chia đã xảy ra xong trước khi gán.

Cách viếtKết quả với 5 và 3Đúng không
5 / 31Sai, mất phần lẻ
(double)5 / 31.666667Đúng
5 / (double)31.666667Đúng
(double)(5 / 3)1.000000Sai, ép kiểu quá muộn
5.0 / 31.666667Đúng, cách viết gọn nhất

Chia cho 0

double trung_binh(const int *a, size_t n)
{
    if (n == 0) return 0.0;         /* bắt buộc phải có */

    return (double)tong(a, n) / (double)n;
}

Chia số nguyên cho 0 làm chương trình chết ngay bằng tín hiệu. Chia số thực cho 0 lại không chết, mà cho ra inf hoặc nan. Cả hai đều không phải điều bạn muốn, nên hãy chặn từ đầu.

terminal
./chia-0
so nguyen: Floating point exception (core dumped)
./chia-0-thuc
so thuc: inf
0.0 / 0.0 = -nan

#Bẫy ba: sai số cộng dồn số thực

sai-so.c
#include <stdio.h>

int main(void)
{
    double s = 0.0;

    for (int i = 0; i < 10000000; ++i)
        s += 0.1;

    printf("%.6f\n", s);
    printf("mong doi: 1000000.000000\n");

    return 0;
}
terminal
./sai-so
999999.999839
mong doi: 1000000.000000

Kết quả lệch 0.000161 sau mười triệu phép cộng, và độ lệch đó lớn dần theo số phép cộng. Nguyên nhân là số 0.1 không biểu diễn chính xác được ở hệ nhị phân, như Bài 2.3 đã nói. Mỗi phép cộng làm tròn một chút, và các lần làm tròn cộng dồn lại.

Thuật toán Kahan

kahan.c
/* Giữ lại phần bị mất do làm tròn trong biến bu, cộng lại ở lượt sau. */
double tong_kahan(const double *a, size_t n)
{
    double s  = 0.0;
    double bu = 0.0;      /* phần đã mất, chờ được bù */

    for (size_t i = 0; i < n; ++i) {
        double y = a[i] - bu;      /* bù phần mất của lượt trước */
        double t = s + y;          /* phép cộng này làm mất phần thấp của y */

        bu = (t - s) - y;          /* tính chính xác phần vừa mất */
        s  = t;
    }

    return s;
}
terminal
./kahan
cong thuong: 999999.999839
cong Kahan:  1000000.000000
Cách cộngSai số tương đốiChi phí
Cộng thẳngTăng theo nn phép cộng
Cộng theo cặp, kiểu chia để trịTăng theo log nn phép cộng, cần đệ quy
KahanGần như không tăng theo n4n phép toán
Sắp xếp tăng dần rồi cộngGiảm đáng kểThêm chi phí sắp xếp

#Phương sai và độ lệch chuẩn

Đây là bài tập tổng hợp cả ba cái bẫy trên, và là hàm bạn sẽ dùng lại nhiều lần.

Công thức một lượt, mất chính xác nghiêm trọng
/* Công thức trung bình bình phương trừ bình phương trung bình.
   Khi phương sai nhỏ so với trung bình, hai số lớn gần bằng nhau bị trừ đi,
   và kết quả mất gần hết chữ số có nghĩa. Có thể ra số âm. */
double phuong_sai_xau(const double *a, size_t n)
{
    double s = 0.0, s2 = 0.0;

    for (size_t i = 0; i < n; ++i) {
        s  += a[i];
        s2 += a[i] * a[i];
    }

    double tb = s / n;

    return s2 / n - tb * tb;
}
Hai lượt, ổn định
/* Lượt một tính trung bình, lượt hai tính tổng bình phương độ lệch. */
double phuong_sai(const double *a, size_t n)
{
    if (n < 2) return 0.0;

    double tb = 0.0;

    for (size_t i = 0; i < n; ++i)
        tb += a[i];

    tb /= (double)n;

    double s = 0.0;

    for (size_t i = 0; i < n; ++i) {
        double d = a[i] - tb;

        s += d * d;
    }

    return s / (double)(n - 1);    /* chia n-1 cho phương sai mẫu */
}

double do_lech_chuan(const double *a, size_t n)
{
    return sqrt(phuong_sai(a, n));
}
terminal
# Với dãy 100000000.0, 100000000.1, 100000000.2
./phuong-sai
cong thuc mot luot: -2.000000   <- am, vo nghia
hai luot:            0.010000   <- dung

Tự làm thử

  1. Tạo mảng một triệu phần tử int giá trị 3000, tính tổng bằng int và bằng long long, so sánh kết quả.
  2. Chạy phiên bản tràn số dưới -fsanitize=undefined và chép lại thông báo.
  3. In kết quả của cả năm cách viết trong bảng phép chia và giải thích từng dòng.
  4. Cộng 0.1 mười triệu lần bằng float, bằng double và bằng Kahan, lập bảng so sánh sai số.
  5. Cài cả hai công thức phương sai, thử với dãy 100000000.0, 100000000.1, 100000000.2 và giải thích vì sao công thức một lượt ra số âm.
  6. Viết hàm trả về trung vị của mảng, xử lý đúng cả trường hợp số phần tử chẵn.

Trình chấm điểm tự động sẽ được bổ sung ở giai đoạn sau. Hiện tại bạn tự chạy thử trên máy.

Tóm tắt

  • Kiểu tích lũy phải rộng hơn kiểu phần tử. Dùng long long khi cộng mảng int.
  • Tràn số nguyên có dấu là hành vi không xác định, phải kiểm tra trước khi cộng chứ không phải sau.
  • Phép chia giữa hai số nguyên cắt bỏ phần lẻ. Ép kiểu phải đặt lên toán hạng, không đặt lên kết quả.
  • Cộng dồn số thực tích lũy sai số theo số phần tử. Thuật toán Kahan giữ lại phần bị làm tròn và bù vào lượt sau.
  • Tính phương sai bằng hai lượt duyệt, không dùng công thức trung bình bình phương trừ bình phương trung bình.