/* integration/tests.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough

 * 

 * 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 3 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

#include <config.h>

#include <math.h>

#include <gsl/gsl_math.h>

#include "tests.h"

/* These are the test functions from table 4.1 of the QUADPACK book */

/* f1(x) = x^alpha * log(1/x) */

/* integ(f1,x,0,1) = 1/(alpha + 1)^2 */

double f1 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(x,alpha) * log(1/x) ;

}

/* f2(x) = 4^-alpha / ((x-pi/4)^2 + 16^-alpha) */

/* integ(f2,x,0,1) = arctan((4-pi)4^(alpha-1)) + arctan(pi 4^(alpha-1)) */

double f2 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(4.0,-alpha) / (pow((x-M_PI/4.0),2.0) + pow(16.0,-alpha)) ;

}

/* f3(x) = cos(2^alpha * sin(x)) */

/* integ(f3,x,0,pi) = pi J_0(2^alpha) */

double f3 (double x, void * params) {

  double alpha = *(double *) params ;

  return cos(pow(2.0,alpha) * sin(x)) ;

}

/* Functions 4, 5 and 6 are duplicates of functions  1, 2 and 3 */

/* ....                                                         */

/* f7(x) = |x - 1/3|^alpha */

/* integ(f7,x,0,1) = ((2/3)^(alpha+1) + (1/3)^(alpha+1))/(alpha + 1) */

double f7 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(fabs(x - (1.0/3.0)),alpha) ;

}

/* f8(x) = |x - pi/4|^alpha */

/* integ(f8,x,0,1) = 

   ((1 - pi/4)^(alpha+1) + (pi/4)^(alpha+1))/(alpha + 1) */

double f8 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(fabs(x - (M_PI/4.0)),alpha) ;

}

/* f9(x) = sqrt(1 - x^2) / (x + 1 + 2^-alpha) */

/* integ(f9,x,-1,1) = pi/sqrt((1+2^-alpha)^2-1) */

double f9 (double x, void * params) {

  double alpha = *(double *) params ;

  return 1 / ((x + 1 + pow(2.0,-alpha)) * sqrt(1-x*x)) ;

}

/* f10(x) = sin(x)^(alpha - 1) */

/* integ(f10,x,0,pi/2) = 2^(alpha-2) ((Gamma(alpha/2))^2)/Gamma(alpha) */

double f10 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(sin(x), alpha-1) ;

}

/* f11(x) = log(1/x)^(alpha - 1) */

/* integ(f11,x,0,1) = Gamma(alpha) */

double f11 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(log(1/x), alpha-1) ;

}

/* f12(x) = exp(20*(x-1)) * sin(2^alpha * x) */

/* integ(f12,x,0,1) = 

   (20 sin(2^alpha) - 2^alpha cos(2^alpha) + 2^alpha exp(-20))

   /(400 + 4^alpha) */

double f12 (double x, void * params) {

  double alpha = *(double *) params ;

  return exp(20*(x-1)) * sin(pow(2.0,alpha) * x) ;

}

/* f13(x) = cos(2^alpha * x)/sqrt(x(1 - x)) */

/* integ(f13,x,0,1) = pi cos(2^(alpha-1)) J_0(2^(alpha-1))  */

double f13 (double x, void * params) {

  double alpha = *(double *) params ;

  return cos(pow(2.0,alpha)*x)/sqrt(x*(1-x)) ;

}

double f14 (double x, void * params) {

  double alpha = *(double *) params ;

  return exp(-pow(2.0,-alpha)*x)*cos(x)/sqrt(x) ;

}

double f15 (double x, void * params) {

  double alpha = *(double *) params ;

  return x*x * exp(-pow(2.0,-alpha)*x) ;

}

double f16 (double x, void * params) {

  double alpha = *(double *) params ;

  if (x==0 && alpha == 1) return 1 ;  /* make the function continuous in x */

  if (x==0 && alpha > 1) return 0 ;   /* avoid problems with pow(0,1) */

  return pow(x,alpha-1)/pow((1+10*x),2.0) ;

}

double f17 (double x, void * params) {

  double alpha = *(double *) params ;

  return pow(2.0,-alpha)/(((x-1)*(x-1)+pow(4.0,-alpha))*(x-2)) ;

}

/* f454(x) = x^3 log|(x^2-1)(x^2-2)| */

/* integ(f454,x,0,inf) = 61 log(2) + (77/4) log(7) - 27 */

double f454 (double x, void * params) {

  double x2 = x * x;

  double x3 = x * x2;

  params = 0 ;

  return x3 * log(fabs((x2 - 1.0) * (x2 - 2.0))) ;

}

/* f455(x) = log(x)/(1+100*x^2) */

/* integ(f455,x,0,inf) = -log(10)/20 */

double f455 (double x, void * params) {

  params = 0 ;

  return log(x) / (1.0 + 100.0 * x * x) ;

}

/* f456(x) = log(x) */

/* integ(f456*sin(10 pi x),x,0,1) = -(gamma + log(10pi) - Ci(10pi))/(10pi) */

double f456 (double x, void * params) {

  params = 0 ;

  if (x == 0.0)

    {

      return 0;

    }

  return log(x) ;

}

/* f457(x) = 1/sqrt(x) */

/* integ(f457*cos(pi x / 2),x,0,+inf) = 1 */

double f457 (double x, void * params) {

  params = 0 ;

  if (x == 0.0)

    {

      return 0;

    }

  return 1/sqrt(x) ;

}

/* f458(x) = 1/(1 + log(x)^2)^2 */

/* integ(log(x) f458(x),x,0,1) = (Ci(1) sin(1) + (pi/2 - Si(1)) cos(1))/pi 

                               = -0.1892752 */

double f458 (double x, void * params) {

  params = 0 ;

  if (x == 0.0) 

    {

      return 0;

    }

  else 

    {

      double u = log(x);

      double v = 1 + u * u;

      return 1.0 / (v * v) ;

    }

}

/* f459(x) = 1/(5 x^3 + 6) */

/* integ(f459/(x-0),x,-1,5) = log(125/631)/18 */

double f459 (double x, void * params) {

  params = 0 ;

  return 1.0 / (5.0 * x * x * x + 6.0) ;

}

/* myfn1(x) = exp(-x - x^2) */

/* integ(myfn1,x,-inf,inf) = sqrt(pi) exp(-1/4) */

double myfn1 (double x, void * params) {

  params = 0;

  return exp(-x - x*x) ;

}

/* myfn2(x) = exp(alpha*x) */

/* integ(myfn2,x,-inf,b) = exp(alpha*b)/alpha */

double myfn2 (double x, void * params) {

  double alpha = *(double *) params ;

  return exp(alpha*x) ;

}

/* f_monomial = constant * x^degree */

double f_monomial(double x, void * params)

{

  struct monomial_params * p = (struct monomial_params *) params;

  return p->constant * gsl_pow_int(x, p->degree);

}

/* integ(f_monomial,x,a b)=constant*(b^(degree+1)-a^(degree+1))/(degree+1) */

double integ_f_monomial(double a, double b, struct monomial_params * p)

{

  const int degreep1 = p->degree + 1;

  const double bnp1 = gsl_pow_int(b, degreep1);

  const double anp1 = gsl_pow_int(a, degreep1);

  return (p->constant / degreep1)*(bnp1 - anp1);

}

/* f(x) = sin(x) */

double f_sin(double x, void * params)

{

    return sin(x);

}

/* integ(f_sin,x,a,b) */

double integ_f_sin(double a, double b)

{

    return -cos(b) + cos(a);

}

/* The test functions. */

double cqf1 ( double x , void *params ) {

  return exp(x);

}

double cqf2 ( double x , void *params ) {

  return x >= 0.3;

}

double cqf3 ( double x , void *params ) {

  return sqrt(x);

}

double cqf4 ( double x , void *params ) {

  return (23.0/25) * cosh(x) - cos(x);

}

double cqf5 ( double x , void *params ) {

  double x2 = x*x;

  return 1.0 / ( x2 * (x2 + 1) + 0.9);

}

double cqf6 ( double x , void *params ) {

  return x * sqrt( x );

}

double cqf7 ( double x , void *params ) {

  return 1.0 / sqrt(x);

}

double cqf8 ( double x , void *params ) {

  double x2 = x*x;

  return 1.0 / (1 + x2*x2);

}

double cqf9 ( double x , void *params ) {

  return 2.0 / (2 + sin(10*M_PI*x));

}

double cqf10 ( double x , void *params ) {

  return 1.0 / (1 + x);

}

double cqf11 ( double x , void *params ) {

  return 1.0 / (1 + exp(x));

}

double cqf12 ( double x , void *params ) {

  return x / (exp(x) - 1.0);

}

double cqf13 ( double x , void *params ) {

  return sin(100 * M_PI * x) / (M_PI * x);

}

double cqf14 ( double x , void *params ) {

  return sqrt(50.0) * exp(-50*M_PI*x*x);

}

double cqf15 ( double x , void *params ) {

  return 25.0 * exp(-25*x);

}

double cqf16 ( double x , void *params ) {

  return 50 / M_PI * (2500 * x*x + 1);

}

double cqf17 ( double x , void *params ) {

  double t1 = 50 * M_PI * x ,t2;

  t2 = sin(t1) / t1;

  return 50 * t2 * t2;

}

double cqf18 ( double x , void *params ) {

  return cos( cos(x) + 3*sin(x) + 2*cos(2*x) + 3*sin(2*x) + 3*cos(3*x) );

}

double cqf19 ( double x , void *params ) {

  return log(x);

}

double cqf20 ( double x , void *params ) {

  return 1 / (x*x + 1.005);

}

double cqf21 ( double x , void *params ) {

  return 1 / cosh( 10 * (x - 0.2) * 2 ) + 

    1 / cosh( 100 * (x - 0.4) * 4 ) + 

    1 / cosh( 1000 * (x - 0.6) * 8 );

}

double cqf22 ( double x , void *params ) {

  return 4 * M_PI*M_PI * x * sin(20*M_PI*x) * cos(2*M_PI*x);

}

double cqf23 ( double x , void *params ) {

  double t = 230*x - 30;

  return 1 / (1 + t*t);

}

double cqf24 ( double x , void *params ) {

  return floor(exp(x));

}

double cqf25 ( double x , void *params ) {

  return (x < 1) * (x + 1) + 

    (1 <= x && x <= 3) * (3 - x) + 

    (x > 3) * 2;

}