Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

solver.cpp 9.0 KiB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359
  1. //
  2. // LolRemez - Remez algorithm implementation
  3. //
  4. // Copyright: (c) 2005-2013 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://www.wtfpl.net/ for more details.
  9. //
  10. #if defined HAVE_CONFIG_H
  11. # include "config.h"
  12. #endif
  13. #include "core.h"
  14. #include <lol/math/real.h>
  15. #include "matrix.h"
  16. #include "solver.h"
  17. using lol::real;
  18. RemezSolver::RemezSolver(int order, int decimals)
  19. : m_order(order),
  20. m_decimals(decimals)
  21. {
  22. }
  23. void RemezSolver::Run(real a, real b,
  24. RemezSolver::RealFunc *func,
  25. RemezSolver::RealFunc *weight)
  26. {
  27. using std::printf;
  28. m_func = func;
  29. m_weight = weight;
  30. m_k1 = (b + a) / 2;
  31. m_k2 = (b - a) / 2;
  32. m_invk2 = re(m_k2);
  33. m_invk1 = -m_k1 * m_invk2;
  34. m_epsilon = pow((real)10, (real)-(m_decimals + 2));
  35. Init();
  36. PrintPoly();
  37. real error = -1;
  38. for (int n = 0; ; n++)
  39. {
  40. real newerror = FindExtrema();
  41. printf("Step %i error: ", n);
  42. newerror.print(m_decimals);
  43. printf("\n");
  44. Step();
  45. if (error >= (real)0 && fabs(newerror - error) < error * m_epsilon)
  46. break;
  47. error = newerror;
  48. PrintPoly();
  49. FindZeroes();
  50. }
  51. PrintPoly();
  52. }
  53. real RemezSolver::EvalCheby(real const &x)
  54. {
  55. real ret = 0.0, xn = 1.0;
  56. for (int i = 0; i < m_order + 1; i++)
  57. {
  58. real mul = 0;
  59. for (int j = 0; j < m_order + 1; j++)
  60. mul += m_coeff[j] * (real)Cheby(j, i);
  61. ret += mul * xn;
  62. xn *= x;
  63. }
  64. return ret;
  65. }
  66. void RemezSolver::Init()
  67. {
  68. /* m_order + 1 Chebyshev coefficients, plus 1 error value */
  69. m_coeff.Resize(m_order + 2);
  70. /* m_order + 1 zeroes of the error function */
  71. m_zeroes.Resize(m_order + 1);
  72. /* m_order + 2 control points */
  73. m_control.Resize(m_order + 2);
  74. /* Pick up x_i where error will be 0 and compute f(x_i) */
  75. real fxn[m_order + 1];
  76. for (int i = 0; i < m_order + 1; i++)
  77. {
  78. m_zeroes[i] = (real)(2 * i - m_order) / (real)(m_order + 1);
  79. fxn[i] = EvalFunc(m_zeroes[i]);
  80. }
  81. /* We build a matrix of Chebishev evaluations: row i contains the
  82. * evaluations of x_i for polynomial order n = 0, 1, ... */
  83. Matrix<real> mat(m_order + 1, m_order + 1);
  84. for (int i = 0; i < m_order + 1; i++)
  85. {
  86. /* Compute the powers of x_i */
  87. real powers[m_order + 1];
  88. powers[0] = 1.0;
  89. for (int n = 1; n < m_order + 1; n++)
  90. powers[n] = powers[n - 1] * m_zeroes[i];
  91. /* Compute the Chebishev evaluations at x_i */
  92. for (int n = 0; n < m_order + 1; n++)
  93. {
  94. real sum = 0.0;
  95. for (int k = 0; k < m_order + 1; k++)
  96. sum += (real)Cheby(n, k) * powers[k];
  97. mat.m(i, n) = sum;
  98. }
  99. }
  100. /* Solve the system */
  101. mat = mat.inv();
  102. /* Compute interpolation coefficients */
  103. for (int j = 0; j < m_order + 1; j++)
  104. {
  105. m_coeff[j] = 0;
  106. for (int i = 0; i < m_order + 1; i++)
  107. m_coeff[j] += mat.m(j, i) * fxn[i];
  108. }
  109. }
  110. void RemezSolver::FindZeroes()
  111. {
  112. /* Find m_order + 1 zeroes of the error function. No need to
  113. * compute the relative error: its zeroes are at the same
  114. * place as the absolute error! */
  115. for (int i = 0; i < m_order + 1; i++)
  116. {
  117. struct { real value, error; } left, right, mid;
  118. left.value = m_control[i];
  119. left.error = EvalCheby(left.value) - EvalFunc(left.value);
  120. right.value = m_control[i + 1];
  121. right.error = EvalCheby(right.value) - EvalFunc(right.value);
  122. static real limit = ldexp((real)1, -500);
  123. static real zero = (real)0;
  124. while (fabs(left.value - right.value) > limit)
  125. {
  126. mid.value = (left.value + right.value) / 2;
  127. mid.error = EvalCheby(mid.value) - EvalFunc(mid.value);
  128. if ((left.error <= zero && mid.error <= zero)
  129. || (left.error >= zero && mid.error >= zero))
  130. left = mid;
  131. else
  132. right = mid;
  133. }
  134. m_zeroes[i] = mid.value;
  135. }
  136. }
  137. real RemezSolver::FindExtrema()
  138. {
  139. using std::printf;
  140. /* Find m_order + 2 extrema of the error function. We need to
  141. * compute the relative error, since its extrema are at slightly
  142. * different locations than the absolute error’s. */
  143. real final = 0;
  144. for (int i = 0; i < m_order + 2; i++)
  145. {
  146. real a = -1, b = 1;
  147. if (i > 0)
  148. a = m_zeroes[i - 1];
  149. if (i < m_order + 1)
  150. b = m_zeroes[i];
  151. for (int round = 0; ; round++)
  152. {
  153. real maxerror = 0, maxweight = 0;
  154. int best = -1;
  155. real c = a, delta = (b - a) / 4;
  156. for (int k = 0; k <= 4; k++)
  157. {
  158. if (round == 0 || (k & 1))
  159. {
  160. real error = fabs(EvalCheby(c) - EvalFunc(c));
  161. real weight = fabs(Weight(c));
  162. /* if error/weight >= maxerror/maxweight */
  163. if (error * maxweight >= maxerror * weight)
  164. {
  165. maxerror = error;
  166. maxweight = weight;
  167. best = k;
  168. }
  169. }
  170. c += delta;
  171. }
  172. switch (best)
  173. {
  174. case 0:
  175. b = a + delta * 2;
  176. break;
  177. case 4:
  178. a = b - delta * 2;
  179. break;
  180. default:
  181. b = a + delta * (best + 1);
  182. a = a + delta * (best - 1);
  183. break;
  184. }
  185. if (delta < m_epsilon)
  186. {
  187. real e = fabs(maxerror / maxweight);
  188. if (e > final)
  189. final = e;
  190. m_control[i] = (a + b) / 2;
  191. break;
  192. }
  193. }
  194. }
  195. return final;
  196. }
  197. void RemezSolver::Step()
  198. {
  199. /* Pick up x_i where error will be 0 and compute f(x_i) */
  200. real fxn[m_order + 2];
  201. for (int i = 0; i < m_order + 2; i++)
  202. fxn[i] = EvalFunc(m_control[i]);
  203. /* We build a matrix of Chebishev evaluations: row i contains the
  204. * evaluations of x_i for polynomial order n = 0, 1, ... */
  205. Matrix<real> mat(m_order + 2, m_order + 2);
  206. for (int i = 0; i < m_order + 2; i++)
  207. {
  208. /* Compute the powers of x_i */
  209. real powers[m_order + 1];
  210. powers[0] = 1.0;
  211. for (int n = 1; n < m_order + 1; n++)
  212. powers[n] = powers[n - 1] * m_control[i];
  213. /* Compute the Chebishev evaluations at x_i */
  214. for (int n = 0; n < m_order + 1; n++)
  215. {
  216. real sum = 0.0;
  217. for (int k = 0; k < m_order + 1; k++)
  218. sum += (real)Cheby(n, k) * powers[k];
  219. mat.m(i, n) = sum;
  220. }
  221. if (i & 1)
  222. mat.m(i, m_order + 1) = fabs(Weight(m_control[i]));
  223. else
  224. mat.m(i, m_order + 1) = -fabs(Weight(m_control[i]));
  225. }
  226. /* Solve the system */
  227. mat = mat.inv();
  228. /* Compute interpolation coefficients */
  229. for (int j = 0; j < m_order + 1; j++)
  230. {
  231. m_coeff[j] = 0;
  232. for (int i = 0; i < m_order + 2; i++)
  233. m_coeff[j] += mat.m(j, i) * fxn[i];
  234. }
  235. /* Compute the error */
  236. real error = 0;
  237. for (int i = 0; i < m_order + 2; i++)
  238. error += mat.m(m_order + 1, i) * fxn[i];
  239. }
  240. int RemezSolver::Cheby(int n, int k)
  241. {
  242. if (k > n || k < 0)
  243. return 0;
  244. if (n <= 1)
  245. return (n ^ k ^ 1) & 1;
  246. return 2 * Cheby(n - 1, k - 1) - Cheby(n - 2, k);
  247. }
  248. int RemezSolver::Comb(int n, int k)
  249. {
  250. if (k == 0 || k == n)
  251. return 1;
  252. return Comb(n - 1, k - 1) + Comb(n - 1, k);
  253. }
  254. void RemezSolver::PrintPoly()
  255. {
  256. using std::printf;
  257. /* Transform Chebyshev polynomial weights into powers of X^i
  258. * in the [-1..1] range. */
  259. real bn[m_order + 1];
  260. for (int i = 0; i < m_order + 1; i++)
  261. {
  262. bn[i] = 0;
  263. for (int j = 0; j < m_order + 1; j++)
  264. bn[i] += m_coeff[j] * (real)Cheby(j, i);
  265. }
  266. /* Transform a polynomial in the [-1..1] range into a polynomial
  267. * in the [a..b] range. */
  268. real k1p[m_order + 1], k2p[m_order + 1];
  269. real an[m_order + 1];
  270. for (int i = 0; i < m_order + 1; i++)
  271. {
  272. k1p[i] = i ? k1p[i - 1] * m_invk1 : (real)1;
  273. k2p[i] = i ? k2p[i - 1] * m_invk2 : (real)1;
  274. }
  275. for (int i = 0; i < m_order + 1; i++)
  276. {
  277. an[i] = 0;
  278. for (int j = i; j < m_order + 1; j++)
  279. an[i] += (real)Comb(j, i) * k1p[j - i] * bn[j];
  280. an[i] *= k2p[i];
  281. }
  282. printf("Polynomial estimate: ");
  283. for (int j = 0; j < m_order + 1; j++)
  284. {
  285. if (j)
  286. printf(" + x**%i * ", j);
  287. an[j].print(m_decimals);
  288. }
  289. printf("\n\n");
  290. }
  291. real RemezSolver::EvalFunc(real const &x)
  292. {
  293. return m_func(x * m_k2 + m_k1);
  294. }
  295. real RemezSolver::Weight(real const &x)
  296. {
  297. if (m_weight)
  298. return m_weight(x * m_k2 + m_k1);
  299. return 1;
  300. }