1 /*
2 * Copyright Nick Thompson, 2018
3 * Use, modification and distribution are subject to the
4 * Boost Software License, Version 1.0. (See accompanying file
5 * LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 */
7 #include <iostream>
8 #include <vector>
9 #include <string>
10 #include <complex>
11 #include <bitset>
12 #include <boost/assert.hpp>
13 #include <boost/multiprecision/cpp_bin_float.hpp>
14 #include <boost/math/constants/constants.hpp>
15 #include <boost/math/tools/polynomial.hpp>
16 #include <boost/math/tools/roots.hpp>
17 #include <boost/math/special_functions/binomial.hpp>
18 #include <boost/multiprecision/cpp_complex.hpp>
19 #include <boost/multiprecision/complex128.hpp>
20 #include <boost/math/quadrature/gauss_kronrod.hpp>
21
22 using std::string;
23 using boost::math::tools::polynomial;
24 using boost::math::binomial_coefficient;
25 using boost::math::tools::schroder_iterate;
26 using boost::math::tools::halley_iterate;
27 using boost::math::tools::newton_raphson_iterate;
28 using boost::math::tools::complex_newton;
29 using boost::math::constants::half;
30 using boost::math::constants::root_two;
31 using boost::math::constants::pi;
32 using boost::math::quadrature::gauss_kronrod;
33 using boost::multiprecision::cpp_bin_float_100;
34 using boost::multiprecision::cpp_complex_100;
35
36 template<class Complex>
find_roots(size_t p)37 std::vector<std::pair<Complex, Complex>> find_roots(size_t p)
38 {
39 // Initialize the polynomial; see Mallat, A Wavelet Tour of Signal Processing, equation 7.96
40 BOOST_ASSERT(p>0);
41 typedef typename Complex::value_type Real;
42 std::vector<Complex> coeffs(p);
43 for (size_t k = 0; k < coeffs.size(); ++k)
44 {
45 coeffs[k] = Complex(binomial_coefficient<Real>(p-1+k, k), 0);
46 }
47
48 polynomial<Complex> P(std::move(coeffs));
49 polynomial<Complex> Pcopy = P;
50 polynomial<Complex> Pcopy_prime = P.prime();
51 auto orig = [&](Complex z) { return std::make_pair<Complex, Complex>(Pcopy(z), Pcopy_prime(z)); };
52
53 polynomial<Complex> P_prime = P.prime();
54
55 // Polynomial is of degree p-1.
56
57 std::vector<Complex> roots(p-1, {std::numeric_limits<Real>::quiet_NaN(),std::numeric_limits<Real>::quiet_NaN()});
58 size_t i = 0;
59 while(P.size() > 1)
60 {
61 Complex guess = {0.0, 1.0};
62 std::cout << std::setprecision(std::numeric_limits<Real>::digits10+3);
63
64 auto f = [&](Complex x)->std::pair<Complex, Complex>
65 {
66 return std::make_pair<Complex, Complex>(P(x), P_prime(x));
67 };
68
69 Complex r = complex_newton(f, guess);
70 using std::isnan;
71 if(isnan(r.real()))
72 {
73 int i = 50;
74 do {
75 // Try a different guess
76 guess *= Complex(1.0,-1.0);
77 r = complex_newton(f, guess);
78 std::cout << "New guess: " << guess << ", result? " << r << std::endl;
79
80 } while (isnan(r.real()) && i-- > 0);
81
82 if (isnan(r.real()))
83 {
84 std::cout << "Polynomial that killed the process: " << P << std::endl;
85 throw std::logic_error("Newton iteration did not converge");
86 }
87 }
88 // Refine r with the original function.
89 // We only use the polynomial division to ensure we don't get the same root over and over.
90 // However, the division induces error which can grow quickly-or slowly! See Numerical Recipes, section 9.5.1.
91 r = complex_newton(orig, r);
92 if (isnan(r.real()))
93 {
94 throw std::logic_error("Found a root for the deflated polynomial which is not a root for the original. Indicative of catastrophic numerical error.");
95 }
96 // Test the root:
97 using std::sqrt;
98 Real tol = sqrt(sqrt(std::numeric_limits<Real>::epsilon()));
99 if (norm(Pcopy(r)) > tol)
100 {
101 std::cout << "This is a bad root: P" << r << " = " << Pcopy(r) << std::endl;
102 std::cout << "Reduced polynomial leading to bad root: " << P << std::endl;
103 throw std::logic_error("Donezo.");
104 }
105
106 BOOST_ASSERT(i < roots.size());
107 roots[i] = r;
108 ++i;
109 polynomial<Complex> q{-r, {1,0}};
110 // This optimization breaks at p = 11. I have no clue why.
111 // Unfortunate, because I expect it to be considerably more stable than
112 // repeatedly dividing by the complex root.
113 /*polynomial<Complex> q;
114 if (r.imag() > sqrt(std::numeric_limits<Real>::epsilon()))
115 {
116 // Then the complex conjugate is also a root:
117 using std::conj;
118 using std::norm;
119 BOOST_ASSERT(i < roots.size());
120 roots[i] = conj(r);
121 ++i;
122 q = polynomial<Complex>({{norm(r), 0}, {-2*r.real(),0}, {1,0}});
123 }
124 else
125 {
126 // The imaginary part is numerical noise:
127 r.imag() = 0;
128 q = polynomial<Complex>({-r, {1,0}});
129 }*/
130
131
132 auto PR = quotient_remainder(P, q);
133 // I should validate that the remainder is small, but . . .
134 //std::cout << "Remainder = " << PR.second<< std::endl;
135
136 P = PR.first;
137 P_prime = P.prime();
138 }
139
140 std::vector<std::pair<Complex, Complex>> Qroots(p-1);
141 for (size_t i = 0; i < Qroots.size(); ++i)
142 {
143 Complex y = roots[i];
144 Complex z1 = static_cast<Complex>(1) - static_cast<Complex>(2)*y + static_cast<Complex>(2)*sqrt(y*(y-static_cast<Complex>(1)));
145 Complex z2 = static_cast<Complex>(1) - static_cast<Complex>(2)*y - static_cast<Complex>(2)*sqrt(y*(y-static_cast<Complex>(1)));
146 Qroots[i] = {z1, z2};
147 }
148
149 return Qroots;
150 }
151
152 template<class Complex>
daubechies_coefficients(std::vector<std::pair<Complex,Complex>> const & Qroots)153 std::vector<typename Complex::value_type> daubechies_coefficients(std::vector<std::pair<Complex, Complex>> const & Qroots)
154 {
155 typedef typename Complex::value_type Real;
156 size_t p = Qroots.size() + 1;
157 // Choose the minimum abs root; see Mallat, discussion just after equation 7.98
158 std::vector<Complex> chosen_roots(p-1);
159 for (size_t i = 0; i < p - 1; ++i)
160 {
161 if(norm(Qroots[i].first) <= 1)
162 {
163 chosen_roots[i] = Qroots[i].first;
164 }
165 else
166 {
167 BOOST_ASSERT(norm(Qroots[i].second) <= 1);
168 chosen_roots[i] = Qroots[i].second;
169 }
170 }
171
172 polynomial<Complex> R{1};
173 for (size_t i = 0; i < p-1; ++i)
174 {
175 Complex ak = chosen_roots[i];
176 R *= polynomial<Complex>({-ak/(static_cast<Complex>(1)-ak), static_cast<Complex>(1)/(static_cast<Complex>(1)-ak)});
177 }
178 polynomial<Complex> a{{half<Real>(), 0}, {half<Real>(),0}};
179 polynomial<Complex> poly = root_two<Real>()*pow(a, p)*R;
180 std::vector<Complex> result = poly.data();
181 // If we reverse, we get the Numerical Recipes and Daubechies convention.
182 // If we don't reverse, we get the Pywavelets and Mallat convention.
183 // I believe this is because of the sign convention on the DFT, which differs between Daubechies and Mallat.
184 // You implement a dot product in Daubechies/NR convention, and a convolution in PyWavelets/Mallat convention.
185 // I won't reverse so I can spot check against Pywavelets: http://wavelets.pybytes.com/wavelet/
186 //std::reverse(result.begin(), result.end());
187 std::vector<Real> h(result.size());
188 for (size_t i = 0; i < result.size(); ++i)
189 {
190 Complex r = result[i];
191 BOOST_ASSERT(r.imag() < sqrt(std::numeric_limits<Real>::epsilon()));
192 h[i] = r.real();
193 }
194
195 // Quick sanity check: We could check all vanishing moments, but that sum is horribly ill-conditioned too!
196 Real sum = 0;
197 Real scale = 0;
198 for (size_t i = 0; i < h.size(); ++i)
199 {
200 sum += h[i];
201 scale += h[i]*h[i];
202 }
203 BOOST_ASSERT(abs(scale -1) < sqrt(std::numeric_limits<Real>::epsilon()));
204 BOOST_ASSERT(abs(sum - root_two<Real>()) < sqrt(std::numeric_limits<Real>::epsilon()));
205 return h;
206 }
207
main()208 int main()
209 {
210 typedef boost::multiprecision::cpp_complex<100> Complex;
211 for(size_t p = 1; p < 200; ++p)
212 {
213 auto roots = find_roots<Complex>(p);
214 auto h = daubechies_coefficients(roots);
215 std::cout << "h_" << p << "[] = {";
216 for (auto& x : h) {
217 std::cout << x << ", ";
218 }
219 std::cout << "} // = h_" << p << "\n\n\n\n";
220 }
221 }
222