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
103
104
105
106
107
108
109
110
111
112
113
114
|
#include "Halide.h"
#include <stdio.h>
using namespace Halide;
float inverse_factorial(int x) {
double y = 1;
for (int i = 2; i <= x; i++) {
y /= i;
}
return (float)y;
}
int main(int argc, char **argv) {
Func f;
Var x;
// Create x scaled down by 256. We're going to intentionally do
// something numerically unstable below, so we prevent folding out
// the /256, or large powers of x will be inf
Expr xf = strict_float(x / 256.0f);
{
// Compute the taylor series approximation for sin
// x - x^3 / 3! + x^5 / 5! + x^7 / 7! ....
Expr y1 = 0.0f;
for (int k = 0; k < 20; k++) {
y1 += Halide::pow(-1, k) * pow(xf, 1 + 2 * k) * inverse_factorial(1 + 2 * k);
}
Func approx_sin_1;
approx_sin_1(x) = y1;
// Try a different way to express the Taylor series that should
// have fewer numerical precision issues. The large inverse
// factorials in the previous version tend to disappear entirely.
// x*(1 - x*x/(2*3) * (1 - x*x/(4*5) * (1 - x*x/(6*7) * (1 - x*x/(8*9) * ( ... )))))
Expr y2 = 1.0f;
for (int k = 20; k > 0; k--) {
y2 = 1 - (y2 * pow(xf, 2)) / (2 * k * (2 * k + 1));
}
y2 *= xf;
Func approx_sin_2;
approx_sin_2(x) = y2;
Func exact_sin;
exact_sin(x) = sin(xf);
// Evaluate from 0 to 5
Buffer<float> approx_result_1 = approx_sin_1.realize({256 * 5});
Buffer<float> approx_result_2 = approx_sin_2.realize({256 * 5});
Buffer<float> exact_result = exact_sin.realize({256 * 5});
Func rms_1, rms_2;
RDom r(exact_result);
rms_1() = sqrt(sum(pow(approx_result_1(r) - exact_result(r), 2), "rms_1_sum"));
rms_2() = sqrt(sum(pow(approx_result_2(r) - exact_result(r), 2), "rms_2_sum"));
Buffer<float> error_1 = rms_1.realize();
Buffer<float> error_2 = rms_2.realize();
if (error_1(0) > 0.0001 || error_2(0) > 0.0001) {
printf("Approximate sin errors too large: %1.20f %1.20f\n", error_1(0), error_2(0));
return 1;
}
}
{
xf = xf + 1;
// Now take negative powers for a spin.
// exp(1/x) = 1 + 1/x + 1/(2*x^2) + 1/(6*x^3) + ...
Func approx_exp_1;
Expr y1 = 0.0f;
for (int k = 0; k < 20; k++) {
y1 += pow(xf, -k) * inverse_factorial(k);
}
approx_exp_1(x) = y1;
// A different factorization.
// exp(1/x) = 1 + (1 + (1 + (1 + ...)/(3x))/(2x))/x;
Func approx_exp_2;
Expr y2 = 0.0f;
for (int k = 20; k > 0; k--) {
y2 = 1 + y2 / (k * xf);
}
approx_exp_2(x) = y2;
Func exact_exp;
exact_exp(x) = exp(1.0f / xf);
// Evaluate from 0 to 5
Buffer<float> approx_result_1 = approx_exp_1.realize({256 * 5});
Buffer<float> approx_result_2 = approx_exp_2.realize({256 * 5});
Buffer<float> exact_result = exact_exp.realize({256 * 5});
Func rms_1, rms_2;
RDom r(exact_result);
rms_1() = sqrt(sum(pow(approx_result_1(r) - exact_result(r), 2), "rms_1_neg_sum"));
rms_2() = sqrt(sum(pow(approx_result_2(r) - exact_result(r), 2), "rms_2_neg_sum"));
Buffer<float> error_1 = rms_1.realize();
Buffer<float> error_2 = rms_2.realize();
if (error_1(0) > 0.0001 || error_2(0) > 0.0001) {
printf("Approximate exp errors too large: %1.20f %1.20f\n", error_1(0), error_2(0));
return 1;
}
}
printf("Success!\n");
return 0;
}
|
