# include <math.h>
# include <stdio.h>
# include <stdlib.h>

# include "sigmoid_derivative.h"

/******************************************************************************/

int nchoosek ( int n, int k )

/******************************************************************************/
/*
  Purpose:

    nchoosek() computes the binomial coefficient C(N,K).

  Discussion:

    The value is calculated in such a way as to avoid overflow and
    roundoff.  The calculation is done in int arithmetic.

    The formula used is:

      C(N,K) = N! / ( K! * (N-K)! )

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    09 December 2013

  Author:

    John Burkardt

  Reference:

    ML Wolfson, HV Wright,
    Algorithm 160:
    Combinatorial of M Things Taken N at a Time,
    Communications of the ACM,
    Volume 6, Number 4, April 1963, page 161.

  Input:

    int N, K, are the values of N and K.

  Output:

    int nchoosek: the number of combinations of N
    things taken K at a time.
*/
{
  int i;
  int mn;
  int mx;
  int value;

  if ( k < n - k )
  {
    mn = k;
    mx = n - k;
  }
  else
  {
    mn = n - k;
    mx = k;
  }

  if ( mn < 0 )
  {
    value = 0;
  }
  else if ( mn == 0 )
  {
    value = 1;
  }
  else
  {
    value = mx + 1;

    for ( i = 2; i <= mn; i++ )
    {
      value = ( value * ( mx + i ) ) / i;
    }
  }

  return value;
}
/******************************************************************************/

double *sigmoid_derivative_coef ( int n )

/******************************************************************************/
/*
  Purpose:

    sigmoid_derivative_coef(): sigmoid derivative expansion coefficients.

  Discussion:

    s(x) = 1 / ( 1 + exp ( - x ) )

    s(x)^(n) = sum ( 1 <= j <= n + 1 ) coef(j) * s(x)^j 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 May 2024

  Author:

    John Burkardt

  Reference:

    Joe McKenna,
    Derivatives of the sigmoid function,
    https://joepatmckenna.github.io/calculus/derivative/sigmoid!20function/linear!20albegra/2018/01/20/sigmoid-derivs/

  Input:

    int n: the order of the derivative.
    0 <= n.

  Output:

    double coef[n+2]: the expansion coefficients.
    coef[0] is always 0, and coef[1] is always 1.
*/
{
  double cnk;
  double *coef;
  int j;
  int k;
  int mop;
  
  coef = ( double * ) malloc ( ( n + 2 ) * sizeof ( double ) );

  for ( k = 0; k < n + 2; k++ )
  {
    coef[k] = 0.0;
  }

  for ( k = 0; k <= n; k++ )
  {
    cnk = 0.0;
    mop = -1;
    for ( j = 0; j <= k; j++ )
    {
      mop = - mop;
      cnk = cnk + mop * pow ( j + 1, n ) * nchoosek ( k, j );
    }
    coef[k+1] = cnk;
  }

  return coef;
}
/******************************************************************************/

double sigmoid_derivative_value ( int n, double x )

/******************************************************************************/
/*
  Purpose:

    sigmoid_derivative_value(): evaluate sigmoid derivative at x.

  Discussion:

    s(x) = 1 / ( 1 + exp ( - x ) )

    s(x)^(n) = sum ( 1 <= j <= n + 1 ) coef(j) * s(x)^j 

  Licensing:

    This code is distributed under the MIT license.

  Modified:

    23 May 2024

  Author:

    John Burkardt

  Reference:

    Joe McKenna,
    Derivatives of the sigmoid function,
    https://joepatmckenna.github.io/calculus/derivative/sigmoid!20function/linear!20albegra/2018/01/20/sigmoid-derivs/

  Input:

    int n: the order of the derivative.

    double x: the evaluation point.

  Output:

    double sigmoid_derivative_value: the value of s^(n)(x).
*/
{
  double cnk;
  int j;
  int k;
  int mop;
  double sig;
  double sigk;
  double value;

  sig = 1.0 / ( 1.0 + exp ( - x ) );
  value = 0.0;
  sigk = 1.0;

  for ( k = 0; k <= n; k++ )
  {
    sigk = sigk * sig;
    cnk = 0.0;
    mop = -1;
    for ( j = 0; j <= k; j++ )
    {
      mop = - mop;
      cnk = cnk + mop * pow ( j + 1, n ) * nchoosek ( k, j );
    }
    value = value + cnk * sigk;
  }

  return value;
}