Search code examples
cnumerical-methods

Golden Section Routine Segmentation Fault


I'm trying to find minimum point of Gamma function by Golden Section method. But when I execute the program I get segmentation fault error. I think since I'm a newbie C user, the problem may be due to calling the function Min_Search_Golden_Section wrong. Here is my complete code. I can't find my mistake. Thanks in advance.

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

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define A 12
#define sqrt5 2.236067977499789696

static int Stopping_Rule(double x0, double x1, double tolerance);


double sp_gamma(double z)
{
  const int a = A;
  static double c_space[A];
  static double *c = NULL;
  int k;
  double accm;

  if ( c == NULL ) {
    double k1_factrl = 1.0; /* (k - 1)!*(-1)^k with 0!==1*/
    c = c_space;
    c[0] = sqrt(2.0*M_PI);
    for(k=1; k < a; k++) {
      c[k] = exp(a-k) * pow(a-k, k-0.5) / k1_factrl;
      k1_factrl *= -k;
    }
  }
  accm = c[0];
  for(k=1; k < a; k++) {
    accm += c[k] / ( z + k );
  }
  accm *= exp(-(z+a)) * pow(z+a, z+0.5); /* Gamma(z+1) */
  return accm/z;
}

void Min_Search_Golden_Section( double (*f)(double), double* a, double *fa,
                                     double* b, double* fb, double tolerance)
{
   static const double lambda = 0.5 * (sqrt5 - 1.0);
   static const double mu = 0.5 * (3.0 - sqrt5);         // = 1 - lambda
   double x1;
   double x2;
   double fx1;
   double fx2;


                // Find first two internal points and evaluate 
                // the function at the two internal points.

   x1 = *b - lambda * (*b - *a);                            
   x2 = *a + lambda * (*b - *a);                         
   fx1 = f(x1);
   fx2 = f(x2);

             // Verify that the tolerance is an acceptable number

   if (tolerance <= 0.0) tolerance = sqrt(DBL_EPSILON) * (*b - *a);

           // Loop by exluding segments from current endpoints a, b
           // to current internal points x1, x2 and then calculating
           // a new internal point until the length of the interval
           // is less than or equal to the tolerance.

   while ( ! Stopping_Rule( *a, *b, tolerance) ) {
      if (fx1 > fx2) {
         *a = x1;
         *fa = fx1;
         if ( Stopping_Rule( *a, *b, tolerance) ) break;
         x1 = x2;
         fx1 = fx2;
         x2 = *b - mu * (*b - *a);
         fx2 = f(x2);
      } else {
         *b = x2;
         *fb = fx2;
         if ( Stopping_Rule( *a, *b, tolerance) ) break;
         x2 = x1;
         fx2 = fx1;
         x1 = *a + mu * (*b - *a);
         fx1 = f(x1);
      }
   }
   return;
}


int main()
{
    double x;
    double a = 0.0, b = 4.0, fa = 0.00001, fb = 6.0;
    double fx = sp_gamma(x);
    Min_Search_Golden_Section( &fx, &a, &fa, &b, &fb, 0.0000001);
    return 0;
}


static int Stopping_Rule(double x0, double x1, double tolerance) 
{
   double xm = 0.5 * fabs( x1 + x0 );

   if ( xm <= 1.0 ) return ( fabs( x1 - x0 ) < tolerance ) ? 1 : 0;
   return ( fabs( x1 - x0 ) < tolerance * xm ) ? 1 : 0;
}

Solution

  • You should be getting a compiler error. The first argument to Min_Search_Golden_Section should be a function pointer, but you pass the address of a variable instead.

    When you get compiler errors, fix them - don't run the program and hope. :)

    I guess you just meant to write:

    Min_Search_Golden_Section( &sp_gamma, &a, &fa, &b, &fb, 0.0000001);