新生 C 语言大赛

• 赛事平台
  • 赛事主页
  • 评测系统
  • 赛事文件
  • 赛事主页(第10届)
• 赛事信息
  • 赛程安排
  • 赛事规则
• 赛事培训
  • 第零次培训
  • 第一次培训
  • 第二次培训
  • 第三次培训
  • 决赛培训

重要通知

初赛截止时间改变

考虑到初赛题目难度可能较大, 决赛入围标准较为严格, 当前已经把初赛截止时间延长到了与决赛相同, 不再要求在 02-15 完成初赛, 但仍要求初赛达到 800 分, 否则无法参与决赛结算.

决赛附加题目

考虑到有超过 1 名选手在原决赛题目中获得了满分, 故增加 1 道附加题目, 与原题目基本一致, 仅参数范围发生变化, 附加总分值为 20 .

赛程安排

赛事规则

初赛规则

初赛在 UOJ 上进行, 一共有 2 组题目, 每组 6 道, 每道题目共 10 个测试点 100 分, 合计 1200 分. 在初赛截止前, 获得 800 分满足决赛结算标准, 不要求完整做出 8 道题目.

初赛超过 800 分的选手, 超出的分数会乘上 $5$ 后计入决赛附加分.

决赛规则

决赛在 UOJ 上进行, 有 1 道主题目, 为 100 分; 在该主题目基础上, 延伸出了 1 道附加题目, 为 20 分.

第零次培训: 安装 C 环境

查看宣讲会与第零次培训的回放与文件.

工物科协为大家提供以下 5 种搭建 C 开发环境的方案, 以下均针对 windows 用户, 如果您为 linux 或 macos 用户, 应该能够更为简单的搭建 C 开发环境.

方案1: 使用 DevC++ 集成开发环境

注意: DevC++ 安装方便, 比较适合刚入门编程的同学, 但开发体验一般.

安装与设置

文件名为 Dev-Cpp 5.11 TDM-GCC 4.9.2 Setup.exe, 访问 sourceforge 上点击 Download 即可下载, 或者可以通过科协文件服务器提供的镜像下载.

下载后双击下载文件, 依次点击 OK, I agree, Next, (选择安装位置) Install 即可完成安装.

安装后点击 Finish (一般会自动打开 Dev-c++, 没有的话自己打开桌面的快捷方式), 随后选择 简体中文, 点击 Next, Next, OK 即可完成设置.

使用测试

在打开的程序左上角点击 文件, 新建, 源代码.

将以下示例代码复制到占较大面积的白色输入框中:

#include <stdio.h>

int main() {
    printf("Hello, World!\n");
    return 0;
}

按键盘上的 F9, 选择一个文件夹保存源代码, 修改文件名后按下 保存, 再按下键盘上的 F10 即可编译运行.

方案2: 使用 Visual Studio 集成开发环境

注意: Visual Studio 操作简单易上手, 功能齐全, 但体积较大.

可参考VS安装教程视频, 视频中的安装链接为 https://visualstudio.microsoft.com/ .

方案3-5: 使用 VSCode 编辑器

注意: Visual Studio Code (编辑器) 与 Visual Studio (集成开发环境) 是两个完全不同的软件.

方案3: 使用 MinGW 作为 windows 上的编译器

可参考VSCode安装教程视频, 视频中的安装链接为 https://code.visualstudio.com/ , 视频中的 MinGW-W64 的 github 下载链接为 https://github.com/niXman/mingw-builds-binaries/releases/ , 清华云盘中的文件已转移到我们的赛事文件中.

方案4: 安装 WSL Debian 发行版在 linux 环境下开发

注意: 使用 WSL 要求了解 Linux 操作.

有关 WSL2 的安装, 可以参考 WSL2 安装教程; 有关 linux 基础, 可以参考 Linux 基础教程.

简单的来说, 你可以在 windows 终端 (cmd 或 PowerShell) 中输入以下命令安装 WSL2 :

wsl --install -d Debian

安装完成后, 需要设置用户名与密码 (注意设置密码时, 输入的密码不会出现在屏幕上), 然后即可进入 WSL. 可以使用终端进行开发, 也可以使用 VSCode 的 WSL 插件来连接到 WSL.

可以使用以下命令来安装 gcc 等编译工具:

sudo apt update
sudo apt install build-essential

对于 gcc 相关命令, 可以参考菜鸟教程. 在初学阶段, 只需要知道可以用以下命令把源代码编译为可执行文件 (可以加 -O2 等优化参数)：

gcc main.c -o main

然后使用如下命令运行:

./main

方案5: 配置远程服务器并使用 ssh 连接

注意: 远程服务器只在网络访问稳定的情况下可用, 且要求了解 Linux 操作.

科协将会为有需要的选手提供服务器账号, 服务器上已经预装了 gcc, 请有需要的选手与比赛负责人私信联系, 负责人将会为您在服务器上创建临时账号 (赛后将收回) 并指导连接.

在远程服务器中可以使用终端进行开发, 也可以使用 VSCode 的 Remote-SSH 插件连接, 获得与本地开发相接近的体验.

第一次培训: C 语言基础语法

查看第一次培训的回放与文件.

第一次培训主要讲解: C语言数据类型与变量、表达式与语句、数组、输入输出、控制流、作用域、函数等.

第二次培训: 编写更复杂的程序

查看第二次培训的回放与文件.

指针

基本知识

定义

指针也就是内存地址, 指针变量是用来存放内存地址的变量, 就像其他变量或常量一样, 必须在使用指针存储其他变量地址之前, 对其进行声明. 指针变量声明的一般形式为:

type *var_name;

在这里, type 是指针的基类型, 它必须是一个有效的 C 数据类型, var_name 是指针变量的名称, 例如:

int* p1;  // int 类型的指针
double *p2, *p3; // 两个 double 类型的指针
char *p4, c5; // 一个 char 类型的指针, 一个 char 类型的变量

所有实际数据类型, 对应指针的值都是一个代表内存地址的长的十六进制数, 对于 64 位的系统, 占 8 字节. 不同数据类型的指针之间唯一的不同是, 指针所指向的变量或常量的数据类型不同. 绝大多数现代通用计算机系统都是小端系统, 即内存地址从低到高排布.

使用

使用指针时会涉及以下几个操作:

• 声明一个指针变量;
• 把变量地址赋值给指针 (通过取地址运算符 & );
• 访问指针变量中可用地址的值.

一个例子是:

// pointer.c
#include <stdio.h>
int main(){
    int a = 20;
    int *p;
    p = &a;
    printf("the address of a: %p\n", &a);
    printf("the value of p: %p\n", p);
    printf("the size of p: %zu\n", sizeof(p));
    printf("the value of *p: %d\n", *p);
    return 0;
}

注: 可以尝试注释掉 p = &a; 一句, 观察未初始化的指针值.

如果要修改变量的值, 既可以直接修改, 也可以修改指针指向的值, 二者会同步变动, 例如:

// pointer_change.c
#include <stdio.h>
int main(){
    int a = 20;
    int *p = &a;
    printf("a, *p:  %d, %d\n", a, *p);
    a = 30;
    printf("a, *p:  %d, %d\n", a, *p);
    *p = 10;
    printf("a, *p:  %d, %d\n", a, *p);
    return 0;
}

多个指针可以指向同一个变量, 例如:

// pointers.c
#include <stdio.h>
int main(){
    short a = 1 << 7;
    short *p1 = &a, *p2 = &a;
    char *pc = &a;
    int *pi = &a;
    printf("*p1, *p2:  %d, %d\n", *p1, *p2);
    printf("*pc:  %d\n", *pc);
    printf("*pi:  %d\n", *pi);
    return 0;
}

思考: *pc 与 *pi 的值各是什么含义.

空指针 NULL

空指针 NULL 是一个定义在 <stddef.h> 库中的值为零的常量:

#define NULL ((void *)0)

通过以下代码可以查看 NULL 的值:

// NULL.c
#include <stdio.h>
int main(){
    int *p = NULL;
    printf("the value of p (%%p): %p\n", p);
    printf("the value of p (%%u): %u\n", p);
    printf("the value of *p: %d\n", *p);
    return 0;
}

NULL 指向的地址的的值是不可访问的, 因此如果一个指针不用了, 或者在其值的修改过程中出现了错误, 可以选择将其值置为 NULL . 可以参考以下方法判断指针是否为空:

// judge_NULL.c
#include <stdio.h>
void judge(int *p){
    if(p) printf("%p\n", p);
    if(!p) printf("NULL\n");
}
int main(){
    int a = 1;
    int *p1 = &a, *p2 = NULL;
    judge(p1);
    judge(p2);
    return 0;
}

运算

指针可以进行加减运算, 即支持以下四种运算符: +, -, ++, --. 不过, 指针的加减运算与一般的 C 数据类型不同, 对于基类型为 type 的指针, 其加减的单位不是 1 , 而是 sizeof(type) . 例如:

// pointer_pm.c
#include <stdio.h>
int main(){
    int a = 1 << 8;
    int *pi = &a;
    char *pc = &a;
    printf("pi:      %p\n", pi);
    printf("pi - 3:  %p\n", pi - 3);
    printf("--pi:    %p\n", --pi);
    printf("pc + 2:  %p\n", pc + 2);
    printf("++pc:    %p\n", ++pc);
    return 0;
}

思考: 这段程序运行完后, *pc 的值是什么? 是随机值还是确定值?

指针作为函数参数

由于函数接收的参数为形式参数, 只接收外部实参的值, 而参数名与外部参数没有关系, 所以通常来讲, 无法在函数内部修改外部参数的值, 如:

// fake_swap.c
#include <stdio.h>
void swap(int a, int b){
    int tmp = a;
    a = b;
    b = tmp;
}
int main(){
    int a = 1, b = 2;
    swap(a, b);
    printf("a, b:  %d, %d\n", a, b);
    return 0;
}

而使用指针, 则可以实现修改外部参数的值的效果, 如:

// swap.c
#include <stdio.h>
void swap(int *a, int *b){
    int tmp = *a;
    *a = *b;
    *b = tmp;
}
int main(){
    int a = 1, b = 2;
    swap(&a, &b);
    printf("a, b:  %d, %d\n", a, b);
    return 0;
}

需要注意的是, 不要把这个当成指针的特性, 而是要从原理层面理解, 对于 a, b 来说, 传入的同样是实参, 指针 a, b 的值 (即其指向的地址) 同样是不可修改的, 但如果在函数内修改了指向地址的值, 则在函数外时可以访问修改后的值的. 以下代码说明了在函数中修改指针本身的值同样是不会影响外部指针的值的:

// pointer_change_func.c
#include <stdio.h>
void change(int *a){
    a++;
}
int main(){
    int *p = (int*)(0x04);
    printf("p (original): %p\n", p);
    change(p);
    printf("p (changed):  %p\n", p);
    return 0;
}

思考: 如何使上述代码可以修改 p 的值?

指针与数组

指针与数组的对应关系

指针与数组有着对应关系:

• 数组名表示数组的地址, 即数组首元素的地址;
• 数组与指针的 [] 运算符与 + 运算符后取值是互通的.

如下所示:

// array.c
#include <stdio.h>
int main(){
    int a[] = {1, 2, 3};
    int *p = &a[0];
    printf("a: %p\n", a);
    printf("p: %p\n", p);
    printf("*(a + 1): %d\n", *(a + 1));
    printf("p[1]:     %d\n", p[1]);
    printf("sizeof(a): %zu\n", sizeof(a));
    printf("sizeof(p): %zu\n", sizeof(p));
    return 0;
}

从上述例子中可以看出 a 与 p 基本可以转换使用, 但指针不包含数组的大小信息.

动态数组

静态数组与动态数组的区别:

• 静态数组在编译时分配内存, 大小固定; 动态数组在运行时手动分配内存, 大小可变.
• 静态数组的内存通常分配在栈上, 随着函数的调用和返回而自动管理; 动态数组的内存空间在运行时通过动态内存分配函数手动分配, 并存储在堆上.
• 静态数组的生命周期始于其定义时, 自动终于在作用域中无法被访问时; 动态数组的生命周期由程序员控制, 需要在使用完数组后手动释放内存, 以避免内存泄漏.

可以使用 <stdlib.h> 库中的函数管理内存:

• void* malloc(size_t size) 函数: 分配所需的内存空间, 并返回一个指向它的指针.
• void* calloc(size_t nitems, size_t size) 函数: 分配所需的内存空间, 设置分配的内存为零, 并返回一个指向它的指针.
• void* realloc(void *ptr, size_t size) 函数: 重新分配内存空间.
• void free(void *ptr) 函数: 释放动态分配函数分配的内存.

当需要一个较大的内存空间时, 通常要使用动态数组, 因为栈的内存是有限的, 可能无法分配足够大的空间, 例如:

// static.c
int main(){
    int N = 100000000;
    int a[N];
    return 0;
}

// dynamic.c
#include <stdlib.h>
int main(){
    int N = 100000000;
    int *a = (int*)malloc(N * sizeof(int));
    if(!a) return 1;
    free(a);
    return 0;
}

其中的 static.c 代码就会出现 Segmentation fault (core dumped) .

以下给出一个二维动态数组的例子:

// pointer2.c
#include <stdio.h>
#include <stdlib.h>
int main(){
    int m = 3, n = 5;
    int **arr = (int**)malloc(m * sizeof(int*));
    for(int i = 0; i < m; i++){
        arr[i] = (int*)calloc(n, sizeof(int));
    }
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            arr[i][j] = 10 * i + j;
        }
    }
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            printf("%02d ", *(arr[i] + j));
        }
        printf("\n");
    }
    for(int i = 0; i < m; i++){
        free(arr[i]);
    }
    free(arr);
    return 0;
}

使用 malloc 等函数分配完的内存, 一定要养成使用 free 释放的习惯! 如果不这样, 可能会出现内存泄露, 即申请的内存被占用但却无法正常访问, 也无法释放, 直到进程结束, 例如 (谨慎运行):

// memory_leak.c
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
void allocate(){
    char *p = malloc(100 * 1024 * 1024);
    memset(p, 0, 100 * 1024 * 1024);
}
int main(){
    for(int i = 0; i < 200; i++){
        allocate();
        printf("%d\n", i);
    }
    return 0;
}

数组作为函数参数或返回值

数组作为函数参数, 实际上是把数组名隐式转化为指针传入函数, 即把数组的首地址传入, 而非把整个数组的值传入, 因此, 在函数内部是可以修改数组的, 例如:

// array_change_func.c
#include <stdio.h>
void change(int a[]){
    printf("%zu\n", sizeof(a));
    a[1] = 4;
}
int main(){
    int a[] = {1, 2, 3};
    printf("%zu\n", sizeof(a));
    change(a);
    for(int i = 0; i < 3; i++){
        printf("%d ", a[i]);
    }
    printf("\n");
    return 0;
}

如果要传递多维数组, 需要把内部的维数写明, 而最外层的维数只做标识用, 如以下代码一般情况下可以正常运行:

// array2_func.c
#include <stdio.h>
void change(int a[2][3]){
    a[2][2] = 1;
}
int main(){
    int a[3][3] = {0};
    change(a);
    printf("a[2][2]: %d\n", a[2][2]);
    return 0;
}

数组不能直接作为函数的返回值类型, 如果想要返回一个在函数中新定义的数组, 可以使用指针返回, 例如:

// array_return.c
#include <stdio.h>
#include <stdlib.h>
int* new_array(int n){
    int *a = (int*)malloc(n * sizeof(int));
    for(int i = 0; i < n; i++){
        a[i] = i;
    }
    return a;
}
int main(){
    int n = 10;
    int *a = new_array(n);
    for(int i = 0; i < n; i++){
        printf("%d\n", a[i]);
    }
    free(a);
    return 0;
}

一个应用是用函数封装一个二维数组内存分配函数:

// allocate.c
#include <stdio.h>
#include <stdlib.h>
int** alloc2(int m, int n){
    int **arr = (int**)malloc(m * sizeof(int*));
    for(int i = 0; i < m; i++){
        arr[i] = (int*)calloc(n, sizeof(int));
    }
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            arr[i][j] = 10 * i + j;
        }
    }
    return arr;
}
void free2(int** arr, int m){
    for(int i = 0; i < m; i++){
        free(arr[i]);
    }
    free(arr);
}
int main(){
    int m = 3, n = 5;
    int **arr = alloc2(m, n);
    for(int i = 0; i < m; i++){
        for(int j = 0; j < n; j++){
            printf("%02d ", arr[i][j]);
        }
        printf("\n");
    }
    free2(arr, m);
    return 0;
}

值得注意的是, 如果在函数中只分配不返回, 函数结束后是无法访问的, 反而会造成内存泄露, 如:

// wrong_allocate.c
#include <stdio.h>
#include <stdlib.h>
void alloc(int *a, int n){
    printf("a (original):   %p\n", a);
    a = (int*)malloc(n * sizeof(int));
    printf("a (allocated):  %p\n", a);
    for(int i = 0; i < n ; i++){
        a[i] = i;
    }
}
int main(){
    int n = 10;
    int *a = &n;
    alloc(a, n);
    printf("a (allocated?): %p\n", a);
    for(int i = 0; i < n; i++){
        printf("%d\n", a[i]);
    }
    free(a);
    return 0;
}

其它数据类型

结构体 struct

结构体是 C 语言中一种用户自定义的可用的数据类型, 它可以存储不同类型的数据项, 其数据成员可以是基本数据类型 (如 int, float, char 等), 也可以是其他结构体类型, 指针类型等.

可以使用 . 运算符访问结构体的成员变量, 也可以使用 -> 运算符访问结构体指针的成员变量 ( a->b 等价于 (*a).b ), 例如:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
struct Person{
    char name[32];
    unsigned age;
};
int main(){
    struct Person p1 = {"dmr", 70};
    struct Person *p2 = malloc(sizeof *p2);
    strcpy(p2->name, "ken");
    p2->age = 83;
    printf("%s, %d\n", p1.name, p1.age);
    printf("%s, %d\n", p2->name, p2->age);
    return 0;
}

结构体的赋值操作为逐成员按值拷贝, 这意味着每个所有成员变量的值会被复制一份, 而非两个结构体只是不同的名称, 但是指针成员变量会指向同一块内存. 例如:

// struct_copy.c
#include <stdio.h>
struct Nums{
    int a;
    int *b;
};
int main(){
    int a = 1, b = 2;
    struct Nums n1 = {a, &b};
    struct Nums n2 = n1;
    n2.a = 3;
    *(n2.b) = 4;
    printf("n1: %d, %d\n", n1.a, *n1.b);
    printf("n2: %d, %d\n", n2.a, *n2.b);
    return 0;
}

共用体 union

共用体是一种特殊的数据类型, 可以在相同的内存位置存储不同的数据类型, 一个共用体可以有多个成员变量, 但任何时候只能有一个成员变量带有有效值. 共用体提供了一种使用相同的内存位置的有效方式, 例如:

// union.c
#include <stdio.h>
union Nums{
    int a;
    short b;
}nums = {1 << 15};
int main(){
    printf("the size of Nums: %zu\n", sizeof(union Nums));
    printf("%d, %d\n", nums.a, nums.b);
    nums.a = 1 << 16;
    printf("%d, %d\n", nums.a, nums.b);
    return 0;
}

枚举 enum

枚举是 C 语言中的一种基本数据类型, 通常用于为程序中的一组相关的常量取名字, 以便于程序的可读性和维护性. 第一个枚举成员的默认值为整型的 0 , 后续枚举成员的值在前一个成员上加 1 , 以下只给出一个简单的例子:

// enum.c
#include <stdio.h>
enum Weekdays{
    MON=1, TUE, WED, THU, FRI, SAT=11, SUN
};
int main(){
    for(int i = MON; i <= FRI; i++){
        printf("%d\n", i);
    }
    printf("%d\n%d\n", SAT, SUN);
    return 0;
}

typedef 关键字

typedef 关键字可用来为类型取一个新的名字, 譬如在 <types.h> 库中, 就定义了 __int8_t, __uint8_t, __int16_t, __uint16_t 等类型, 分别对应于 char, unsigned char, short, unsigned short 等; 在 <stddef.h> 库中, 定义了 size_t 类型, 是用来表示对象大小, 内存大小, 元素个数的专用无符号类型.

typedef 的一个常用用途是给结构体取别名, 使得在使用时不需要再在前面加上 struct 关键字, 例如:

// typedef.c
#include <stdio.h>
#include <stdlib.h>
typedef struct Vector{
    size_t n;
    int *data;
}Vector;
int main(){
    size_t n = 5;
    Vector a;
    a.n = n;
    a.data = (int*)malloc(n * sizeof(int));
    for(int i = 0; i < a.n; i++){
        a.data[i] = i;
    }
    for(size_t i = 0; i < a.n; i++){
        printf("%d\n", a.data[i]);
    }
    free(a.data);
    return 0;
}

补充知识

main 函数参数

与其它函数类似, main 函数也可以有参数, 对于 main.c 编译好的程序 main , 可以通过

./main args1 args2 ...

运行, 在 main 函数中:

• int argc 表示参数的数目;
• char **argv 表示参数列表, 为字符串类型;

需要注意的是, 程序会默认把程序名作为第一个参数 (argv[0]), 对应上例即 ./main .

以下给出一个可以接收并打印所有参数的代码示例:

// main_args.c
#include <stdio.h>
int main(int argc, char **argv){
    printf("argc: %d\n", argc);
    for(int i = 0; i < argc; i++){
        printf("argv[%d] = %s\n", i, argv[i]);
    }
    return 0;
}

预处理器

预处理器是编译过程中的独立阶段, 在实际编译前对源代码进行文本处理, 主要功能包括:

• 宏展开;
• 头文件包含;
• 条件编译;
• 特殊指令处理.

可以用 #define 定义宏, 在编译过程中, 会直接进行文本替换; 可以用 #undef 删除对某个宏的定义 (需要注意不要随意 #undef 标准库中定义的宏, 否则很有可能发生编译错误). 比如:

#define MY_INT_MAX 0x7fffffff
#undef MY_INT_MAX

可以用 #include 引入头文件.

可以用 #if, #ifdef, #ifndef 等来控制条件编译, 一个可能的应用场景是头文件保护, 即防止同一个头文件被多次包含, 甚至两个头文件互相包含导致错误的问题. 标准库中已经做好了头文件保护, 如果自己写的头文件, 可以按照如下格式做头文件保护:

#ifndef MY_HEADER_H
#define MY_HEADER_H

/* 头文件代码 */

#endif

在现代主流编译器, 对于头文件保护有更简单的写法:

#pragma once

内存布局

一个运行中的 C 程序, 其内存通常可以分为五大区:

• 代码区
• 常量 / 只读数据区
• 全局 / 静态区
• 堆 (Heap)
• 栈 (Stack)

栈

栈是一块连续的内存区域, 专门用于函数调用管理和局部数据存储, 其特点为:

• 先进后出;
• 自动分配和释放, 不需要手动管理;
• 内存有限制, 通常为 MB 量级.

在栈中存放的有:

• 局部变量;
• 函数参数;
• 函数返回地址与某些寄存器的状态.

每次函数调用都会在栈上创建一个栈帧(Stack Frame), 也叫活动记录(Activation Record). 在大部分现代系统中, 栈帧从高地址向低地址增长. 函数调用结束时会弹出栈帧, 局部变量会自动销毁, 即生命周期结束.

栈有两个常见的问题: 一个为栈溢出(Stack Overflow), 通常由递归太深或局部数组过大导致; 另一个为野指针/悬空指针, 即指向已释放内存区域的指针, 例如:

// dangling_pointer.c
#include <stdio.h>
int* foo(){
    int x = 1;
    printf("&x: %p\n", &x);
    return &x;
}
int main(){
    int *p = foo();
    printf("p:  %p\n", p);
    printf("*p: %d\n", *p);
    return 0;
}

堆

堆是程序运行时用于动态分配内存的一块区域, 其与栈相比的特点为:

• 变量的生命周期由程序员控制, 在程序结束前不会自动回收.
• 内存较大, 基本与系统的内存相当.
• 堆的地址一般比栈低, 栈和堆从不同方向增长, 即堆是从低地址向高地址增长.
• 堆的内存不一定连续, 可能会因多次 malloc 与 free 产生碎片化空间.
• 堆的内存分配比栈涉及更多流程, 故分配的速度比栈慢.

运行错误分析

以下只给出一些常见的运行时报错, 由于现代的 IDE 或编辑器插件都具有代码分析功能, 所以编译时报错的问题通常可被直接指出 (如下划红曲线).

根据具体环境不同, 有些未定义行为 (UB) 在某些平台会被认为是运行错误.

算数错误

算术错误常发生于除以 0 的情形, 典型的报错信息为:

Floating point exception (core dumped)

段错误

段错误一般发生于违规访问的情形, 如数组越界访问, 空指针解引用等, 典型的报错信息为:

Segmentation fault (core dumped)

栈溢出也会导致报错:

Segmentation fault (stack overflow)

第三次培训: 数据结构

查看第三次培训的回放与文件.

第三次培训主要讲解: 数据结构基本概念、时间复杂度、向量、列表、栈与队列、散列表、树、优先级队列等.

对于 C++ 版本的介绍, 可参考 第十届比赛的资料.

决赛培训: 数值计算求解中子扩散方程

查看决赛培训的回放与文件.

有限差分法

问题背景

考虑一个常系数稳态扩散方程

$$- \lambda \nabla^2 T = f$$

式中的 $T$ 为待求的场 (如温度场), $f$ 是已知的源项, $\lambda$ 是扩散系数, 在这里设常数, $\nabla^2$ 是 Laplace 算子, 在 $n$ 维空间中, 其含义为

$$\nabla^2 T = \sum_{i=1}^{n} \frac{\partial^2 T}{\partial x_{i}^2}$$

先考虑一维的情形, 则上述方程可写为

$$- \lambda \frac{\mathrm{d}^2}{\mathrm{d} x^2} T(x) = f(x)$$

对于比较复杂的方程形式或者比较复杂的源项 $f(x)$ , 可能无法简单地求出解析解, 在此情况下, 可以尝试求上述方程的数值解.

所谓有限差分法, 就是将求解域划分为差分网格, 用有限个网格节点代替连续的求解域, 把导数用网格节点上的函数值的差商代替进行离散, 建立以网格节点上的值为未知数的代数方程组, 从而将微分问题变为代数问题, 从而可利用计算机求解.

一维问题的网格划分与导数表示

有限差分法所研究的节点位于网格顶点 (与此相对, 有限体积法的节点位于网格中心), 下图即为一个一维的网格划分示意图:

假设划分的网格是均匀的 (在决赛题目中亦遵从此假设), 每两个节点的间距为 $h$ , 以 $T_i$ 点为研究对象, 其一阶导数可表示为

$$T'$$i \approx \frac{T{i+1} - T_{i}}{h} \ \ \ T'i \approx \frac{T{i} - T_{i-1}}{h}

这两种表达可分别被称为前向导数与后向导数. 二阶导数可以表示为

$$T''$$i \approx \frac{T{i+1} - 2 T_{i} + T_{i-1}}{h^2}

这样, 上述一维常系数稳态扩散方程则可写为离散化的形式

$$- \frac{\lambda}{h^2} T_{i+1} + \frac{2 \lambda}{h^2} T_{i} - \frac{\lambda}{h^2} T_{i-1} = f_i$$

但是, 对于边界处的点, 一般需要额外处理, 即边界条件. 在决赛中, 我们所用的是第一类边界条件 (Dirichlet 边界条件), 即边界处的值为 $0$ , 即

$$T_{0} = T_{N_x} = 0$$

这样, 由上述 $N_{x} - 1$ 个内部方程, 与 $2$ 个边界方程, 即组成了一个由 $N_{x} + 1$ 个方程组成的线性方程组, 求解该线性方程组, 得到的解向量即为所有节点处的 $T$ 值.

有限差分法的精度分析

有限差分法的精度, 可以通过泰勒展开来分析, 以上述问题为例, 由泰勒展开可得

$$T_{i+1} = T_{i} + \frac{1}{1!} \frac{\mathrm{d} T_{i}}{\mathrm{d} x} h + \frac{1}{2!} \frac{\mathrm{d}^2 T_{i}}{\mathrm{d} x^2} h^2 + \frac{1}{3!} \frac{\mathrm{d}^3 T_{i}}{\mathrm{d} x^3} h^3 + \frac{1}{4!} \frac{\mathrm{d}^4 T_{i}}{\mathrm{d} x^4} h^4 + O(h^5) \ \ \ T_{i-1} = T_{i} - \frac{1}{1!} \frac{\mathrm{d} T_{i}}{\mathrm{d} x} h + \frac{1}{2!} \frac{\mathrm{d}^2 T_{i}}{\mathrm{d} x^2} h^2 - \frac{1}{3!} \frac{\mathrm{d}^3 T_{i}}{\mathrm{d} x^3} h^3 + \frac{1}{4!} \frac{\mathrm{d}^4 T_{i}}{\mathrm{d} x^4} h^4 + O(h^5)$$

将其代入差分的表达式, 可得

$$\frac{T_{i+1} - 2 T_{i} + T_{i-1}}{h^2} = \frac{\mathrm{d}^2 T_{i}}{\mathrm{d} x^2} + \frac{1}{12} \frac{\mathrm{d}^4 T_{i}}{\mathrm{d} x^4} h^2 + O(h^3)$$

将其代入原始方程, 可得

$$- \lambda \frac{\mathrm{d}^2 T_i}{\mathrm{d} x^2} - f_i = - \lambda \frac{\mathrm{d}^2 T_i}{\mathrm{d} x^2} + \lambda \frac{T_{i+1} - 2 T_{i} + T_{i-1}}{h^2} = \frac{\lambda}{12} \frac{\mathrm{d}^4 T_{i}}{\mathrm{d} x^4} h^2 + O(h^3)$$

由此可见, 该离散格式在空间上是二阶精度的.

通过选取更多的节点来表示导数, 可以获得更高精度的离散格式, 在决赛题目中, 由于限定了离散格式, 并以限定的离散格式来评分, 所以请不要更改离散格式, 否则即使更接近真实情况, 分数也会更低.

二维问题的网格划分与拉普拉斯算符表示

二维问题与一维问题类似, 首先, 我们需要把一个二维的平面用一维的向量来表示, 采用下图所示的方式为节点编号 (从 0 开始):

二维的拉普拉斯算符为

$$\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$$

采用五点差分格式进行离散, 使用 $[i, j]$ 表示编号为 $i + j(N+1)$ 的格点, 记 $\mathrm{d} x = \mathrm{d} y = h$ , 则对于内部节点:

$$\left(\nabla^2 T\right)$${[i, j]} = \frac{T{[i-1, j]} + T_{[i+1, j]} + T_{[i, j-1]} + T_{[i, j+1]} -4 T_{[i, j]}}{h^2}

对于第一类边界条件, 有

$$T_{[0, j]} = T_{[N, j]} = T_{[i, 0]} = T_{[i, N]} = 0$$

中子扩散方程

单群中子扩散方程

常系数单群中子扩散方程为:

$$\frac{\partial n}{\partial t} = \frac{1}{v} \frac{\partial \phi(\vec{r}, t)}{\partial t} = S(\vec{r}, t) + D \nabla^2 \phi(\vec{r}, t) - \Sigma_{\text{a}} \phi(\vec{r}, t)$$

式中 $n$ 为中子数密度, $v$ 为中子速度, $\phi = n v$ 为中子通量密度, $S$ 为中子源项, $D$ 为扩散系数, $\Sigma_{\text{a}}$ 为中子宏观吸收截面.

对于上述中子扩散方程:

• 左侧表示中子数密度的变化.
• 右侧第一项 $S(\vec{r}, t)$ 为中子的产生项, 在具体反应堆中, 由裂变反应产生, 与中子通量密度 $\phi$ 有关, 这里假定研究的区域是一个无裂变的区域, 中子源呈与中子通量密度相独立的分布.
• 右侧第二项 $D \nabla^2 \phi(\vec{r}, t)$ 为中子的泄露项, 就其推导而言, 其原始形式为 $- \nabla \cdot \vec{J}(\vec{r}, t)$ , 根据菲克定律, 中子流 $\vec{J} = - D \nabla \phi$ , 当扩散系数 $D$ 为常数时, 即得到该项表达式.
• 右侧第三项 $- \Sigma_{\text{a}} \phi(\vec{r}, t)$ 为中子的吸收项, 同样假定空间物质分布是均匀的, 物质的宏观吸收截面为常数.

当上述方程左侧为 $0$ 时, 即变为稳态常系数单群中子扩散方程:

$$- D \nabla^2 \phi(\vec{r}) + \Sigma_{\text{a}} \phi(\vec{r}) = S(\vec{r})$$

该方程的含义即为 $\text{泄露项} + \text{吸收项} = \text{产生项}$ . 决赛的前一半分数, 即为该方程的离散求解.

双群中子扩散方程

在实际的反应堆问题中, 中子的能量是有一定的分布的, 且扩散系数与宏观吸收截面等参数是与能量有关的. 为了更细致的刻画中子通量密度的分布, 常采用的方法是把整个能量区间划分为若干区, 对每个区写出一个中子扩散方程, 上述的与能量有关的参数则取区间内的平均值, 称为群常数. 最简单的划分方式即划分为双群, 即高能的快群与低能的热群.

划分为双群后, 中子扩散方程也会发生一定的变化, 需要增加中子在两群中转移的项, 用两个散射转移截面 $\Sigma_{1\to2}, \Sigma_{2\to1}$ 来表征. 在理论上, 反应堆中的中子发生散射后能量会降低, 也就是说通常只有快群向热群的散射, 而热群向快群的散射是难以发生的, 即 $\Sigma_{2\to1}$ 是接近于 $0$ 的, 但是为了使两群的中子扩散方程耦合 (在实际的反应堆问题中, 热中子发生链式反应后会产生快中子), 在决赛中假定 $\Sigma_{2\to1}$ 是非零的.

假定中子源产生的都是快中子, 则稳态常系数双群中子扩散方程为:

$$- D_1 \nabla^2 \phi_1(\vec{r}) + \Sigma_{\text{a}1} \phi_1(\vec{r}) + \Sigma_{1\to2} \phi_1(\vec{r}) - \Sigma_{2\to1} \phi_2(\vec{r}) = S(\vec{r}) \ \ \ - D_2 \nabla^2 \phi_2(\vec{r}) + \Sigma_{\text{a}2} \phi_2(\vec{r}) + \Sigma_{2\to1} \phi_2(\vec{r}) - \Sigma_{1\to2} \phi_1(\vec{r}) = 0$$

对于该方程组, 不难发现任何一个方程都无法独立求解, 为了求解该方程组, 有以下两种思路:

• 分离求解: 使用 Picard 迭代方法, 求解一个方程时, 将一个物理场视为未知量, 其他物理场视为常数, 不同方程之间不断迭代, 并传递耦合参数, 直到收敛.
• 耦合求解: 将每个方程离散得到矩阵组合起来, 形成的分块矩阵每一块表示物理场之间的耦合, 对这个矩阵进行整体求解.

稀疏矩阵

稀疏矩阵的定义

稀疏矩阵 (sparse matrix) 是一种元素大部分为零的矩阵, 反之, 如果大部分元素都非零, 则称矩阵是稠密 (dense) 的.

对于一个 $n$ 阶的方阵, 假设每行的平均非零元素数目为 $d$ , 则使用稠密矩阵的方式存储所消耗的空间为 $O(n^2)$ , 而稀疏矩阵只存储非零元素的信息, 所消耗的空间为 $O(nd)$ . 一般来讲, $d \ll n$ , 这样稀疏矩阵可以大大减少存储所用的空间, 同时也在一定程度上加快了相对应的运算速度.

对于上述的 PDE 离散后形成的线性方程组, 非耦合的条件下每行至多有 $5$ 个非零元素, 对于双群理论的耦合方程组, 每行至多有 $6$ 个非零元素, 明显是稀疏矩阵.

稀疏矩阵的存储格式

COO 格式

COO 格式 (Coordinate list) 存储一个 (行, 列, 值) ((row, column, value)) 元组的列表. 在具体的 C/C++ 实现中: 可以创建一个含有 r, c, v 三个元素的结构体, 然后用一个结构体数组来存储; 也可以直接创建三个数组 row, col, val 来存储矩阵信息, 这时需要保证元素的对应关系, 即 row[i], col[i], val[i] 对应的是同一个矩阵元的信息.

在理想情况下, 条目应先按行索引排序, 再按列索引排序, 这样可以方便矩阵元素的修改维护, 也可以方便计算. 假设 COO 中存储了 n 个矩阵元素的信息, 则随机的添加/删除/修改一个元素, 所需的平均时间复杂度为 $O(n)$ , 对于工程化的实现, 如果想要尽可能降低随机修改的时间消耗, 可以使用 DOK 格式 (Dictionary of keys), 即使用词典 ( 散列表 ) 来存储矩阵元素, 可以使用该格式进行矩阵的构建, 然后再转化为便于计算的格式.

一个稀疏矩阵的例子如下:

用 COO 格式表示为:

• row: 0, 0, 1, 1, 2, 4, 4;
• col: 1, 4, 0, 3, 2, 2, 4;
• val: 1, 2, 6, 5, 4, 9, 3.

CSR/CSC 格式

CSR 格式 (compressed sparse row, 压缩行储存) 和 CSC 格式(compressed sparse column, 压缩列储存) 是另外两种常用的稀疏矩阵储存方法. 以 CSR 为例, 矩阵被分为三个数组: 行指针 row , 列索引 col 和值 val, 行指针数组记录每行的起始位置 (row[i+1] - row[i] 即为第 i 行的非零元素数目) , 列索引数组存储每个非零元素的列索引, 值数组包含相应的非零元素值; CSC 格式与此类似.

CSR/CSC 格式与 COO/DOK 格式相比, 对顺序性有着较强的要求, 不适用于矩阵的构建 (当然, 如果已知矩阵的内容, 也可以直接按顺序强行构建); CSR/CSC 格式对于快速访问与计算有优势.

对于上述的稀疏矩阵, 用 CSR 格式表示为:

• row: 0, 2, 4, 5, 5, 7;
• col: 1, 4, 0, 3, 2, 2, 4;
• val: 1, 2, 6, 5, 4, 9, 3.

用 CSC 格式表示为:

• row: 0, 0, 1, 1, 2, 4, 4;
• col: 0, 1, 2, 4, 5, 7;
• val: 1, 2, 6, 5, 4, 9, 3.

稀疏矩阵的运算

稀疏矩阵涉及的元素, 主要是稀疏矩阵向量乘 (Sparse matrix–vector multiplication, SpMV). 以 COO 格式为例, 假设使用 row, col, val 三个数组来存储非零元素信息, 用 nnz 来表示非零元素数目, 则 $\vec{y} = A \vec{x}$ 的代码如下 (假定大小匹配且 y 各元素初值为0) :

for(int i = 0; i < nnz; i++){
    y[row[i]] += val[i] * x[col[i]];
}

线性方程组迭代求解器

决赛要求使用以迭代法为基础的算法求解线性方程组.

迭代求解的基本思想

对于一个稠密矩阵构成的小规模线性方程组, 我们在线性代数课程中学习过高斯消元法或 LU 分解法, 这种直接法在确定的时间内能够求出准确的解 (忽略计算机的舍入误差因素), 对于一个 $n$ 阶的方阵, 其时间复杂度为 $O(n^3)$ .

但是, 对于稀疏矩阵构成的大规模线性方程组, 如果仍然采用上述的直接法求解, 则会破坏掉矩阵的稀疏性, 在内存有限的条件下无法满足需求, 且计算效率与准确度均无法得到保证, 这时使用迭代法能够更好的解决问题.

迭代法的基本思想, 就是对于方程组 $A \vec{x} = \vec{b}$ , 其中 $A \in \mathbb{R}^{n \times n}$ 且可逆, $\vec{b} \in \mathbb{R}^{n}$ , 可以把矩阵 $A$ 分裂为 $A = M - N$ , 其中 $M$ 是一个方便求逆的矩阵, 则

$$M \vec{x} = \vec{b} + N \vec{x} \Longleftrightarrow \vec{x} = M^{-1} ( \vec{b} + N \vec{x} )$$

这样, 原方程可以等价的表示为一个迭代过程

$$\vec{x}^{(k+1)} = B \vec{x}^{(k)} + \vec{f}$$

其中 $B = M^{-1} N = I - M^{-1} A$ , $\vec{f} = M^{-1} \vec{b}$ . 如果迭代矩阵 $B$ 保持不变, 则称为定常迭代法.

可以证明, 若该迭代序列收敛, 则 $\lim_{k \to \infty} x^{(k)} = x^{*}$ 是原方程组的解. 迭代收敛的充要条件是, 迭代矩阵的谱半径满足

$$\rho(B) = \max_{\lambda \in \sigma(B)} |\lambda| < 1$$

充分条件是, 迭代矩阵的谱范数 (矩阵的谱半径不大于任意矩阵范数) 满足

$$|B|$$2 = \max{|\vec{x}| = 1} |B \vec{x}| = \sqrt{\rho(B^T B)} < 1

简单的迭代求解器

Jacobi 迭代法

若系数矩阵 $A$ 的对角线元素不全为 $0$ , 则可分解为

$$A = D - (L + U)$$

其中 $D$ 为 $A$ 的对角元素所构成的对角矩阵, $L$ 为严格下三角矩阵, $U$ 为严格上三角矩阵, 则 Jacobi 迭代法的迭代公式为

$$\vec{x}^{(k+1)} = (I - D^{-1} A) \vec{x}^{(k)} + D^{-1} \vec{b}$$

用矩阵元的形式可写为

$$x_{i}^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1, j \ne i}^{n} a_{ij} x_{j}^{(k)} \right)$$

每步迭代中 $i$ 从 $1$ 到 $n$ 可并行计算.

如果 $A$ 严格对角占优, 则可以证明 Jacobi 迭代法收敛.

Gauss-Seidel 迭代法

若系数矩阵 $A$ 的对角线元素不全为 $0$ , 则可分解为

$$A = (D - L) - U$$

其中 $D$ 为 $A$ 的对角元素所构成的对角矩阵, $L$ 为严格下三角矩阵, $U$ 为严格上三角矩阵, 则 Gauss-Seidel 迭代法的迭代公式为

$$\vec{x}^{(k+1)} = (D - L)^{-1} U \vec{x}^{(k)} + (D - L)^{-1} \vec{b}$$

用矩阵元的形式可写为

$$x_{i}^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j = 1}^{i-1} a_{ij} x_{j}^{(k+1)} - \sum_{j = i+1}^{n} a_{ij} x_{j}^{(k)} \right)$$

每步迭代中 $i$ 从 $1$ 到 $n$ 需按序计算.

如果 $A$ 严格对角占优, 则可以证明 Gauss-Seidel 迭代法收敛.

SOR 迭代法

若系数矩阵 $A$ 的对角线元素不全为 $0$ , 引入松弛因子 $\omega \in (0, 2)$ , 则 $A$ 可分解为

$$A = \left( \frac{1}{\omega} D - L \right) - \left( \frac{1 - \omega}{\omega} D + U \right)$$

其中 $D$ 为 $A$ 的对角元素所构成的对角矩阵, $L$ 为严格下三角矩阵, $U$ 为严格上三角矩阵, 则 SOR 迭代法 (松弛法) 的迭代公式为

$$\vec{x}^{(k+1)} = (D - \omega L)^{-1} [ (1 - \omega) D + \omega U ] \vec{x}^{(k)} + \omega (D - \omega L)^{-1} \vec{b}$$

用矩阵元的形式可写为

$$x_{i}^{(k+1)} = (1 - \omega) x_{i}^{(k)} + \frac{\omega}{a_{ii}} \left( b_i - \sum_{j = 1}^{i-1} a_{ij} x_{j}^{(k+1)} - \sum_{j = i+1}^{n} a_{ij} x_{j}^{(k)} \right)$$

每步迭代中 $i$ 从 $1$ 到 $n$ 需按序计算.

收敛速度分析

迭代法第 $k$ 步的收敛误差为

$$\vec{\varepsilon}^{,(k)} = \vec{x}^{(k)} - \vec{x}^{$$} = B^{k} (\vec{x}^{(0)} - \vec{x}^{}) = B^{k} \vec{\varepsilon}^{,(0)}

所以

$$\frac{|\vec{\varepsilon}^{,(k)}|}{|\vec{\varepsilon}^{,(0)}|} < |B^{k}|$$

因此平均每次迭代后误差的压缩率为 $|B^{k}|^{1/k}$ , 定义平均收敛速度为:

$$R_{k}(B) = - \ln |B^{k}|^{1/k}$$

渐进收敛速度定义为

$$R(B) = \lim_{k \to \infty} R_{k}(B) = - \ln \rho(B)$$

如果 $0 < \rho(B) < 1$ , 则线性收敛; 且一般来说, $\rho(B)$ 越小, 收敛速度越快.

一般情况下 Gauss-Seidel 比 Jacobi 的 $\rho(B)$ 更小, 收敛速度更快; 对于 SOR 迭代法, 当 $\omega = 1$ 时退化为 Gauss-Seidel 迭代, 为使迭代速度最快, 可以选取最佳松弛因子

$$\omega_{\text{opt}} = \frac{2}{1 + \sqrt{1 - [\rho(I - D^{-1} A)]^2}}$$

其它的迭代求解方法

• 代数多重网格法 (Algebraic Multigrid, AMG): 先用简单迭代 (平滑器) 快速消掉解中的高频误差, 再把问题投影到粗网格上去处理低频误差, 通过多层网格循环, 使所有误差都被快速消除, 从而高效求解大型稀疏线性方程组.
• Krylov 子空间法: 猜测一个初始解 $\vec{x}_0$ , 从初始残差 $\vec{r}$0 = \vec{b} - A \vec{x}{0} 出发，构造由 $\vec{r}_0, A \vec{r}_0, A^2 \vec{r}_0 \cdots$ 生成的 Krylov 子空间, 在这个子空间里逐步寻找最优解, 通过迭代不断减小残差, 直至收敛.