//+------------------------------------------------------------------+
//|                                                       polevl.mqh |
//|                                  Copyright 2025, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include"config.mqh"

namespace special
{
namespace cephes
{
inline double polevl(double x, double &coef[], int N)
  {
   double ans;
   int i;
   int p;

   p = 0;
   ans = coef[p++];
   i = N;

   do
     {
      ans = ans * x + coef[p++];
     }
   while(--i != 0);

   return (ans);
  }

/*                                                     p1evl() */
/*                                          N
 * Evaluate polynomial when coefficient of x  is 1.0.
 * That is, C_{N} is assumed to be 1, and that coefficient
 * is not included in the input array coef.
 * coef must have length N and contain the polynomial coefficients
 * stored as
 *     coef[0] = C_{N-1}
 *     coef[1] = C_{N-2}
 *          ...
 *     coef[N-2] = C_1
 *     coef[N-1] = C_0
 * Otherwise same as polevl.
 */

inline double p1evl(double x, double &coef[], int N)
  {
   double ans;
   int p;
   int i;

   p = 0;
   ans = x + coef[p++];
   i = N - 1;

   do
      ans = ans * x + coef[p++];
   while(--i != 0);

   return (ans);
  }

/* Evaluate a rational function. See [1]. */

/* The function ratevl is only used once in cephes/lanczos.h. */
inline double ratevl(double x, double &num[], int M, double &denom[], int N)
  {
   int i, dir;
   double y, num_ans, denom_ans;
   double absx = abs(x);
   int p;

   if(absx > 1)
     {
      /* Evaluate as a polynomial in 1/x. */
      dir = -1;
      p = M;
      y = 1 / x;
     }
   else
     {
      dir = 1;
      p = 0;
      y = x;
     }

   /* Evaluate the numerator */
   num_ans = num[p];
   p += dir;
   for(i = 1; i <= M; i++)
     {
      num_ans = num_ans * y + num[p];
      p += dir;
     }

   /* Evaluate the denominator */
   if(absx > 1)
     {
      p = N;
     }
   else
     {
      p = 0;
     }

   denom_ans = denom[p];
   p += dir;
   for(i = 1; i <= N; i++)
     {
      denom_ans = denom_ans * y + denom[p];
      p += dir;
     }

   if(absx > 1)
     {
      i = M - N;
      return pow(x, i) * num_ans / denom_ans;
     }
   else
     {
      return num_ans / denom_ans;
     }
  }
} // namespace cephes
} // namespace special