Skip to main content

never say never? ——简介现代类型系统

· 17 min read

注意:本文之探讨现代类型系统中的一些实践,供大家参考。多有不严谨之处,望学术大佬海涵。

类型与类型系统

我们先介绍一下什么是类型。

这里我们简单的把类型(type)定义为值(value)的集合。

例如有两个取值的布尔类型:

bool := {false, true}

一些类型的取值范围是有限的(例如大部分的基本类型),一些是无限的(如字符串)

字面类型

在不引起混淆的情况下,值本身也可以作为一种类型的记号,只包含这个值本身,即

false := {false}
true := {true}

这样的类型成为字面类型(Literal Type)。

联合类型

接着我们引入联合类型(Union Type)或和类型(Sum Type)的概念:

AB 的联合类型记作 A | B,该类型的值可以取 A 的值,也可以取 B 的值,即在 ABA\cup B 中取值。

那么我们可以把布尔类型定义为:

bool := false | true

例如在 Python 中

from typing import Literal

Mode = Literal['r', 'rb', 'w', 'wb']
TrueType = Literal[True]
def open(file: str, mode: Mode) -> TrueType:
return True

参数 mode 取某几个固定的字面量取值,而返回值一定是 True,这都可以用字面类型很好的表示。用上面的写法,即

Mode := 'r' | 'rb' | 'w' | 'wb'

交类型

既然有 Union Type 就有 Intersect Type。

AB 的交类型(Intersect Type)记作 A & B,该类型的值必须在 ABA\cap B 中取。例如

ReadMode := 'r' | 'rb'
BinaryMode := 'rb' | 'wb'

那么这两个类型的交类型就是 'rb'

积类型

既然有 Sum Type 就有 Product Type。所谓 Product,指的就是笛卡尔积。

AB 的交类型(Product Type)记作 A x B,该类型的值必须在 A×BA\times B 中取。例如

A = 1 | 2
B = a | b | c

那么

A x B = (1, a) | (1, b) | (1, c) | (2, a) | (2, b) | (2, c) 

可以看出,这其实就是元组(Tuple),因此积类型也可以记作 (A, B)

二元组还可以被扩展为 n 元组,这里不再赘述。

一元组的取值与原类型相同,故 (A) = A

零元组非常特殊,只有 () 一个取值,被称为单位类型 Unit

子类型

如果 BAB \subset A,则称 BA 的子类型(Subtyping)。

子类型是一个偏序关系。即满足下面三个性质:

  • 自反性:AAA \subset A
  • 反对称性:ABBA    A=BA\subset B \land B \subset A \implies A = B
  • 传递性:ABBC    ACA\subset B \land B\subset C \implies A\subset C

从数学本质上来说,子类型似乎只是子集在类型系统里的新名字。但由于面向对象编程(OOP)的发扬光大,子类型已经成为了非常重要的一个概念。

面向对象

若一个类型继承自另外一个类型,类图如下所示

根据 OOP 的规则,所有 B 的对象都是一个(is-a)A 的对象。即 bB,bA\forall b \in B, b \in A,也就是说 BAB \subset A。因此我们可以说 BA 的子类型。

一般来说,子类型的对象可以隐式转换为父类型,毕竟子类型的取值范围更小,这样的转换天经地义。这个概念可以推广到一般的类型系统中,如我们上面提到的 bool 类型,类图如下所示

因此,当一个类的所有子类确定,且本身不会有额外的对象(在面向对象里,就是一个抽象类或是接口),那么它的派生关系实际上可以写作所有子类型的联合类型。

显然,根据定义,任意两个类型 A B 都有公共父类型 A | B,同时有公共子类型 A & B,类图如下所示

A & B 常用于表示同时实现两个接口类型。根据定义,它能作为其中任意一个类型的对象被使用。

如果我们要尝试把所有类型的图画出来,你可能会本能的追问两个问题:这个图的顶端是谁?这个图的底端是谁?

其实很多 OOP 语言都实现成了“单根模型”,即所有的类都有一个共同的父类。如 java.lang.Object kotlin.Any System.Object 等等。基于这样的模型,类图的顶端自然就是这个单根。在之后的叙述中,我们把这个类型记作 any,表示所有类型的联合、所有类型的父类型。

null

现在我们来探究一个问题,空指针 null 是什么类型?

你可能会下意识的认为,既然所有类型都可以赋值为 null,那它理所应当是 any 类型的。这忽略了一个事实:只有子类型才可以直接转换为父类型。any 不是除了本身之外任何类型的父类型,所以 null 若是 any 则几乎不能赋值给任何类型。这显然是矛盾的。

稍加思索你会注意到,由于 null 可以赋值给任何类型,所以它实际上是任何类型的子类。即

实际上,还有一种非常容易思考的方法得出这个结论:既然 null 可以赋值给任何类型,那么所有类型的交集,自然而然有且仅有这一个 null 的取值。

“一切类型的子类型”,这一点可能相当的反直觉,但却符合逻辑与实践。许多带 null 语言的类型系统就是这样实现的。

根据类图,我们可以把把 any 称为 Top Type,null 称为 Bottom Type。用 null 作为 Bottom Type 并非是唯一的选择,我们在后面还会提到别的考虑。

函数类型

我们先来考虑一种特殊的类型——函数类型。

所谓函数类型,即用于描述接受某些参数 A,返回某些返回值 B 的类型,记作 A -> B

这里只有一个参数,有多个参数的情况,可以用各个参数类型的笛卡尔积来表示。同理,多个返回值的情况也可以轻松表示。即 (A1, A2, ..., An) -> (B1, B2, ..., Bn)

那么 void 怎么办?考虑到它什么也没有接受/返回,可以用 ()Unit 来表示。这个方案把多个参数和返回值用元组表示的规律相统一,成为了很多现代编程语言的实践。

函数型变

下面我们来讨论函数类型之间的类型关系,尤其是其形变(variance)。假设 BA 的子类型。

协变

考虑下面两个函数及其类型

f: () -> A
g: () -> B

由于 BA 的子类型,所以 g() 的结果是 B 也是 A。那么我们就能说,g 的类型同时也可以是 () -> A。既然所有 () -> B 都是 () -> A,那么前者就是后者的子类型,即

这样函数派生关系与类型派生关系方向相同的型变成为协变(covariance)

逆变

考虑下面两个函数和两个对象,及其类型

f: A -> ()
g: B -> ()
a: A
b: B

由于 BA 的子类型,所以不仅 f(a) 调用是合法的,而且 f(b) 调用也是合法的。那么我们就能说,f 的类型同时也可以是 B -> ()。既然所有 A -> () 都是 B -> (),那么前者就是后者的子类型,即

这样函数派生关系与类型派生关系方向相反的型变成为逆变(contravariance)

不变

考虑下面两个函数及其类型

f: A -> A
g: B -> B

显然对于参数而言,可以构成逆变,对于返回值而言,可以构成协变。但是合在一起,派生关系的方向就相互矛盾了。因此这样的函数类型称为不变(invariance)。而

f: A -> B
g: B -> A

则不同,虽然同时构成协变和逆变,但是两者的派生关系是相互平行的。故最终可以得出 A -> BB -> A 的子类型。

never

永不返回

当我们说,一个函数不返回任何东西,和一个函数不返回的时候,实际上含义是不一样的。前者会返回,只是没有返回值,后者则是不会回到调用函数的地方。一个典型的例子便是 C/C++ 的 exit 函数,声明如下所示:

[[noreturn]] void exit(int status);

在这里,[[noreturn]] 用来标注函数永远不会返回。你可能已经发觉,[[noreturn]] void 这两个在一起似乎有些多余。难道不会返回的函数还能声明为非 void 的情况吗?许多现代语言引入了 never 类型(一些语言称为 nothing),专门用于表示这种不会返回的函数类型。例如上面的函数就可以表示为

exit: int -> never

如果我们想要认为构造出一个返回 never 的函数,最简单的方法莫过于

  • 调用其它返回 never 的函数
  • 死循环
  • 抛出异常

抛出异常看似会返回到上一层栈帧,但实际上不会返回到函数的调用点,故同样是永不返回。

分支语句

考虑下面的语句:

a: A
b: B
value = cond ? a : b;

value 是什么类型?为了保证可以接受两个分支的值,显然应该是 A | B

假如把代码改为

value = cond ? a : throw Exception();

显然两个分支的类型分别是 Anever。由于第二个分支不可能返回,所以 value 的类型必然取一个分支的类型 A。根据前面的公式,即 A | never = A。这个恒等式实际上就说明 neverA\text{never} \subset A

稍加思索你会发现,对于任何类型,这个结果都满足。换句话,任何类型都是 never 的父类型。这也与实践相吻合:既然不会返回,那么实际上可以 never 当做任意类型的值来看待。这或许也说明,前面的 [[noreturn]] void 并列,并非是完全多余的。如果你希望它参与运算,或许换成其它类型也是合理的。

由于 never 实际上不可能有任何的取值,我们可以得出 never=\text{never} = \varnothing,这也符合上面的数学逻辑。

很显然,我们得出了一种新的 Bottom Type 的方案,那就是 never

现在还有个小问题没有解决:如果一个类型系统里面同时有 nevernull,怎么会出现两个 Bottom Type 呢?难道它们两个都是所有类型的子类型?

never?

空语义

要先回答这个问题,我们不妨来看看 null 究竟是什么意思。

null 最初用来表示空指针,引申为,一个不存在的对象。由于人们常常忘记检查它的存在,所以带来的问题比较多。现代的编程语言有两个实际上等价的解决方案:Option<T> 和 Nullability。

一种方法,Option<T> = T | (),即要么有值,要么是 ()。这样在使用的时候就必须先排除是 () 的情形。

另一种方法称为 Nullability,即默认普通的类型 T 不能包含 null,只有显式添加一个问号 T? 的时候,才可以是接受 null。显然 T? = T | null,是 T 的父类型,不能直接转换为 T,进而就需要额外的检查才能使用。

前一种方案往往用于没有 null 的情况。而我们如果要讨论 nevernull 并存,故讨论后一种方案更为合理。

结论

讲到这里其实答案已经呼之欲出,根据 T? = T | null

T==neverT = \varnothing = \text{never},那么 never?=null=null\text{never}? = \varnothing | \text{null} = \text{null}

也就是说,never? = null

根据定义,TT? 的子类型,所以最终我们得出了结论:never 才是真正的 Bottom Type。

我们也可以画出一个完整的,带 Nullability 的类图:

推荐 BGM:Never Enough

表连接算法概览

· 5 min read

一般连接

形式语言与 SQL

连接两个表

RθSR \bowtie_\theta S
select *
from R, S
where theta
select *
from R
join S on theta

连接多个表

θ(R1,R2,,Rn)\Large \bowtie_\theta \normalsize(R_1, R_2, \cdots, R_n)
select *
from R1, R2, ... Rn
where theta
R1θ2R2θnRnR_1 \bowtie_{\theta_2} R_2 \cdots \bowtie_{\theta_{n}} R_n
select *
from R1
join R2 on theta2
...
join Rn on thetan
θ=i=2nθi\theta=\bigwedge_{i=2}^n\theta_i

Nested Loop Join

对于一般情况,多重循环是最朴素的算法

for r in R.tuples():
for s in S.tuples():
if predicate(r, s):
yield join(r, s)

可以通过分块遍历降低 IO 的性能损失

for rb in R.blocks():
for sb in S.blocks():
for r in rb:
for s in sb:
if predicate(r, s):
yield join(r, s)

等值连接

形式语言与 SQL

一种常见的连接是等值连接

关于一个属性的等值连接

R(id,),S(id,)RR.id=S.idSR(id, \cdots), S(id, \cdots)\\ R \bowtie_{R.id=S.id} S
select *
from R, S
where R.id = S.id
select *
from R
join S on R.id = S.id
select *
from R
join S using (id)
select *
from R
natural join S

关于多个属性的等值连接

R(id1,id2,),S(id1,id2,)RR.id1=S.id1R.id2=S.id2SR(id1, id2, \cdots), S(id1, id2, \cdots)\\ R \bowtie_{R.id1=S.id1 \land R.id2=S.id2 \land \cdots} S
select *
from R, S
where R.id1 = S.id1 and R.id2 = S.id2 and ...
select *
from R
join S on R.id1 = S.id1 and R.id2 = S.id2 and ...
select *
from R
join S using (id1, id2, ...)
select *
from R
natural join S

涉及等值查询问题,呼之欲出的就是两种方法:

  • 排序方法
  • 哈希方法

这两种方法在关联容器(即有序和无序的 setdict)的构造中相当有用,在表连接中也是如此。

如果其中一个关系在连接键上有索引,则可以直接采取 Index Nested Loop Join,否则则需要进行下面两种额外的工作。

Sort-Merge Join

RRSS 按照等值连接所需键排序,随后分别遍历两个关系;当遇到左右等值时,将等值的区间取出求笛卡尔积。

import functools
import itertools


def join(p):
r, s = p
return *r, *s


def smj(R, S, cmp):
R.sort(key=functools.cmp_to_key(cmp))
S.sort(key=functools.cmp_to_key(cmp))
ri = iter(R)
si = iter(S)
r = next(ri, None)
s = next(si, None)
while True:
if r is None or s is None:
return
while (c := cmp(r, s)) != 0:
if c < 0:
if (r := next(ri, None)) is None:
return
else:
if (s := next(si, None)) is None:
return

rc = [r]
sc = [s]

while (r := next(ri, None)) is not None and cmp(rc[0], r) == 0:
rc.append(r)

while (s := next(si, None)) is not None and cmp(sc[0], s) == 0:
sc.append(s)

yield from map(join, itertools.product(rc, sc))

如果连接的键是唯一的,则可以更进一步省去笛卡尔积的部分:

def join(r, s):
return *r, *s


def smj(R, S, cmp):
ri = iter(R)
si = iter(S)
while True:
if (r := next(ri, None)) is None or (s := next(si, None)) is None:
return
while (c := cmp(r, s)) != 0:
if c < 0:
if (r := next(ri, None)) is None:
return
else:
if (s := next(si, None)) is None:
return

yield join(r, s)

如果不想编写笛卡尔积的代码,又需要考虑到键不唯一的情况,可以考虑使用迭代器标记恢复之前的过程,但需要谨慎处理其中的边界条件:

struct iterator {
int *p = nullptr, *q = nullptr;
operator bool() const {
return p < q;
}
int operator*() const {
return *p;
}
iterator& operator++() {
++p;
return *this;
}
};

void smj(iterator r, iterator s) {
iterator mark;
while (r and s) {
if (!mark) {
while (*r != *s) {
if (*r < *s) {
++r;
if (!r) return;
} else {
++s;
if (!s) return;
}
}
mark = s;
}
if (*r == *s) {
emit(*r, *s);
++s;
if (s) {
continue;
}
}
s = mark;
++r;
mark = {};
}
}

Hash Join

用等值连接所需的键的哈希值建立 RR 的哈希表,再用 SS 等值连接的键去查询这个哈希表,得到的相同的行进行连接。

def join(r, s):
return *r, *s


def hj(R, S, hsh, eq):
build = {}
for s in S:
build.setdefault(hsh(s), []).append(s)
for r in R:
for s in build.get(hsh(r), []):
if eq(r, s):
yield join(r, s)

光线—三角形相交算法

· 3 min read

背景

在图形学编程中,经常需要求玩家准心对准的物体。该问题可转化为求与玩家视线相交且离玩家最近的三角形面片所对应的物体。因此需要一种光线-三角形相交算法来求出视线是否与某三角形相交、若相交,距离又是多少。下面介绍一种简单的基于参数和矩阵求解的算法。

原理

已知三角形的三个顶点 v0,v1,v2\pmb v_0,\pmb v_1,\pmb v_2,可以给出三角形内点的参数方程:

h=v0+α(v1v0)+β(v2v0)\pmb h=\pmb v_0+\alpha(\pmb v_1- \pmb v_0)+\beta(\pmb v_2- \pmb v_0)

其中参数 α,β\alpha,\beta 满足下面的条件:

α0,β0,α+β1\alpha\ge0,\beta\ge0,\alpha+\beta\le1

e1=v1v0e2=v2v0\pmb e_1=\pmb v_1-\pmb v_0\\ \pmb e_2=\pmb v_2-\pmb v_0\\

则可简记为:

h=v0+αe1+βe2\pmb h=\pmb v_0+\alpha\pmb e_1+\beta\pmb e_2

玩家从 p\pmb p 点往 r\pmb r 方向看去,与三角形所在平面的交点为 h\pmb h 的参数方程:

h=p+tr\pmb h=\pmb p+t\pmb r

其中参数 tt 满足:

t0t\ge0

如图所示:

联立两个方程

v0+αe1+βe2=p+tr\pmb v_0+\alpha\pmb e_1+\beta\pmb e_2=\pmb p+t\pmb r

移项

αe1+βe2tr=pv0\alpha\pmb e_1+\beta\pmb e_2-t\pmb r=\pmb p-\pmb v_0

(e1,e2,r)(αβt)=pv0(\pmb e_1,\pmb e_2,-\pmb r) \left(\begin{array}{r} \alpha\\ \beta\\ t\\ \end{array}\right)=\pmb p-\pmb v_0

A=(e1,e2,r)x=(αβt)b=pv0\pmb A=(\pmb e_1,\pmb e_2,-\pmb r)\\ \pmb x=\left(\begin{array}{r} \alpha\\ \beta\\ t\\ \end{array}\right)\\ \pmb b=\pmb p-\pmb v_0

也就是说

Ax=b\pmb A\pmb x=\pmb b

根据克拉默法则

xj=detAj(b)detAx_j=\dfrac{\det \pmb A_j(\pmb b)}{\det \pmb A}

α=det(b,e2,r)det(e1,e2,r)β=det(e1,b,r)det(e1,e2,r)t=det(e1,e2,b)det(e1,e2,r)\begin{array}{lll} \alpha&=&\dfrac{\det(\pmb b,\pmb e_2,-\pmb r)}{\det(\pmb e_1,\pmb e_2,-\pmb r)}\\ \beta&=&\dfrac{\det(\pmb e_1,\pmb b,-\pmb r)}{\det(\pmb e_1,\pmb e_2,-\pmb r)}\\ t&=&\dfrac{\det(\pmb e_1,\pmb e_2,\pmb b)}{\det(\pmb e_1,\pmb e_2,-\pmb r)} \end{array}

如果计算得出的参数在其限制范围内,则可以认为玩家的视线与三角形相交。解出的参数 tt 即为三角形面片与玩家的距离。

实现

该实现使用了 GLM 库。这里构造的矩阵实际上是所需矩阵的转置,但由于 detAT=detA\det A^T=\det A,因此不影响结果。

函数返回 -1 表示光线与三角形不相交,返回非负数表示三角形面片与玩家的距离。

#include <glm/glm.hpp>

float raytrace(glm::vec3 v0, glm::vec3 v1, glm::vec3 v2, glm::vec3 p, glm::vec3 r) {
glm::vec3 e1 = v1 - v0;
glm::vec3 e2 = v2 - v0;
glm::vec3 b = p - v0;
float den = glm::determinant(glm::mat3(e1, e2, -r));
if (den == 0) return -1;
float nom_u = glm::determinant(glm::mat3(b, e2, -r));
float nom_v = glm::determinant(glm::mat3(e1, b, -r));
float nom_t = glm::determinant(glm::mat3(e1, e2, b));
float u = nom_u / den;
float v = nom_v / den;
float t = nom_t / den;
if (u >= 0.0 && v >= 0.0 && u + v <= 1.0 && t >= 0) {
return t;
}
return -1;
}

<experimental/simd> 初体验

· 3 min read

之前用非标准的 simd 写过矩阵乘法加速,近期注意到 C++ 新的 Technical specifications——Parallelism library extensions v2 加入了 <experimental/simd>。于是我尝试用它重写了一下。

#include <iostream>
#include <string>
#include <functional>
#include <chrono>
#include <vector>
#include <cstring>
using namespace std;


double timeit(std::function<void()> test) {
auto start = std::chrono::system_clock::now();
test();
auto stop = std::chrono::system_clock::now();
std::chrono::duration<double, std::milli> time = stop - start;
return time.count();
}

#include <experimental/simd>
namespace stdx = std::experimental;

struct vec {
constexpr static int N = 256;

stdx::fixed_size_simd<double, 4> a[N];
double& operator[](int x) {
return ((double*) a)[x];
}
const double& operator[](int x) const {
return ((double*) a)[x];
}
double dot(const vec& x) const {
stdx::fixed_size_simd<double, 4> sum = 0;
for (int i = 0; i < N; i++)
sum += a[i] * x.a[i];
return stdx::reduce(sum);
}
};

class matrix {
constexpr static size_t m = 1024, n = 1024;
vector<double> e;
public:
explicit matrix(): e(m * n) {}
void random() {
for (size_t i = 0; i < m; ++i)
for (size_t j = 0; j < n; ++j)
at(i, j) = rand();
}
matrix(matrix const& that) = default;
matrix(matrix&&) = default;
matrix& operator=(matrix const& that) = default;
matrix& operator=(matrix&& that) = default;
double& at(size_t i, size_t j) { return e[i * m + j]; }
double const& at(size_t i, size_t j) const { return e[i * m + j]; }

matrix mul_plain(matrix const& that) const {
size_t p = that.m;
matrix product;
for (size_t i = 0; i < m; ++i)
for (size_t j = 0; j < n; ++j) {
for (size_t k = 0; k < p; ++k)
product.at(i, j) += at(i, k) * that.at(k, j);
}
return product;
}

matrix mul_simd(matrix const& that) const {
vector<vec> lines(m), columns(n);
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
lines[i][j] = at(i, j);
columns[j][i] = that.at(i, j);
}
}
matrix r;
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
r.at(i, j) = lines[i].dot(columns[j]);
}
}
return r;
}
};

double timeit(matrix (matrix::* test)(matrix const&) const, int times, const char* impl, double base = 0) {
double time = timeit([=] {
srand(1);
for (int i = 0; i < times; ++i) {
matrix m1, m2;
m1.random(); m2.random();
(m1.*test)(m2);
}
});
printf("'%5s' took %.2f ms to complete %d times 1024 x 1024 matrices multiplication", impl, time, times);
if (base && base != time) {
printf(", %.2f%% ", abs(base - time) / base * 100);
if (base > time) printf("faster");
else printf("slower");
}
puts("");
return time;
}


int main() {
double base = timeit(&matrix::mul_plain, 1, "plain");
timeit(&matrix::mul_simd, 1, "simd", base);
}

输出:

'plain' took 5458.48 ms to complete 1 times 1024 x 1024 matrices multiplication
' simd' took 232.27 ms to complete 1 times 1024 x 1024 matrices multiplication, 95.74% faster

速通统计学

· 2 min read

本文为 UESTC 《概率论与数理统计》课程统计学部分(6-10章)期末考试复习速通笔记,旨在背多分。

四个统计学分布

正态分布

XN(μ,σ2)E(X)=μ,D(X)=σ2X\sim N(\mu, \sigma^2)\\ E(X)=\mu,D(X)=\sigma^2

卡方分布

χ2=i=1nXi2χ2(n) 其中 XN(0,1)E(χ2)=n,D(χ2)=2n\chi^2=\sum_{i=1}^n X_i^2\sim\chi^2(n)\\ \text{ 其中 } X \sim N(0, 1)\\ E(\chi^2)=n,D(\chi^2)=2n

t 分布

又叫学生分布

T=XY/nt(n) 其中 XN(0,1),Yχ2(n)T=\dfrac{X}{\sqrt{Y/n}} \sim t(n)\\ \text{ 其中 } X \sim N(0, 1), Y\sim\chi^2(n)

F 分布

又叫费希尔(Fisher)分布

F=X/n1Y/n2F(n1,n2) 其中 Xχ2(n1),Yχ2(n2)F=\dfrac{X/n_1}{Y/n_2}\sim F(n_1,n_2)\\ \text{ 其中 } X\sim\chi^2(n_1), Y\sim\chi^2(n_2)

参数估计与假设检验

距估计

μ^=Xσ^2=n1nS2\hat{\mu}=\overline{X}\\ \hat{\sigma}^2=\dfrac{n-1}nS^2

极大似然估计

核心思想是使似然函数 LL 取极大值。即求出一个 θ\theta,使取出这个样本的概率最大。

L=i=1nP(Xi=xi)lnL=i=1nlnP(Xi=xi)lnLθ=0L=\prod_{i=1}^nP(X_i=x_i)\\ \ln L=\sum_{i=1}^n \ln P(X_i=x_i)\\ \dfrac{\partial \ln L}{\partial\theta}=0

区间估计与假设检验

  • 估计或检验 μ\mu,当 σ\sigma 已知时
U=Xμσ/nN(0,1)U=\dfrac{\overline{X}-\mu}{\sigma/\sqrt{n}}\sim N(0,1)
  • 估计或检验 μ\mu,当 σ\sigma 未知时
T=XμS/nt(n1)T=\dfrac{\overline{X}-\mu}{S/\sqrt{n}}\sim t(n-1)
  • 估计或检验 σ\sigma,当 μ\mu 未知时
χ2=n1σ2S2χ2(n1)\chi^2=\dfrac{n-1}{\sigma^2}S^2\sim\chi^2(n-1)

第一类错误:弃真(犯错概率为 α\alpha

第二类错误:纳伪

线性回归

lxy=i=1nxiyinxˉyˉlxx=i=1nxi2nxˉ2lyy=i=1nyi2nyˉ2b^=lxylxx,a^=yˉb^xˉσ^2=lyyb^2lxxn2R=lxylxxlyy>Rα(n2)l_{xy}=\sum_{i=1}^nx_iy_i-n\bar{x}\bar{y}\\ l_{xx}=\sum_{i=1}^nx_i^2-n\bar{x}^2\\ l_{yy}=\sum_{i=1}^ny_i^2-n\bar{y}^2\\ \hat{b}=\dfrac{l_{xy}}{l_{xx}}, \hat{a}=\bar{y}-\hat{b}\bar{x}\\ \hat{\sigma}^2=\dfrac{l_{yy}-\hat{b}^2l_{xx}}{n-2}\\ R=\dfrac{l_{xy}}{\sqrt{l_{xx}l_{yy}}}>R_\alpha(n-2)

SSO? SOS!

· 8 min read

这是我在编写编译器的过程中遇到的一个让人拍案叫绝的陷阱。特此记录下来。

问题

假设你正在编写源代码处理系统。你设计了这样一个函数和类:

std::vector<std::string_view> splitLines(std::string_view view) {
std::vector<std::string_view> lines;
const char *p = view.begin(), *q = p, *r = view.end();
while (q != r) {
if (*q++ == '\n') {
lines.emplace_back(p, q - 1);
p = q;
}
}
if (p != q) lines.emplace_back(p, q);
return lines;
}

struct Source {
std::vector<std::string> snippets;
std::vector<std::string_view> lines;

void append(std::string code) {
snippets.push_back(std::move(code));
for (auto line : splitLines(snippets.back())) {
lines.push_back(line);
}
}
};

为了效率,你没有真的去分割字符串——而是利用 string_view,引用字符串中的一段作为一行,并把它们存入 lines。但是 string_view 并没有字符串的所有权,为了保证它引用的字符串不被回收,再使用 snippets 来保存这些字符串。如下图所示:

下面我们来试一试:

int main() {
Source source;
source.append("ab\ncd\nef");
for (auto line : source.lines) {
std::cout << line << std::endl;
}
}

输出:

ab
cd
ef

很好!

可是,当我们试图多次调用 append 的时候

int main() {
Source source;
source.append("line 1");
source.append("line 2");
source.append("line 3");
for (auto line : source.lines) {
std::cout << line << std::endl;
}
}

大跌眼镜的事情发生了:输出了一团乱码!

思考

输出乱码原因有哪些呢?你肯定已经想到了以下两种可能:

  • string 管理的 char[] 被覆写或被释放了
  • string_view 被修改了

先来测一下第一种可能

int main() {
Source source;
source.append("line 1");
source.append("line 2");
source.append("line 3");
for (auto const& snippet : source.snippets) {
std::cout << snippet << std::endl;
}
}

输出:

line 1
line 2
line 3

看上去字符串完好无损。

再来测测第二种可能

int main() {
Source source;
source.append("line 1");
for (auto line : source.lines) {
std::cout << (void*) line.data() << ' ';
}
std::cout << std::endl;
source.append("line 2");
for (auto line : source.lines) {
std::cout << (void*) line.data() << ' ';
}
std::cout << std::endl;
source.append("line 3");
for (auto line : source.lines) {
std::cout << (void*) line.data() << ' ';
}
}

输出:

0x26a40fc18e0
0x26a40fc18e0 0x26a40fc1c90
0x26a40fc18e0 0x26a40fc1c90 0x26a40fc1d00

string_view 也没有变化!

规律

你可能百思不得其解。先别急,我们再来做一个实验。

之前我们就发现,如果只 append 一次是没问题的。那么不妨测试一下乱码与 append 次数有没有什么关系。

int main() {
Source source;
for (int i = 0; i < 10; ++i) {
source.append("OOO");
for (auto line: source.lines) {
if (line == "OOO") {
std::cout << line << ' ';
} else {
std::cout << "***" << ' ';
}
}
std::cout << std::endl;
}
}

输出:

OOO
*** OOO
*** *** OOO
*** *** OOO OOO
*** *** *** *** OOO
*** *** *** *** OOO OOO
*** *** *** *** OOO OOO OOO
*** *** *** *** OOO OOO OOO OOO
*** *** *** *** *** *** *** *** OOO
*** *** *** *** *** *** *** *** OOO OOO

这太有规律了!每次 append 乱码的 string_view 数量要么保持不变,要么就翻一倍,这似乎暗合了某个容器的某个特性。是的,vector。于是我们对前面的 Source 作出一点修改:

struct Source {

- std::vector<std::string> snippets;
+ std::deque<std::string> snippets;

};

再执行上面的代码,输出:

OOO
OOO OOO
OOO OOO OOO
OOO OOO OOO OOO
OOO OOO OOO OOO OOO
OOO OOO OOO OOO OOO OOO
OOO OOO OOO OOO OOO OOO OOO
OOO OOO OOO OOO OOO OOO OOO OOO
OOO OOO OOO OOO OOO OOO OOO OOO OOO
OOO OOO OOO OOO OOO OOO OOO OOO OOO OOO

修复了!太棒了!

但,这应该不是 vector 的问题吧?

第三种可能

vectorsize > capacity 的时候,会重新分配两倍的空间,并把原来的元素全部移动过去。但是移动的是 string,它管理的 char[] 难道也跟着发生了移动吗?还是说,发生的不是移动,而是复制?

带着这个疑问,我们再做一个实验:

struct Observer {
int i;
Observer(int i): i(i) {
std::cout << " make " << i;
}
Observer(Observer&& other) noexcept: i(other.i) {
std::cout << " move " << i;
}
Observer(const Observer& other): i(other.i) {
std::cout << " copy " << i;
}
};

int main() {
std::vector<Observer> observers;
std::cout << "line 1: ";
observers.emplace_back(1);
std::cout << std::endl << "line 2: ";
observers.emplace_back(2);
std::cout << std::endl << "line 3: ";
observers.emplace_back(3);
std::cout << std::endl << "line 4: ";
observers.emplace_back(4);
}

输出:

line 1:  make 1
line 2: make 2 move 1
line 3: make 3 move 1 move 2
line 4: make 4

需要注意的是,Observer 移动构造器的 noexcept 是不可以省略的,否则它会优先选择复制构造器而不是可能抛出异常的移动构造器。

经过查询,我们可以确认 string 的移动构造器是 noexcept 的。可见发生的是移动而不是复制。

也就是说:我们移动了 string,而它管理的 char[] 也跟着移动了!为什么?

在引出今天的主角之前,我们再测试一个例子。Source 回退到未修复的版本,运行下面的代码:

int main() {
Source source;
source.append("a very long line 1");
source.append("a very long line 2");
source.append("a very long line 3");
for (auto line : source.lines) {
std::cout << line << std::endl;
}
}

输出:

a very long line 1
a very long line 2
a very long line 3

问题奇迹般的消失了。

SSO

没错,这一切的罪魁祸首便是:短字符串优化(Short String Optimization,SSO)。

这是一种常见的字符串实现,被大部分的标准库(如 GCC libstdc++、Clang libc++、MSVC STL 等)采用。当字符串比较长的时候,字符串会被存在别处,并用 string 对象中的指针引用,如我们预料的一样;当字符串比较短的时候,字符串将直接存在 string 对象内部!

例如 libstdc++ 的实现如下图所示:

std::string s1 = "line 1";
std::string s2 = "a very long line 1";

在这个实现中,短字符串优化的阈值是 15。末尾多余的一个字符用来存 0,来保证它同时是一个 C-Style String。此时 capacity() 永远返回 15。

综上所述,当字符串较短的时候,由它构造的 string_view 实际上是指向 string 自己(的一部分)的!那么当 string 本身被移走,它所在的位置被释放,string_view 自然而然就引用了一段非法的内存。deque 在插入新元素的时候不会移动旧元素,因此避免了这个问题。

后记

我实际上很早就了解过 SSO,但是实在没料到会以这样的形式遭遇,过程曲折离奇,叹为观止。只能说如果没有这个前置知识我可能还得埋头苦苦寻找原因很长时间。

高中圆锥曲线结论

· One min read

若椭圆(b2>0b^2>0)或双曲线(b2<0b^2<0

x2a2+y2b2=1\frac{x^2}{a^2}+\frac{y^2}{b^2}=1

与直线

y=kx+my=kx+m

交于A(x1,y1),B(x2,y2)A(x_1,y_1),B(x_2,y_2),联立消 yy

(b2+a2k2)x2+2kma2x+a2(m2b2)=0(b^2+a^2k^2)x^2+2kma^2x+a^2(m^2-b^2)=0

Δ=4a2b2(b2+a2k2m2)0\Delta = 4a^2b^2(b^2+a^2k^2 - m^2) \geq 0 x1+x2=2kma2b2+a2k2,x1x2=a2(m2b2)b2+a2k2x_1+x_2=\frac{-2kma^2}{b^2+a^2k^2}, x_1x_2=\frac{a^2(m^2-b^2)}{b^2+a^2k^2} y1+y2=2mb2b2+a2k2,y1y2=b2(m2a2k2)b2+a2k2y_1+y_2=\frac{2mb^2}{b^2+a^2k^2}, y_1y_2=\frac{b^2(m^2-a^2k^2)}{b^2+a^2k^2} x1y2+x2y1=2ka2b2b2+a2k2x_1y_2+x_2y_1=\frac{-2ka^2b^2}{b^2+a^2k^2} x1y2x2y1=mΔb2+a2k2|x_1y_2-x_2y_1|=|m|\frac{\sqrt{\Delta}}{b^2+a^2k^2} AB=1+k2Δb2+a2k2|AB|=\sqrt{1+k^2}\frac{\sqrt{\Delta}}{b^2+a^2k^2}