子非鱼
0%

浮点数精度判定

发表于 更新于 分类于

在测试Bjerksund and Stensland American Option Approximation模型的时候发现如下这段计算二维高斯分布的函数数值溢出了。开始发现是输入的数值过于极端导致的，但是后面仔细研究发现，C++内浮点数的影响要比预期的大。

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
using flow::model::norm_dist::cdf;

// Sign function
double sign(const double &a) {
    if (a < 0.)
        return -1.0;
    else if (a > 0.)
        return 1.0;
    else
        return 0.0;
}

// Support function
double f(const double &x, const double &y, const double &aprime, const double &bprime, const double &rho) {
    double r = aprime * (2.0 * x - aprime) + bprime * (2.0 * y - bprime) + 2.0 * rho * (x - aprime) * (y - bprime);
    return exp(r);
}

double bi_cdf(const double &a, const double &b, const double &rho) {
    double pi = 3.14159265359;
    if ((a <= 0.0) && (b <= 0.0) && (rho <= 0.0)) {
        double aprime = a / sqrt(2.0 * (1.0 - rho * rho));
        double bprime = b / sqrt(2.0 * (1.0 - rho * rho));
        double A[] = {0.3253030, 0.4211071, 0.1334425, 0.006374323};
        double B[] = {0.1337764, 0.6243247, 1.3425378, 2.2626645};
        double sum = 0.0;
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
                sum = sum + A[i] * A[j] * f(B[i], B[j], aprime, bprime, rho);
        sum = sum * (sqrt(1.0 - rho * rho) / pi);
        return sum;
    } else if (a * b * rho <= 0.0) {
        if ((a <= 0.0) && (b >= 0.0) && (rho >= 0.0))
            return cdf(a) - bi_cdf(a, -1.0 * b, -1.0 * rho);
        else if ((a >= 0.0) && (b <= 0.0) && (rho >= 0.0))
            return cdf(b) - bi_cdf(-1.0 * a, b, -1.0 * rho);
        else if ((a >= 0.0) && (b >= 0.0) && (rho <= 0.0))
            return cdf(a) + cdf(b) - 1.0 + bi_cdf(-1.0 * a, -1.0 * b, rho);
    } else if (a * b * rho >= 0.0) {
        double denum = sqrt(a * a - 2.0 * rho * a * b + b * b);
        double rho1 = ((rho * a - b) * sign(a)) / denum;
        double rho2 = ((rho * b - a) * sign(b)) / denum;
        double delta = (1.0 - sign(a) * sign(b)) / 4.0;
        return bi_cdf(a, 0.0, rho1) + bi_cdf(b, 0.0, rho2) - delta;
    } else {
        throw std::out_of_range("This part of the code should never be reached.");
    }
    throw std::out_of_range("This part of the code should never be reached.");
}

为了进一步增强模型实现代码“皮实耐操”的属性，添加了浮点数比较功能；

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
// equal
inline bool eq(const double &a, const double &b) {
return fabs(a - b) < 1e-6;
}

// not equal
inline bool neq(const double &a, const double &b) {
return !eq(a, b);
}

// great than
inline bool gt(const double &a, const double &b) {
return a > b && neq(a, b);
}

// great equal than
inline bool geq(const double &a, const double &b) {
return gt(a, b) || eq(a, b);
}

// less than
inline bool lt(const double &a, const double &b) {
return a < b && neq(a, b);
}

//less equal than
inline bool leq(const double &a, const double &b) {
return lt(a, b) || eq(a, b);
}

代码修改为

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
using flow::model::utils::geq;
using flow::model::utils::leq;
using flow::model::utils::lt;
using flow::model::utils::gt;
using flow::model::norm_dist::cdf;

// Sign function
double sign(const double &a) {
    if (lt(a , 0.))
        return -1.0;
    else if (gt(a , 0.))
        return 1.0;
    else
        return 0.0;
}

// Support function
double f(const double &x, const double &y, const double &a_prime, const double &b_prime, const double &rho) {
    double r = a_prime * (2.0 * x - a_prime) + b_prime * (2.0 * y - b_prime) + 2.0 * rho * (x - a_prime) * (y - b_prime);
    return exp(r);
}

double bi_cdf(const double &a, const double &b, const double &rho) {
    double pi = 3.14159265359;
    if (leq(a, 0.0) && leq(b , 0.0) && leq(rho , 0.0)) {
        double a_prime = a / sqrt(2.0 * (1.0 - rho * rho));
        double b_prime = b / sqrt(2.0 * (1.0 - rho * rho));
        double A[] = {0.3253030, 0.4211071, 0.1334425, 0.006374323};
        double B[] = {0.1337764, 0.6243247, 1.3425378, 2.2626645};
        double sum = 0.0;
        for (int i = 0; i < 4; i++)
            for (int j = 0; j < 4; j++)
                sum = sum + A[i] * A[j] * f(B[i], B[j], a_prime, b_prime, rho);
        sum = sum * (sqrt(1.0 - rho * rho) / pi);
        return sum;
    } else if (leq(a * b * rho , 0.0)) {
        if (leq(a , 0.0) && geq(b , 0.0) && geq(rho , 0.0))
            return cdf(a) - bi_cdf(a, -1.0 * b, -1.0 * rho);
        else if (geq(a , 0.0) && leq(b , 0.0) && geq(rho , 0.0))
            return cdf(b) - bi_cdf(-1.0 * a, b, -1.0 * rho);
        else if (geq(a , 0.0) && geq(b , 0.0) && leq(rho , 0.0))
            return cdf(a) + cdf(b) - 1.0 + bi_cdf(-1.0 * a, -1.0 * b, rho);
    } else if (geq(a * b * rho , 0.0)) {
        double de_num = sqrt(a * a - 2.0 * rho * a * b + b * b);
        double rho1 = ((rho * a - b) * sign(a)) / de_num;
        double rho2 = ((rho * b - a) * sign(b)) / de_num;
        double delta = (1.0 - sign(a) * sign(b)) / 4.0;
        return bi_cdf(a, 0.0, rho1) + bi_cdf(b, 0.0, rho2) - delta;
    } else {
        throw std::out_of_range("This part of the code should never be reached.");
    }
    throw std::out_of_range("This part of the code should never be reached.");
}

再进行极限测试发现数值溢出的情况消失。