test_sf.c

math library from gnu

/* specfunc/test_sf.c
 * 
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2004 Gerard Jungman
 * 
 * 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"


double
test_sf_frac_diff(double x1, double x2)
{
  if(x1 == 0.0 && x2 == 0.0) {
    return 0.0;
  }
  else if (x1 == 0.0) {
    return fabs(x2);
  } else if(x1 <= DBL_MAX && x2 <= DBL_MAX && (x1 + x2 != 0.0)) {
    return fabs((x1-x2)/(x1+x2));
  }
  else {
    return 1.0;
  }
}


/* Check a result against a given value and tolerance.
 */
int
test_sf_check_result(char * message_buff, gsl_sf_result r, double val, double tol)
{
  int    s = 0;
  double f = 0, d = 0;

  if (gsl_isnan(r.val) || gsl_isnan(val)) 
    {
      s = (gsl_isnan(r.val) != gsl_isnan(val)) ? TEST_SF_INCONS : s; 
    }
  else if (gsl_isinf(r.val) || gsl_isinf(val)) 
    {
      s = (gsl_isinf(r.val) != gsl_isinf(val)) ? TEST_SF_INCONS : s; 
    }
  else
    {
      f = test_sf_frac_diff(val, r.val);
      d = fabs(val - r.val);

      if(d > 2.0 * TEST_SIGMA * r.err) s |= TEST_SF_INCONS;
      if(r.err < 0.0)                                  s |= TEST_SF_ERRNEG;
      if(gsl_isinf(r.err))                             s |= TEST_SF_ERRBAD;
#if TEST_EXCESSIVE_ERROR
      if(d > 0 && r.err > 1e4 * fabs(val)*tol)         s |= TEST_SF_ERRBIG;
#endif
      if(f > TEST_FACTOR * tol)                        s |= TEST_SF_TOLBAD;
    }

  if(s != 0) {
    char buff[2048];
    sprintf(buff, "  expected: %20.16e\n", val);
    strcat(message_buff, buff);
    sprintf(buff, "  obtained: %20.16e +/- %.16e (rel=%g)\n", r.val, r.err, r.err/(fabs(r.val) + r.err));
    strcat(message_buff, buff);
    sprintf(buff, "  fracdiff: %20.16e\n", f);
    strcat(message_buff, buff);
    sprintf(buff, " tolerance: %20.16e\n", tol);
    strcat(message_buff, buff);
  }

  if(s & TEST_SF_INCONS) {
    strcat(message_buff, "  value/expected not consistent within reported error\n");
  }
  if(s & TEST_SF_ERRNEG) {
    strcat(message_buff, "  reported error negative\n");
  }
  if(s & TEST_SF_ERRBAD) {
    strcat(message_buff, "  reported error is bad\n");
  }
  if(s & TEST_SF_ERRBIG) {
    strcat(message_buff, "  reported error is much too big\n");
  }
  if(s & TEST_SF_TOLBAD) {
    strcat(message_buff, "  value not within tolerance of expected value\n");
  }

  return s;
}

/* Check a result against a given value and tolerance.
 */
int
test_sf_check_e10(char * message_buff, int e10, int e10_in)
{
  int    s = 0;

  if (e10 != e10_in) 
    {
      s = TEST_SF_EXPBAD;
    }

  if(s != 0) {
    char buff[2048];
    sprintf(buff, "  expected exponent: 10^%d\n", e10_in);
    strcat(message_buff, buff);
    sprintf(buff, "  obtained exponent: 10^%d", e10);
    strcat(message_buff, buff);
  }

  if(s & TEST_SF_EXPBAD) {
    strcat(message_buff, "  exponent is incorrect\n");
  }

  return s;
}

int
test_sf_check_val(char * message_buff, double rval, double val, double tol)
{
  int    s = 0;
  double f = test_sf_frac_diff(val, rval);

  if(f > TEST_FACTOR * tol)                       s |= TEST_SF_TOLBAD;

  if(s != 0) {
    char buff[2048];
    sprintf(buff, "  expected: %20.16e\n", val);
    strcat(message_buff, buff);
    sprintf(buff, "  obtained: %20.16e\n", rval);
    strcat(message_buff, buff);
    sprintf(buff, "  fracdiff: %20.16e\n", f);
    strcat(message_buff, buff);
  }

  if(s & TEST_SF_TOLBAD) {
    strcat(message_buff, "  value not within tolerance of expected value\n");
  }

  return s;
}


/* Relax the condition on the agreement and on the usefulness
 * of the error estimate.
 */
int
test_sf_check_result_relax(char * message_buff, gsl_sf_result r, double val, double tol)
{
  int    s = 0;
  double f = test_sf_frac_diff(val, r.val);

  if(f > GSL_MAX_DBL(TEST_SNGL, TEST_FACTOR * tol))   s |= TEST_SF_INCONS;
  if(r.err < 0.0)     s |= TEST_SF_ERRNEG;
  if(gsl_isinf(r.err))              s |= TEST_SF_ERRBAD;
  if(f > TEST_FACTOR * tol)         s |= TEST_SF_TOLBAD;

  if(s != 0) {
    char buff[2048];
    sprintf(buff, "  expected: %20.16e\n", val);
    strcat(message_buff, buff);
    sprintf(buff, "  obtained: %20.16e +/- %.16e  (rel=%g)\n", r.val, r.err, r.err/(fabs(r.val) + r.err));
    strcat(message_buff, buff);
    sprintf(buff, "  fracdiff: %20.16e\n", f);
    strcat(message_buff, buff);
  }

  if(s & TEST_SF_INCONS) {
    strcat(message_buff, "  value/expected not consistent MAX(tol,SINGLE_PREC)\n");
  }
  if(s & TEST_SF_ERRNEG) {
    strcat(message_buff, "  reported error negative\n");
  }
  if(s & TEST_SF_ERRBAD) {
    strcat(message_buff, "  reported error is bad\n");
  }
  if(s & TEST_SF_TOLBAD) {
    strcat(message_buff, "  value not within tolerance of expected value\n");
  }

  return s;
}


/* Check a return value.
 */
int
test_sf_check_return(char * message_buff, int val_return, int expected_return)
{
  if(val_return != expected_return) {
    char buff[256];
    sprintf(buff, "  unexpected return code: %d\n", val_return);
    strcat(message_buff, buff);
    return TEST_SF_RETBAD;
  }
  else {
    return 0;
  }
}


int
test_sf (gsl_sf_result r, double val_in, double tol, int status,
         int expect_return, const char * desc)
{
  char message_buff[4096];
  int local_s = 0;

  message_buff[0] = '\0';