Line data    Source code

       1             : #include <cmath>
       2             : #include <limits>
       3             : #include <string>
       4             : #include <vector>
       5             : 
       6             : extern "C"
       7             : {
       8             : #include "dist.h"
       9             : }
      10             : #include "data.hpp"
      11             : #include "exceptions.hpp"
      12             : #include "matrix.hpp"
      13             : #include "qrfactorization.hpp"
      14             : #include "vectorf.hpp"
      15             : 
      16           3 : void constructLinearLeastSquaresSystem(const std::vector<double> &coeff,
      17             :                                        const int kstart, matrix<double> &W,
      18             :                                        vectorf<double> &b)
      19             : {
      20             :   // Storage required to construct Top-Line system
      21           3 :   if ((int)coeff.size() - kstart > (int)W.get_rows())
      22             :   {
      23             :     std::string message{
      24           2 :         "Insufficient storage to construct Top-Line system. Have [" +
      25           4 :         std::to_string(W.get_rows()) + "] Need [" +
      26           3 :         std::to_string((int)coeff.size() - kstart) + "]\n"};
      27           1 :     throw sayMessage(message);
      28             :   }
      29             : 
      30          47 :   for (int k{kstart}; k < (int)coeff.size(); k++)
      31             :   {
      32          45 :     int row = k - kstart + 1;
      33          45 :     W(row, 1) = 1.0;
      34          45 :     W(row, 2) = (double)k;
      35          45 :     b(row) = log10(fabs(coeff[ k ]));
      36             :   }
      37           2 : }
      38             : 
      39             : // To fit the tail of the Taylor series, ignore the first kstart terms
      40             : // of the TCs.
      41             : //
      42             : // If Taylor coefficients TC(i) are unscaled, then the scaled TCs are
      43             : // T(i) = (1/scale)^i TC(i).
      44             : //
      45             : // Multiply Rc by scale to get Rc when computing with scaled TCs.
      46             : //
      47             : // See this from ratio test as Rc is 1/L where
      48             : //
      49             : //   L = (1/scale)*limit of abs(TC(i+1)/TC(i)) = limit of abs(T(i+1)/T(i)).
      50             : //
      51             : // To justify the formula for Rc, consider the problem of fitting m and b
      52             : //
      53             : // log10(abs(T(i))) = y(i) = b + m*i = [1 i](b,m)
      54             : //
      55             : // by linear least squares. Then
      56             : //
      57             : // Log10(L) = Log10(abs(T(i+1))/abs(T(i))) =
      58             : // Log10(abs(T(i+1)))-Log10(abs(T(i)))
      59             : //          = m*(i+1)+b - (m*i+b) = m
      60             : //
      61             : // It follows that
      62             : //
      63             : // rc = 1/L = scale/pow(10, m).
      64           2 : double topline(const std::vector<double> &coeff, const double &scale,
      65             :                double &rc, double &order)
      66             : {
      67           2 :   int m{(int)coeff.size() - TOPLINE_KSTART}, n{2};
      68           4 :   matrix<double> W(m, n);
      69           4 :   vectorf<double> beta(n);
      70           4 :   vectorf<double> b(m);
      71             : 
      72             :   // Less than TOPLINE_NUSE coefficients to estimate Rc
      73           2 :   if (m < TOPLINE_NUSE)
      74             :   {
      75             :     std::string message{
      76           2 :         "Need at least [" + std::to_string(TOPLINE_NUSE) +
      77           2 :         "] Taylor coefficients to estimate Radius of Convergence\n"};
      78           1 :     throw sayMessage(message);
      79             :   }
      80             : 
      81             :   // Construct the least squares linear system
      82           1 :   constructLinearLeastSquaresSystem(coeff, TOPLINE_KSTART, W, b);
      83             : 
      84             :   // Save a copy of linear system for computing residual
      85           2 :   matrix<double> WSaved{W};
      86           1 :   vectorf<double> bSaved{b};
      87             : 
      88             :   // Solve W beta = b for beta
      89           1 :   qr(m, n, W, b, beta);
      90             : 
      91           1 :   rc = scale / pow(10, beta(2));
      92           1 :   order = std::numeric_limits<double>::quiet_NaN();
      93             : 
      94             :   // Compute and return 2-norm of residuals
      95          21 :   for (int i{1}; i <= m; i++)
      96          20 :     bSaved(i) -= WSaved(i, 1) * beta(1) + WSaved(i, 2) * beta(2);
      97           2 :   return dnrm2(m, &bSaved(1), 1);
      98             : }