用数学浅谈浮点数

浮点数很神奇,Prof. Kahan是个鬼才。可以说浮点数的诞生很大程度上推动了计算机的发展。

浮点数的概况

​ 在计算机中,浮点数主要以科学计数法存储,其科学计数法底数为2,于是一般而言:浮点数分为三个部分,符号域(Sign),指数域(Exponent),尾数域(Fraction)。 分别对应科学计数法的有效数字的符号,有效数字的指数,有效数字的小数部分。

浮点数的指数域

​ 为了让浮点数的比较能够由整形的比较器完成,其指数采取移码来表示,(移码和补码只有符 号位相反,其余都一样。对于正数而言,原码、反码和补码都一样;对于负数而言,补码就是其绝对值的原码全部取反,然后加1(不包括符号位)) 形象地来来说移码是原码的移动,在数值上移码等于原码减去2^(k-1)-1,k是位数。但是并不是所有的浮点数都采取此种表示方法。有一小部分浮点数的指数域, 与之不同,我们在后面说。

浮点数的尾数域

​ 为什么叫尾数域,而不叫有效数字域,其实很容易想到二进制的科学计数法的整数部分必然是1,于是这个1被省略掉了,计算机不存储他。

浮点数中的特例

​ 如果我们真的像刚刚分析的那样,结束了浮点数的定义,会有一个问题,那就是0不见了,并且靠近0的一小部分有一个巨大的gap,我们后面在证明这个gap存在,于是我们的浮点数设计失败, 但是Prof. Kahan想到了一个解决方法,他发现当指数域的指是一个最小的负数的时候,当前浮点数所表示的数非常小,与0很接近,并且那一部分的精度非常高,精度高,但是靠近0的地 方有一个巨大的gap,怎么办?他考虑到,用一部分高精度的丢失换取gap的填充。也就是说,当指数部分是最小的负数的时候,我们不采取科学计数法了,当指数域为最小的负数 的时候,我们把尾数域看作一个无符号整形,他会可以表示一个整形区间,我们hash掉这个区间,让他来表示gap即可。后面来证明这个做法的可靠性。

​ 既然花了一部分指数域做了一个0,为什么不来一个无穷大呢?于是Prof. Kahan创建了无穷大的表示方法,当指数部分是最大的正数的时候,此数字表示无穷大怎么样?若符号域为正,我们表示 正无穷,若符号位为负,我们表示负无穷怎么样?于是我们发现我们思维的漏洞了,刚刚的0,岂不是有+0与-0?是的就是有,等下我们证明他的合理性。 再回到无穷上,我们与0的表示做对比,发现了指数固定的时候可是一个区间啊,有很多的无穷大,是啊,我们要不了这么多无穷大,那怎么办呢?于是Prof. Kahan说多余的他们都叫做 nan,意味着not a number。我们依然在后面证明这个做法的优点很大很大。

浮点数的具体表示

​ 这些都是理论指导,实际上浮点数到底怎么存的呢?见图 (照片丢了,我也没找到)

浮点数的相关证明:gap存在性

​ 给出特例,对于单精度浮点而言,我们先考虑浮点数的表示精度问题,如果放开指数不谈,令指数为0 ,那么我们可以表示出大致区间[1,2),如果指数为1,我们又可以表示区间[2,4) 等等。。。 我们仔细分析,这里的区间长度是变化的,但是表示这段区间的数的数量是固定的。显然一个问题出现了,我们可以大胆猜测,指数越大,精度越低,我们把几个区间都写一下 当指数为-2 -1 0 1 2 分别表示了区间[1/4,1/2)[1/2,1)[1,2)[2,4)[4,8),是的,区间长度递增,猜测正确。

​ 然后我们还发现,精度是有规律的,同一个区间的精度固定,因为尾数是平分区间的,不同区间的精度怎么样呢?我们来看看,[1,2)的精度是[2,4)的精度的两倍, 于是我们得到精度的详细情况,离0越远的区间,精度越低,且是他的前一个区间的精度的一半,(定义一个区间的前一个区间为与之相邻的离0更近的区间)

​ 我们再来考虑最小的正浮点数,在什么地方,对的,就是$2^{-127}$,我们假设这个数为x,他右边部分的精度达到了$\frac{2^{-126}-2^{-127}}{2^{23}}=2^{-150}$在x的左侧呢,哈哈一个gap,他与0之间相隔$2^{-127}$,$2^{-150}$ 和$2^{-127}$ 区别可大了,负数那边也是样的。

浮点数的相关证明:gap填充的可靠性

​ 如果考虑不要精度为$2^{-150}$的区间的,用它来填补gap,代价是什么?会不会导致浮点数优秀的精度递增模式被打破呢,很遗憾不会。我们来计算,如果考虑丢弃此区间 换来的最小的正科学计数法所表示的浮点数的值为$2^{-126}$ ,因为区间$[2^{-127} ,2^{-126})$ 拿走了,这时候他右边的精度为$\frac{2^{-125}-2^{-126}}{2^{23}}=2^{-149}$,emm,还行,我们来计算那个究极大gap的精度,$\frac{2^{-126}}{2^{23}}=2^{-149}$奇迹出现了,精度一摸一样,精度的优秀性质基本得到了保留,这种做法使得从0到正无穷的过程中没有变小。很神奇。

浮点数的相关证明:+0与-0的优点

​ 为什么要搞+0和-0,这两者不是相同的吗?是的他们是相同的,+0.0==-0.0 返回值是true ,这个时候我们开始思考0的意义,0到底是什么,我们参考无穷大,重新定义0,+0.0代表正无穷小 -0.0代表负无穷小,这才是他们本质上的意义,sorry,我们又把0给弄没了。这次我们不把它找回来了,+0.0和-0.0共同组成了0。为什么要这样做,还有一个额外的点,我们的数域里面可是有正负 无穷大的,我们,要搞在浮点数里面搞一套新的,特别的运算法则,极限的运算。这是他的另一个点。

浮点数的相关证明:nan的实用性

​ 为什么要搞nan,这是给程序员用的,哈哈哈哈哈哈嗝,哈哈嗝,用来debug,但笔者不太懂一点,为什么要搞那么多nan呢,指数最大值的时候,除了正负无穷,其他的数字都是nan,为什么要这么复杂呢? 限于水平,笔者大致猜到了,极有可能,不正确的运算直接会算出nan,这会加速浮点数的运算,(也就是说,不需要我们自己去做判断是否运算合法)(也就是说浮点数的运算不是封闭的,错误 或者不合法的运算会直接算出nan,而不是计算机去判断运算是否合法)这只是笔者的一个猜想。

浮点数细节:

​ 若x为nan ,那么x==x为假 x不管之后怎么运算,得到的永远是nan +0.0 == -0.0 为真

多项式倍增

让你用lg的时间复杂度求下面这东西(n<1e18)

$$
\begin{aligned}
&\prod_{i=1…n}{(2i-1)} %2^{64}\&
suppose\quad that\quad f(x,n)=(2x+1)(2x+3)(2x+5)…(2x+2n-1)%2^{64}\&
then \quad the \quad 0th \quad item \quad of \quad f(x,n) \quad is \quad answer\&
we \quad can \quad try \quad to\quad calculate \quad f(x,n) \quad by \quad f(x,\left \lfloor \frac{n}{2} \right \rfloor) \&
let \quad y=x+n\&
then \quad f(y,n)=(2y+1)(2y+3)(2y+5)…(2y+2n-1)%2^{64}\&
so \quad f(x+n,n)=(2x+2n+1)(2x+2n+3)(2x+2n+5)…(2y+4n-1)%2^{64}\&
Surprisedly \quad we \quad find \quad that \quad f(x,2n)=f(x,n)*f(x+n,n)\&
which \quad means \quad we \quad can \quad calculate\quad f(x,2n)\quad by\quad f(x,n) \quad easy\&
\&
becase \quad we \quad can \quad calculate \quad f(x+n,n)\quad by \quad f(x,n) \&
but \quad we \quad can’t \quad calculate \quad it \quad faster \quad because\quad of\quad the\quad huge\quad numbers \quad of \quad item\&
considering \quad that \quad we \quad need \quad mod \quad 2^{64} , we\quad can \quad reserve \quad the \quad items \quad with \quad index \quad less\quad than\quad 64\&
because \quad the \quad useful \quad information \quad is \quad the \quad 0th \quad item \&
and \quad in \quad f(x,n)\quad when \quad the\quad index\quad of\quad x \quad is \quad larger \quad than\quad 63 ,the \quad coefficient \quad of \quad it \quad must\quad divisible\quad 2^{64}\&
it \quad has \quad no \quad contribution \quad to \quad answer \quad and \quad f(x+n,n)\&
so \quad we \quad can \quad solve \quad it \quad by \quad this \quad way
\&
for \quad every \quad polynomial \quad we \quad reserve\quad the \quad first \quad 64\quad items \quad of\quad x
\& and \quad then \quad calculate \quad f(x,n) \quad by
\left{\begin{matrix}
\quad f(x,\left \lfloor \frac{n}{2} \right \rfloor)*f(x+\left \lfloor \frac{n}{2} \right \rfloor,\left \lfloor \frac{n}{2} \right \rfloor)\quad if \quad n%2=0\&
\quad f(x,\left \lfloor \frac{n}{2} \right \rfloor)f(x+\left \lfloor \frac{n}{2} \right \rfloor,\left \lfloor \frac{n}{2} \right \rfloor)(2x+2n-1) \quad if \quad n%2=1
\end{matrix}\right.
\&
and \quad we \quad can \quad get \quad the \quad answer \quad no \quad more \quad than \quad lg \quad times \quad recursion\&
because \quad of \quad the \quad 64 \quad items \quad only \quad and \quad we \quad don’t\quad care \quad the \quad higher \quad items ,\&
it \quad is \quad very \quad fast \quad to \quad get \quad f(x,n)\quad by \quad f(x+n,n)
\&
so\quad the\quad total\quad time\quad complexity \quad is \quad O(lgn)

\end{aligned}
$$

假设f(x,n)=(2x+1)(2x+3)(2x+5)…(2x+2n-1)%64

然后这东西的0次项系数就是答案

我们尝试通过f(x,n/2)来求f(x,n)

令y=x+n

则f(y,n)=(2y+1)(2y+3)(2y+5)…(2y+2n-1)%64

所以f(x+n,n)=(2x+2n+1)(2x+2n+3)(2x+2n+5)…(2x+4n-1)%64

我们惊讶地发现了f(x,2n)=f(x,n)*(fx+n,n)

这意味着我们可以通过f(x,n)来求f(x,2n)因为我们可以通过f(x,n)求出f(x+n,n)

很遗憾的是这些东西项数太多了

考虑到我们要的是模上2^64的答案,我们可以只保留前64项

因为有用的只有0次项,但是在f(x,n)转移到f(x+n,n)的时候也只有前64项有效,因为大于指数64的项,他们前面的系数一定整除2^64次方,

于是我们就有了做法了

我们保留前64项

时间复杂度为lg级别

贪心加暴力

大范围贪心,小范围暴力

这一类做法,通常不是很好想,但是很多题目都能这样做。

例题

有一个超级大的背包,物品的价值等于容量,但是物品只有8种容量分别为1-8 给你每个物品的数量ai,和背包总容量w,问能背走的最大价值是多少 w<1e18;

朴素背包

dp[i][j]为前i类物品,背满j是否可行,复杂度$O(8\times n)$

复杂度为什么这么大?

状态过多,状态冗余,可能存在某些能合并的状态

分解背包

​ 840是1-8的数的lcm, 于是我们可以把背包W分为 $840k + (W-840k)$ 两个部分。 且后一部分小于$840\times8$

为什么要这么分?

先谈谈唯一分解

我们对每一种放法 ,进行唯一分解:

先把每一类容量相等的物品唯一分解为两部分, 第一部分的总容量为840的一个倍数,第二部分总容量小于840

比方说某种方法选择了1000个1,3000个2,5000个3…

那么我们将其唯一分解为:

$1\times840+160个1$

$7\times 420+ 60个2$

$17\times280+240个3$
然后把每类容量为840的倍数的那一部分合在一起:
于是成了$(1\times840+7\times420\times2+17\times280\times3) + 160\times1+60\times2+240\times3+0\times4+0\times5+…$
=$840\times25 + 160\times1+60\times2+240\times3+0\times4+0\times5+…$ 这就是唯一分解。

根据唯一分解优化dp

然后我们考虑,为什么这一题不能使用,背包,因为容量大?物品多?是的,

但是这些都只是表面上的。我们深入思考,能不能把某些没有意义的方案合并到一个状态里面呢?

我们不要设dp[i][j]表示什么前i类物品装满容量j是否可行这样的状态。因为这是$8\times n$级别的状态
我们根据唯一分解,设dp[k][i][j]代表背包的第一个部分容量为840*k,第二部分为前i类物品装满容量j
若值为1,代表可行,否则不可行。

状态的个数

显然k<n/840,i<8, j<840*8

现在的状态数为(n/840)*8*(840*8)

状态数目变多了。转移也更加复杂,看似此状态还不如朴素做法。

单调性

这一点应该是很容易想到的。这个dp[k][i][j],具有单调性,一定存在某个值t, 使得dp[k][i][j]的值在i和j固定,k<=t的时候全为1,在k>t的时候全为0

显然这个t是最优的

优化

我们换一种状态的设法,设dp[i][j]为只取前i类物品的方案的唯一分解下,不考虑背包容量上限,第二部分容量为j,第一部分的k能取到的大值。

转移方程

dp[i][j] —> dp[i+1][j+t*i] t是选取的数量,j+t*i<8*840

这样的做法就已经很快了。

数位dp

常见数位dp

如果一个在数字上的计数问题只与数位和数的大小有关的时候 我们可以尝试用数位dp来解决。最经典的就像不要62那道题。

数位dp状态

我们常常设dp[val][len][limit][lead]来表示以val开头 数位长度剩余len(包含val),limit表示数有没有上限,后面发现这一维度没有作用 。lead表示val及以前是否含有前导数,来特判某些跟前导0有关的题目

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
ll dfs(int*num, int n, int pre, int pos, bool limit, bool lead){//try to fill the pos bit
if(pos==n)return 1;
else{
if(!limit&&dp[pre==6][num[pos]][n-pos][lead]!=-1)
return dp[pre==6][num[pos]][n-pos][lead];
else {
int upper=limit ? num[pos] : 9;
ll ans=0;
for(int i=0;i<=upper;i++){
if(i==4)continue;
if(pre==6&&i==2)continue;
ans+=dfs(num,n,i,pos+1,limit&&i==upper,lead||i!=0);
}
if(!limit) dp[pre==6][num[pos]][n-pos][lead]=ans;
return ans;
}
}
}

莫队算法

例题1

http://codeforces.com/contest/617/problem/E

给长度为n的序列,给数字k,给q次询问,每次询问给出[L,R],求这个区间内有多少个连续区间的异或和等于k。

题解

预处理前缀异或和,化简区间为两个点异或值为k,再来莫队。

例题2

hdu5213,给你序列a,其长度为n,给你q个询问,每次询问给两个不相交的区间,求a[i]+a[j]=k的方案数,i属于第一个区间,j属于第二个区间。

(1≤N≤30000) (2≤K≤2*N) (1≤ai≤N)(1≤M≤30000) (1≤Li≤Ri<Ui≤Vi≤N)

题解

化简多区间为单区间询问,再来莫队

例题3

fzu2226, 给你长度为n的序列,定义第i个元素和第j个元素间的距离dis(i,j)=abs(i-j) 。给q个询问,每次询问一个区间[l,r],要求你求出一对数字(i,j)(l<=i<=j<=r),使得a[i]=a[j]并且dis(i,j)最大,由于这样的数对可能有多个,因此答案只要输出dis。

N<=10^5 Q<=10^4 1<=a[i]<=10^3 1<=l<=r<=n

题解

莫队的时候,维护区间中每个值出现的最左下标和最右下标。

例题4

hdu4638, 给你一个长度为n的序列,m个询问,每次查询给一个区间,我们来划分这个区间,对于划分出的任意一个子区间,他的所有元素构成的集合所包含的数必须是一段连续的区间,然后定义区间的代价为区间长度的平方,定义划分的代价为划分结束后所有划分出的子区间的代价的和,询问此代价最高时的划分出的子区间的个数。

题解

可以证明,划分唯一,我们把划分操作看成切割原区间,定义划分的操作为切割点的位置的集合,若两个划分不一样,我们对这两个划分操作切割点集取交集,交集构成新的划分,此划分优于前两个划分。于是按照能合并就合并的贪心,我们可以找的最优划分。考虑莫队,当区间变化的时候:添加一个数,如果他左右的数字都存在,显然段数-1 若左右存在某一个,段数不变,若左右均不存在,则段数+1 ,删除同理。

例题5

hdu4676,给你一个排列,q个区间询问,求$\sum_{L\leq i<j\leq R}gcd(a_i,a_j) $

题解

反演出结果后,直接莫队,对新加进来的数,枚举他的因子,计算莫队区间的变化。
$$
\begin{aligned}
&n=id(n)=(1*\varphi)(n)=\sum_{d|n}\varphi(d)
\&\sum_{L\leq i<j\leq R}gcd(a_i,a_j)
\=&\sum_{L\leq i<j\leq R}\sum_{d|gcd(a_i,a_j)}\varphi(d)
\=&\sum_{L\leq i<j\leq R}\sum_{d|a_i}\sum_{d|a_j}\varphi(d)
\=&\sum_{d}\varphi(d)*C_{\sum_{i=L}^{R}[d|a_i]}^2
\end{aligned}
$$

最大团

DFS计算最大团

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include<bits/stdc++.h>
using namespace std;

const int maxn=105;

bool edge[maxn][maxn]; int vertn; //

int dfs_ans,found,mcp[maxn],suf[maxn][maxn]; // dfs
void dfs(int d){ // begin with d =1 which means choose one point
if(suf[d][0]==0){
if(dfs_ans<d) {
dfs_ans=d;
found=1;
}
return ;
}

for(int i=1;i<=suf[d][0]&&!found;++i) {
if(d+suf[d][0]-i+1<=dfs_ans) break;// cut
if(d+mcp[suf[d][i]]<=dfs_ans) break;// cut
for(int j=i+1;j<=suf[d][0];++j){
if(edge[suf[d][i]][suf[d][j]]) suf[d+1][++suf[d+1][0]]=suf[d][j];
}
dfs(d+1);
suf[d+1][0]=0;
}
}

int max_cluster(){
mcp[vertn+1]=0;
for(int i=vertn;i>=1;i--) {
for(int j=i+1;j<=vertn;++j) if(edge[i][j]) suf[1][++suf[1][0]]=j;
dfs(1);
found=0;
suf[1][0]=0;
mcp[i]=dfs_ans;
}
return mcp[1];
}

bool check(int x){
int s=sqrt(x)+0.5;
return s*s==x;
}
int main(){
vertn=100;
for(int i=1;i<=vertn;++i){
for(int j=1;j<=vertn;++j){
if(check(i)||check(j)||check(i+j)) edge[i][j]=0;
else edge[i][j]=1;
}
}
cout<<max_cluster()<<endl;
}

bzoj2121

题意

BX正在进行一个字符串游戏,他手上有一个字符串L,以及其他一些字符串的集合S,然后他可以进行以下操作:对于一个在集合S中的字符串p,如果p在L中出现,BX就可以选择是否将其删除,如果删除,则将删除后L分裂成的左右 两部分合并。举个例子,L=’abcdefg’ , S={‘de’},如果BX选择将’de’从L中删去,则删后的L=’abcfg’。现在BX可 以进行任意多次操作(删的次数,顺序都随意),他想知道最后L串的最短长度是多少。

限制条件

|L|<=150

|S[i]|<=20

|S|<=30

Read More

后缀自动机

前言

本文力争用理性分析的手段,来推测此算法发明者的思维过程, 尝试感受其在设计此算法的时所展现出的思维方式, 力求用数学证明的手段,尽可能多的为读者证明相关结论,建议有其他自动机学习的基础,最好已经学会AC自动机和回文自动机,后缀自动机很难,他 和其他自动机不一样,它的状态更加复杂,一个算法的创作过程很 复杂,学起来当然会感到很难。强烈建议看陈立杰的ppt,看一遍肯定看不懂,仔细看,一遍看不懂看多遍,第一次可能只能看懂几面,第二次可 能就能看懂到十几面了,慢慢的就全懂了。

Read More

Educational Codeforces Round 60 (Rated for Div. 2) - D

链接

http://codeforces.com/contest/1117/problem/D

题意

你有两个数字:1和m,

你需要构造一个和为n的序列,问你能构造出多少种序列。答案对$10^9+7$取模。

范围

$N\le10^{18}$

$M\le100$

题解

ans(n,m) = ans(n-1,m) + ans(n-m,m)

由于m一样,所以dp[n] = dp[n-1] + dp[n-m]

矩阵快速幂

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include<bits/stdc++.h>
using namespace std;

//函数参数输入,返回值输出
//特别注意lenth一定要改,不然每次都遍历到最大的矩阵会tle
typedef long long ll;
ll lenth=200;
ll mod=1e9+7;
struct Sarray{
ll data[200][200];
//看着改
Sarray operator *(const Sarray&a){
Sarray tem;
for(ll i=0;i<lenth;i++){
for(ll j=0;j<lenth;j++){
for(ll k=0;k<lenth;k++){
tem.data[i][j]=(tem.data[i][j]+data[i][k]*a.data[k][j])%mod;
}
}
}
return tem;
}

Sarray operator +(const Sarray&a){
Sarray tem;
for(ll i=0;i<lenth;i++){
for(ll j=0;j<lenth;j++){
tem.data[i][j]=(data[i][j]+a.data[i][j])%mod;
}
}
return tem;
}

Sarray(const Sarray&a){
for(ll i=0;i<lenth;i++){
for(ll j=0;j<lenth;j++){
data[i][j]=a.data[i][j];
}
}
}

Sarray(){
memset(data,0,sizeof(data));
}

};

Sarray E;
void ini(){
for(ll i=0;i<lenth;i++)
for(ll j=0;j<lenth;j++)
if(i==j)E.data[i][j]=1;
else E.data[i][j]=0;
}

Sarray qpow(Sarray a,ll b){
Sarray tem=E;
while(b){
if(b&1)tem=a*tem;
a=a*a;
b>>=1;
}
return tem;
}

//sigma(p^i) from 0 to n [0,n]
Sarray sigma(Sarray&p,ll n){
if(n==0)return E;
if(n==1)return E+p;
if(n&1) return (E+qpow(p,n/2+1))*sigma(p,n>>1);
else return (E+qpow(p,n/2+1))*sigma(p,n/2-1)+qpow(p,n>>1);
}
///////////////////////////////////////


int main(){
ini();

ll n,m;
cin>>n>>m;
lenth=m;
Sarray p,r;

p.data[0][0]=1;
p.data[0][m-1]=1;
for(ll i=1;i<m;i++){
p.data[i][i-1]=1;
}

r.data[0][0]=2;
for(ll i=1;i<m;i++){
r.data[i][0]=1;
}

if(n>=m){
Sarray l=qpow(p,n-m)*r;
cout<<l.data[0][0]<<endl;
}
else{
cout<<1<<endl;
}
}