module Math
  def self.factorial(n)
    return if n < 0

    n = n.to_i # Only integers.

    return 1 if n == 0 || n == 1
    Math.gamma(n + 1) # Math.gamma(x) == (n - 1)! for integer values
  end

  def self.combination(n, r)
    self.factorial(n)/(self.factorial(r) * self.factorial(n - r)).to_r # n!/(r! * [n - r]!)
  end

  def self.permutation(n, k)
    self.factorial(n)/self.factorial(n - k).to_r
  end

  # Function adapted from the python implementation that exists in https://en.wikipedia.org/wiki/Simpson%27s_rule#Sample_implementation
  # Finite integral in the interval [a, b] split up in n-intervals
  def self.simpson_rule(a, b, n, &block)
    unless n.even?
      puts "The composite simpson's rule needs even intervals!"
      return
    end

    h = (b - a)/n.to_r

    resA = yield(a)
    resB = yield(b)

    sum = resA + resB

    (1..n).step(2).each do |number|
      res = yield(a + number * h)
      sum += 4 * res
    end

    (1..(n-1)).step(2).each do |number|
      res = yield(a + number * h)
      sum += 2 * res
    end

    return sum * h / 3.0
  end

  def self.lower_incomplete_gamma_function(s, x)
    base_iterator = x.round(1)
    base_iterator += 1 if x < 1.0 && !x.zero?

    # The greater the iterations, the better. That's why we are iterating 10_000 * x times
    iterator = (10_000 * base_iterator).round
    iterator = 100_000 if iterator.zero?

    self.simpson_rule(0, x.to_r, iterator) do |t|
      (t ** (s - 1)) * Math.exp(-t)
    end
  end

  # Algorithm implementation translated from the ASA147 C++ version https://people.sc.fsu.edu/~jburkardt/cpp_src/asa147/asa147.html
  # translated from FORTRAN by John Burkardt. Original algorithm written by Chi Leung Lau.
  # It contains a modification on the error and underflow parameters to use maximum available float number
  # and it performs the series using `Rational` objects to avoid memory exhaustion and reducing precision errors.
  #
  # This algorithm is licensed with MIT license.
  #
  # Reference:
  #
  #    Chi Leung Lau,
  #    Algorithm AS 147:
  #    A Simple Series for the Incomplete Gamma Integral,
  #    Applied Statistics,
  #    Volume 29, Number 1, 1980, pages 113-114.
  def self.normalised_lower_incomplete_gamma_function(s, x)
    return 0.0 if s.negative? || x.zero? || x.negative?

    # e = 1.0e-308
    # uflo = 1.0e-47
    e = Float::MIN
    uflo = Float::MIN

    lgamma, sign = Math.lgamma(s + 1.0)
    arg = s * Math.log(x) - (sign * lgamma) - x

    return 0.0 if arg < Math.log(uflo)

    f = Math.exp(arg).to_r

    return 0.0 if f.zero?

    c = 1r
    value = 1r
    a = s.to_r
    rational_x = x.to_r

    loop do
      a += 1r
      c = c * (rational_x / a)
      value += c

      break if c <= (e * value)
    end

    (value * f).to_f
  end

  def self.beta_function(x, y)
    return 1 if x == 1 && y == 1

    (Math.gamma(x) * Math.gamma(y))/Math.gamma(x + y)
  end

  ### This implementation is an adaptation of the incomplete beta function made in C by
  ### Lewis Van Winkle, which released the code under the zlib license.
  ### The whole math behind this code is described in the following post: https://codeplea.com/incomplete-beta-function-c
  def self.incomplete_beta_function(x, alp, bet)
    return if x < 0.0
    return 1.0 if x > 1.0

    tiny = 1.0E-50

    if x > ((alp + 1.0)/(alp + bet + 2.0))
      return 1.0 - self.incomplete_beta_function(1.0 - x, bet, alp)
    end

    # To avoid overflow problems, the implementation applies the logarithm properties
    # to calculate in a faster and safer way the values.
    lbet_ab = (Math.lgamma(alp)[0] + Math.lgamma(bet)[0] - Math.lgamma(alp + bet)[0]).freeze
    front = (Math.exp(Math.log(x) * alp + Math.log(1.0 - x) * bet - lbet_ab) / alp.to_r).freeze

    # This is the non-log version of the left part of the formula (before the continuous fraction)
    # down_left = alp * self.beta_function(alp, bet)
    # upper_left = (x ** alp) * ((1.0 - x) ** bet)
    # front = upper_left/down_left

    f, c, d = 1.0, 1.0, 0.0

    returned_value = nil

    # Let's do more iterations than the proposed implementation (200 iters)
    (0..500).each do |number|
      m = number/2

      numerator = if number == 0
                    1.0
                  elsif number % 2 == 0
                    (m * (bet - m) * x)/((alp + 2.0 * m - 1.0)* (alp + 2.0 * m))
                  else
                    top = -((alp + m) * (alp + bet + m) * x)
                    down = ((alp + 2.0 * m) * (alp + 2.0 * m + 1.0))

                    top/down
                  end

      d = 1.0 + numerator * d
      d = tiny if d.abs < tiny
      d = 1.0 / d.to_r

      c = 1.0 + numerator / c.to_r
      c = tiny if c.abs < tiny

      cd = (c*d).freeze
      f = f * cd


      if (1.0 - cd).abs < 1.0E-10
        returned_value = front * (f - 1.0)
        break
      end
    end

    returned_value
  end
end