The R Project SVN R

Rev

Rev 88128 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
42303 ripley 1 
/*
2 
 *  Mathlib : A C Library of Special Functions
87903 ripley 3 
 *  Copyright (C) 1998-2025 The R Core Team
42303 ripley 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, a copy is available at
68947 ripley 17 
 *  https://www.R-project.org/Licenses/
42303 ripley 18 
 */
19 
 
20 
#ifdef HAVE_CONFIG_H
21 
# include <config.h>
22 
# undef fprintf
23 
#endif
24 
#include "nmath.h"
25 
 
26 
#ifdef MATHLIB_STANDALONE
27 
/*
28 
 *  based on code in ../main/arithmetic.c
29 
 *  used only in standalone Rmath lib.
30 
 */
31 
 
32 
int R_finite(double x)
33 
{
34 
#ifdef HAVE_WORKING_ISFINITE
35 
    return isfinite(x);
36 
# else
37 
    return (!isnan(x) & (x != ML_POSINF) & (x != ML_NEGINF));
38 
#endif
39 
}
40 
 
41 
/* C++ math header undefines any isnan macro. This file
42 
   doesn't get C++ headers and so is safe. */
43 
int R_isnancpp(double x)
44 
{
42422 ripley 45 
    return (isnan(x) != 0);
42303 ripley 46 
}
47 
 
48 
static double myfmod(double x1, double x2)
49 
{
50 
    double q = x1 / x2;
51 
    return x1 - floor(q) * x2;
52 
}
53 
 
54 
double R_pow(double x, double y) /* = x ^ y */
55 
{
56 
    if(x == 1. || y == 0.)
57 
	return(1.);
58 
    if(x == 0.) {
59 
	if(y > 0.) return(0.);
86145 maechler 60 
	else if(y < 0) return(ML_POSINF);
61 
	else return(y); /* y is NA or NaN, we assert */
42303 ripley 62 
    }
63 
    if (R_FINITE(x) && R_FINITE(y))
64 
	return(pow(x,y));
65 
    if (ISNAN(x) || ISNAN(y)) {
66 
#ifdef IEEE_754
67 
	return(x + y);
68 
#else
69 
	return(ML_NAN);
70 
#endif
71 
    }
72 
    if(!R_FINITE(x)) {
73 
	if(x > 0)		/* Inf ^ y */
74 
	    return((y < 0.)? 0. : ML_POSINF);
75 
	else {			/* (-Inf) ^ y */
76 
	    if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */
86145 maechler 77 
		return (y < 0.) ? 0. : (myfmod(y,2.) != 0 ? x  : -x);
42303 ripley 78 
	}
79 
    }
80 
    if(!R_FINITE(y)) {
81 
	if(x >= 0) {
82 
	    if(y > 0)		/* y == +Inf */
83 
		return((x >= 1)? ML_POSINF : 0.);
84 
	    else		/* y == -Inf */
85 
		return((x < 1) ? ML_POSINF : 0.);
86 
	}
87 
    }
88 
    return(ML_NAN);		/* all other cases: (-Inf)^{+-Inf,
89 
				   non-int}; (neg)^{+-Inf} */
90 
}
91 
 
92 
double R_pow_di(double x, int n)
93 
{
94 
    double pow = 1.0;
95 
 
96 
    if (ISNAN(x)) return x;
97 
    if (n != 0) {
98 
	if (!R_FINITE(x)) return R_pow(x, (double)n);
99 
	if (n < 0) { n = -n; x = 1/x; }
100 
	for(;;) {
101 
	    if(n & 01) pow *= x;
102 
	    if(n >>= 1) x *= x; else break;
103 
	}
104 
    }
105 
    return pow;
106 
}
107 
 
88128 ripley 108 
/* It is not clear why these are being defined in standalone nmath:
109 
 * but that they are is stated in the R-admin manual.
110 
 *
88135 ripley 111 
 * In R NA_REAL is a specific NaN computed during initialization.
88128 ripley 112 
 */
113 
#if defined(__clang__) && defined(NAN)
114 
// C99 (optionally) has NAN, which is a float but will coerce to double.
115 
double NA_REAL = NAN;
116 
#else
117 
// ML_NAN is defined as (0.0/0.0) in nmath.h
118 
// Fails to compile in Intel ics 2025.0, Apple clang 17, LLVM clang 20
119 
double NA_REAL = ML_NAN;
120 
#endif
121 
 
42303 ripley 122 
double R_PosInf = ML_POSINF, R_NegInf = ML_NEGINF;
123 
 
124 
#include <stdio.h>
125 
#include <stdarg.h>
87903 ripley 126 
attribute_hidden void REprintf(const char *format, ...)
42303 ripley 127 
{
128 
    va_list(ap);
129 
    va_start(ap, format);
130 
    fprintf(stderr, format, ap);
131 
    va_end(ap);
132 
}
133 
 
134 
#endif /* MATHLIB_STANDALONE */