coe_db.cxx

#include <algorithm>
#include <cmath>
#include <complex>
#include <iostream>
#include <ios>
#include <iomanip>
#include <limits>
#include <numeric>
#include <vector>
#include <boost/math/tools/polynomial.hpp>
#include <boost/multiprecision/gmp.hpp>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/mpc.hpp>
#define EIGEN_NO_DEBUG
#define EIGEN_NO_STATIC_ASSERT
#include <Eigen/Dense>
namespace mt = boost::math::tools;
namespace mp = boost::multiprecision;
namespace e  = Eigen;
constexpr size_t MD10f  = std::numeric_limits<     float >::max_digits10;
constexpr size_t MD10d  = std::numeric_limits<     double>::max_digits10;
constexpr size_t MD10ld = std::numeric_limits<long double>::max_digits10;
constexpr size_t MD10   = MD10f * 16; // NOTE: For exact a[.] in cxx_float for N <= 38.
typedef mp::mpz_int cxx_int;
typedef mp::mpq_rational cxx_rational;
typedef mp::number<mp::mpfr_float_backend<MD10>> cxx_float;
typedef mp::number<mp::mpc_complex_backend<MD10>> cxx_complex;
int main(int argc, char *argv[]) {
	constexpr size_t Nmin =  1;
	constexpr size_t Nmax = 38;
	for (size_t N = Nmin; N <= Nmax; N++) {
		std::clog << "N = " << N << std::endl;
		const size_t Nsup = 2 * N;
		std::vector<cxx_int> n(N);      std::iota(n.begin(), n.end(), 1);
		std::vector<cxx_int> sup(Nsup); std::iota(sup.begin(), sup.end(), -static_cast<ssize_t>(N - 1));
		std::vector<cxx_rational> a(N); std::fill(a.begin(), a.end(), 1);
		for (size_t i = 0; i < N; i++) {
			for (size_t j = 0; j < Nsup; j++) {
				if (sup[j] != n[i]) {
					a[i] *= 1 - sup[j] * 2; a[i] /= 2;
					a[i] /= n[i] - sup[j];
				}
			}
			std::clog << "a[" << i << "]"
			          << " = " << a[i]
			          << " = " << std::scientific << std::setprecision(MD10) << a[i].convert_to<cxx_float>() << std::endl;
		}
		const size_t Np = 4 * N - 1;
		std::vector<cxx_rational> p(Np); std::fill(p.begin(), p.end(), 0);
		for (size_t i = 0; i < N; i++) {
			p[2 * i] = a[N - 1 - i];
			p[2 * i + 2 * N] = a[i];
		}
		p[2 * N - 1] = 1;
		for (size_t i = 0; i < Np; i++) {
			std::clog << "p[" << i << "]"
			          << " = " << p[i]
			          << " = " << std::scientific << std::setprecision(MD10) << p[i].convert_to<cxx_float>() << std::endl;
		}
		// FIXME: Following implementation using eigen3 doesn't produce exact eigen values
		//        b/c maybe internal e::NumTraits<...>::dummy_precision() and wrong cmath functions are used?
		//        So just use eigen3 to guess r.
		const size_t NA = Np - 1;
		e::Matrix<cxx_float, e::Dynamic, e::Dynamic> A = e::Matrix<cxx_float, e::Dynamic, e::Dynamic>::Zero(NA, NA);
		for (size_t i = 0; i < NA - 1; i++) {
			A(i + 1, i) = 1;
		}
		for (size_t i = 0; i < NA; i++) {
			const cxx_rational pi = -p[i + 1] / p[0];
			A(0, i) = pi.convert_to<cxx_float>();
		}
		e::EigenSolver<e::Matrix<cxx_float, e::Dynamic, e::Dynamic>> es(A, false);
		const e::Matrix<std::complex<cxx_float>, 1, e::Dynamic> ev = es.eigenvalues(); // NOTE: Do not use cxx_complex, b/c eigen3 expects std::complex<cxx_float> for ev.
		std::clog //<< "A = "                                                                   << std::endl
		          //<< std::scientific << std::setprecision(MD10) << A                          << std::endl
		          //<< "A.diagonal(-1).transpose() = "                                          << std::endl
		          //<< std::scientific << std::setprecision(MD10) << A.diagonal(-1).transpose() << std::endl
		            << "A.row(0).transpose() = "                                                << std::endl
		            << std::scientific << std::setprecision(MD10) << A.row(0).transpose()       << std::endl
		            << "ev.transpose() = "                                                      << std::endl
		            << std::scientific << std::setprecision(MD10) << ev.transpose()             << std::endl;
		std::vector<cxx_complex> r(0);
		for (size_t i = 0; i < NA; i++) {
			const cxx_complex evi(ev(0, i).real(), ev(0, i).imag()); // NOTE: Do not use ev directly, b/c eigen3 expects std::complex<cxx_float> for ev.
			if (mp::abs(evi) < 1 && evi.real() > 0 && evi.imag() >= 0) { // FIXME: Consider error in imaginary part.
				r.push_back(evi);
			}
		}
		sort(
			r.begin(), r.end(),
			[] (const auto &l, const auto &r) { return (mp::abs(l) > mp::abs(r) || mp::abs(mp::arg(l)) > mp::abs(mp::arg(r))); } // NOTE: MATLAB/Octave sort(..., 'descend') compatible.
		);
		const size_t nr = r.size();
		for (size_t i = 0; i < nr; i++) {
			std::clog << "r[" << i << "] = " << std::scientific << std::setprecision(MD10) << r[i] << std::endl;
		}
		assert(nr == (N - 1) / 2 + (N - 1) % 2);
		// NOTE: Polish r by Newton-Raphson method.
		// FIXME: Or, tune eigen3 that can produces in desired precision?
		//        Or, use Hirano method that can produce r in ascending order?
		const size_t nb = static_cast<size_t>(std::ceil(std::log2(          // FIXME: Be constexpr?
			std::numeric_limits<cxx_float>::digits /
			floor(-std::log2(e::NumTraits<long double>::dummy_precision())) // NOTE: Estimated bit accuracy of initial r.
		)));
		const size_t Npd = Np - 1;
		for (size_t i = 0; i < nr; i++) {
			for (size_t j = 0; j < nb; j++) {
				// FIXME: Use std::accumulate()?
				cxx_complex pv = p[0].convert_to<cxx_float>(); // FIXME: ceil(log10(Np)) + 2 * MD10 digits10 capable?
				for (size_t k = 1; k < Np; k++) {
					pv = pv * r[i] + p[k].convert_to<cxx_float>(); // FIXME: Use FMA?
				}
				const cxx_rational pd0 = Npd * p[0];
				cxx_complex pdv = pd0.convert_to<cxx_float>(); // FIXME: ceil(log10(Np - 1)) + 2 * MD10 digits10 capable?
				for (size_t k = 1; k < Npd; k++) {
					const cxx_rational pdk = (Npd - k) * p[k];
					pdv = pdv * r[i] + pdk.convert_to<cxx_float>();
				}
				r[i] -= pv / pdv;
				std::clog << "abs(polyval(p,  r[" << i << "])) = " << std::scientific << std::setprecision(MD10) << mp::abs(pv ) << std::endl
				          << "abs(polyval(pd, r[" << i << "])) = " << std::scientific << std::setprecision(MD10) << mp::abs(pdv) << std::endl
				          << "r[" << i << "] <- "                  << std::scientific << std::setprecision(MD10) << r[i]         << std::endl;
			}
		}
		// NOTE: Boost polynomial storage is reversed order of MATLAB/Octave's.
		const std::vector<cxx_float> cc = { 1 };
		      mt::polynomial<cxx_float> c(cc.rbegin(), cc.rend());
		const std::vector<cxx_float> cc1 = { 1, 1 };
		const mt::polynomial<cxx_float> c1(cc1.rbegin(), cc1.rend());
		for (size_t i = 0; i < N; i++) {
			c *= c1;
		}
		for (size_t i = 0; i < nr; i++) {
			if (imag(r[i]) == 0) {
				const std::vector<cxx_float> ccz = { 1, -r[i].real() };
				const mt::polynomial<cxx_float> cz(ccz.rbegin(), ccz.rend());
				c *= cz;
			} else {
				const std::vector<cxx_float> ccc = { 1, -2 * r[i].real(), mp::norm(r[i]) };
				const mt::polynomial<cxx_float> cc(ccc.rbegin(), ccc.rend());
				c *= cc;
			}
		}
		const size_t nc = c.size();
		for (size_t i = 0; i < nc; i++) {
			std::clog << "c[" << i << "] = " << std::scientific << std::setprecision(MD10) << c[i] << std::endl;
		}
		const cxx_float rssc = mp::sqrt(std::accumulate(
			c.data().begin(), c.data().end(), cxx_float(0), // FIXME: ceil(log10(nc)) + 2 * MD10 digits10 capable?
			//[] (const auto &acc, const auto &arg) { return acc + mp::pow(arg, 2);  } // FIXME: mp::pow(.., int) is missing?
			//[] (const auto &acc, const auto &arg) { return acc + arg * arg;        }
			  [] (const auto &acc, const auto &arg) { return mp::fma(arg, arg, acc); }
		));
		const size_t ng = nc;
		std::vector<cxx_float> g(ng);
		std::transform(
			c.data().rbegin(), c.data().rend(), g.begin(),
			[&rssc] (const auto &arg) { return arg / rssc; }
		);
		for (size_t i = 0; i < ng; i++) {
			std::clog << "g[" << i << "] = " << std::scientific << std::setprecision(MD10) << g[i] << std::endl;
		}
		std::cout << "case " << N << std::endl;
			std::cout << "\tg = [" << std::endl;
				for (size_t i = 0; i < ng; i++) {
					std::cout << "\t\t" << ((g[i] >= 0) ? " " : "")
					          << std::scientific << std::setprecision(MD10d) << g[i]
					          << ((i == ng - 1) ? "" : " ...") << std::endl;
				}
			std::cout << "\t];" << std::endl;
		std::vector<size_t> dist(ng); std::iota(dist.begin(), dist.end(), 0);
		std::vector<cxx_float> ccsl(ng); std::transform(g.begin(),  g.end(),  ccsl.begin(),                 [] (const auto &arg) { return arg; });
		std::vector<cxx_float> ccal(ng); std::transform(g.rbegin(), g.rend(), ccal.begin(),                 [] (const auto &arg) { return arg; });
		std::vector<cxx_float> ccsh(ng); std::transform(g.rbegin(), g.rend(), dist.rbegin(), ccsh.begin(),  [] (const auto &arg1, const auto &arg2) { return (arg2 % 2 == 0) ? arg1 : -arg1; });
		std::vector<cxx_float> ccah(ng); std::transform(g.rbegin(), g.rend(), dist.rbegin(), ccah.rbegin(), [] (const auto &arg1, const auto &arg2) { return (arg2 % 2 == 0) ? arg1 : -arg1; });
		const mt::polynomial<cxx_float> csl(ccsl.rbegin(), ccsl.rend());
		const mt::polynomial<cxx_float> cal(ccal.rbegin(), ccal.rend());
		const mt::polynomial<cxx_float> csh(ccsh.rbegin(), ccsh.rend());
		const mt::polynomial<cxx_float> cah(ccah.rbegin(), ccah.rend());
		// NOTE: Boost polynomial reduces leading zero coefficients.
		std::vector<cxx_float> cc2ref(2 * ng - 1); std::fill(cc2ref.begin(), cc2ref.end(), 0); cc2ref[ng - 1] = 2;
		const mt::polynomial<cxx_float> c2ref(cc2ref.rbegin(), cc2ref.rend()); const size_t nc2ref = c2ref.size(); assert(nc2ref <= 2 * ng - 1);
		const mt::polynomial<cxx_float> c2 = csl * cal + csh * cah;            const size_t nc2    = c2.size();    assert(nc2    <= 2 * ng - 1);
		const mt::polynomial<cxx_float> e2 = c2 - c2ref;                       const size_t ne2    = e2.size();    assert(ne2    <= std::max(nc2ref, nc2));
		for (size_t i = 0; i < ne2; i++) {
			std::clog << "e2[" << i << "] = " << std::scientific << std::setprecision(MD10) << e2[i] << std::endl;
		}
		const cxx_float b2 = -mp::log2(mp::abs(*max_element(
			e2.data().begin(), e2.data().end(),
			[] (const auto &l, const auto &r) { return (mp::abs(l) < mp::abs(r)); }
		)));
		std::clog << "b2 = " << std::fixed << std::setprecision(1) << b2 << std::endl;
	}
	return 0;
}