math

总览

这篇博客用于记录数学库的实现,

项目地址

链接

先介绍我们的数学函数库 math_function

这里记录了很多数学函数,
math_function代码
treeview raw
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
#pragma once
namespace math {
template <class T>
T abs(const T& a) {
return a > 0 ? a : -a;
}
template <class T>
T max(const T& a, const T& b) {
if (a < b)
return b;
else
return a;
}
template <class T>
T min(const T& a, const T& b) {
if (a < b)
return a;
else
return b;
}
template <class T>
T newton_sqrt(const T& x, const double& eps) {
return sqrt(x);
if (x == 0) return 0;
assert(x > 0);
T res = x;
while (true) {
T d = abs(res * res - x);
if (d / x < eps) break;
res = (res * res + x) / 2 / res;
}
return res;
}
const double PI = 3.1415926535897;
const double E = 2.71828182845904553488;

template <class T>
T sin(T x, const int& cnt) {
if (x < 0) x = -x + PI;
x = x - int(x / 2 / PI) * PI * 2;
if (x < 0) x += 2 * PI;
T x2 = x * x;
T res = 0, d = x;
for (int i = 1; i <= cnt; i++) {
res += d;
d *= -x2 / (2 * i) / (2 * i + 1);
}
return res;
}
template <class T>
T cos(T x, const int& cnt) {
return sin(x + PI / 2, cnt);
}

// e^x = 1 + x + x^2/2 +x^3/2/3 + x^4/2/3/4
double taylor_exp(const double& x) {
if (x < 0) return 1 / taylor_exp(-x);
if (x > 0.0001) { //快速幂
double res = taylor_exp(x / 2);
return res * res;
}
double res = 1, d = 1;
for (int i = 1; i <= 10; i++) {
d = d * x / i;
res += d;
}
return res;
}

// 用牛顿迭代法求解lg(c)=x
// -> e^x=c
// -> f(x)=e^x-c
// -> f'(x) = e^x
// x-(x)/f'(x) = x-1+c/e^x
double newton_log(double c) {
double res1 = 0, res2 = 0;
while (c < E) c *= E, res1--;
while (c > E) c /= E, res1++;
for (int i = 0; i < 10; i++) res2 -= 1 - c / taylor_exp(res2);
return res1 + res2;
}

} // namespace math

重大问题

这些算法的复杂度感觉都是\(lgc\)级别的,应该是可以通过倍增来达到\(lg(lgc)\)的复杂度,我下次再仔细思考思考。

介绍我们的牛顿迭代法求\(\sqrt{c}\)

\[ \begin{aligned} &x^2=c\\ &f(x) = x^2-c\\ &f'(x) = 2x \\ &g(x) = x-\frac{f(x)}{f'(x)} = x-\frac{x^2-c}{2x} =\frac{x^2+c}{2x} \end{aligned} \] 按照\(x_i=g(x_{i-1})\)进行迭代即可求出结果。 更新: 我下次找个机会用下面求\(e^c\)的方法实现一下先让c变小,看看能不能加速。

介绍我们的泰勒展开求\(e^c\)

首先根据公式\(e^{-t}=\frac{1}{e^t}\),可以递归为c恒非负 然后根据公式\(e^t=(e^\frac{t}{2})^2\), 可以递归为c在范围\([0,0.001]\)上 最后使用泰勒展开,\(e^x=1+x+\frac{x^2}{2!}+...\),这里我们取前10项就能够达到很高的精度了。 为什么要用第二步将c保证在范围\([0,0.001]\)上? 因为如果c过大,我们的第三部需要展开更多的项才够,这在c达到10000这种,你至少要展开10000项,这不现实。

介绍我们的牛顿迭代法求\(ln(c)\)

\[ \begin{aligned} &e^x=c \\&f(x)=e^x-c \\&f'(x) = e^x \\&g(x)=x-\frac{f(x)}{f'(x)} = x-1+\frac{c}{e^x} \end{aligned} \] 还是一样的,为了减少迭代次数,我们先对c进行变小,根据公式\(ln(x)=ln(\frac{x}{e})+1\),我们可以保证c的值在e附近, 最后使用迭代,\(x_i=g(x_{i-1})\), 更新: 我刚刚突然想到如果第二步使用泰勒展开而不是牛顿迭代,可能会快很多,考虑到这一点,我们有时间去实现一下泰勒展开求对数函数。