math

发表于 分类于 阅读次数： 本文字数： 2.5k 阅读时长 ≈ 2 分钟

总览

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

项目地址

链接

先介绍我们的数学函数库 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})\), 更新： 我刚刚突然想到如果第二步使用泰勒展开而不是牛顿迭代，可能会快很多，考虑到这一点，我们有时间去实现一下泰勒展开求对数函数。