The R Project SVN R

Rev

Rev 18834 | Rev 29441 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
574 ihaka 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
12976 pd 3 
 *  Copyright (C) 1998-2001 Ross Ihaka and the R Development Core Team
574 ihaka 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 2 of the License, or
8 
 *  (at your option) any later version.
9 
 *
10 
 *  This program is distributed in the hope that it will be useful,
11 
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 
 *  GNU 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
5458 ripley 17 
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA.
574 ihaka 18 
 */
19 
 
8431 ripley 20 
#ifdef HAVE_CONFIG_H
18834 hornik 21 
# include <config.h>
574 ihaka 22 
#endif
8431 ripley 23 
#include "nmath.h"
574 ihaka 24 
 
25 
#ifndef IEEE_754
26 
 
27 
void ml_error(int n)
28 
{
29 
    switch(n) {
30 
 
31 
    case ME_NONE:
32 
	errno = 0;
33 
	break;
34 
 
35 
    case ME_DOMAIN:
36 
    case ME_NOCONV:
37 
	errno = EDOM;
38 
	break;
39 
 
40 
    case ME_RANGE:
41 
	errno = ERANGE;
42 
	break;
43 
 
44 
    default:
45 
	break;
46 
    }
47 
}
48 
 
49 
#endif
8431 ripley 50 
 
51 
#ifdef MATHLIB_STANDALONE
52 
/*
53 
 *  based on code in ../main/arithmetic.c
54 
 */
55 
 
56 
 
57 
#ifdef IEEE_754
58 
 
59 
int R_IsNaNorNA(double x)
60 
{
61 
/* NOTE: some systems do not return 1 for TRUE. */
62 
    return (isnan(x) != 0);
63 
}
64 
 
65 
/* Include the header file defining finite() */
66 
#ifdef HAVE_IEEE754_H
67 
# include <ieee754.h>		/* newer Linuxen */
68 
#else
69 
# ifdef HAVE_IEEEFP_H
70 
#  include <ieeefp.h>		/* others [Solaris 2.5.x], .. */
71 
# endif
72 
#endif
73 
#if defined(Win32) && defined(_MSC_VER)
15252 hornik 74 
# include <float.h>
8431 ripley 75 
#endif
76 
 
77 
int R_finite(double x)
78 
{
15252 hornik 79 
#ifdef HAVE_WORKING_FINITE
8431 ripley 80 
    return finite(x);
15252 hornik 81 
#else
82 
# ifdef _AIX
83 
#  include <fp.h>
84 
     return FINITE(x);
8431 ripley 85 
# else
86 
    return (!isnan(x) & (x != ML_POSINF) & (x != ML_NEGINF));
15252 hornik 87 
# endif
8431 ripley 88 
#endif
89 
}
90 
 
91 
#else /* not IEEE_754 */
92 
 
93 
int R_IsNaNorNA(double x)
94 
{
95 
# ifndef HAVE_ISNAN
96 
    return (x == ML_NAN);
97 
# else
98 
    return (isnan(x) != 0 || x == ML_NAN);
99 
# endif
100 
}
101 
 
102 
int R_finite(double x)
103 
{
104 
# ifndef HAVE_FINITE
105 
    return (x !=  ML_NAN && x < ML_POSINF && x > ML_NEGINF);
106 
# else
107 
    int finite(double);
108 
    return finite(x);
109 
# endif
110 
}
111 
#endif /* IEEE_754 */
112 
 
9363 ripley 113 
static double myfmod(double x1, double x2)
114 
{
115 
    double q = x1 / x2;
116 
    return x1 - floor(q) * x2;
117 
}
118 
 
15252 hornik 119 
#ifdef HAVE_WORKING_LOG
120 
# define R_log	log
121 
#else
9363 ripley 122 
double R_log(double x) { return(x > 0 ? log(x) : x < 0 ? ML_NAN : ML_NEGINF); }
123 
#endif
124 
 
125 
double R_pow(double x, double y) /* = x ^ y */
126 
{
127 
    if(x == 1. || y == 0.)
128 
	return(1.);
129 
    if(x == 0.) {
130 
	if(y > 0.) return(0.);
131 
	/* y < 0 */return(ML_POSINF);
132 
    }
133 
    if (R_FINITE(x) && R_FINITE(y))
134 
	return(pow(x,y));
135 
    if (ISNAN(x) || ISNAN(y)) {
136 
#ifdef IEEE_754
137 
	return(x + y);
138 
#else
139 
	return(NA_REAL);
140 
#endif
141 
    }
142 
    if(!R_FINITE(x)) {
143 
	if(x > 0)		/* Inf ^ y */
144 
	    return((y < 0.)? 0. : ML_POSINF);
145 
	else {			/* (-Inf) ^ y */
146 
	    if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */
147 
		return((y < 0.) ? 0. : (myfmod(y,2.) ? x  : -x));
148 
	}
149 
    }
150 
    if(!R_FINITE(y)) {
151 
	if(x >= 0) {
152 
	    if(y > 0)		/* y == +Inf */
153 
		return((x >= 1)? ML_POSINF : 0.);
154 
	    else		/* y == -Inf */
155 
		return((x < 1) ? ML_POSINF : 0.);
156 
	}
157 
    }
158 
    return(ML_NAN);		/* all other cases: (-Inf)^{+-Inf,
159 
				   non-int}; (neg)^{+-Inf} */
160 
}
161 
 
162 
double R_pow_di(double x, int n)
163 
{
164 
    double pow = 1.0;
165 
 
166 
    if (ISNAN(x)) return x;
167 
    if (n != 0) {
168 
	if (!R_FINITE(x)) return R_pow(x, (double)n);
169 
	if (n < 0) { n = -n; x = 1/x; }
170 
	for(;;) {
171 
	    if(n & 01) pow *= x;
172 
	    if(n >>= 1) x *= x; else break;
173 
	}
174 
    }
175 
    return pow;
176 
}
177 
 
178 
double NA_REAL = ML_NAN;
179 
double R_PosInf = ML_POSINF, R_NegInf = ML_NEGINF;
180 
 
8431 ripley 181 
#endif /* MATHLIB_STANDALONE */