123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225 |
- #pragma once
- #include <cmath>
- #include <functional>
- namespace WwiseGTE
- {
- template <typename Real>
- class RootsBrentsMethod
- {
- public:
-
-
-
-
-
- static bool Find(std::function<Real(Real)> const& F, Real t0, Real t1,
- unsigned int maxIterations, Real negFTolerance, Real posFTolerance,
- Real stepTTolerance, Real convTTolerance, Real& root)
- {
-
- if (t1 <= t0
- || maxIterations == 0
- || negFTolerance > (Real)0
- || posFTolerance < (Real)0
- || stepTTolerance < (Real)0
- || convTTolerance < (Real)0)
- {
-
- return false;
- }
- Real f0 = F(t0);
- if (negFTolerance <= f0 && f0 <= posFTolerance)
- {
-
-
- root = t0;
- return true;
- }
- Real f1 = F(t1);
- if (negFTolerance <= f1 && f1 <= posFTolerance)
- {
-
-
- root = t1;
- return true;
- }
- if (f0 * f1 > (Real)0)
- {
-
- return false;
- }
- if (std::fabs(f0) < std::fabs(f1))
- {
-
-
- std::swap(t0, t1);
- std::swap(f0, f1);
- }
-
- Real t2 = t0, t3 = t0, f2 = f0;
- bool prevBisected = true;
-
- for (unsigned int i = 0; i < maxIterations; ++i)
- {
- Real fDiff01 = f0 - f1, fDiff02 = f0 - f2, fDiff12 = f1 - f2;
- Real invFDiff01 = ((Real)1) / fDiff01;
- Real s;
- if (fDiff02 != (Real)0 && fDiff12 != (Real)0)
- {
-
- Real infFDiff02 = ((Real)1) / fDiff02;
- Real invFDiff12 = ((Real)1) / fDiff12;
- s =
- t0 * f1 * f2 * invFDiff01 * infFDiff02 -
- t1 * f0 * f2 * invFDiff01 * invFDiff12 +
- t2 * f0 * f1 * infFDiff02 * invFDiff12;
- }
- else
- {
-
- s = (t1 * f0 - t0 * f1) * invFDiff01;
- }
-
- Real tDiffSAvr = s - ((Real)0.75) * t0 - ((Real)0.25) * t1;
- Real tDiffS1 = s - t1;
- Real absTDiffS1 = std::fabs(tDiffS1);
- Real absTDiff12 = std::fabs(t1 - t2);
- Real absTDiff23 = std::fabs(t2 - t3);
- bool currBisected = false;
- if (tDiffSAvr * tDiffS1 > (Real)0)
- {
-
-
-
-
- currBisected = true;
- }
- else if (prevBisected)
- {
-
-
- currBisected =
- (absTDiffS1 >= ((Real)0.5) * absTDiff12) ||
- (absTDiff12 <= stepTTolerance);
- }
- else
- {
-
-
- currBisected =
- (absTDiffS1 >= ((Real)0.5) * absTDiff23) ||
- (absTDiff23 <= stepTTolerance);
- }
- if (currBisected)
- {
-
-
- s = ((Real)0.5) * (t0 + t1);
- if (s == t0 || s == t1)
- {
-
-
- root = s;
- return true;
- }
- prevBisected = true;
- }
- else
- {
- prevBisected = false;
- }
-
-
- Real fs = F(s);
- if (negFTolerance <= fs && fs <= posFTolerance)
- {
- root = s;
- return true;
- }
-
-
- t3 = t2;
- t2 = t1;
- f2 = f1;
- if (f0 * fs < (Real)0)
- {
- t1 = s;
- f1 = fs;
- }
- else
- {
- t0 = s;
- f0 = fs;
- }
-
-
- if (std::fabs(t1 - t0) <= convTTolerance)
- {
- root = t1;
- return true;
- }
-
-
- if (std::fabs(f0) < std::fabs(f1))
- {
- std::swap(t0, t1);
- std::swap(f0, f1);
- }
- }
-
- return false;
- }
- };
- }
|