Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

326 lignes
8.9 KiB

  1. //
  2. // Lol Engine - Sample math program: Chebyshev polynomials
  3. //
  4. // Copyright: (c) 2005-2011 Sam Hocevar <sam@hocevar.net>
  5. // This program is free software; you can redistribute it and/or
  6. // modify it under the terms of the Do What The Fuck You Want To
  7. // Public License, Version 2, as published by Sam Hocevar. See
  8. // http://sam.zoy.org/projects/COPYING.WTFPL for more details.
  9. //
  10. #if !defined __REMEZ_SOLVER_H__
  11. #define __REMEZ_SOLVER_H__
  12. template<int ORDER> class RemezSolver
  13. {
  14. public:
  15. typedef real RealFunc(real const &x);
  16. RemezSolver()
  17. {
  18. }
  19. void Run(real a, real b, RealFunc *func, RealFunc *weight, int steps)
  20. {
  21. m_func = func;
  22. m_weight = weight;
  23. m_k1 = (b + a) >> 1;
  24. m_k2 = (b - a) >> 1;
  25. m_invk2 = re(m_k2);
  26. m_invk1 = -m_k1 * m_invk2;
  27. Init();
  28. PrintPoly();
  29. for (int n = 0; n < steps; n++)
  30. {
  31. FindExtrema();
  32. Step();
  33. PrintPoly();
  34. FindZeroes();
  35. }
  36. FindExtrema();
  37. Step();
  38. PrintPoly();
  39. }
  40. real ChebyEval(real const &x)
  41. {
  42. real ret = 0.0, xn = 1.0;
  43. for (int i = 0; i < ORDER + 1; i++)
  44. {
  45. real mul = 0;
  46. for (int j = 0; j < ORDER + 1; j++)
  47. mul += coeff[j] * (real)Cheby(j, i);
  48. ret += mul * xn;
  49. xn *= x;
  50. }
  51. return ret;
  52. }
  53. void Init()
  54. {
  55. /* Pick up x_i where error will be 0 and compute f(x_i) */
  56. real fxn[ORDER + 1];
  57. for (int i = 0; i < ORDER + 1; i++)
  58. {
  59. zeroes[i] = (real)(2 * i - ORDER) / (real)(ORDER + 1);
  60. fxn[i] = Value(zeroes[i]);
  61. }
  62. /* We build a matrix of Chebishev evaluations: row i contains the
  63. * evaluations of x_i for polynomial order n = 0, 1, ... */
  64. Matrix<ORDER + 1> mat;
  65. for (int i = 0; i < ORDER + 1; i++)
  66. {
  67. /* Compute the powers of x_i */
  68. real powers[ORDER + 1];
  69. powers[0] = 1.0;
  70. for (int n = 1; n < ORDER + 1; n++)
  71. powers[n] = powers[n - 1] * zeroes[i];
  72. /* Compute the Chebishev evaluations at x_i */
  73. for (int n = 0; n < ORDER + 1; n++)
  74. {
  75. real sum = 0.0;
  76. for (int k = 0; k < ORDER + 1; k++)
  77. sum += (real)Cheby(n, k) * powers[k];
  78. mat.m[i][n] = sum;
  79. }
  80. }
  81. /* Solve the system */
  82. mat = mat.inv();
  83. /* Compute interpolation coefficients */
  84. for (int j = 0; j < ORDER + 1; j++)
  85. {
  86. coeff[j] = 0;
  87. for (int i = 0; i < ORDER + 1; i++)
  88. coeff[j] += mat.m[j][i] * fxn[i];
  89. }
  90. }
  91. void FindZeroes()
  92. {
  93. /* Find ORDER + 1 zeroes of the error function. No need to
  94. * compute the relative error: its zeroes are at the same
  95. * place as the absolute error! */
  96. for (int i = 0; i < ORDER + 1; i++)
  97. {
  98. struct { real value, error; } left, right, mid;
  99. left.value = control[i];
  100. left.error = ChebyEval(left.value) - Value(left.value);
  101. right.value = control[i + 1];
  102. right.error = ChebyEval(right.value) - Value(right.value);
  103. static real limit = real::R_1 >> 500;
  104. while (fabs(left.value - right.value) > limit)
  105. {
  106. mid.value = (left.value + right.value) >> 1;
  107. mid.error = ChebyEval(mid.value) - Value(mid.value);
  108. if ((left.error < real::R_0 && mid.error < real::R_0)
  109. || (left.error > real::R_0 && mid.error > real::R_0))
  110. left = mid;
  111. else
  112. right = mid;
  113. }
  114. zeroes[i] = mid.value;
  115. }
  116. }
  117. void FindExtrema()
  118. {
  119. /* Find ORDER + 2 extrema of the error function. We need to
  120. * compute the relative error, since its extrema are at slightly
  121. * different locations than the absolute error’s. */
  122. real final = 0;
  123. for (int i = 0; i < ORDER + 2; i++)
  124. {
  125. real a = -1, b = 1;
  126. if (i > 0)
  127. a = zeroes[i - 1];
  128. if (i < ORDER + 1)
  129. b = zeroes[i];
  130. for (;;)
  131. {
  132. real c = a, delta = (b - a) >> 3;
  133. real maxerror = 0;
  134. real maxweight = 0;
  135. int best = -1;
  136. for (int k = 1; k <= 7; k++)
  137. {
  138. real error = ChebyEval(c) - Value(c);
  139. real weight = Weight(c);
  140. if (fabs(error * maxweight) >= fabs(maxerror * weight))
  141. {
  142. maxerror = error;
  143. maxweight = weight;
  144. best = k;
  145. }
  146. c += delta;
  147. }
  148. b = a + (real)(best + 1) * delta;
  149. a = a + (real)(best - 1) * delta;
  150. if (b - a < (real)1e-18)
  151. {
  152. real e = maxerror / maxweight;
  153. if (e > final)
  154. final = e;
  155. control[i] = (a + b) >> 1;
  156. break;
  157. }
  158. }
  159. }
  160. printf("Final error: ");
  161. final.print(40);
  162. }
  163. void Step()
  164. {
  165. /* Pick up x_i where error will be 0 and compute f(x_i) */
  166. real fxn[ORDER + 2];
  167. for (int i = 0; i < ORDER + 2; i++)
  168. fxn[i] = Value(control[i]);
  169. /* We build a matrix of Chebishev evaluations: row i contains the
  170. * evaluations of x_i for polynomial order n = 0, 1, ... */
  171. Matrix<ORDER + 2> mat;
  172. for (int i = 0; i < ORDER + 2; i++)
  173. {
  174. /* Compute the powers of x_i */
  175. real powers[ORDER + 1];
  176. powers[0] = 1.0;
  177. for (int n = 1; n < ORDER + 1; n++)
  178. powers[n] = powers[n - 1] * control[i];
  179. /* Compute the Chebishev evaluations at x_i */
  180. for (int n = 0; n < ORDER + 1; n++)
  181. {
  182. real sum = 0.0;
  183. for (int k = 0; k < ORDER + 1; k++)
  184. sum += (real)Cheby(n, k) * powers[k];
  185. mat.m[i][n] = sum;
  186. }
  187. if (i & 1)
  188. mat.m[i][ORDER + 1] = fabs(Weight(control[i]));
  189. else
  190. mat.m[i][ORDER + 1] = -fabs(Weight(control[i]));
  191. }
  192. /* Solve the system */
  193. mat = mat.inv();
  194. /* Compute interpolation coefficients */
  195. for (int j = 0; j < ORDER + 1; j++)
  196. {
  197. coeff[j] = 0;
  198. for (int i = 0; i < ORDER + 2; i++)
  199. coeff[j] += mat.m[j][i] * fxn[i];
  200. }
  201. /* Compute the error */
  202. real error = 0;
  203. for (int i = 0; i < ORDER + 2; i++)
  204. error += mat.m[ORDER + 1][i] * fxn[i];
  205. }
  206. int Cheby(int n, int k)
  207. {
  208. if (k > n || k < 0)
  209. return 0;
  210. if (n <= 1)
  211. return (n ^ k ^ 1) & 1;
  212. return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
  213. }
  214. int Comb(int n, int k)
  215. {
  216. if (k == 0 || k == n)
  217. return 1;
  218. return Comb(n - 1, k - 1) + Comb(n - 1, k);
  219. }
  220. void PrintPoly()
  221. {
  222. /* Transform Chebyshev polynomial weights into powers of X^i
  223. * in the [-1..1] range. */
  224. real bn[ORDER + 1];
  225. for (int i = 0; i < ORDER + 1; i++)
  226. {
  227. bn[i] = 0;
  228. for (int j = 0; j < ORDER + 1; j++)
  229. bn[i] += coeff[j] * (real)Cheby(j, i);
  230. }
  231. /* Transform a polynomial in the [-1..1] range into a polynomial
  232. * in the [a..b] range. */
  233. real k1p[ORDER + 1], k2p[ORDER + 1];
  234. real an[ORDER + 1];
  235. for (int i = 0; i < ORDER + 1; i++)
  236. {
  237. k1p[i] = i ? k1p[i - 1] * m_invk1 : real::R_1;
  238. k2p[i] = i ? k2p[i - 1] * m_invk2 : real::R_1;
  239. }
  240. for (int i = 0; i < ORDER + 1; i++)
  241. {
  242. an[i] = 0;
  243. for (int j = i; j < ORDER + 1; j++)
  244. an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
  245. an[i] *= k2p[i];
  246. }
  247. printf("Polynomial estimate:\n");
  248. for (int j = 0; j < ORDER + 1; j++)
  249. {
  250. if (j)
  251. printf("+");
  252. printf("x^%i*", j);
  253. an[j].print(40);
  254. }
  255. printf("\n");
  256. }
  257. real Value(real const &x)
  258. {
  259. return m_func(x * m_k2 + m_k1);
  260. }
  261. real Weight(real const &x)
  262. {
  263. return m_weight(x * m_k2 + m_k1);
  264. }
  265. /* ORDER + 1 Chebyshev coefficients and 1 error value */
  266. real coeff[ORDER + 2];
  267. /* ORDER + 1 zeroes of the error function */
  268. real zeroes[ORDER + 1];
  269. /* ORDER + 2 control points */
  270. real control[ORDER + 2];
  271. private:
  272. RealFunc *m_func, *m_weight;
  273. real m_k1, m_k2, m_invk1, m_invk2;
  274. };
  275. #endif /* __REMEZ_SOLVER_H__ */