/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  //  (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)
  
  package org.apache.commons.numbers.gamma;
  
  import java.util.function.DoubleSupplier;
  
  /**
   * Utility tools used by the Boost functions.
   *
   * <p>This code has been adapted from the <a href="https://www.boost.org/">Boost</a>
   * {@code c++} implementations in {@code <boost/math/tools/>}.
   * All work is copyright John Maddock 2006 and subject to the Boost Software License.
   */
  final class BoostTools {
      /**
       * The minimum epsilon value for relative error in the summation.
       * Equal to Math.ulp(1.0) or 2^-52.
       *
       * <h2>Note</h2>
       *
       * <p>The summation will terminate when any additional terms are too small to
       * change the sum. Assuming additional terms are reducing in magnitude this
       * occurs when the term is 1 ULP from the sum:
       * <pre>
       * ulp(sum) >= term
       * </pre>
       *
       * <p>The epsilon is used to set a configurable threshold using:
       * <pre>
       * sum * eps >= term
       * </pre>
       * <p>The minimum epsilon is the smallest value that will create a value
       * {@code >= 1} ULP and {@code < 2} ULP of the sum. For any normal number the ULP
       * of all values with the same exponent b is scalb(1.0, b - 52). This can
       * be achieved by multiplication by 2^-52.
       */
      private static final double EPSILON = 0x1.0p-52;
      /**
       * The minimum epsilon value for relative error in the Kahan summation.
       * This can be lower than {@link #EPSILON}. Set to 2^-62.
       *
       * <h2>Note</h2>
       *
       * <p>The Kahan summation uses a carry term to extend the precision of the sum.
       * This extends the 53-bit mantissa by adding more bits to hold round-off error.
       * Thus the term may effect the sum when it has a magnitude smaller than 1 ULP
       * of the current sum. The epsilon can be lowered from 2^-52 to include the
       * extra bits in the convergence criteria. The lower limit for the epsilon is
       * 2^-104. Boost uses an epsilon specified as a number of bits of accuracy. Each
       * lowering of the epsilon by a factor of 2 adds a guard digit to the sum.
       * Slower converging series will benefit from a lower epsilon. This uses 2^-62
       * to add 10 guard digits and allow testing of different thresholds when the
       * Kahan summation is used, for example in the gamma function. Lowering the
       * epsilon further is only of benefit if the terms can be computed exactly.
       * Otherwise the rounding errors of the early terms affect the final result as
       * much as the inclusion of extra guard digits.
       */
      private static final double KAHAN_EPSILON = 0x1.0p-62;
      /** Message for failure to converge. */
      private static final String MSG_FAILED_TO_CONVERGE = "Failed to converge within %d iterations";
  
      /** Private constructor. */
      private BoostTools() {
          // intentionally empty.
      }
  
      /**
       * Sum the series.
       *
       * <p>Adapted from {@code boost/math/tools/series.hpp}.
       *
       * @param func Series generator
       * @param epsilon Maximum relative error allowed
       * @param maxTerms Maximum number of terms
       * @return result
       */
      static double sumSeries(DoubleSupplier func, double epsilon, int maxTerms) {
          return sumSeries(func, epsilon, maxTerms, 0);
      }
 
     /**
      * Sum the series.
      *
      * <p>Adapted from {@code boost/math/tools/series.hpp}.
      *
      * @param func Series generator
      * @param epsilon Maximum relative error allowed
      * @param maxTerms Maximum number of terms
      * @param initValue Initial value
      * @return result
      */
     static double sumSeries(DoubleSupplier func, double epsilon, int maxTerms, double initValue) {
         // Note:
         // The Boost code requires eps to be non-zero. It is created in the
         // <boost/math/policies/policy.hpp> as a non-zero relative error term.
         // An alternative termination condition with a divide is:
         // (eps < Math.abs(nextTerm / result))
         //
         // Here the argument is checked against the minimum epsilon for a double
         // to provide functional equivalence with the Boost policy.
         // In the min eps case the loop terminates if the most recently added term is
         // 0 or 1 ulp of the result. This condition is acceptable if the next
         // computed term will be at most half of the most recent term (thus
         // cannot be added to the current result).
 
         final double eps = getEpsilon(epsilon, EPSILON);
 
         int counter = maxTerms;
 
         double result = initValue;
         double nextTerm;
         do {
             nextTerm = func.getAsDouble();
             result += nextTerm;
         } while (Math.abs(eps * result) < Math.abs(nextTerm) && --counter > 0);
 
         if (counter <= 0) {
             throw new ArithmeticException(
                String.format(MSG_FAILED_TO_CONVERGE, maxTerms));
         }
 
         return result;
     }
 
     /**
      * Sum the series using Kahan summation.
      *
      * <p>Adapted from {@code boost/math/tools/series.hpp}.
      *
      * @param func Series generator
      * @param epsilon Maximum relative error allowed
      * @param maxTerms Maximum number of terms
      * @return result
      */
     static double kahanSumSeries(DoubleSupplier func, double epsilon, int maxTerms) {
         return kahanSumSeries(func, epsilon, maxTerms, 0);
     }
 
     /**
      * Sum the series using Kahan summation.
      *
      * <p>Adapted from {@code boost/math/tools/series.hpp}.
      *
      * @param func Series generator
      * @param epsilon Maximum relative error allowed
      * @param maxTerms Maximum number of terms
      * @param initValue Initial value
      * @return result
      */
     static double kahanSumSeries(DoubleSupplier func, double epsilon, int maxTerms, double initValue) {
         final double eps = getEpsilon(epsilon, KAHAN_EPSILON);
 
         int counter = maxTerms;
 
         // Kahan summation:
         // https://en.wikipedia.org/wiki/Kahan_summation_algorithm
         // This summation is accurate if the term is smaller in magnitude
         // than the current sum. This is a condition required for the
         // series termination thus the extended precision sum need not
         // check magnitudes of terms to compute the carry.
 
         double result = initValue;
         double carry = 0;
         double nextTerm;
         do {
             nextTerm = func.getAsDouble();
             final double y = nextTerm - carry;
             final double t = result + y;
             carry = t - result;
             carry -= y;
             result = t;
         } while (Math.abs(eps * result) < Math.abs(nextTerm) && --counter > 0);
 
         if (counter <= 0) {
             throw new ArithmeticException(
                String.format(MSG_FAILED_TO_CONVERGE, maxTerms));
         }
 
         return result;
     }
 
     /**
      * Gets the epsilon ensuring it satisfies the minimum allowed value.
      *
      * <p>This is returning the maximum of the two arguments.
      * Do not use Math.max as it returns NaN if either value is NaN.
      * In this case the desired result in the default minEpsilon, not NaN.
      * Math.max will also check ordering when terms are equal to support
      * -0.0 and 0.0. This does not apply here and a single conditional
      * returns the desired result.
      *
      * @param epsilon Configured epsilon
      * @param minEpsilon Minimum allowed epsilon
      * @return the epsilon
      */
     private static double getEpsilon(double epsilon, double minEpsilon) {
         return epsilon > minEpsilon ? epsilon : minEpsilon;
     }
 
     /**
      * Evaluate the polynomial using Horner's method.
      * The coefficients are used in descending order, for example a polynomial of order
      * 3 requires 4 coefficients:
      * <pre>
      * f(x) = c[3] * x^3 + c[2] * x^2 + c[1] * x + c[0]
      * </pre>
      *
      * @param c Polynomial coefficients (must have length > 0)
      * @param x Argument x
      * @return polynomial value
      */
     static double evaluatePolynomial(double[] c, double x) {
         final int count = c.length;
         double sum = c[count - 1];
         for (int i = count - 2; i >= 0; --i) {
             sum *= x;
             sum += c[i];
         }
         return sum;
     }
 }