Skip to main content

· 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;
}

· 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)

· 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}