The R Project SVN R

Rev

Rev 68947 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 68947 Rev 74312
Line 1... Line 1...
1 
/*
 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
 2 
 *  Mathlib : A C Library of Special Functions
3 
 *  Copyright (C) 1998 Ross Ihaka
 3 
 *  Copyright (C) 1998 Ross Ihaka
4 
 *  Copyright (C) 2000, 2003, 2011 The R Core Team
 4 
 *  Copyright (C) 2000-2018 The R Core Team
5 
 *
 5 
 *
6 
 *  This program is free software; you can redistribute it and/or modify
 6 
 *  This program is free software; you can redistribute it and/or modify
7 
 *  it under the terms of the GNU General Public License as published by
 7 
 *  it under the terms of the GNU General Public License as published by
8 
 *  the Free Software Foundation; either version 2 of the License, or
 8 
 *  the Free Software Foundation; either version 2 of the License, or
9 
 *  (at your option) any later version.
 9 
 *  (at your option) any later version.
Line 39... Line 39...
39 
#ifdef HAVE_CONFIG_H
 39 
#ifdef HAVE_CONFIG_H
40 
# include <config.h>
 40 
# include <config.h>
41 
#endif
 41 
#endif
42 
#include "nmath.h"
 42 
#include "nmath.h"
43 
 
 43 
 
44 
/* want to compile log1p as Rlog1p if HAVE_LOG1P && !HAVE_WORKING_LOG1P */
 - 
 
45 
#if defined(HAVE_LOG1P) && !defined(HAVE_WORKING_LOG1P)
 - 
 
46 
#undef HAVE_LOG1P
 - 
 
47 
#endif
 - 
 
48 
 
 - 
 
49 
#ifndef HAVE_LOG1P
 44 
#ifndef HAVE_WORKING_LOG1P
50 
double log1p(double x)
 45 
double Rlog1p(double x)
51 
{
 46 
{
52 
    /* series for log1p on the interval -.375 to .375
 47 
    /* series for log1p on the interval -.375 to .375
53 
     *				     with weighted error   6.35e-32
 48 
     *				     with weighted error   6.35e-32
54 
     *				      log weighted error  31.20
 49 
     *				      log weighted error  31.20
55 
     *			    significant figures required  30.93
 50 
     *			    significant figures required  30.93
Line 135... Line 130...
135 
	ML_ERROR(ME_PRECISION, "log1p");
 130 
	ML_ERROR(ME_PRECISION, "log1p");
136 
    }
 131 
    }
137 
    return log(1 + x);
 132 
    return log(1 + x);
138 
}
 133 
}
139 
#endif
 134 
#endif
140 
 
 - 
 
141 
 
 - 
 
142 
 
 - 
 
143 
#ifndef HAVE_HYPOT
 - 
 
144 
/* Used as a substitute for the C99 function hypot, which all currently
 - 
 
145 
   known platforms have */
 - 
 
146 
 
 - 
 
147 
/* hypot(a,b)	finds sqrt(a^2 + b^2)
 - 
 
148 
 *		without overflow or destructive underflow.
 - 
 
149 
 */
 - 
 
150 
 
 - 
 
151 
double hypot(double a, double b)
 - 
 
152 
{
 - 
 
153 
    double p, r, s, t, tmp, u;
 - 
 
154 
 
 - 
 
155 
    if(ISNAN(a) || ISNAN(b)) /* propagate Na(N)s: */
 - 
 
156 
        return
 - 
 
157 
#ifdef IEEE_754
 - 
 
158 
	  a + b;
 - 
 
159 
#else
 - 
 
160 
          ML_NAN;
 - 
 
161 
#endif
 - 
 
162 
    if (!R_FINITE(a) || !R_FINITE(b)) {
 - 
 
163 
        return ML_POSINF;
 - 
 
164 
    }
 - 
 
165 
    p = fmax2(fabs(a), fabs(b));
 - 
 
166 
    if (p != 0.0) {
 - 
 
167 
 
 - 
 
168 
	/* r = (min(|a|,|b|) / p) ^2 */
 - 
 
169 
	tmp = fmin2(fabs(a), fabs(b))/p;
 - 
 
170 
	r = tmp * tmp;
 - 
 
171 
	for(;;) {
 - 
 
172 
	    t = 4.0 + r;
 - 
 
173 
	    /* This was a test of 4.0 + r == 4.0, but optimizing
 - 
 
174 
		compilers nowadays infinite loop on that. */
 - 
 
175 
	    if(fabs(r) < 2*DBL_EPSILON) break;
 - 
 
176 
	    s = r / t;
 - 
 
177 
	    u = 1. + 2. * s;
 - 
 
178 
	    p *= u ;
 - 
 
179 
 
 - 
 
180 
	    /* r = (s / u)^2 * r */
 - 
 
181 
	    tmp = s / u;
 - 
 
182 
	    r *= tmp * tmp;
 - 
 
183 
	}
 - 
 
184 
    }
 - 
 
185 
    return p;
 - 
 
186 
}
 - 
 
187 
#endif
 -