htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/result.c

   1 /* specfunc/result.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* Author:  G. Jungman */
  21
  22 #include <config.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_sf_exp.h>
  26 #include <gsl/gsl_sf_result.h>
  27
  28
  29 int
  30 gsl_sf_result_smash_e(const gsl_sf_result_e10 * re, gsl_sf_result * r)
  31 {
  32   if(re->e10 == 0) {
  33     /* nothing to smash */
  34     r->val = re->val;
  35     r->err = re->err;
  36     return GSL_SUCCESS;
  37   }
  38   else {
  39     const double av = fabs(re->val);
  40     const double ae = fabs(re->err);
  41
  42     if(   GSL_SQRT_DBL_MIN < av && av < GSL_SQRT_DBL_MAX
  43        && GSL_SQRT_DBL_MIN < ae && ae < GSL_SQRT_DBL_MAX
  44        && 0.49*GSL_LOG_DBL_MIN  < re->e10 && re->e10 < 0.49*GSL_LOG_DBL_MAX
  45        ) {
  46       const double scale = exp(re->e10 * M_LN10);
  47       r->val = re->val * scale;
  48       r->err = re->err * scale;
  49       return GSL_SUCCESS;
  50     }
  51     else {
  52       return gsl_sf_exp_mult_err_e(re->e10*M_LN10, 0.0, re->val, re->err, r);
  53     }
  54   }
  55 /*
  56   int stat_v;
  57   int stat_e;
  58
  59   if(re->val == 0.0) {
  60     r->val = 0.0;
  61     stat_v = GSL_SUCCESS;
  62   }
  63   else {
  64     gsl_sf_result r_val;
  65     const double s = GSL_SIGN(re->val);
  66     const double x_v = re->e10*M_LN10 + log(fabs(re->val));
  67     stat_v = gsl_sf_exp_e(x_v, &r_val);
  68     r->val = s * r_val.val;
  69   }
  70
  71   if(re->err == 0.0) {
  72     r->err = 0.0;
  73     stat_e = GSL_SUCCESS;
  74   }
  75   else if(re->val != 0.0) {
  76     r->err = fabs(r->val * re->err/re->val);
  77     stat_e = GSL_SUCCESS;
  78   }
  79   else {
  80     gsl_sf_result r_err;
  81     const double x_e = re->e10*M_LN10 + log(fabs(re->err));
  82     stat_e = gsl_sf_exp_e(x_e, &r_err);
  83     r->err = r_err.val;
  84   }
  85
  86   return GSL_ERROR_SELECT_2(stat_v, stat_e);
  87 */
  88 }
  89