#include <iostream>
#include <iomanip>
using namespace std;

#include "MpIeee.hh"
#include "CMpIeee.hh"
#include "ArithmosIO.hh"

/* deriv/test.c
 * 
 * Copyright (C) 2000 David Morrison
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_deriv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

MpIeee f1(MpIeee x, void *params)
{
  return exp (x);
}

MpIeee df1(MpIeee x, void *params)
{
  return exp (x);
}

MpIeee f2(MpIeee x, void *params)
{
  if (x >= MpIeee( "0.0" ))
    {
      return x * sqrt (x);
    }
  else
    {
      return MpIeee( "0.0" );
    }
}

MpIeee df2(MpIeee x, void *params)
{
  if (x >= MpIeee( "0.0" ))
    {
      return MpIeee( "1.5" ) * sqrt (x);
    }
  else
    {
      return MpIeee( "0.0" );
    }
}

MpIeee f3(MpIeee x, void *params)
{
  if (x != MpIeee( "0.0" ))
    {
      return sin (MpIeee( "1" ) / x);
    }
  else
    {
      return MpIeee( "0.0" );
    }
}

MpIeee df3(MpIeee x, void *params)
{
  if (x != MpIeee( "0.0" ))
    {
      return -cos (MpIeee( "1" ) / x) / (x * x);
    }
  else
    {
      return MpIeee( "0.0" );
    }
}

MpIeee f4(MpIeee x, void *params)
{
  return exp (-x * x);
}

MpIeee df4(MpIeee x, void *params)
{
  return -MpIeee( "2.0" ) * x * exp (-x * x);
}

MpIeee f5(MpIeee x, void *params)
{
  return x * x;
}

MpIeee df5(MpIeee x, void *params)
{
  return MpIeee( "2.0" ) * x;
}

MpIeee f6(MpIeee x, void *params)
{
  return MpIeee( "1.0" ) / x;
}

MpIeee df6(MpIeee x, void *params)
{
  return -MpIeee( "1.0" ) / (x * x);
}

typedef int (deriv_fn) (const gsl_function * f, MpIeee x, MpIeee h, MpIeee * res, MpIeee *abserr);

void
test (deriv_fn * deriv, gsl_function * f, gsl_function * df, MpIeee x, 
      const char * desc)
{
  MpIeee result;MpIeee  abserr;
  MpIeee expected=  GSL_FN_EVAL (df, x);
  (*deriv) (f, x, 1e-4, &result, &abserr);

  gsl_test_abs (result, expected, GSL_MIN(1e-4,fabs(expected)) + GSL_DBL_EPSILON, desc);

  if (abserr < fabs(result-expected)) 
    {
      gsl_test_factor (abserr, fabs(result-expected), 2, "%s error estimate", desc);
    }
  else if (result == expected || expected == MpIeee( "0.0" ))
    {
      gsl_test_abs (abserr, 0.0, 1e-6, "%s abserr", desc);
    }
  else
    {
      MpIeee d=  fabs(result - expected);
      gsl_test_abs (abserr, fabs(result-expected), 1e6*d, "%s abserr", desc);
    }
}

int
main ()
{
  gsl_function F1, DF1, F2, DF2, F3, DF3, F4, DF4, F5, DF5, F6, DF6;

  gsl_ieee_env_setup ();

  F1.function = &f1;
  DF1.function = &df1;

  F2.function = &f2;
  DF2.function = &df2;

  F3.function = &f3;
  DF3.function = &df3;

  F4.function = &f4;
  DF4.function = &df4;

  F5.function = &f5;
  DF5.function = &df5;

  F6.function = &f6;
  DF6.function = &df6;
  
  test (&gsl_deriv_central, &F1, &DF1, 1.0, "exp(x), x=1, central deriv");
  test (&gsl_deriv_forward, &F1, &DF1, 1.0, "exp(x), x=1, forward deriv");
  test (&gsl_deriv_backward, &F1, &DF1, 1.0, "exp(x), x=1, backward deriv");

  test (&gsl_deriv_central, &F2, &DF2, 0.1, "x^(3/2), x=0.1, central deriv");
  test (&gsl_deriv_forward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, forward deriv");
  test (&gsl_deriv_backward, &F2, &DF2, 0.1, "x^(3/2), x=0.1, backward deriv");

  test (&gsl_deriv_central, &F3, &DF3, 0.45, "sin(1/x), x=0.45, central deriv");
  test (&gsl_deriv_forward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, forward deriv");
  test (&gsl_deriv_backward, &F3, &DF3, 0.45, "sin(1/x), x=0.45, backward deriv");

  test (&gsl_deriv_central, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, central deriv");
  test (&gsl_deriv_forward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, forward deriv");
  test (&gsl_deriv_backward, &F4, &DF4, 0.5, "exp(-x^2), x=0.5, backward deriv");

  test (&gsl_deriv_central, &F5, &DF5, 0.0, "x^2, x=0, central deriv");
  test (&gsl_deriv_forward, &F5, &DF5, 0.0, "x^2, x=0, forward deriv");
  test (&gsl_deriv_backward, &F5, &DF5, 0.0, "x^2, x=0, backward deriv");

  test (&gsl_deriv_central, &F6, &DF6, 10.0, "1/x, x=10, central deriv");
  test (&gsl_deriv_forward, &F6, &DF6, 10.0, "1/x, x=10, forward deriv");
  test (&gsl_deriv_backward, &F6, &DF6, 10.0, "1/x, x=10, backward deriv");

  exit (gsl_test_summary ());
}


