The R Project SVN – R – Blame – /trunk/src/main/arithmetic.c

Rev	Author	Line No.	Line
2	r	1	`/*`
834	maechler	2	`* R : A Computer Language for Statistical Data Analysis`
2	r	3	`* Copyright (C) 1995, 1996, 1997 Robert Gentleman and Ross Ihaka`
24516	ripley	4	`* Copyright (C) 1998--2003 The R Development Core Team.`
27167	maechler	5	`* Copyright (C) 2003 The R Foundation`
2	r	6	`*`
		7	`* This program is free software; you can redistribute it and/or modify`
		8	`* it under the terms of the GNU General Public License as published by`
		9	`* the Free Software Foundation; either version 2 of the License, or`
		10	`* (at your option) any later version.`
		11	`*`
		12	`* This program is distributed in the hope that it will be useful,`
		13	`* but WITHOUT ANY WARRANTY; without even the implied warranty of`
		14	`* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the`
		15	`* GNU General Public License for more details.`
		16	`*`
27167	maechler	17	`* A copy of the GNU General Public License is available via WWW at`
		18	`* http://www.gnu.org/copyleft/gpl.html. You can also obtain it by`
		19	`* writing to the Free Software Foundation, Inc., 59 Temple Place,`
		20	`* Suite 330, Boston, MA 02111-1307 USA.`
2	r	21	`*/`
		22
5187	hornik	23	`#ifdef HAVE_CONFIG_H`
7701	hornik	24	`#include <config.h>`
5187	hornik	25	`#endif`
		26
14701	ripley	27	`#ifdef __OpenBSD__`
		28	`/* for definition of "struct exception" in math.h */`
		29	`# define __LIBM_PRIVATE`
		30	`#endif`
		31	`#include "Defn.h" /-> Arith.h -> math.h /`
		32	`#ifdef __OpenBSD__`
		33	`# undef __LIBM_PRIVATE`
		34	`#endif`
		35
11499	ripley	36	`#define MATHLIB_PRIVATE`
		37	`#include <Rmath.h>`
		38	`#undef MATHLIB_PRIVATE`
2632	maechler	39	`#include "arithmetic.h"`
2	r	40
1881	ihaka	41	`/* Error Handling for Floating Point Errors */`
2	r	42
580	ihaka	43	`#ifndef IEEE_754`
7658	ripley	44	`#ifdef Unix`
2	r	45	`#include <signal.h>`
		46
		47	`static RETSIGTYPE handle_fperror(int dummy)`
		48	`{`
571	ihaka	49	`errno = ERANGE;`
		50	`signal(SIGFPE, handle_fperror);`
7658	ripley	51	`}`
13997	duncan	52	`#endif /* Unix */`
7655	ripley	53	`#endif /* not IEEE_754 */`
2	r	54
		55	`#ifdef HAVE_MATHERR`
		56
1881	ihaka	57	`/* Override the SVID matherr function */`
2	r	58
		59	`int matherr(struct exception *exc)`
		60	`{`
571	ihaka	61	`switch (exc->type) {`
		62	`case DOMAIN:`
		63	`case SING:`
		64	`errno = EDOM;`
		65	`break;`
		66	`case OVERFLOW:`
		67	`errno = ERANGE;`
		68	`break;`
		69	`case UNDERFLOW:`
		70	`exc->retval = 0.0;`
		71	`break;`
		72	`}`
		73	`return 1;`
2	r	74	`}`
		75	`#endif`
		76
580	ihaka	77	`#ifdef IEEE_754`
17162	tlumley	78	`#ifndef _AIX`
13997	duncan	79	`const double R_Zero_Hack = 0.0; /* Silence the Sun compiler */`
17162	tlumley	80	`#else`
		81	`double R_Zero_Hack = 0.0;`
		82	`#endif`
580	ihaka	83	`typedef union`
		84	`{`
16704	ripley	85	`double value;`
		86	`unsigned int word[2];`
580	ihaka	87	`} ieee_double;`
		88
13997	duncan	89	`/* These variables hw and lw are only used if IEEE_754 is defined.`
		90	`The value of each is fixed once we determine the endiannes`
		91	`of the machine, and this cna be done via WORDS_BIGENDIAN.`
580	ihaka	92
13997	duncan	93	`Earlier code used to use establish_endianness()`
		94	`to compute these, but this is uncessary and makes them`
		95	`non-constant, but initialize once. This would involve`
14701	ripley	96	`a song and dance in a threaded application (e.g pthread_once(),`
13997	duncan	97	`etc.).`
		98	`*/`
21057	ripley	99
		100	`/* gcc has problems with static const on AIX and Solaris`
		101	`Solaris is for gcc 3.1 and 3.2 under -O2 32-bit on 64-bit kernel */`
		102	`#ifdef _AIX`
		103	`#define CONST`
		104	`#elif defined(sparc) && defined (__GNUC__) && __GNUC__ == 3`
		105	`#define CONST`
14632	tlumley	106	`#else`
21057	ripley	107	`#define CONST const`
14632	tlumley	108	`#endif`
580	ihaka	109
21057	ripley	110	`#ifdef WORDS_BIGENDIAN`
		111	`static CONST int hw = 0;`
		112	`static CONST int lw = 1;`
		113	`#else /* !WORDS_BIGENDIAN */`
		114	`static CONST int hw = 1;`
		115	`static CONST int lw = 0;`
		116	`#endif /* WORDS_BIGENDIAN */`
13997	duncan	117
21057	ripley	118
580	ihaka	119	`static double R_ValueOfNA(void)`
		120	`{`
23886	pd	121	`/* The gcc shipping with RedHat 9 gets this wrong without`
23894	pd	122	`* the volatile declaration. Thanks to Marc Schwartz. */`
23886	pd	123	`volatile ieee_double x;`
1881	ihaka	124	`x.word[hw] = 0x7ff00000;`
		125	`x.word[lw] = 1954;`
		126	`return x.value;`
580	ihaka	127	`}`
		128
		129	`int R_IsNA(double x)`
		130	`{`
3888	hornik	131	`if (isnan(x)) {`
1881	ihaka	132	`ieee_double y;`
		133	`y.value = x;`
		134	`return (y.word[lw] == 1954);`
		135	`}`
		136	`return 0;`
580	ihaka	137	`}`
1096	ihaka	138
		139	`int R_IsNaN(double x)`
		140	`{`
3888	hornik	141	`if (isnan(x)) {`
1881	ihaka	142	`ieee_double y;`
		143	`y.value = x;`
		144	`return (y.word[lw] != 1954);`
		145	`}`
		146	`return 0;`
1096	ihaka	147	`}`
7655	ripley	148
		149	`int R_IsNaNorNA(double x)`
		150	`{`
7671	maechler	151	`/* True for both NA and NaN.`
7655	ripley	152	`NOTE: some systems do not return 1 for TRUE. */`
		153	`return (isnan(x) != 0);`
		154	`}`
		155
		156	`/* Include the header file defining finite() */`
		157	`#ifdef HAVE_IEEE754_H`
		158	`# include <ieee754.h> /* newer Linuxen */`
		159	`#else`
		160	`# ifdef HAVE_IEEEFP_H`
		161	`# include <ieeefp.h> /* others [Solaris 2.5.x], .. */`
		162	`# endif`
580	ihaka	163	`#endif`
7655	ripley	164	`#if defined(Win32) && defined(_MSC_VER)`
15252	hornik	165	`# include <float.h>`
7671	maechler	166	`#endif`
580	ihaka	167
7655	ripley	168	`int R_finite(double x)`
		169	`{`
15252	hornik	170	`#ifdef HAVE_WORKING_FINITE`
7655	ripley	171	`return finite(x);`
15252	hornik	172	`#else`
18510	maechler	173	`# if defined(HAVE_DECL_ISFINITE) && HAVE_DECL_ISFINITE`
17382	luke	174	`return isfinite(x);`
		175	`# else`
		176	`# ifdef _AIX`
15252	hornik	177	`# include <fp.h>`
		178	`return FINITE(x);`
17382	luke	179	`# else`
7655	ripley	180	`return (!isnan(x) & (x != R_PosInf) & (x != R_NegInf));`
17382	luke	181	`# endif`
15252	hornik	182	`# endif`
7655	ripley	183	`#endif`
		184	`}`
		185
		186	`#else /* not IEEE_754 */`
		187
		188	`int R_IsNA(double x)`
		189	`{`
		190	`return (x == R_NaReal);`
		191	`}`
		192
		193	`/* NaN but not NA: never true */`
		194	`int R_IsNaN(double x)`
		195	`{`
		196	`return 0;`
		197	`}`
		198
		199	`int R_IsNaNorNA(double x)`
		200	`{`
		201	`# ifndef HAVE_ISNAN`
		202	`return (x == R_NaReal);`
		203	`# else`
		204	`return (isnan(x) != 0 \|\| x == R_NaReal);`
		205	`# endif`
		206	`}`
		207
		208	`int R_finite(double x)`
		209	`{`
		210	`# ifndef HAVE_FINITE`
7873	ripley	211	`return (x != R_NaReal && x < R_PosInf && x > R_NegInf);`
7655	ripley	212	`# else`
		213	`int finite(double);`
		214	`return finite(x);`
		215	`# endif`
		216	`}`
		217	`#endif /* IEEE_754 */`
		218
1881	ihaka	219	`/* Arithmetic Initialization */`
580	ihaka	220
2	r	221	`void InitArithmetic()`
		222	`{`
571	ihaka	223	`R_NaInt = INT_MIN;`
2	r	224
580	ihaka	225	`#ifdef IEEE_754`
7655	ripley	226	`/* establish_endianness(); */`
571	ihaka	227	`R_NaN = 0.0/R_Zero_Hack;`
580	ihaka	228	`R_NaReal = R_ValueOfNA();`
571	ihaka	229	`R_PosInf = 1.0/R_Zero_Hack;`
		230	`R_NegInf = -1.0/R_Zero_Hack;`
2	r	231	`#else`
7658	ripley	232	`R_NaN = -DBL_MAX*(1-1e-15);`
571	ihaka	233	`R_NaReal = R_NaN;`
7134	ripley	234	`R_PosInf = DBL_MAX;`
		235	`R_NegInf = -DBL_MAX;`
2	r	236	`#ifdef Unix`
571	ihaka	237	`signal(SIGFPE, handle_fperror);`
2	r	238	`#endif`
		239	`#endif`
		240	`}`
		241
		242
		243
4847	maechler	244	`static double myfmod(double x1, double x2)`
		245	`{`
		246	`double q = x1 / x2;`
		247	`return x1 - floor(q) * x2;`
		248	`}`
2	r	249
4878	maechler	250
15252	hornik	251	`#ifdef HAVE_WORKING_LOG`
		252	`# define R_log log`
		253	`#else`
4878	maechler	254	`double R_log(double x) { return(x > 0 ? log(x) : x < 0 ? R_NaN : R_NegInf); }`
4834	maechler	255	`#endif`
		256
4878	maechler	257	`#ifdef POW_DIRTY`
		258
		259	`# define R_pow pow`
		260
		261	`#else`
4839	maechler	262	`double R_pow(double x, double y) /* = x ^ y */`
4834	maechler	263	`{`
4839	maechler	264	`if(x == 1. \|\| y == 0.)`
		265	`return(1.);`
7449	maechler	266	`if(x == 0.) {`
		267	`if(y > 0.) return(0.);`
		268	`/* y < 0 */return(R_PosInf);`
		269	`}`
5107	maechler	270	`if (R_FINITE(x) && R_FINITE(y))`
4859	maechler	271	`return(pow(x,y));`
4839	maechler	272	`if (ISNAN(x) \|\| ISNAN(y)) {`
		273	`#ifdef IEEE_754`
		274	`return(x + y);`
		275	`#else`
		276	`return(NA_REAL);`
		277	`#endif`
		278	`}`
5107	maechler	279	`if(!R_FINITE(x)) {`
6098	pd	280	`if(x > 0) /* Inf ^ y */`
4839	maechler	281	`return((y < 0.)? 0. : R_PosInf);`
6098	pd	282	`else { /* (-Inf) ^ y */`
5107	maechler	283	`if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */`
4839	maechler	284	`return((y < 0.) ? 0. : (myfmod(y,2.) ? x : -x));`
		285	`}`
		286	`}`
5107	maechler	287	`if(!R_FINITE(y)) {`
4859	maechler	288	`if(x >= 0) {`
6098	pd	289	`if(y > 0) /* y == +Inf */`
4859	maechler	290	`return((x >= 1)? R_PosInf : 0.);`
6098	pd	291	`else /* y == -Inf */`
4859	maechler	292	`return((x < 1) ? R_PosInf : 0.);`
		293	`}`
4839	maechler	294	`}`
6098	pd	295	`return(R_NaN); /* all other cases: (-Inf)^{+-Inf,`
		296	`non-int}; (neg)^{+-Inf} */`
4834	maechler	297	`}`
4878	maechler	298	`#endif`
4834	maechler	299
7867	ripley	300	`double R_pow_di(double x, int n)`
		301	`{`
10866	maechler	302	`double xn = 1.0;`
7867	ripley	303
		304	`if (ISNAN(x)) return x;`
		305	`if (n == NA_INTEGER) return NA_REAL;`
		306	`if (n != 0) {`
		307	`if (!R_FINITE(x)) return R_pow(x, (double)n);`
		308	`if (n < 0) { n = -n; x = 1/x; }`
		309	`for(;;) {`
10866	maechler	310	`if(n & 01) xn *= x;`
7867	ripley	311	`if(n >>= 1) x *= x; else break;`
		312	`}`
		313	`}`
10866	maechler	314	`return xn;`
7867	ripley	315	`}`
		316
		317
5451	ihaka	318	`/* General Base Logarithms */`
2	r	319
7824	ripley	320	`static double logbase(double x, double base)`
834	maechler	321	`{`
4834	maechler	322	`return R_log(x) / log(base);`
2	r	323	`}`
		324
10777	luke	325	`SEXP R_unary(SEXP, SEXP, SEXP);`
		326	`SEXP R_binary(SEXP, SEXP, SEXP, SEXP);`
10718	maechler	327	`static SEXP integer_unary(ARITHOP_TYPE, SEXP);`
13997	duncan	328	`static SEXP real_unary(ARITHOP_TYPE, SEXP, SEXP);`
10718	maechler	329	`static SEXP real_binary(ARITHOP_TYPE, SEXP, SEXP);`
16033	luke	330	`static SEXP integer_binary(ARITHOP_TYPE, SEXP, SEXP, SEXP);`
2	r	331
13997	duncan	332	`#if 0`
2	r	333	`static int naflag;`
		334	`static SEXP lcall;`
13997	duncan	335	`#endif`
2	r	336
		337
1881	ihaka	338	`/* Unary and Binary Operators */`
2	r	339
		340	`SEXP do_arith(SEXP call, SEXP op, SEXP args, SEXP env)`
		341	`{`
571	ihaka	342	`SEXP ans;`
2	r	343
1881	ihaka	344	`if (DispatchGroup("Ops", call, op, args, env, &ans))`
571	ihaka	345	`return ans;`
2	r	346
571	ihaka	347	`switch (length(args)) {`
		348	`case 1:`
10777	luke	349	`return R_unary(call, op, CAR(args));`
571	ihaka	350	`case 2:`
10777	luke	351	`return R_binary(call, op, CAR(args), CADR(args));`
571	ihaka	352	`default:`
6203	maechler	353	`error("operator needs one or two arguments");`
571	ihaka	354	`}`
6098	pd	355	`return ans; /* never used; to keep -Wall happy */`
2	r	356	`}`
		357
24538	luke	358	`#define COERCE_IF_NEEDED(v, tp, vpi) do { \`
		359	`if (TYPEOF(v) != (tp)) { \`
		360	`int __vo__ = OBJECT(v); \`
		361	`REPROTECT(v = coerceVector(v, (tp)), vpi); \`
		362	`if (__vo__) SET_OBJECT(v, 1); \`
		363	`} \`
		364	`} while (0)`
2	r	365
24538	luke	366	`#define FIXUP_NULL_AND_CHECK_TYPES(v, vpi) do { \`
		367	`switch (TYPEOF(v)) { \`
		368	`case NILSXP: REPROTECT(v = allocVector(REALSXP,0), vpi); break; \`
		369	`case CPLXSXP: case REALSXP: case INTSXP: case LGLSXP: break; \`
		370	`default: errorcall(lcall, "non-numeric argument to binary operator"); \`
		371	`} \`
		372	`} while (0)`
		373
10777	luke	374	`SEXP R_binary(SEXP call, SEXP op, SEXP x, SEXP y)`
2	r	375	`{`
24538	luke	376	`SEXP class, dims, tsp, xnames, ynames, val;`
571	ihaka	377	`int mismatch, nx, ny, xarray, yarray, xts, yts;`
24538	luke	378	`int xattr, yattr;`
13997	duncan	379	`SEXP lcall = call;`
10777	luke	380	`PROTECT_INDEX xpi, ypi;`
24538	luke	381	`int nprotect = 2; /* x and y */`
2	r	382
		383
10777	luke	384	`PROTECT_WITH_INDEX(x, &xpi);`
		385	`PROTECT_WITH_INDEX(y, &ypi);`
		386
24538	luke	387	`FIXUP_NULL_AND_CHECK_TYPES(x, xpi);`
		388	`FIXUP_NULL_AND_CHECK_TYPES(y, ypi);`
2	r	389
24538	luke	390	`nx = LENGTH(x);`
		391	`if (ATTRIB(x) != R_NilValue) {`
		392	`xattr = TRUE;`
		393	`xarray = isArray(x);`
		394	`xts = isTs(x);`
1155	maechler	395	`}`
24538	luke	396	`else xarray = xts = xattr = FALSE;`
		397	`ny = LENGTH(y);`
		398	`if (ATTRIB(y) != R_NilValue) {`
		399	`yattr = TRUE;`
		400	`yarray = isArray(y);`
		401	`yts = isTs(y);`
		402	`}`
		403	`else yarray = yts = yattr = FALSE;`
2	r	404
3865	pd	405	`/* If either x or y is a matrix with length 1 and the other is a`
		406	`vector, we want to coerce the matrix to be a vector. */`
1010	ihaka	407
1881	ihaka	408	`/* FIXME: Danger Will Robinson.`
1155	maechler	409	`* ----- We might be trashing arguments here.`
		410	`* If we have NAMED(x) or NAMED(y) we should duplicate!`
		411	`*/`
1881	ihaka	412	`if (xarray != yarray) {`
24538	luke	413	`if (xarray && nx==1 && ny!=1) {`
10777	luke	414	`REPROTECT(x = duplicate(x), xpi);`
1010	ihaka	415	`setAttrib(x, R_DimSymbol, R_NilValue);`
834	maechler	416	`}`
24538	luke	417	`if (yarray && ny==1 && nx!=1) {`
10777	luke	418	`REPROTECT(y = duplicate(y), ypi);`
834	maechler	419	`setAttrib(y, R_DimSymbol, R_NilValue);`
785	pd	420	`}`
		421	`}`
		422
24538	luke	423	`mismatch = 0;`
571	ihaka	424	`if (xarray \|\| yarray) {`
		425	`if (xarray && yarray) {`
		426	`if (!conformable(x, y))`
5731	ripley	427	`errorcall(lcall, "non-conformable arrays");`
571	ihaka	428	`PROTECT(dims = getAttrib(x, R_DimSymbol));`
2	r	429	`}`
571	ihaka	430	`else if (xarray) {`
		431	`PROTECT(dims = getAttrib(x, R_DimSymbol));`
2	r	432	`}`
6098	pd	433	`else { /* (yarray) */`
571	ihaka	434	`PROTECT(dims = getAttrib(y, R_DimSymbol));`
2	r	435	`}`
24538	luke	436	`nprotect++;`
		437	`if (xattr) {`
		438	`PROTECT(xnames = getAttrib(x, R_DimNamesSymbol));`
		439	`nprotect++;`
		440	`}`
		441	`else xnames = R_NilValue;`
		442	`if (yattr) {`
		443	`PROTECT(ynames = getAttrib(y, R_DimNamesSymbol));`
		444	`nprotect++;`
		445	`}`
		446	`else ynames = R_NilValue;`
571	ihaka	447	`}`
		448	`else {`
24538	luke	449	`if (nx == ny \|\| nx == 1 \|\| ny == 1) mismatch = 0;`
		450	`else if (nx > 0 && ny > 0) {`
1881	ihaka	451	`if (nx > ny) mismatch = nx % ny;`
571	ihaka	452	`else mismatch = ny % nx;`
		453	`}`
24538	luke	454	`dims = R_NilValue;`
		455	`if (xattr) {`
		456	`PROTECT(xnames = getAttrib(x, R_NamesSymbol));`
		457	`nprotect++;`
		458	`}`
		459	`else xnames = R_NilValue;`
		460	`if (yattr) {`
		461	`PROTECT(ynames = getAttrib(y, R_NamesSymbol));`
		462	`nprotect++;`
		463	`}`
		464	`else ynames = R_NilValue;`
571	ihaka	465	`}`
2	r	466
571	ihaka	467	`if (xts \|\| yts) {`
		468	`if (xts && yts) {`
		469	`if (!tsConform(x, y))`
5731	ripley	470	`errorcall(lcall, "Non-conformable time-series");`
571	ihaka	471	`PROTECT(tsp = getAttrib(x, R_TspSymbol));`
		472	`PROTECT(class = getAttrib(x, R_ClassSymbol));`
2	r	473	`}`
571	ihaka	474	`else if (xts) {`
24538	luke	475	`if (nx < ny)`
571	ihaka	476	`ErrorMessage(lcall, ERROR_TSVEC_MISMATCH);`
		477	`PROTECT(tsp = getAttrib(x, R_TspSymbol));`
		478	`PROTECT(class = getAttrib(x, R_ClassSymbol));`
2	r	479	`}`
6203	maechler	480	`else { /* (yts) */`
24538	luke	481	`if (ny < nx)`
571	ihaka	482	`ErrorMessage(lcall, ERROR_TSVEC_MISMATCH);`
		483	`PROTECT(tsp = getAttrib(y, R_TspSymbol));`
		484	`PROTECT(class = getAttrib(y, R_ClassSymbol));`
2	r	485	`}`
24538	luke	486	`nprotect += 2;`
571	ihaka	487	`}`
24538	luke	488	`else class = tsp = NULL; /* -Wall */`
4719	maechler	489
3865	pd	490	`if (mismatch)`
4179	rgentlem	491	`warningcall(lcall, "longer object length\n\tis not a multiple of shorter object length");`
2	r	492
24538	luke	493	`/* need to preserve object here, as _binary copies class attributes /`
571	ihaka	494	`if (TYPEOF(x) == CPLXSXP \|\| TYPEOF(y) == CPLXSXP) {`
24538	luke	495	`COERCE_IF_NEEDED(x, CPLXSXP, xpi);`
		496	`COERCE_IF_NEEDED(y, CPLXSXP, ypi);`
		497	`val = complex_binary(PRIMVAL(op), x, y);`
571	ihaka	498	`}`
24538	luke	499	`else if (TYPEOF(x) == REALSXP \|\| TYPEOF(y) == REALSXP) {`
		500	`COERCE_IF_NEEDED(x, REALSXP, xpi);`
		501	`COERCE_IF_NEEDED(y, REALSXP, ypi);`
		502	`val = real_binary(PRIMVAL(op), x, y);`
		503	`}`
		504	`else val = integer_binary(PRIMVAL(op), x, y, lcall);`
2	r	505
24539	luke	506	`/* quick return if there are no attributes */`
24538	luke	507	`if (! xattr && ! yattr) {`
		508	`UNPROTECT(nprotect);`
		509	`return val;`
		510	`}`
		511
		512	`PROTECT(val);`
		513	`nprotect++;`
		514
6098	pd	515	`/* Don't set the dims if one argument is an array of size 0 and the`
6203	maechler	516	`other isn't of size zero, cos they're wrong */`
21144	tlumley	517	`/* Not if the other argument is a scalar (PR#1979) */`
571	ihaka	518	`if (dims != R_NilValue) {`
21144	tlumley	519	`if (!((xarray && (nx == 0) && (ny > 1)) \|\|`
		520	`(yarray && (ny == 0) && (nx > 1)))){`
24538	luke	521	`setAttrib(val, R_DimSymbol, dims);`
1881	ihaka	522	`if (xnames != R_NilValue)`
24538	luke	523	`setAttrib(val, R_DimNamesSymbol, xnames);`
1881	ihaka	524	`else if (ynames != R_NilValue)`
24538	luke	525	`setAttrib(val, R_DimNamesSymbol, ynames);`
1881	ihaka	526	`}`
571	ihaka	527	`}`
		528	`else {`
24538	luke	529	`if (length(val) == length(xnames))`
		530	`setAttrib(val, R_NamesSymbol, xnames);`
		531	`else if (length(val) == length(ynames))`
		532	`setAttrib(val, R_NamesSymbol, ynames);`
571	ihaka	533	`}`
		534
6098	pd	535	`if (xts \|\| yts) { /* must set after dims! */`
24538	luke	536	`setAttrib(val, R_TspSymbol, tsp);`
		537	`setAttrib(val, R_ClassSymbol, class);`
4562	pd	538	`}`
		539
24538	luke	540	`UNPROTECT(nprotect);`
		541	`return val;`
2	r	542	`}`
		543
10777	luke	544	`SEXP R_unary(SEXP call, SEXP op, SEXP s1)`
2	r	545	`{`
571	ihaka	546	`switch (TYPEOF(s1)) {`
		547	`case LGLSXP:`
		548	`case INTSXP:`
		549	`return integer_unary(PRIMVAL(op), s1);`
		550	`case REALSXP:`
13997	duncan	551	`return real_unary(PRIMVAL(op), s1, call);`
571	ihaka	552	`case CPLXSXP:`
		553	`return complex_unary(PRIMVAL(op), s1);`
		554	`default:`
13997	duncan	555	`errorcall(call, "Invalid argument to unary operator");`
571	ihaka	556	`}`
6098	pd	557	`return s1; /* never used; to keep -Wall happy */`
2	r	558	`}`
		559
10718	maechler	560	`static SEXP integer_unary(ARITHOP_TYPE code, SEXP s1)`
2	r	561	`{`
571	ihaka	562	`int i, n, x;`
		563	`SEXP ans;`
2	r	564
1881	ihaka	565	`switch (code) {`
571	ihaka	566	`case PLUSOP:`
		567	`return s1;`
		568	`case MINUSOP:`
		569	`ans = duplicate(s1);`
16336	rgentlem	570	`SET_TYPEOF(ans, INTSXP);`
571	ihaka	571	`n = LENGTH(s1);`
		572	`for (i = 0; i < n; i++) {`
		573	`x = INTEGER(s1)[i];`
		574	`INTEGER(ans)[i] = (x == NA_INTEGER) ?`
		575	`NA_INTEGER : ((x == 0.0) ? 0 : -x);`
2	r	576	`}`
571	ihaka	577	`return ans;`
		578	`default:`
5731	ripley	579	`error("illegal unary operator");`
571	ihaka	580	`}`
6098	pd	581	`return s1; /* never used; to keep -Wall happy */`
2	r	582	`}`
		583
13997	duncan	584	`static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall)`
2	r	585	`{`
571	ihaka	586	`int i, n;`
		587	`SEXP ans;`
2	r	588
1881	ihaka	589	`switch (code) {`
571	ihaka	590	`case PLUSOP: return s1;`
		591	`case MINUSOP:`
		592	`ans = duplicate(s1);`
		593	`n = LENGTH(s1);`
		594	`for (i = 0; i < n; i++) {`
		595	`#ifdef IEEE_754`
		596	`REAL(ans)[i] = -REAL(s1)[i];`
		597	`#else`
1400	maechler	598	`double x;`
571	ihaka	599	`x = REAL(s1)[i];`
580	ihaka	600	`REAL(ans)[i] = ISNA(x) ? NA_REAL :`
571	ihaka	601	`((x == 0.0) ? 0.0 : -x);`
		602	`#endif`
2	r	603	`}`
571	ihaka	604	`return ans;`
		605	`default:`
5731	ripley	606	`errorcall(lcall, "illegal unary operator");`
571	ihaka	607	`}`
6098	pd	608	`return s1; /* never used; to keep -Wall happy */`
2	r	609	`}`
988	maechler	610
90	ihaka	611	`/* i1 = i % n1; i2 = i % n2;`
834	maechler	612	`* this macro is quite a bit faster than having real modulo calls`
		613	`* in the loop (tested on Intel and Sparc)`
90	ihaka	614	`*/`
		615	`#define mod_iterate(n1,n2,i1,i2) for (i=i1=i2=0; i<n; \`
6218	pd	616	`i1 = (++i1 == n1) ? 0 : i1,\`
		617	`i2 = (++i2 == n2) ? 0 : i2,\`
		618	`++i)`
2	r	619
16033	luke	620
		621
		622	`/* The tests using integer comparisons are a bit faster than the tests`
		623	`using doubles, but they depend on a two's complement representation`
		624	`(but that is almost universal). The tests that compare results to`
		625	`double's depend on being able to accurately represent all int's as`
		626	`double's. Since int's are almost universally 32 bit that should be`
		627	`OK. */`
		628
16216	luke	629	`#ifndef INT_32_BITS`
		630	`/* configure checks whehter int is 32 bits. If not this code will`
		631	`need to be rewritten. Since 32 bit ints are pretty much universal,`
		632	`we can worry about writing alternate code when the need arises.`
		633	`To be safe, we signal a compiler error if int is not 32 bits. */`
		634	`# error code requires that int have 32 bits`
		635	`#else`
16365	luke	636	`/* Just to be on the safe side, configure ought to check that the`
		637	`mashine uses two's complement. A define like`
16033	luke	638	`#define USES_TWOS_COMPLEMENT (~0 == (unsigned) -1)`
16365	luke	639	`might work, but at least one compiler (CodeWarrior 6) chokes on it.`
		640	`So for now just assume it is true.`
		641	`*/`
		642	`#define USES_TWOS_COMPLEMENT 1`
16033	luke	643
		644	`#if USES_TWOS_COMPLEMENT`
		645	`# define OPPOSITE_SIGNS(x, y) ((x < 0) ^ (y < 0))`
		646	`# define GOODISUM(x, y, z) (((x) > 0) ? ((y) < (z)) : ! ((y) < (z)))`
		647	`# define GOODIDIFF(x, y, z) (!(OPPOSITE_SIGNS(x, y) && OPPOSITE_SIGNS(x, z)))`
		648	`#else`
		649	`# define GOODISUM(x, y, z) ((double) (x) + (double) (y) == (z))`
		650	`# define GOODIDIFF(x, y, z) ((double) (x) - (double) (y) == (z))`
		651	`#endif`
		652	`#define GOODIPROD(x, y, z) ((double) (x) * (double) (y) == (z))`
		653	`#define INTEGER_OVERFLOW_WARNING "NAs produced by integer overflow"`
16216	luke	654	`#endif`
16033	luke	655
		656	`static SEXP integer_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2, SEXP lcall)`
2	r	657	`{`
571	ihaka	658	`int i, i1, i2, n, n1, n2;`
		659	`int x1, x2;`
		660	`SEXP ans;`
16033	luke	661	`Rboolean naflag = FALSE;`
2	r	662
571	ihaka	663	`n1 = LENGTH(s1);`
		664	`n2 = LENGTH(s2);`
16613	ripley	665	`/* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */`
		666	`if (n1 == 0 \|\| n2 == 0) n = 0; else n = (n1 > n2) ? n1 : n2;`
2	r	667
571	ihaka	668	`if (code == DIVOP \|\| code == POWOP)`
		669	`ans = allocVector(REALSXP, n);`
		670	`else`
		671	`ans = allocVector(INTSXP, n);`
16613	ripley	672	`if (n1 == 0 \|\| n2 == 0) return(ans);`
		673	`/* note: code below was surely wrong in DIVOP and POWOP cases,`
		674	`since ans was a REALSXP.`
		675	`*/`
2	r	676
16613	ripley	677	`/* if (n1 < 1 \|\| n2 < 1) {`
571	ihaka	678	`for (i = 0; i < n; i++)`
		679	`INTEGER(ans)[i] = NA_INTEGER;`
		680	`return ans;`
16613	ripley	681	`} */`
2	r	682
571	ihaka	683	`switch (code) {`
		684	`case PLUSOP:`
		685	`mod_iterate(n1, n2, i1, i2) {`
		686	`x1 = INTEGER(s1)[i1];`
		687	`x2 = INTEGER(s2)[i2];`
		688	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		689	`INTEGER(ans)[i] = NA_INTEGER;`
16033	luke	690	`else {`
		691	`int val = x1 + x2;`
		692	`if (val != NA_INTEGER && GOODISUM(x1, x2, val))`
		693	`INTEGER(ans)[i] = val;`
		694	`else {`
		695	`INTEGER(ans)[i] = NA_INTEGER;`
		696	`naflag = TRUE;`
		697	`}`
		698	`}`
2	r	699	`}`
16033	luke	700	`if (naflag)`
		701	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
571	ihaka	702	`break;`
		703	`case MINUSOP:`
834	maechler	704	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	705	`x1 = INTEGER(s1)[i1];`
		706	`x2 = INTEGER(s2)[i2];`
		707	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		708	`INTEGER(ans)[i] = NA_INTEGER;`
16033	luke	709	`else {`
		710	`int val = x1 - x2;`
		711	`if (val != NA_INTEGER && GOODIDIFF(x1, x2, val))`
		712	`INTEGER(ans)[i] = val;`
		713	`else {`
		714	`naflag = TRUE;`
		715	`INTEGER(ans)[i] = NA_INTEGER;`
		716	`}`
		717	`}`
571	ihaka	718	`}`
16033	luke	719	`if (naflag)`
		720	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
571	ihaka	721	`break;`
		722	`case TIMESOP:`
834	maechler	723	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	724	`x1 = INTEGER(s1)[i1];`
		725	`x2 = INTEGER(s2)[i2];`
		726	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		727	`INTEGER(ans)[i] = NA_INTEGER;`
16033	luke	728	`else {`
		729	`int val = x1 * x2;`
		730	`if (val != NA_INTEGER && GOODIPROD(x1, x2, val))`
		731	`INTEGER(ans)[i] = val;`
		732	`else {`
		733	`naflag = TRUE;`
		734	`INTEGER(ans)[i] = NA_INTEGER;`
		735	`}`
		736	`}`
571	ihaka	737	`}`
16033	luke	738	`if (naflag)`
		739	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
571	ihaka	740	`break;`
		741	`case DIVOP:`
834	maechler	742	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	743	`x1 = INTEGER(s1)[i1];`
		744	`x2 = INTEGER(s2)[i2];`
575	ihaka	745	`#ifdef IEEE_754`
		746	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		747	`#else`
16613	ripley	748	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER \|\| x2 == 0)`
575	ihaka	749	`#endif`
16613	ripley	750	`REAL(ans)[i] = NA_REAL;`
		751	`else`
		752	`REAL(ans)[i] = (double) x1 / (double) x2;`
571	ihaka	753	`}`
		754	`break;`
		755	`case POWOP:`
834	maechler	756	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	757	`x1 = INTEGER(s1)[i1];`
		758	`x2 = INTEGER(s2)[i2];`
		759	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		760	`REAL(ans)[i] = NA_REAL;`
		761	`else {`
4834	maechler	762	`REAL(ans)[i] = MATH_CHECK(R_pow((double) x1, (double) x2));`
571	ihaka	763	`}`
		764	`}`
		765	`break;`
		766	`case MODOP:`
834	maechler	767	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	768	`x1 = INTEGER(s1)[i1];`
834	maechler	769	`x2 = INTEGER(s2)[i2];`
571	ihaka	770	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER \|\| x2 == 0)`
		771	`INTEGER(ans)[i] = NA_INTEGER;`
		772	`else {`
4719	maechler	773	`INTEGER(ans)[i] = /* till 0.63.2: x1 % x2 */`
		774	`(x1 >= 0 && x2 > 0) ? x1 % x2 :`
		775	`(int)myfmod((double)x1,(double)x2);`
571	ihaka	776	`}`
		777	`}`
		778	`break;`
		779	`case IDIVOP:`
834	maechler	780	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	781	`x1 = INTEGER(s1)[i1];`
834	maechler	782	`x2 = INTEGER(s2)[i2];`
571	ihaka	783	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		784	`INTEGER(ans)[i] = NA_INTEGER;`
1881	ihaka	785	`else if (x2 == 0)`
571	ihaka	786	`INTEGER(ans)[i] = 0;`
		787	`else`
		788	`INTEGER(ans)[i] = floor((double)x1 / (double)x2);`
		789	`}`
		790	`break;`
		791	`}`
1010	ihaka	792
18510	maechler	793
24539	luke	794	`/* quick return if there are no attributes */`
		795	`if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)`
		796	`return ans;`
		797
16613	ripley	798	`/* Copy attributes from longer argument. */`
1010	ihaka	799
		800	`if (n1 > n2)`
1155	maechler	801	`copyMostAttrib(s1, ans);`
1010	ihaka	802	`else if (n1 == n2) {`
		803	`copyMostAttrib(s2, ans);`
		804	`copyMostAttrib(s1, ans);`
		805	`}`
		806	`else`
		807	`copyMostAttrib(s2, ans);`
		808
571	ihaka	809	`return ans;`
2	r	810	`}`
		811
10718	maechler	812	`static SEXP real_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2)`
2	r	813	`{`
571	ihaka	814	`int i, i1, i2, n, n1, n2;`
		815	`SEXP ans;`
1400	maechler	816	`#ifndef IEEE_754`
		817	`double x1, x2;`
		818	`#endif`
2	r	819
1010	ihaka	820	`/* Note: "s1" and "s2" are protected above. */`
571	ihaka	821	`n1 = LENGTH(s1);`
		822	`n2 = LENGTH(s2);`
16613	ripley	823
		824	`/* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */`
		825	`if (n1 == 0 \|\| n2 == 0) return(allocVector(REALSXP, 0));`
		826
18510	maechler	827	`n = (n1 > n2) ? n1 : n2;`
571	ihaka	828	`ans = allocVector(REALSXP, n);`
2	r	829
16613	ripley	830	`/* if (n1 < 1 \|\| n2 < 1) {`
		831	`for (i = 0; i < n; i++)`
		832	`REAL(ans)[i] = NA_REAL;`
		833	`return ans;`
		834	`} */`
2	r	835
571	ihaka	836	`switch (code) {`
		837	`case PLUSOP:`
834	maechler	838	`mod_iterate(n1, n2, i1, i2) {`
571	ihaka	839	`#ifdef IEEE_754`
		840	`REAL(ans)[i] = REAL(s1)[i1] + REAL(s2)[i2];`
		841	`#else`
		842	`x1 = REAL(s1)[i1];`
		843	`x2 = REAL(s2)[i2];`
580	ihaka	844	`if (ISNA(x1) \|\| ISNA(x2))`
571	ihaka	845	`REAL(ans)[i] = NA_REAL;`
		846	`else`
		847	`REAL(ans)[i] = MATH_CHECK(x1 + x2);`
		848	`#endif`
2	r	849	`}`
571	ihaka	850	`break;`
		851	`case MINUSOP:`
		852	`mod_iterate(n1, n2, i1, i2) {`
		853	`#ifdef IEEE_754`
		854	`REAL(ans)[i] = REAL(s1)[i1] - REAL(s2)[i2];`
		855	`#else`
834	maechler	856	`x1 = REAL(s1)[i1];`
571	ihaka	857	`x2 = REAL(s2)[i2];`
580	ihaka	858	`if (ISNA(x1) \|\| ISNA(x2))`
571	ihaka	859	`REAL(ans)[i] = NA_REAL;`
		860	`else`
		861	`REAL(ans)[i] = MATH_CHECK(x1 - x2);`
		862	`#endif`
		863	`}`
		864	`break;`
		865	`case TIMESOP:`
		866	`mod_iterate(n1, n2, i1, i2) {`
		867	`#ifdef IEEE_754`
		868	`REAL(ans)[i] = REAL(s1)[i1] * REAL(s2)[i2];`
		869	`#else`
		870	`x1 = REAL(s1)[i1];`
		871	`x2 = REAL(s2)[i2];`
580	ihaka	872	`if (ISNA(x1) && ISNA(x2))`
571	ihaka	873	`REAL(ans)[i] = NA_REAL;`
		874	`else`
		875	`REAL(ans)[i] = MATH_CHECK(x1 * x2);`
		876	`#endif`
		877	`}`
		878	`break;`
		879	`case DIVOP:`
		880	`mod_iterate(n1, n2, i1, i2) {`
		881	`#ifdef IEEE_754`
		882	`REAL(ans)[i] = REAL(s1)[i1] / REAL(s2)[i2];`
		883	`#else`
		884	`x1 = REAL(s1)[i1];`
		885	`x2 = REAL(s2)[i2];`
580	ihaka	886	`if (ISNA(x1) \|\| ISNA(x2) \|\| x2 == 0)`
571	ihaka	887	`REAL(ans)[i] = NA_REAL;`
		888	`else`
		889	`REAL(ans)[i] = MATH_CHECK(x1 / x2);`
		890	`#endif`
		891	`}`
		892	`break;`
		893	`case POWOP:`
		894	`mod_iterate(n1, n2, i1, i2) {`
		895	`#ifdef IEEE_754`
4834	maechler	896	`REAL(ans)[i] = R_pow(REAL(s1)[i1], REAL(s2)[i2]);`
571	ihaka	897	`#else`
		898	`x1 = REAL(s1)[i1];`
		899	`x2 = REAL(s2)[i2];`
580	ihaka	900	`if (ISNA(x1) \|\| ISNA(x2))`
571	ihaka	901	`REAL(ans)[i] = NA_REAL;`
		902	`else`
4834	maechler	903	`REAL(ans)[i] = MATH_CHECK(R_pow(x1, x2));`
571	ihaka	904	`#endif`
		905	`}`
		906	`break;`
		907	`case MODOP:`
		908	`mod_iterate(n1, n2, i1, i2) {`
		909	`#ifdef IEEE_754`
		910	`REAL(ans)[i] = myfmod(REAL(s1)[i1], REAL(s2)[i2]);`
		911	`#else`
		912	`x1 = REAL(s1)[i1];`
		913	`x2 = REAL(s2)[i2];`
580	ihaka	914	`if (ISNA(x1) \|\| ISNA(x2) \|\| x2 == 0)`
571	ihaka	915	`REAL(ans)[i] = NA_REAL;`
		916	`else`
		917	`REAL(ans)[i] = MATH_CHECK(myfmod(x1, x2));`
		918	`#endif`
		919	`}`
		920	`break;`
		921	`case IDIVOP:`
		922	`mod_iterate(n1, n2, i1, i2) {`
		923	`#ifdef IEEE_754`
		924	`REAL(ans)[i] = floor(REAL(s1)[i1] / REAL(s2)[i2]);`
		925	`#else`
		926	`x1 = REAL(s1)[i1];`
		927	`x2 = REAL(s2)[i2];`
580	ihaka	928	`if (ISNA(x1) \|\| ISNA(x2))`
571	ihaka	929	`REAL(ans)[i] = NA_REAL;`
		930	`else {`
1881	ihaka	931	`if (x2 == 0)`
16613	ripley	932	`REAL(ans)[i] = 0;`
571	ihaka	933	`else`
		934	`REAL(ans)[i] = MATH_CHECK(floor(x1 / x2));`
		935	`}`
		936	`#endif`
		937	`}`
		938	`break;`
		939	`}`
1010	ihaka	940
24539	luke	941	`/* quick return if there are no attributes */`
		942	`if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)`
		943	`return ans;`
18510	maechler	944
16613	ripley	945	`/* Copy attributes from longer argument. */`
1010	ihaka	946
		947	`if (n1 > n2)`
		948	`copyMostAttrib(s1, ans);`
		949	`else if (n1 == n2) {`
		950	`copyMostAttrib(s2, ans);`
		951	`copyMostAttrib(s1, ans);`
		952	`}`
		953	`else`
		954	`copyMostAttrib(s2, ans);`
		955
571	ihaka	956	`return ans;`
2	r	957	`}`
		958
7615	maechler	959
1881	ihaka	960	`/* Mathematical Functions of One Argument */`
2	r	961
13997	duncan	962	`static SEXP math1(SEXP sa, double(*f)(), SEXP lcall)`
2	r	963	`{`
571	ihaka	964	`SEXP sy;`
		965	`double y, a;`
21984	ripley	966	`int i, n, sao = OBJECT(sa);`
13997	duncan	967	`int naflag;`
2	r	968
7615	maechler	969	`if (!isNumeric(sa))`
9147	maechler	970	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	971
		972	`n = length(sa);`
21984	ripley	973	`/* coercion can lose the object bit */`
7615	maechler	974	`PROTECT(sa = coerceVector(sa, REALSXP));`
		975	`PROTECT(sy = allocVector(REALSXP, n));`
		976	`a = REAL(sa);`
		977	`y = REAL(sy);`
		978	`naflag = 0;`
		979	`for (i = 0; i < n; i++) {`
		980	`if (ISNAN(a[i]))`
		981	`y[i] = a[i];`
		982	`else {`
		983	`y[i] = MATH_CHECK(f(a[i]));`
		984	`if (ISNAN(y[i])) naflag = 1;`
2	r	985	`}`
571	ihaka	986	`}`
7615	maechler	987	`if(naflag)`
9147	maechler	988	`warningcall(lcall, R_MSG_NA);`
7615	maechler	989
10172	luke	990	`SET_ATTRIB(sy, duplicate(ATTRIB(sa)));`
21984	ripley	991	`SET_OBJECT(sy, sao);`
7615	maechler	992	`UNPROTECT(2);`
		993	`return sy;`
2	r	994	`}`
		995
7615	maechler	996
2	r	997	`SEXP do_math1(SEXP call, SEXP op, SEXP args, SEXP env)`
		998	`{`
571	ihaka	999	`SEXP s;`
2	r	1000
1881	ihaka	1001	`if (DispatchGroup("Math", call, op, args, env, &s))`
571	ihaka	1002	`return s;`
2	r	1003
22644	ripley	1004	`checkArity(op, args);`
		1005
571	ihaka	1006	`if (isComplex(CAR(args)))`
		1007	`return complex_math1(call, op, args, env);`
2	r	1008
13997	duncan	1009	`#define MATH1(x) math1(CAR(args), x, call);`
571	ihaka	1010	`switch (PRIMVAL(op)) {`
13997	duncan	1011	`case 1: return MATH1(floor);`
		1012	`case 2: return MATH1(ceil);`
		1013	`case 3: return MATH1(sqrt);`
		1014	`case 4: return MATH1(sign);`
		1015	`case 5: return MATH1(trunc);`
2	r	1016
13997	duncan	1017	`case 10: return MATH1(exp);`
18510	maechler	1018	`case 11: return MATH1(expm1);`
13997	duncan	1019	`case 12: return MATH1(log1p);`
		1020	`case 20: return MATH1(cos);`
		1021	`case 21: return MATH1(sin);`
		1022	`case 22: return MATH1(tan);`
		1023	`case 23: return MATH1(acos);`
		1024	`case 24: return MATH1(asin);`
14701	ripley	1025
13997	duncan	1026	`case 30: return MATH1(cosh);`
		1027	`case 31: return MATH1(sinh);`
		1028	`case 32: return MATH1(tanh);`
		1029	`case 33: return MATH1(acosh);`
		1030	`case 34: return MATH1(asinh);`
		1031	`case 35: return MATH1(atanh);`
14701	ripley	1032
13997	duncan	1033	`case 40: return MATH1(lgammafn);`
		1034	`case 41: return MATH1(gammafn);`
14701	ripley	1035
13997	duncan	1036	`case 42: return MATH1(digamma);`
		1037	`case 43: return MATH1(trigamma);`
		1038	`case 44: return MATH1(tetragamma);`
		1039	`case 45: return MATH1(pentagamma);`
14701	ripley	1040
13997	duncan	1041	`case 46: return MATH1(gamma_cody);`
2	r	1042
571	ihaka	1043	`default:`
7615	maechler	1044	`errorcall(call, "unimplemented real function (of 1 arg.)");`
571	ihaka	1045	`}`
6098	pd	1046	`return s; /* never used; to keep -Wall happy */`
2	r	1047	`}`
		1048
18454	ripley	1049	`SEXP do_abs(SEXP call, SEXP op, SEXP args, SEXP env)`
		1050	`{`
		1051	`SEXP s;`
		1052	`if (DispatchGroup("Math", call, op, args, env, &s))`
		1053	`return s;`
18466	ripley	1054	`return do_cmathfuns(call, op, args, env);`
18454	ripley	1055	`}`
		1056
7615	maechler	1057	`/* Mathematical Functions of Two Numeric Arguments (plus 1 int) */`
		1058
		1059	`#define if_NA_Math2_set(y,a,b) \`
		1060	`if (ISNA (a) \|\| ISNA (b)) y = NA_REAL; \`
		1061	`else if (ISNAN(a) \|\| ISNAN(b)) y = R_NaN;`
		1062
13997	duncan	1063	`static SEXP math2(SEXP sa, SEXP sb, double (*f)(), SEXP lcall)`
2	r	1064	`{`
571	ihaka	1065	`SEXP sy;`
21984	ripley	1066	`int i, ia, ib, n, na, nb, sao = OBJECT(sa), sbo = OBJECT(sb);`
571	ihaka	1067	`double ai, bi, a, b, *y;`
13997	duncan	1068	`int naflag;`
2	r	1069
571	ihaka	1070	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1071	`errorcall(lcall, R_MSG_NONNUM_MATH);`
2	r	1072
27167	maechler	1073	`/* for 0-length a we want the attributes of a, not those of b`
27095	ripley	1074	`as no recycling will occur */`
7615	maechler	1075	`#define SETUP_Math2 \`
		1076	`na = LENGTH(sa); \`
		1077	`nb = LENGTH(sb); \`
27095	ripley	1078	`if ((na == 0) \|\| (nb == 0)) { \`
		1079	`PROTECT(sy = allocVector(REALSXP, 0)); \`
		1080	`if (na == 0) { \`
		1081	`SET_ATTRIB(sy, duplicate(ATTRIB(sa)));\`
		1082	`SET_OBJECT(sy, sao); \`
		1083	`} \`
		1084	`UNPROTECT(1); \`
		1085	`return(sy); \`
		1086	`} \`
7615	maechler	1087	`n = (na < nb) ? nb : na; \`
		1088	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1089	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1090	`PROTECT(sy = allocVector(REALSXP, n)); \`
		1091	`a = REAL(sa); \`
		1092	`b = REAL(sb); \`
		1093	`y = REAL(sy); \`
		1094	`naflag = 0`
		1095
		1096	`SETUP_Math2;`
		1097
6098	pd	1098	`mod_iterate(na, nb, ia, ib) {`
		1099	`ai = a[ia];`
		1100	`bi = b[ib];`
7615	maechler	1101	`if_NA_Math2_set(y[i], ai, bi)`
6098	pd	1102	`else {`
		1103	`y[i] = MATH_CHECK(f(ai, bi));`
7615	maechler	1104	`if (ISNAN(y[i])) naflag = 1;`
2	r	1105	`}`
571	ihaka	1106	`}`
7615	maechler	1107
		1108	`#define FINISH_Math2 \`
		1109	`if(naflag) \`
9147	maechler	1110	`warningcall(lcall, R_MSG_NA); \`
7615	maechler	1111	`\`
		1112	`if (n == na) { \`
10172	luke	1113	`SET_ATTRIB(sy, duplicate(ATTRIB(sa))); \`
21984	ripley	1114	`SET_OBJECT(sy, sao); \`
7615	maechler	1115	`} \`
		1116	`else if (n == nb) { \`
10172	luke	1117	`SET_ATTRIB(sy, duplicate(ATTRIB(sb))); \`
21984	ripley	1118	`SET_OBJECT(sy, sbo); \`
7615	maechler	1119	`} \`
		1120	`UNPROTECT(3)`
		1121
		1122	`FINISH_Math2;`
		1123
		1124	`return sy;`
		1125	`} /* math2() */`
		1126
13997	duncan	1127	`static SEXP math2_1(SEXP sa, SEXP sb, SEXP sI, double (*f)(), SEXP lcall)`
7615	maechler	1128	`{`
		1129	`SEXP sy;`
21984	ripley	1130	`int i, ia, ib, n, na, nb, sao = OBJECT(sa), sbo = OBJECT(sb);`
7615	maechler	1131	`double ai, bi, a, b, *y;`
		1132	`int m_opt;`
13997	duncan	1133	`int naflag;`
7615	maechler	1134
		1135	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1136	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	1137
		1138	`SETUP_Math2;`
		1139	`m_opt = asInteger(sI);`
		1140
		1141	`mod_iterate(na, nb, ia, ib) {`
		1142	`ai = a[ia];`
		1143	`bi = b[ib];`
		1144	`if_NA_Math2_set(y[i], ai, bi)`
		1145	`else {`
		1146	`y[i] = MATH_CHECK(f(ai, bi, m_opt));`
		1147	`if (ISNAN(y[i])) naflag = 1;`
		1148	`}`
571	ihaka	1149	`}`
7615	maechler	1150	`FINISH_Math2;`
		1151	`return sy;`
		1152	`} /* math2_1() */`
		1153
13997	duncan	1154	`static SEXP math2_2(SEXP sa, SEXP sb, SEXP sI1, SEXP sI2, double (*f)(), SEXP lcall)`
7615	maechler	1155	`{`
		1156	`SEXP sy;`
21984	ripley	1157	`int i, ia, ib, n, na, nb, sao = OBJECT(sa), sbo = OBJECT(sb);`
7615	maechler	1158	`double ai, bi, a, b, *y;`
		1159	`int i_1, i_2;`
13997	duncan	1160	`int naflag;`
7615	maechler	1161	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1162	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	1163
		1164	`SETUP_Math2;`
		1165	`i_1 = asInteger(sI1);`
		1166	`i_2 = asInteger(sI2);`
		1167
		1168	`mod_iterate(na, nb, ia, ib) {`
		1169	`ai = a[ia];`
		1170	`bi = b[ib];`
		1171	`if_NA_Math2_set(y[i], ai, bi)`
		1172	`else {`
		1173	`y[i] = MATH_CHECK(f(ai, bi, i_1, i_2));`
		1174	`if (ISNAN(y[i])) naflag = 1;`
		1175	`}`
571	ihaka	1176	`}`
7615	maechler	1177	`FINISH_Math2;`
571	ihaka	1178	`return sy;`
7615	maechler	1179	`} /* math2_2() */`
2	r	1180
13997	duncan	1181	`#define Math2(A, FUN) math2(CAR(A), CADR(A), FUN, call);`
		1182	`#define Math2_1(A, FUN) math2_1(CAR(A), CADR(A), CADDR(A), FUN, call);`
		1183	`#define Math2_2(A, FUN) math2_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call)`
2	r	1184
		1185	`SEXP do_math2(SEXP call, SEXP op, SEXP args, SEXP env)`
		1186	`{`
571	ihaka	1187	`checkArity(op, args);`
2	r	1188
571	ihaka	1189	`if (isComplex(CAR(args)))`
		1190	`return complex_math2(call, op, args, env);`
2	r	1191
14701	ripley	1192
571	ihaka	1193	`switch (PRIMVAL(op)) {`
13997	duncan	1194
7615	maechler	1195	`case 0: return Math2(args, atan2);`
22644	ripley	1196	`case 10001: return Math2(args, rround);`
		1197	`case 10004: return Math2(args, prec);`
2	r	1198
7615	maechler	1199	`case 2: return Math2(args, lbeta);`
		1200	`case 3: return Math2(args, beta);`
		1201	`case 4: return Math2(args, lchoose);`
		1202	`case 5: return Math2(args, choose);`
2	r	1203
7671	maechler	1204	`case 6: return Math2_1(args, dchisq);`
		1205	`case 7: return Math2_2(args, pchisq);`
		1206	`case 8: return Math2_2(args, qchisq);`
14701	ripley	1207
7671	maechler	1208	`case 9: return Math2_1(args, dexp);`
		1209	`case 10: return Math2_2(args, pexp);`
		1210	`case 11: return Math2_2(args, qexp);`
14701	ripley	1211
7671	maechler	1212	`case 12: return Math2_1(args, dgeom);`
		1213	`case 13: return Math2_2(args, pgeom);`
		1214	`case 14: return Math2_2(args, qgeom);`
14701	ripley	1215
7615	maechler	1216	`case 15: return Math2_1(args, dpois);`
		1217	`case 16: return Math2_2(args, ppois);`
		1218	`case 17: return Math2_2(args, qpois);`
14701	ripley	1219
7671	maechler	1220	`case 18: return Math2_1(args, dt);`
		1221	`case 19: return Math2_2(args, pt);`
		1222	`case 20: return Math2_2(args, qt);`
14701	ripley	1223
7671	maechler	1224	`case 21: return Math2_1(args, dsignrank);`
		1225	`case 22: return Math2_2(args, psignrank);`
		1226	`case 23: return Math2_2(args, qsignrank);`
2239	hornik	1227
7615	maechler	1228	`case 24: return Math2(args, bessel_j);`
		1229	`case 25: return Math2(args, bessel_y);`
2632	maechler	1230
571	ihaka	1231	`default:`
7615	maechler	1232	`errorcall(call, "unimplemented real function of 2 numeric arg.s");`
571	ihaka	1233	`}`
6098	pd	1234	`return op; /* never used; to keep -Wall happy */`
2	r	1235	`}`
		1236
7615	maechler	1237
2	r	1238	`SEXP do_atan(SEXP call, SEXP op, SEXP args, SEXP env)`
		1239	`{`
571	ihaka	1240	`SEXP s;`
		1241	`int n;`
1881	ihaka	1242	`if (DispatchGroup("Math", call, op, args, env, &s))`
571	ihaka	1243	`return s;`
13997	duncan	1244
1881	ihaka	1245	`switch (n = length(args)) {`
571	ihaka	1246	`case 1:`
		1247	`if (isComplex(CAR(args)))`
		1248	`return complex_math1(call, op, args, env);`
		1249	`else`
13997	duncan	1250	`return math1(CAR(args), atan, call);`
571	ihaka	1251	`case 2:`
		1252	`if (isComplex(CAR(args)) \|\| isComplex(CDR(args)))`
		1253	`return complex_math2(call, op, args, env);`
		1254	`else`
13997	duncan	1255	`return math2(CAR(args), CADR(args), atan2, call);`
571	ihaka	1256	`default:`
5731	ripley	1257	`error("%d arguments passed to \"atan\" which requires 1 or 2", n);`
571	ihaka	1258	`}`
6098	pd	1259	`return s; /* never used; to keep -Wall happy */`
2	r	1260	`}`
		1261
22644	ripley	1262
		1263	`/* The S4 Math2 group, round and signif */`
		1264	`SEXP do_Math2(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1265	`{`
22644	ripley	1266	`SEXP a;`
1820	ihaka	1267	`if (DispatchGroup("Math", call, op, args, env, &a))`
		1268	`return a;`
13997	duncan	1269
22644	ripley	1270	`checkArity(op, args);`
		1271	`if (length(CADR(args)) == 0)`
		1272	`errorcall(call, "illegal 2nd arg of length 0");`
		1273	`return do_math2(call, op, args, env);`
2	r	1274	`}`
		1275
		1276	`SEXP do_log(SEXP call, SEXP op, SEXP args, SEXP env)`
		1277	`{`
571	ihaka	1278	`SEXP s;`
		1279	`int n;`
1881	ihaka	1280	`if (DispatchGroup("Math", call, op, args, env, &s))`
571	ihaka	1281	`return s;`
13997	duncan	1282
1881	ihaka	1283	`switch (n = length(args)) {`
571	ihaka	1284	`case 1:`
		1285	`if (isComplex(CAR(args)))`
		1286	`return complex_math1(call, op, args, env);`
		1287	`else`
13997	duncan	1288	`return math1(CAR(args), R_log, call);`
571	ihaka	1289	`case 2:`
6098	pd	1290	`if (length(CADR(args)) == 0)`
6203	maechler	1291	`errorcall(call, "illegal 2nd arg of length 0");`
571	ihaka	1292	`if (isComplex(CAR(args)) \|\| isComplex(CDR(args)))`
		1293	`return complex_math2(call, op, args, env);`
		1294	`else`
13997	duncan	1295	`return math2(CAR(args), CADR(args), logbase, call);`
571	ihaka	1296	`default:`
5731	ripley	1297	`error("%d arguments passed to \"log\" which requires 1 or 2", n);`
571	ihaka	1298	`}`
6098	pd	1299	`return s; /* never used; to keep -Wall happy */`
2	r	1300	`}`
		1301
13997	duncan	1302
7615	maechler	1303	`/* Mathematical Functions of Three (Real) Arguments */`
		1304
		1305	`#define if_NA_Math3_set(y,a,b,c) \`
		1306	`if (ISNA (a) \|\| ISNA (b)\|\| ISNA (c)) y = NA_REAL; \`
		1307	`else if (ISNAN(a) \|\| ISNAN(b)\|\| ISNAN(c)) y = R_NaN;`
		1308
90	ihaka	1309	`#define mod_iterate3(n1,n2,n3,i1,i2,i3) for (i=i1=i2=i3=0; i<n; \`
7615	maechler	1310	`i1 = (++i1==n1) ? 0 : i1, \`
		1311	`i2 = (++i2==n2) ? 0 : i2, \`
		1312	`i3 = (++i3==n3) ? 0 : i3, \`
		1313	`++i)`
90	ihaka	1314
13997	duncan	1315	`static SEXP math3(SEXP sa, SEXP sb, SEXP sc, double (*f)(), SEXP lcall)`
2	r	1316	`{`
571	ihaka	1317	`SEXP sy;`
25788	maechler	1318	`int i, ia, ib, ic, n, na, nb, nc,`
21984	ripley	1319	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc);`
571	ihaka	1320	`double ai, bi, ci, a, b, c, y;`
13997	duncan	1321	`int naflag;`
2	r	1322
7615	maechler	1323	`#define SETUP_Math3 \`
		1324	`if (!isNumeric(sa) \|\| !isNumeric(sb) \|\| !isNumeric(sc)) \`
9147	maechler	1325	`errorcall(lcall, R_MSG_NONNUM_MATH); \`
7615	maechler	1326	`\`
		1327	`na = LENGTH(sa); \`
		1328	`nb = LENGTH(sb); \`
		1329	`nc = LENGTH(sc); \`
		1330	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0)) \`
		1331	`return(allocVector(REALSXP, 0)); \`
		1332	`n = na; \`
		1333	`if (n < nb) n = nb; \`
		1334	`if (n < nc) n = nc; \`
		1335	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1336	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1337	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		1338	`PROTECT(sy = allocVector(REALSXP, n)); \`
		1339	`a = REAL(sa); \`
		1340	`b = REAL(sb); \`
		1341	`c = REAL(sc); \`
		1342	`y = REAL(sy); \`
		1343	`naflag = 0`
834	maechler	1344
7615	maechler	1345	`SETUP_Math3;`
		1346
6098	pd	1347	`mod_iterate3 (na, nb, nc, ia, ib, ic) {`
		1348	`ai = a[ia];`
		1349	`bi = b[ib];`
		1350	`ci = c[ic];`
7615	maechler	1351	`if_NA_Math3_set(y[i], ai,bi,ci)`
6098	pd	1352	`else {`
		1353	`y[i] = MATH_CHECK(f(ai, bi, ci));`
7615	maechler	1354	`if (ISNAN(y[i])) naflag = 1;`
2	r	1355	`}`
571	ihaka	1356	`}`
7615	maechler	1357
		1358	`#define FINISH_Math3 \`
		1359	`if(naflag) \`
9147	maechler	1360	`warningcall(lcall, R_MSG_NA); \`
7615	maechler	1361	`\`
		1362	`if (n == na) { \`
10172	luke	1363	`SET_ATTRIB(sy, duplicate(ATTRIB(sa))); \`
21984	ripley	1364	`SET_OBJECT(sy, sao); \`
7615	maechler	1365	`} \`
		1366	`else if (n == nb) { \`
10172	luke	1367	`SET_ATTRIB(sy, duplicate(ATTRIB(sb))); \`
21984	ripley	1368	`SET_OBJECT(sy, sbo); \`
7615	maechler	1369	`} \`
		1370	`else if (n == nc) { \`
10172	luke	1371	`SET_ATTRIB(sy, duplicate(ATTRIB(sc))); \`
21984	ripley	1372	`SET_OBJECT(sy, sco); \`
7615	maechler	1373	`} \`
		1374	`UNPROTECT(4)`
		1375
		1376	`FINISH_Math3;`
		1377
		1378	`return sy;`
		1379	`} /* math3 */`
		1380
13997	duncan	1381	`static SEXP math3_1(SEXP sa, SEXP sb, SEXP sc, SEXP sI, double (*f)(), SEXP lcall)`
7615	maechler	1382	`{`
		1383	`SEXP sy;`
25788	maechler	1384	`int i, ia, ib, ic, n, na, nb, nc,`
21984	ripley	1385	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc);`
7615	maechler	1386	`double ai, bi, ci, a, b, c, y;`
		1387	`int i_1;`
13997	duncan	1388	`int naflag;`
7615	maechler	1389
		1390	`SETUP_Math3;`
		1391	`i_1 = asInteger(sI);`
		1392
		1393	`mod_iterate3 (na, nb, nc, ia, ib, ic) {`
		1394	`ai = a[ia];`
		1395	`bi = b[ib];`
		1396	`ci = c[ic];`
		1397	`if_NA_Math3_set(y[i], ai,bi,ci)`
		1398	`else {`
		1399	`y[i] = MATH_CHECK(f(ai, bi, ci, i_1));`
		1400	`if (ISNAN(y[i])) naflag = 1;`
		1401	`}`
571	ihaka	1402	`}`
7615	maechler	1403
		1404	`FINISH_Math3;`
		1405	`return sy;`
		1406	`} /* math3_1 */`
		1407
13997	duncan	1408	`static SEXP math3_2(SEXP sa, SEXP sb, SEXP sc, SEXP sI, SEXP sJ, double (*f)(), SEXP lcall)`
7615	maechler	1409	`{`
		1410	`SEXP sy;`
25788	maechler	1411	`int i, ia, ib, ic, n, na, nb, nc,`
21984	ripley	1412	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc);`
7615	maechler	1413	`double ai, bi, ci, a, b, c, y;`
		1414	`int i_1,i_2;`
13997	duncan	1415	`int naflag;`
7615	maechler	1416
		1417	`SETUP_Math3;`
		1418	`i_1 = asInteger(sI);`
		1419	`i_2 = asInteger(sJ);`
		1420
		1421	`mod_iterate3 (na, nb, nc, ia, ib, ic) {`
		1422	`ai = a[ia];`
		1423	`bi = b[ib];`
		1424	`ci = c[ic];`
		1425	`if_NA_Math3_set(y[i], ai,bi,ci)`
		1426	`else {`
		1427	`y[i] = MATH_CHECK(f(ai, bi, ci, i_1, i_2));`
		1428	`if (ISNAN(y[i])) naflag = 1;`
		1429	`}`
571	ihaka	1430	`}`
7615	maechler	1431
		1432	`FINISH_Math3;`
571	ihaka	1433	`return sy;`
7615	maechler	1434	`} /* math3_2 */`
2	r	1435
13997	duncan	1436	`#define Math3(A, FUN) math3 (CAR(A), CADR(A), CADDR(A), FUN, call);`
		1437	`#define Math3_1(A, FUN) math3_1(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call);`
		1438	`#define Math3_2(A, FUN) math3_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), CAD4R(A), FUN, call)`
2	r	1439
		1440	`SEXP do_math3(SEXP call, SEXP op, SEXP args, SEXP env)`
		1441	`{`
571	ihaka	1442	`checkArity(op, args);`
2	r	1443
571	ihaka	1444	`switch (PRIMVAL(op)) {`
2	r	1445
7671	maechler	1446	`case 1: return Math3_1(args, dbeta);`
		1447	`case 2: return Math3_2(args, pbeta);`
		1448	`case 3: return Math3_2(args, qbeta);`
2	r	1449
7671	maechler	1450	`case 4: return Math3_1(args, dbinom);`
		1451	`case 5: return Math3_2(args, pbinom);`
		1452	`case 6: return Math3_2(args, qbinom);`
2	r	1453
7671	maechler	1454	`case 7: return Math3_1(args, dcauchy);`
		1455	`case 8: return Math3_2(args, pcauchy);`
		1456	`case 9: return Math3_2(args, qcauchy);`
2	r	1457
7671	maechler	1458	`case 10: return Math3_1(args, df);`
		1459	`case 11: return Math3_2(args, pf);`
		1460	`case 12: return Math3_2(args, qf);`
2	r	1461
7671	maechler	1462	`case 13: return Math3_1(args, dgamma);`
		1463	`case 14: return Math3_2(args, pgamma);`
		1464	`case 15: return Math3_2(args, qgamma);`
2	r	1465
7639	maechler	1466	`case 16: return Math3_1(args, dlnorm);`
		1467	`case 17: return Math3_2(args, plnorm);`
		1468	`case 18: return Math3_2(args, qlnorm);`
2	r	1469
7671	maechler	1470	`case 19: return Math3_1(args, dlogis);`
		1471	`case 20: return Math3_2(args, plogis);`
		1472	`case 21: return Math3_2(args, qlogis);`
2	r	1473
7671	maechler	1474	`case 22: return Math3_1(args, dnbinom);`
		1475	`case 23: return Math3_2(args, pnbinom);`
		1476	`case 24: return Math3_2(args, qnbinom);`
2	r	1477
7639	maechler	1478	`case 25: return Math3_1(args, dnorm);`
		1479	`case 26: return Math3_2(args, pnorm);`
		1480	`case 27: return Math3_2(args, qnorm);`
2	r	1481
7671	maechler	1482	`case 28: return Math3_1(args, dunif);`
		1483	`case 29: return Math3_2(args, punif);`
		1484	`case 30: return Math3_2(args, qunif);`
2	r	1485
7671	maechler	1486	`case 31: return Math3_1(args, dweibull);`
		1487	`case 32: return Math3_2(args, pweibull);`
		1488	`case 33: return Math3_2(args, qweibull);`
2	r	1489
7671	maechler	1490	`case 34: return Math3_1(args, dnchisq);`
		1491	`case 35: return Math3_2(args, pnchisq);`
		1492	`case 36: return Math3_2(args, qnchisq);`
2	r	1493
7671	maechler	1494	`case 37: return Math3_1(args, dnt);`
		1495	`case 38: return Math3_2(args, pnt);`
602	ihaka	1496	`#ifdef UNIMP`
7671	maechler	1497	`case 39: return Math3_2(args, qnt);`
602	ihaka	1498	`#endif`
		1499
7671	maechler	1500	`case 40: return Math3_1(args, dwilcox);`
		1501	`case 41: return Math3_2(args, pwilcox);`
		1502	`case 42: return Math3_2(args, qwilcox);`
1276	hornik	1503
7615	maechler	1504	`case 43: return Math3(args, bessel_i);`
		1505	`case 44: return Math3(args, bessel_k);`
2632	maechler	1506
		1507
571	ihaka	1508	`default:`
7615	maechler	1509	`errorcall(call, "unimplemented real function of 3 numeric arg.s");`
571	ihaka	1510	`}`
6098	pd	1511	`return op; /* never used; to keep -Wall happy */`
4719	maechler	1512	`} /* do_math3() */`
988	maechler	1513
7615	maechler	1514	`/* Mathematical Functions of Four (Real) Arguments */`
		1515
		1516	`#define if_NA_Math4_set(y,a,b,c,d) \`
		1517	`if (ISNA (a)\|\| ISNA (b)\|\| ISNA (c)\|\| ISNA (d)) y = NA_REAL;\`
		1518	`else if (ISNAN(a)\|\| ISNAN(b)\|\| ISNAN(c)\|\| ISNAN(d)) y = R_NaN;`
		1519
90	ihaka	1520	`#define mod_iterate4(n1,n2,n3,n4,i1,i2,i3,i4) for (i=i1=i2=i3=i4=0; i<n; \`
7615	maechler	1521	`i1 = (++i1==n1) ? 0 : i1, \`
		1522	`i2 = (++i2==n2) ? 0 : i2, \`
		1523	`i3 = (++i3==n3) ? 0 : i3, \`
		1524	`i4 = (++i4==n4) ? 0 : i4, \`
		1525	`++i)`
2	r	1526
13997	duncan	1527	`static SEXP math4(SEXP sa, SEXP sb, SEXP sc, SEXP sd, double (*f)(), SEXP lcall)`
2	r	1528	`{`
571	ihaka	1529	`SEXP sy;`
25788	maechler	1530	`int i, ia, ib, ic, id, n, na, nb, nc, nd,`
21984	ripley	1531	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc), sdo = OBJECT(sd);`
571	ihaka	1532	`double ai, bi, ci, di, a, b, c, d, *y;`
13997	duncan	1533	`int naflag;`
2	r	1534
7615	maechler	1535	`#define SETUP_Math4 \`
		1536	`if(!isNumeric(sa)\|\| !isNumeric(sb)\|\| !isNumeric(sc)\|\| !isNumeric(sd))\`
9147	maechler	1537	`errorcall(lcall, R_MSG_NONNUM_MATH); \`
7615	maechler	1538	`\`
		1539	`na = LENGTH(sa); \`
		1540	`nb = LENGTH(sb); \`
		1541	`nc = LENGTH(sc); \`
		1542	`nd = LENGTH(sd); \`
		1543	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0) \|\| (nd == 0)) \`
		1544	`return(allocVector(REALSXP, 0)); \`
		1545	`n = na; \`
		1546	`if (n < nb) n = nb; \`
		1547	`if (n < nc) n = nc; \`
		1548	`if (n < nd) n = nd; \`
		1549	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1550	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1551	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		1552	`PROTECT(sd = coerceVector(sd, REALSXP)); \`
		1553	`PROTECT(sy = allocVector(REALSXP, n)); \`
		1554	`a = REAL(sa); \`
		1555	`b = REAL(sb); \`
		1556	`c = REAL(sc); \`
		1557	`d = REAL(sd); \`
		1558	`y = REAL(sy); \`
		1559	`naflag = 0`
834	maechler	1560
7615	maechler	1561	`SETUP_Math4;`
		1562
6098	pd	1563	`mod_iterate4 (na, nb, nc, nd, ia, ib, ic, id) {`
		1564	`ai = a[ia];`
		1565	`bi = b[ib];`
		1566	`ci = c[ic];`
		1567	`di = d[id];`
7615	maechler	1568	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
6098	pd	1569	`else {`
		1570	`y[i] = MATH_CHECK(f(ai, bi, ci, di));`
7615	maechler	1571	`if (ISNAN(y[i])) naflag = 1;`
2	r	1572	`}`
571	ihaka	1573	`}`
7615	maechler	1574
		1575	`#define FINISH_Math4 \`
		1576	`if(naflag) \`
9147	maechler	1577	`warningcall(lcall, R_MSG_NA); \`
7615	maechler	1578	`\`
		1579	`if (n == na) { \`
10172	luke	1580	`SET_ATTRIB(sy, duplicate(ATTRIB(sa))); \`
21984	ripley	1581	`SET_OBJECT(sy, sao); \`
7615	maechler	1582	`} \`
		1583	`else if (n == nb) { \`
10172	luke	1584	`SET_ATTRIB(sy, duplicate(ATTRIB(sb))); \`
21984	ripley	1585	`SET_OBJECT(sy, sbo); \`
7615	maechler	1586	`} \`
		1587	`else if (n == nc) { \`
10172	luke	1588	`SET_ATTRIB(sy, duplicate(ATTRIB(sc))); \`
21984	ripley	1589	`SET_OBJECT(sy, sco); \`
7615	maechler	1590	`} \`
		1591	`else if (n == nd) { \`
10172	luke	1592	`SET_ATTRIB(sy, duplicate(ATTRIB(sd))); \`
21984	ripley	1593	`SET_OBJECT(sy, sdo); \`
7615	maechler	1594	`} \`
		1595	`UNPROTECT(5)`
		1596
		1597	`FINISH_Math4;`
		1598
		1599	`return sy;`
		1600	`} /* math4() */`
		1601
13997	duncan	1602	`static SEXP math4_1(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, double (*f)(), SEXP lcall)`
7615	maechler	1603	`{`
		1604	`SEXP sy;`
25788	maechler	1605	`int i, ia, ib, ic, id, n, na, nb, nc, nd,`
21984	ripley	1606	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc), sdo = OBJECT(sd);`
7615	maechler	1607	`double ai, bi, ci, di, a, b, c, d, *y;`
		1608	`int i_1;`
13997	duncan	1609	`int naflag;`
7615	maechler	1610
		1611	`SETUP_Math4;`
		1612	`i_1 = asInteger(sI);`
		1613
		1614	`mod_iterate4 (na, nb, nc, nd, ia, ib, ic, id) {`
		1615	`ai = a[ia];`
		1616	`bi = b[ib];`
		1617	`ci = c[ic];`
		1618	`di = d[id];`
		1619	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
		1620	`else {`
		1621	`y[i] = MATH_CHECK(f(ai, bi, ci, di, i_1));`
		1622	`if (ISNAN(y[i])) naflag = 1;`
		1623	`}`
571	ihaka	1624	`}`
7615	maechler	1625	`FINISH_Math4;`
		1626	`return sy;`
		1627	`} /* math4_1() */`
		1628
		1629	`static SEXP math4_2(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, SEXP sJ,`
13997	duncan	1630	`double (*f)(), SEXP lcall)`
7615	maechler	1631	`{`
		1632	`SEXP sy;`
25788	maechler	1633	`int i, ia, ib, ic, id, n, na, nb, nc, nd,`
21984	ripley	1634	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc), sdo = OBJECT(sd);`
7615	maechler	1635	`double ai, bi, ci, di, a, b, c, d, *y;`
		1636	`int i_1, i_2;`
13997	duncan	1637	`int naflag;`
7615	maechler	1638
		1639	`SETUP_Math4;`
		1640	`i_1 = asInteger(sI);`
		1641	`i_2 = asInteger(sJ);`
		1642
		1643	`mod_iterate4 (na, nb, nc, nd, ia, ib, ic, id) {`
		1644	`ai = a[ia];`
		1645	`bi = b[ib];`
		1646	`ci = c[ic];`
		1647	`di = d[id];`
		1648	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
		1649	`else {`
		1650	`y[i] = MATH_CHECK(f(ai, bi, ci, di, i_1, i_2));`
		1651	`if (ISNAN(y[i])) naflag = 1;`
		1652	`}`
571	ihaka	1653	`}`
7615	maechler	1654	`FINISH_Math4;`
571	ihaka	1655	`return sy;`
7615	maechler	1656	`} /* math4_2() */`
2	r	1657
		1658
7615	maechler	1659	`#define CAD3R CADDDR`
		1660	`/* This is not (yet) in Rinternals.h : */`
7671	maechler	1661	`#define CAD5R(e) CAR(CDR(CDR(CDR(CDR(CDR(e))))))`
4719	maechler	1662
13997	duncan	1663	`#define Math4(A, FUN) math4 (CAR(A), CADR(A), CADDR(A), CAD3R(A), FUN, call)`
7671	maechler	1664	`#define Math4_1(A, FUN) math4_1(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \`
13997	duncan	1665	`FUN, call)`
7671	maechler	1666	`#define Math4_2(A, FUN) math4_2(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \`
13997	duncan	1667	`CAD5R(A), FUN, call)`
7615	maechler	1668
7671	maechler	1669
2	r	1670	`SEXP do_math4(SEXP call, SEXP op, SEXP args, SEXP env)`
		1671	`{`
571	ihaka	1672	`checkArity(op, args);`
2	r	1673
13997	duncan	1674
571	ihaka	1675	`switch (PRIMVAL(op)) {`
7671	maechler	1676
		1677	`/* Completely dummy for -Wall -- math4() at all! : */`
		1678	`case -99: return Math4(args, dhyper);`
		1679
		1680	`case 1: return Math4_1(args, dhyper);`
		1681	`case 2: return Math4_2(args, phyper);`
		1682	`case 3: return Math4_2(args, qhyper);`
		1683
		1684	`case 4: return Math4_1(args, dnbeta);`
		1685	`case 5: return Math4_2(args, pnbeta);`
602	ihaka	1686	`#ifdef UNIMP`
7671	maechler	1687	`case 6: return Math4_2(args, qnbeta);`
602	ihaka	1688	`#endif`
		1689	`#ifdef UNIMP`
7671	maechler	1690	`case 7: return Math4_1(args, dnf);`
602	ihaka	1691	`#endif`
7671	maechler	1692	`case 8: return Math4_2(args, pnf);`
602	ihaka	1693	`#ifdef UNIMP`
7671	maechler	1694	`case 9: return Math4_2(args, qnf);`
602	ihaka	1695	`#endif`
623	ihaka	1696	`#ifdef UNIMP`
7671	maechler	1697	`case 10: return Math4_1(args, dtukey);`
1350	ihaka	1698	`#endif`
7671	maechler	1699	`case 11: return Math4_2(args, ptukey);`
		1700	`case 12: return Math4_2(args, qtukey);`
571	ihaka	1701	`default:`
7615	maechler	1702	`errorcall(call, "unimplemented real function of 4 numeric arg.s");`
571	ihaka	1703	`}`
6098	pd	1704	`return op; /* never used; to keep -Wall happy */`
2	r	1705	`}`
4719	maechler	1706
		1707
10866	maechler	1708	`#ifdef WHEN_MATH5_IS_THERE/* as in ./arithmetic.h */`
4719	maechler	1709
7615	maechler	1710	`/* Mathematical Functions of Five (Real) Arguments */`
		1711
		1712	`#define if_NA_Math5_set(y,a,b,c,d,e) \`
		1713	`if (ISNA (a)\|\| ISNA (b)\|\| ISNA (c)\|\| ISNA (d)\|\| ISNA (e)) \`
		1714	`y = NA_REAL; \`
		1715	`else if(ISNAN(a)\|\| ISNAN(b)\|\| ISNAN(c)\|\| ISNAN(d)\|\| ISNAN(e)) \`
		1716	`y = R_NaN;`
		1717
		1718	`#define mod_iterate5(n1,n2,n3,n4,n5, i1,i2,i3,i4,i5) \`
		1719	`for (i=i1=i2=i3=i4=i5=0; i<n; \`
		1720	`i1 = (++i1==n1) ? 0 : i1, \`
		1721	`i2 = (++i2==n2) ? 0 : i2, \`
		1722	`i3 = (++i3==n3) ? 0 : i3, \`
		1723	`i4 = (++i4==n4) ? 0 : i4, \`
		1724	`i5 = (++i5==n5) ? 0 : i5, \`
		1725	`++i)`
		1726
		1727	`static SEXP math5(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP se, double (*f)())`
4719	maechler	1728	`{`
		1729	`SEXP sy;`
25788	maechler	1730	`int i, ia, ib, ic, id, ie, n, na, nb, nc, nd, ne,`
		1731	`sao = OBJECT(sa), sbo = OBJECT(sb), sco = OBJECT(sc),`
21984	ripley	1732	`sdo = OBJECT(sd), seo = OBJECT(se);`
4719	maechler	1733	`double ai, bi, ci, di, ei, a, b, c, d, e, y;`
		1734
7615	maechler	1735	`#define SETUP_Math5 \`
		1736	`if (!isNumeric(sa) \|\| !isNumeric(sb) \|\| !isNumeric(sc) \|\| \`
		1737	`!isNumeric(sd) \|\| !isNumeric(se)) \`
9147	maechler	1738	`errorcall(lcall, R_MSG_NONNUM_MATH); \`
7615	maechler	1739	`\`
		1740	`na = LENGTH(sa); \`
		1741	`nb = LENGTH(sb); \`
		1742	`nc = LENGTH(sc); \`
		1743	`nd = LENGTH(sd); \`
		1744	`ne = LENGTH(se); \`
		1745	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0) \|\| (nd == 0) \|\| (ne == 0)) \`
		1746	`return(allocVector(REALSXP, 0)); \`
		1747	`n = na; \`
		1748	`if (n < nb) n = nb; \`
		1749	`if (n < nc) n = nc; \`
		1750	`if (n < nd) n = nd; \`
		1751	`if (n < ne) n = ne; /* n = max(na,nb,nc,nd,ne) */ \`
		1752	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1753	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1754	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		1755	`PROTECT(sd = coerceVector(sd, REALSXP)); \`
		1756	`PROTECT(se = coerceVector(se, REALSXP)); \`
		1757	`PROTECT(sy = allocVector(REALSXP, n)); \`
		1758	`a = REAL(sa); \`
		1759	`b = REAL(sb); \`
		1760	`c = REAL(sc); \`
		1761	`d = REAL(sd); \`
		1762	`e = REAL(se); \`
		1763	`y = REAL(sy); \`
		1764	`naflag = 0`
4719	maechler	1765
7615	maechler	1766	`SETUP_Math5;`
		1767
6098	pd	1768	`mod_iterate5 (na, nb, nc, nd, ne,`
		1769	`ia, ib, ic, id, ie) {`
		1770	`ai = a[ia];`
		1771	`bi = b[ib];`
		1772	`ci = c[ic];`
		1773	`di = d[id];`
		1774	`ei = e[ie];`
7615	maechler	1775	`if_NA_Math5_set(y[i], ai,bi,ci,di,ei)`
6098	pd	1776	`else {`
		1777	`y[i] = MATH_CHECK(f(ai, bi, ci, di, ei));`
7615	maechler	1778	`if (ISNAN(y[i])) naflag = 1;`
4719	maechler	1779	`}`
		1780	`}`
7615	maechler	1781
		1782	`#define FINISH_Math5 \`
		1783	`if(naflag) \`
9147	maechler	1784	`warningcall(lcall, R_MSG_NA); \`
7615	maechler	1785	`\`
		1786	`if (n == na) { \`
10172	luke	1787	`SET_ATTRIB(sy, duplicate(ATTRIB(sa))); \`
21984	ripley	1788	`SET_OBJECT(sy, sao); \`
7615	maechler	1789	`} \`
		1790	`else if (n == nb) { \`
10172	luke	1791	`SET_ATTRIB(sy, duplicate(ATTRIB(sb))); \`
21984	ripley	1792	`SET_OBJECT(sy, sbo); \`
7615	maechler	1793	`} \`
		1794	`else if (n == nc) { \`
10172	luke	1795	`SET_ATTRIB(sy, duplicate(ATTRIB(sc))); \`
21984	ripley	1796	`SET_OBJECT(sy, sco); \`
7615	maechler	1797	`} \`
		1798	`else if (n == nd) { \`
10172	luke	1799	`SET_ATTRIB(sy, duplicate(ATTRIB(sd))); \`
21984	ripley	1800	`SET_OBJECT(sy, sdo); \`
7615	maechler	1801	`} \`
		1802	`else if (n == ne) { \`
10172	luke	1803	`SET_ATTRIB(sy, duplicate(ATTRIB(se))); \`
21984	ripley	1804	`SET_OBJECT(sy, seo); \`
7615	maechler	1805	`} \`
		1806	`UNPROTECT(6)`
		1807
		1808	`FINISH_Math5;`
		1809
4719	maechler	1810	`return sy;`
		1811	`} /* math5() */`
		1812
7615	maechler	1813	`#define Math5(A, FUN) \`
		1814	`math5(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), FUN);`
4719	maechler	1815
		1816	`SEXP do_math5(SEXP call, SEXP op, SEXP args, SEXP env)`
		1817	`{`
		1818	`checkArity(op, args);`
4859	maechler	1819	`lcall = call;`
4719	maechler	1820
		1821	`switch (PRIMVAL(op)) {`
7671	maechler	1822
		1823	`/* Completely dummy for -Wall -- use math5() at all! : */`
		1824	`case -99: return Math5(args, dhyper);`
4719	maechler	1825	`#ifdef UNIMP`
7671	maechler	1826	`case 2: return Math5(args, p...);`
		1827	`case 3: return Math5(args, q...);`
4719	maechler	1828	`#endif`
		1829	`default:`
7615	maechler	1830	`errorcall(call, "unimplemented real function of 5 numeric arg.s");`
4719	maechler	1831	`}`
6098	pd	1832	`return op; /* never used; to keep -Wall happy */`
4719	maechler	1833	`} /* do_math5() */`
		1834
6098	pd	1835	`#endif /* Math5 is there */`

The R Project SVN R

(root)/trunk/src/main/arithmetic.c – Rev 28362