mirror of
https://github.com/saitohirga/WSJT-X.git
synced 2024-08-24 23:26:36 -04:00
207 lines
5.1 KiB
C++
207 lines
5.1 KiB
C++
|
// Copyright John Maddock 2006.
|
||
|
// Use, modification and distribution are subject to the
|
||
|
// Boost Software License, Version 1.0. (See accompanying file
|
||
|
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
|
||
|
|
||
|
#include <boost/math/special_functions/log1p.hpp>
|
||
|
#include <boost/math/special_functions/erf.hpp>
|
||
|
#include <boost/math/constants/constants.hpp>
|
||
|
#include <map>
|
||
|
#include <iostream>
|
||
|
#include <iomanip>
|
||
|
#include "mp_t.hpp"
|
||
|
|
||
|
using namespace std;
|
||
|
using namespace boost::math;
|
||
|
|
||
|
//
|
||
|
// This program calculates the coefficients of the polynomials
|
||
|
// used for the regularized incomplete gamma functions gamma_p
|
||
|
// and gamma_q when parameter a is large, and sigma is small
|
||
|
// (where sigma = fabs(1 - x/a) ).
|
||
|
//
|
||
|
// See "The Asymptotic Expansion of the Incomplete Gamma Functions"
|
||
|
// N. M. Temme.
|
||
|
// Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
|
||
|
// Coeffient calculation is described from Eq 3.8 (p762) onwards.
|
||
|
//
|
||
|
|
||
|
//
|
||
|
// Alpha:
|
||
|
//
|
||
|
mp_t alpha(unsigned k)
|
||
|
{
|
||
|
static map<unsigned, mp_t> data;
|
||
|
if(data.empty())
|
||
|
{
|
||
|
data[1] = 1;
|
||
|
}
|
||
|
|
||
|
map<unsigned, mp_t>::const_iterator pos = data.find(k);
|
||
|
if(pos != data.end())
|
||
|
return (*pos).second;
|
||
|
//
|
||
|
// OK try and calculate the value:
|
||
|
//
|
||
|
mp_t result = alpha(k-1);
|
||
|
for(unsigned j = 2; j <= k-1; ++j)
|
||
|
{
|
||
|
result -= j * alpha(j) * alpha(k-j+1);
|
||
|
}
|
||
|
result /= (k+1);
|
||
|
data[k] = result;
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
mp_t gamma(unsigned k)
|
||
|
{
|
||
|
static map<unsigned, mp_t> data;
|
||
|
|
||
|
map<unsigned, mp_t>::const_iterator pos = data.find(k);
|
||
|
if(pos != data.end())
|
||
|
return (*pos).second;
|
||
|
|
||
|
mp_t result = (k&1) ? -1 : 1;
|
||
|
|
||
|
for(unsigned i = 1; i <= (2 * k + 1); i += 2)
|
||
|
result *= i;
|
||
|
result *= alpha(2 * k + 1);
|
||
|
data[k] = result;
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
mp_t Coeff(unsigned n, unsigned k)
|
||
|
{
|
||
|
map<unsigned, map<unsigned, mp_t> > data;
|
||
|
if(data.empty())
|
||
|
data[0][0] = mp_t(-1) / 3;
|
||
|
|
||
|
map<unsigned, map<unsigned, mp_t> >::const_iterator p1 = data.find(n);
|
||
|
if(p1 != data.end())
|
||
|
{
|
||
|
map<unsigned, mp_t>::const_iterator p2 = p1->second.find(k);
|
||
|
if(p2 != p1->second.end())
|
||
|
{
|
||
|
return p2->second;
|
||
|
}
|
||
|
}
|
||
|
|
||
|
//
|
||
|
// If we don't have the value, calculate it:
|
||
|
//
|
||
|
if(k == 0)
|
||
|
{
|
||
|
// special case:
|
||
|
mp_t result = (n+2) * alpha(n+2);
|
||
|
data[n][k] = result;
|
||
|
return result;
|
||
|
}
|
||
|
// general case:
|
||
|
mp_t result = gamma(k) * Coeff(n, 0) + (n+2) * Coeff(n+2, k-1);
|
||
|
data[n][k] = result;
|
||
|
return result;
|
||
|
}
|
||
|
|
||
|
void calculate_terms(double sigma, double a, unsigned bits)
|
||
|
{
|
||
|
cout << endl << endl;
|
||
|
cout << "Sigma: " << sigma << endl;
|
||
|
cout << "A: " << a << endl;
|
||
|
double lambda = 1 - sigma;
|
||
|
cout << "Lambda: " << lambda << endl;
|
||
|
double y = a * (-sigma - log1p(-sigma));
|
||
|
cout << "Y: " << y << endl;
|
||
|
double z = -sqrt(2 * (-sigma - log1p(-sigma)));
|
||
|
cout << "Z: " << z << endl;
|
||
|
double dom = erfc(sqrt(y)) / 2;
|
||
|
cout << "Erfc term: " << dom << endl;
|
||
|
double lead = exp(-y) / sqrt(2 * constants::pi<double>() * a);
|
||
|
cout << "Remainder factor: " << lead << endl;
|
||
|
double eps = ldexp(1.0, 1 - static_cast<int>(bits));
|
||
|
double target = dom * eps / lead;
|
||
|
cout << "Target smallest term: " << target << endl;
|
||
|
|
||
|
unsigned max_n = 0;
|
||
|
|
||
|
for(unsigned n = 0; n < 10000; ++n)
|
||
|
{
|
||
|
double term = tools::real_cast<double>(Coeff(n, 0) * pow(z, (double)n));
|
||
|
if(fabs(term) < target)
|
||
|
{
|
||
|
max_n = n-1;
|
||
|
break;
|
||
|
}
|
||
|
}
|
||
|
cout << "Max n required: " << max_n << endl;
|
||
|
|
||
|
unsigned max_k;
|
||
|
for(unsigned k = 1; k < 10000; ++k)
|
||
|
{
|
||
|
double term = tools::real_cast<double>(Coeff(0, k) * pow(a, -((double)k)));
|
||
|
if(fabs(term) < target)
|
||
|
{
|
||
|
max_k = k-1;
|
||
|
break;
|
||
|
}
|
||
|
}
|
||
|
cout << "Max k required: " << max_k << endl << endl;
|
||
|
|
||
|
bool code = false;
|
||
|
cout << "Print code [0|1]? ";
|
||
|
cin >> code;
|
||
|
|
||
|
int prec = 2 + (static_cast<double>(bits) * 3010LL)/10000;
|
||
|
std::cout << std::scientific << std::setprecision(40);
|
||
|
|
||
|
if(code)
|
||
|
{
|
||
|
cout << " T workspace[" << max_k+1 << "];\n\n";
|
||
|
for(unsigned k = 0; k <= max_k; ++k)
|
||
|
{
|
||
|
cout <<
|
||
|
" static const T C" << k << "[] = {\n";
|
||
|
for(unsigned n = 0; n < 10000; ++n)
|
||
|
{
|
||
|
double term = tools::real_cast<double>(Coeff(n, k) * pow(a, -((double)k)) * pow(z, (double)n));
|
||
|
if(fabs(term) < target)
|
||
|
{
|
||
|
break;
|
||
|
}
|
||
|
cout << " " << Coeff(n, k) << "L,\n";
|
||
|
}
|
||
|
cout <<
|
||
|
" };\n"
|
||
|
" workspace[" << k << "] = tools::evaluate_polynomial(C" << k << ", z);\n\n";
|
||
|
}
|
||
|
cout << " T result = tools::evaluate_polynomial(workspace, 1/a);\n\n";
|
||
|
}
|
||
|
}
|
||
|
|
||
|
|
||
|
int main()
|
||
|
{
|
||
|
bool cont;
|
||
|
do{
|
||
|
cont = false;
|
||
|
double sigma;
|
||
|
cout << "Enter max value for sigma (sigma = |1 - x/a|): ";
|
||
|
cin >> sigma;
|
||
|
double a;
|
||
|
cout << "Enter min value for a: ";
|
||
|
cin >> a;
|
||
|
unsigned precision;
|
||
|
cout << "Enter number of bits precision required: ";
|
||
|
cin >> precision;
|
||
|
|
||
|
calculate_terms(sigma, a, precision);
|
||
|
|
||
|
cout << "Try again[0|1]: ";
|
||
|
cin >> cont;
|
||
|
|
||
|
}while(cont);
|
||
|
|
||
|
|
||
|
return 0;
|
||
|
}
|
||
|
|