高精度误差

关于高精度误差

高精度问题

高精度问题总是带给竞赛选手噩梦般的体验，每个oier选手都或多或少遇到过。有必要对高精度做一个说明，避免掉入高精度的大坑。

对于oier选手而言，float基本可以弃用了。表示实数，最少也要用double的。

数据类型	字节数	有效位数/数据范围
double	8	16~17位
long double	16	19位
int	4	10位 （-2147483648 ~2147483647）
long long int	8	18位 （-9223372036854775808 ~ 9223372036854775807 ）

一. 类型转换问题

1. 整数转成double类型(看有效位数)

对于double类型和long double类型，有效位数决定了该类型能够表示数值的精度。

如果一个 long long int 转成 double，因为double的有效位数只有16~17位，所以如果long long int 的数位超过17位，则转化为double后会出错。

long long int a = 123456789123456789;
double b = a;
printf("%.1f", b);
// 输出为:123456789123456784.0, 最后一位已经无法表示了。

而int最大只有10位，所以int类型转成double类型不会丢失有效位，没有任何问题。

总结：int 转成 double 是毫无问题， long long int 只有当数据特别大(超过16位)时，转成double会产生问题，丢失有效位； 不管是int，还是long long int， 转成long double 是没有任何问题的。因为long double的有效位数超过了int和long long int。

2. double类型转成整数

（1）double直接转成整数，是会去掉小数部分，保留整数部分。

而如果double的整数部分超过int范围，转成int，会发生溢出。

 	double a = 12345678901234.567;
    int d = (int) a;
    printf("%d", d);//此处发生溢出，输出的是-2147483648. double转int，只要发生溢出，都会得到这个值-2147483648。

（2）放大转换：

有时double的小数位数不太多，我们可能会将double放大若干倍，将它变成整数再处理，以避免精度问题。这种放大转换也有问题。

比如：

	double a = 19.40;
    int c = (int) (a * 10);
	printf("%d\n", c); //此处输出194
    int d =(int) (a * 100);
    printf("%d\n", d); //此处输出1939
    //原因是19.40无法用二进制精确表示，采用的是近似值：19.3999999999。让人崩溃的是第一个输出的为何不是193。

放大转换时，如果放大的倍数已经保证没有小数位数了，可以采用round，这样得到的结果就没有问题。

    double a = 19.40;
    int c = round (a * 10);
    int d =round (a * 100);
    printf("%d\n%d\n", c, d); // 输出194和1940

总结：double的小数部分，并不能精确的表达所有数值，许多时候，都是用一些近似值表示的，显示出来的结果是四舍五入的结果。

比如上面的19.40，实际上使用19.39999999近似表示的。所以，double总是会有一些莫名其妙的精度问题。

3. floor、ceil

floor(a)： 取小于等于a的最大整数值。

ceil(a)：取大于等于a的最小整数值。

round(a): 对小数部分四舍五入。

floor(-3.9) : 值为-4
floor(3.9): 值为3
ceil(-3.9): 值为-3
ceil(3.9): 值为4
round(3.9): 值为4
round(-3.9)：值为-4

因为double许多时候是以近似值来表示的，比如$(0.5^2-0.3^2-0.4^2)$, 本来应该是0，但是因为$0.3,0.4$都是近似值的关系，最后的结果实际上是接近0的一个近似值.看下面的例子:

 	double a = 0.3, b = 0.4, c = 0.5;
    printf("%.0f %.0f", floor(a - sqrt(c * c - b * b)), ceil(a - sqrt(c * c - b * b)));
    //输出为：0 1

采用long double也会有相似的问题。

    long double a = 0.0003, b = 0.0004, c = 0.0005;
    printf("%.0Lf %.0Lf", floor(a - sqrt(c * c - b * b)), ceil(a - sqrt(c * c - b * b)));
    //这下可好，输出的是 -1 0

为了解决floor，ceil的问题，我们需要用一个很小的量去修正它。

比如floor(x)，因为$x$可能是一个整点附近的值。我们应当给$x$加上一个很小的量$eps$,比如$eps$取值为$10^{-13}$

而求ceil(x), 我们需要让$x$减去一个很小的量$eps$.

例1：表格点

题意描述

有一个半径为$R$的圆，圆心为$(X,Y)$.
求圆内及圆周上的所有表格点的数量，所谓表格点，即整数坐标点。

输入格式

第一行三个数，表示$X,Y,R$, 每个数最多包含4位小数。

输出格式

一个整数，表示答案。

数据范围

$|X|\le 10^5,|Y|\le 10^5,0<R\le 10^5$

分析：本题可以枚举所有的与圆相交或相切的$y$值为整数的水平线，求出它与圆的两个交点$(x_1,y),(x_2,y)$,设$x_1\le x_2$,则该水平线上一共有$floor(x_2)-ceil(x_1)+1$ 个表格点。

$x_1,x_2$怎么求呢？

可以通过勾股定理来求,如下图所示：

$x_1=X-\sqrt{(R^2-(Y-y)\ast (Y-y))}$

$x_2=X+\sqrt{(R^2-(Y-y)\ast (Y-y))}$

但本题有两个细节要格外注意，否则可能因精度问题出错。

1. 需要用long double存储，因为double精度已经不够了。

比如$R=10^5,Y-y=10^{-5}$时，$\sqrt{(R^2-(Y-y)\ast (Y-y))}$的数值已经超过了double的有效位数，导致double不能准确的表示平方根，会带来误差。所以应该用long double，或者将所有坐标换成整数，用long long int。

2. 使用floor和ceil时，会有精度的问题。

比如$x_1=0.0003-\sqrt{(0.0005^2-0.0004^2)}$.

$x_1$应该为$0$，但是因为精度问题，$x_1$里面存储的是0的近似值，可能比$x_1$大，比如$0.000000000000000101$, 它与$0$是很接近的。此时,$ceil(x_1)$就到了$1$了,而$floor(x)$没有问题。

也有可能$x_1$存储的是小于$0$的一个近似值，比如$-0.000000000000001$之类，此时$ceil(x_1)$没有问题， 但$floor(x_1)$就等于$-1$了。这样，我们在统计答案$ans=floor(x_2)-ceil(x_1)+1$时，就可能出现错误。

可以设一个极小的值：$EPS=1E-14$，将答案修正为：

$ans=floor(x_2+EPS)-ceil(x_1-EPS)+1$

参考代码如下：

#include <iostream>
#include <cmath>
using namespace std;
#define rep(i, n) for (int i = 0; i < (n); i++)
using ll = long long;
using ld = long double;
using P = pair<int, int>;

ld X, Y, R;
const ld EPS = 1e-14;

int main() {
    cin >> X >> Y >> R;
    R += EPS;
    ll ans = 0;
    for (ll y = ceil(Y - R); y <= floor(Y + R); y++) {
        ll xl = ceil(X - sqrt(R * R - ((ld)y - Y) * ((ld)y - Y)) - EPS); //此处用EPS进行了修正
        ll xr = floor(X + sqrt(R * R - ((ld)y - Y) * ((ld)y - Y)) + EPS); // 此处用EPS进行了修正
        if (xr < xl)
            continue;
        ans += xr - xl + 1;
    }
    cout << ans << endl;
    return 0;
}

第二种方法：放大10000倍，将所有坐标点变成整数。

因为输入数据限定带4位小数，所以，放大10000倍，所有输入都变成了整数。但是，有些double数据，乘以10000倍，得到的整数却不是我们期望的值。

 double a = 0.0003;

 int b = a * 10000;

 printf("%d", b); //得到的是2

我们在转换时，要用round四舍五入，这样才能正确。

参考代码：

#include <bits/stdc++.h>
using namespace std;
#define LL long long int
#define ld long double
double X, Y, R;
#define EPS 1e-12
LL ymin, ymax, x1, x2, xd, yd, rd;
LL ans;
int main(){
    scanf("%lf %lf %lf", &X, &Y, &R);
    ymin = ceil(Y - R) * 10000;
    ymax = floor(Y + R) * 10000;
    xd = round(X * 10000), yd = round(Y * 10000), rd = round(R * 10000);
    //cerr << xd << ' ' << yd << ' ' << rd << ' '<< endl;
    for(LL i = ymin; i <= ymax; i = i + 10000){
        LL tmp = sqrt((ld) rd * rd - (ld) (yd - i) * (yd - i)); // 此处一定要转成ld，否则会有精度误差
        x1 = xd - tmp;
        x2 = xd + tmp;
        x1 = ceil(1.0 * x1 / 10000);
        x2 = floor(1.0 * x2 / 10000);
 
        ans += (x2 - x1) + 1;
    }
    printf("%lld", ans);
    return 0;
}