2018年7月 – 第2页

网络流讲稿

【省选难度】

[pdf-embedder url=”https://mhy12345.xyz/wp-content/uploads/securepdfs/2018/07/网络流.pdf” title=”网络流 v1.0″]

数论部分讲稿

[pdf-embedder url=”https://mhy12345.xyz/wp-content/uploads/securepdfs/2018/07/数论-讲义.pdf” title=”数论讲义v1.0″]

[pdf-embedder url=”https://mhy12345.xyz/wp-content/uploads/securepdfs/2018/07/数论-课件.pdf” title=”数论课件 v1.0″]

bool Miller_Rabbin(int n,int a)//a属于[2,n）
{
        if (n<2)return false;
        if (!(n%a))return false;
        int r=0,d=n-1;
        while (!(d&1))
        {
                d>>=1;
                r++;
        }
        int k=pow_mod(a,d,n);
        if (k==1)return true;
        for (int i=0;i<r;i++)
        {
                if (k==n-1)return true;
                k=k*k%n;
        }
        return false;
}

typedef long long qword;
qword ext_gcd(qword a,qword b,qword &x,qword &y)
{
        if (a%b==0)
        {
                //a*x+b*y==b
                x=0;y=1;
                return y;
        }
        qword ret=ext_gcd(b,a%b,x,y);
        qword tx=x,ty=y;
        x=ty;
        y=tx-a/b*ty;
        return ret;
}

#include<iostream>
using namespace std;
#define MAXN 1010
bool pflag[MAXN];
int prime[MAXN],topp=-1;
void init() {
	for (int i=2;i<MAXN;i++)
	{
		if (!pflag[i])
			prime[++topp] = i;
		for (int j=0;j<=topp && prime[j]*i<MAXN;j++)
		{
			pflag[i*prime[j]] = true;
			if (i%prime[j] == 0)break;
		}
	}
	for (int i=0;i<=topp;i++)
		cout<<prime[i]<<endl;
}

int main()
{
	init();
}

#include<iostream>
using namespace std;
#define MAXN 100
#define MOD 1000000007
bool pflag[MAXN];
int prime[MAXN],topp=-1;
int inv[MAXN];
int phi[MAXN];
int h[MAXN];//存最小的那个质数所在的括号
int hh[MAXN];//存其他括号的乘积
int g[MAXN];//存最小那个质数所在的括号有多少个数
int gg[MAXN];//存其他括号的个数的乘积

long long pow_mod(long long a,long long b) {
	long long res = 1;
	while (b) {
		if (b&1)
			res = res * a%MOD;
		a = a*a%MOD;
		b>>=1;
	}
	return res;
}
//(1+5+25)*(1+3)   *   (1+2)
//(1+5+25)*(1+3)   *   (1+2+4)
//    hh[i]              h[i]
//
//(1+2)*2+1 = (1+2+4)
//(1+2+4)*2+1 = (1+2+4+8)
//
//(1+5+25)         *   (1+3)
//(1+5+25)*(1+3)   *   (1+2)

void init() {
	inv[1] = 1;
	for (int i=2;i<MAXN;i++) {
		if (!pflag[i]) {
			prime[++topp] = i;
			inv[i] = pow_mod(i,MOD-2);
			phi[i] = i-1;
			h[i] = i+1;
			hh[i] = 1;
			g[i] = 2;
			gg[i] = 1;
		}
		for (int j=0;j<=topp && i*prime[j]<MAXN;j++) {
			pflag[i*prime[j]] = true;
			inv[i*prime[j]] = (long long)inv[i] * inv[prime[j]] %MOD;
			if (i%prime[j] == 0) {
				phi[i*prime[j]] = phi[i] * prime[j];
				h[i*prime[j]] = h[i] * prime[j] + 1;
				hh[i*prime[j]] = hh[i];
				g[i*prime[j]] = g[i]+1;
				gg[i*prime[j]] = gg[i];
				break;
			}else {
				phi[i*prime[j]] = phi[i] * (prime[j]-1);
				h[i*prime[j]] = prime[j] + 1;
				hh[i*prime[j]] = h[i] * hh[i];
				g[i*prime[j]] = 2;
				gg[i*prime[j]] = g[i] * gg[i];
			}
		}
	}
}

int main() {
	init();
	/*
	for (int i=0;i<=topp;i++)
		cout<<prime[i]<<"\t";*/
	for (int i=2;i<MAXN;i++) {
		cout<<i<<" INV="<<inv[i]<<" PHI="<<phi[i]<<" H="<<h[i]*hh[i]<<" G="<<g[i]*gg[i]<<endl;
	}
}

网站访问速度是搜索引擎衡量一个网站的核心因素，因为这在另一方面意味着用户点选这个词条的用户体验。因此，我们十分有必要优化Wordpress的速度。这里，我使用了WP Super Cache和Jetpack插件，WP Super Cache可以缓存你的博客页面，特别是对于匿名用户，不论如何，网站对匿名用户返回的始终是同一个页面，因此将他们缓存下来可以达到很好的效果（从原来的2s变成了几乎看不到延迟）。同时Jetpack插件提供了免费的CDN服务。也可以优化图片以及视频的加载速度。

我是如何和搜索引擎斗智斗勇的.要有一个固定的域名

要有一个固定的域名

之前曾经提到百度将我的博客列为“用户不友好网站”，追溯其原因。回想之前对网站做的危险操作，不外乎变更子域名了。最初我的博客使用的是blog.mhy12345.xyz，提交了百度收录之后，我又心血来潮将站点的域名改变成了mhy12345.xyz。当时的想仅仅是，也许我改变了域名，搜索引擎自己会直接迁移过去blablabla，但是可能并不是这样。

从搜索引擎视角来看，首先出现了一个blog.mhy12345.xyz网站，接着有出现了一个冒充blog.mhy12345.xyz网站的另一个网站，内容完全相似。直接判定为抄袭站点。接下来发生的事情就更好解释了，百度只收录blog.mhy12345.xyz，我一怒之下，删去了对于blog.mhy12345.xyz的解析，剩下了一个mhy12345.xyz，而这个网站已经被搜索引擎判定为了恶意站点。然后……就没有然后了……（等等，为啥google就没有判断出错呢？唔，毕竟百度不是世界一流互联网公司吧）

我的网站

自己写博客，特别是技术类博客，当然希望能够有其他人能够看到。我的第一篇正式文章Ubuntu 16.04 & Docker 搭建 WordPress 写于2017年4月，而到现在已经一年多了。大大小小杂文有七十于篇（怎么会那么多！？）自己认为很有用的也有大约10篇了。

有了足够的干货，当然希望更多人能够看到，这就依赖于搜索引擎了。然而搜索引擎却貌似并不买我的账。怎么说呢，google是个好公司啊，因为至少在我一点一点写博客这个时间段，google索引量在逐渐增加的，google给我带来了95%的搜索引擎索引量。而百度，甚至将我的博客列为了“用户不友好网站”==。想一想其实不无道理，google面向国外，有非常多的优质个人站点，所以理所当然赋予了个人站点更高的权值，反观百度的用户主体，国内，随随便便点开一个小站，各种广告网站满天飞，是我也想一枪毙掉所有的个人网站了。

真实感图像渲染系列：小结

这是真实感图像渲染系列的第十篇文章。

这里写一写个人在编写过程中发现的一些坑，分享出来以免之后再踩吧。

CMake好啊，比Makefile好，可以很容易整合第三方库，但是……编译超级慢，别人的工程3s编译完成我的需要半分钟……
eps设多少确实是个坑，我设的是1e-7，但是得保证牛顿迭代的精度严格高于eps，否则焦点可能存在大于eps的误差，eps的设置就没有意义了
MAC的错误报告挺有用的 == 比如下面这个直接秒调

真实感图像渲染系列：相机(Camera)与景深效果(Field of Depth)

这是真实感图像渲染系列的第九篇文章。

普通相机模型

不同的相机模型中，有一个人眼位置和一个屏幕，对于屏幕中每一个像素，连接人眼和该像素位置得到一条射出光线。其2D示意图如下：

首先，一直屏幕像素坐标，以及屏幕本身在场景中的位置，我们可以计算出该像素的位置，如图中红点所述。接着，连接人眼和像素的光线就可以用于计算该光线的颜色。这一部分挺容易理解的。

带景深的相机模型

在普通相机模型的屏幕和人眼位置之外，带景深的相机模型额外存在一个焦平面。如下图所示

通过普通的相机模型，我们可以求得普通相机的光线的方向。

普通相机的光线一定会和焦平面在某个位置碰撞，为了达到景深效果，一方面，我们需要保证同一个像素对应的光线一定会击打到焦平面上的同一个点，另一方面，在焦平面之外的其他位置，同一个像素对应的光线应该击打到不同的点。

解决方案如下：我们求得在普通相机模型中，光线和焦平面的焦点，将该焦点景深相机模型的焦点。接下来，我们对光线的起始位置做一个微小扰动，并且修改光线方向使得光线恒定过焦平面上的点。这样，无论初始光线如何扰动，焦平面上的点都可以最清晰的表现出来。而不在焦平面上的点，会产生随机的模糊效果。

真实感图像渲染系列：贝塞尔曲线旋转体与直线求交

这是真实感图像渲染系列的第八篇文章。

概述

贝塞尔曲线是一个由多个控制点组合而成的曲线。在本文探讨的贝塞尔曲线，是由多条简单的4阶贝塞尔曲线拼接而成。绕z轴旋转贝塞尔曲线，我们可以得到一个旋转体曲面，该曲面就是待渲染的参数曲面。本文将探讨如何使用牛顿迭代法与直线求交。求交过程中

关于迭代初值的求法
迭代中的一些优化技巧
关于浮点数精度误差eps的求法。

贝塞尔曲线的定义

贝塞尔曲线是一个由多个控制点组合而成的曲线，n阶贝塞尔曲线定义为

(1) $\begin{equation*} P(u) = \left[ \begin{matrix}P^{(x)}(u) \\ P^{(y)}(u) \end{matrix} \right] = \left[ \begin{matrix} \sum _{i=0}^n P_i^{(x)} C(n,i) (1-u)^{n-i}u^i \\ \sum _{i=0}^n P_i^{(y)}C(n,i) (1-u)^{n-i}u^i \end{matrix} \right] \end{equation*}$

式(1)定义了该曲线 $P(u)$ 的关于参数 $u \in [0,1]$ 的参数方程，其中 $P_i$ 为n个曲线的控制点。当我看到这个公式的时候，整个人都已经懵掉了。实际上理解起来并不复杂。

函数 $P(u)$ 由参数 $u$ 定义了一个二维向量，其中 $P^{(x)}$ 与 $P^{(y)}$ 相互独立。只需要理解其中一个的产生式即可。考虑函数 $f_i(u) =C(n,i) u^i (1-u)^{n-i-1}$ ，发现这 $f_0(u),f_1(i),\cdots,f_{n-1}(u)$ 实际上是一个和为1的数列，数列在杨辉三角形中定义。换句话来说， $f_i(u)$ 实际上可以得到一个权值的分配，使用 $f_i(u)$ 对于所有的控制点 $P_i$ 加权平均，就是贝塞尔曲线尝试完成的事情。当 $f_0(u) = C(n,0) u^0 (1-u)^{n-i} = 1$ 权值全部分配给了 $P_0$ ，因此贝塞尔曲线的起点( $u=0$ )和 $P_0$ 重合。

贝塞尔曲线旋转体

在三维空间中，我们通过旋转贝塞尔曲线得到一个曲面。曲面强制对z轴旋转对称。（注：在代码中是对x轴旋转对称，不过z轴对称相对好理解一些）。

(2) $\begin{equation*} S(u,\theta) = \left[ \begin{matrix} Q_x + sin\theta P^{(x)}(u) \\ Q_y + cos\theta P^{(x)}(u) \\ Q_z + P^{(y)}(u) \end{matrix} \right] \end{equation*}$

一个贝塞尔曲线旋转体要能够成功渲染，其核心就是需要支持与光线的求交。光线我们使用 $(rayO,rayD)$ 表示， $rayO$ 表示光线起点， $rayD$ 表示光线方向，是一个单位向量。则光线上一点可表示为 $C(t)$ ，满足：

(3) $\begin{equation*} C(t) = rayO + t \times rayD \end{equation*}$

倘若直线 $C(t)$ 与 $S(u,\theta)$ 相交，则有

(4) $\begin{equation*} F(t,u,\theta) = C(t) - S(u,\theta) = 0 \end{equation*}$

其中, $F(t,u,\theta)$ 为 $t,u,\theta$ 的函数，要求 $t>0$ 且 $0\leq u \leq 1$ ，原问题可转化为求函数零点问题。函数零点问题使用牛顿迭代解决。

(5) $\begin{equation*} x_{i+1} = x_i - [F'(x_i)]^{-1} \cdot F(x_i) \end{equation*}$

其中 $F'(x_i)$ 为函数 $F(x_i)$ 的Jacobian矩阵，定义为：

(6) $\begin{equation*} \frac{\partial F}{\partial t} = \left[ \begin{matrix} \frac{\partial F_1}{\partial t} & \frac{\partial F_1}{\partial u} & \frac{\partial F_1}{\partial \theta} \\ \frac{\partial F_2}{\partial t} & \frac{\partial F_2}{\partial u} & \frac{\partial F_2}{\partial \theta} \\ \frac{\partial F_3}{\partial t} & \frac{\partial F_3}{\partial u} & \frac{\partial F_3}{\partial \theta} \\ \end{matrix} \right] \end{equation*}$

要试图计算出式(\ref(eq:jacobian))，我们需要分别求出 $F(t,u,\theta)$ 对于 $t,u,\theta$ 偏导数的每一个分量的值（共九个）。结果为：

(7) $\begin{equation*} \frac{\partial F}{\partial args} = \left[ \begin{matrix} rayD_x & -sin\theta \frac{dP_x}{du} & -cos\theta P_x(u)\\ rayD_y & -cos\theta \frac{dP_x}{du} & +sin\theta P_x(u)\\ rayD_z & -\frac{dP_y}{du} & 0\\ \end{matrix} \right] \end{equation*}$

其中 $\frac{dP}{du}$ 是对于式子(1)按照最原始的求导法则求导的结果。

贝塞尔曲面物体设计

如图

共6条4阶贝塞尔曲线，每一条由绿蓝红绿四个控制点定义。

注意事项

牛顿迭代需要保证答案中 $u \in [0,1]$ ，而牛顿迭代本身无法附加条件，因此需要在答案的邻域中查找合适的 $u$ , $t$ , $\theta$ ，我用到的方法是，对于给定曲线，生成一个圆柱体包含框，随机该圆柱体里面高度 $h$ 的一个圆形面片，将光线与该面片的交点的 $u$ , $t$ , $\theta$ 作为迭代的初值。随机 $h$ 约30次即可得解。
由于牛顿迭代不太精确，可能是的结果产生微小偏移，这种偏移会对判断点在平面哪一侧产生影响，我们可以通过调整eps的取值，不同部分赋予不同的eps，使得这些误差不至于相互影响。
贝塞尔曲线求交本身费时间，但是判断是否一定没有交点却相对简单。可以求出贝塞尔曲线旋转体的“圆柱体包围盒”，（更精确的话，可以是空心圆柱体包围盒）。然后判断直线是否与包围盒有交点。判断直线与圆柱交点 $P = rayO + k \times rayD$ 的位置，可以先讲所有东西映射到 $z=0$ 屏幕，这是变成2维直线和圆求交，算出 $k$ ，回带到三维情况验证。
实际上，像本文采用的多条低阶贝塞尔曲线求交并不明智，倘若将本文那个碗的6条贝塞尔曲线改为1条简单一些的曲线，然后做一个玻璃杯什么的，渲染效率将提高6倍！