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`
77193	maechler	3	`* Copyright (C) 1998--2019 The R Core Team.`
		4	`* Copyright (C) 2003--2019 The R Foundation`
71206	maechler	5	`* Copyright (C) 1995--1997 Robert Gentleman and Ross Ihaka`
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	`*`
42307	ripley	17	`* You should have received a copy of the GNU General Public License`
		18	`* along with this program; if not, a copy is available at`
68947	ripley	19	`* https://www.R-project.org/Licenses/`
2	r	20	`*/`
		21
5187	hornik	22	`#ifdef HAVE_CONFIG_H`
7701	hornik	23	`#include <config.h>`
5187	hornik	24	`#endif`
		25
77193	maechler	26	`// LDBL_EPSILON`
		27	`#include <float.h>`
60166	ripley	28
		29	`/* interval at which to check interrupts, a guess */`
		30	`#define NINTERRUPT 10000000`
		31
		32
14701	ripley	33	`#ifdef __OpenBSD__`
		34	`/* for definition of "struct exception" in math.h */`
		35	`# define __LIBM_PRIVATE`
		36	`#endif`
60667	ripley	37	`#include <Defn.h> /-> Arith.h -> math.h /`
14701	ripley	38	`#ifdef __OpenBSD__`
		39	`# undef __LIBM_PRIVATE`
		40	`#endif`
		41
60667	ripley	42	`#include <Internal.h>`
		43
49076	ripley	44	`#define R_MSG_NA _("NaNs produced")`
60844	ripley	45	`#define R_MSG_NONNUM_MATH _("non-numeric argument to mathematical function")`
49076	ripley	46
11499	ripley	47	`#include <Rmath.h>`
42422	ripley	48
63918	luke	49	`#include <R_ext/Itermacros.h>`
		50
2632	maechler	51	`#include "arithmetic.h"`
2	r	52
42398	ripley	53	`#include <errno.h>`
		54
2	r	55	`#ifdef HAVE_MATHERR`
		56
36491	ripley	57	`/* Override the SVID matherr function:`
		58	`the main difference here is not to print warnings.`
		59	`*/`
39866	duncan	60	`#ifndef __cplusplus`
2	r	61	`int matherr(struct exception *exc)`
		62	`{`
571	ihaka	63	`switch (exc->type) {`
		64	`case DOMAIN:`
		65	`case SING:`
		66	`errno = EDOM;`
		67	`break;`
		68	`case OVERFLOW:`
		69	`errno = ERANGE;`
		70	`break;`
		71	`case UNDERFLOW:`
		72	`exc->retval = 0.0;`
		73	`break;`
37954	maechler	74	`/*`
36491	ripley	75	`There are cases TLOSS and PLOSS which are ignored here.`
		76	`According to the Solaris man page, there are for`
		77	`trigonometric algorithms and not needed for good ones.`
		78	`*/`
571	ihaka	79	`}`
		80	`return 1;`
2	r	81	`}`
		82	`#endif`
39866	duncan	83	`#endif`
2	r	84
580	ihaka	85	`typedef union`
		86	`{`
16704	ripley	87	`double value;`
		88	`unsigned int word[2];`
580	ihaka	89	`} ieee_double;`
		90
42416	ripley	91	`/* gcc had problems with static const on AIX and Solaris`
		92	`Solaris was for gcc 3.1 and 3.2 under -O2 32-bit on 64-bit kernel */`
21057	ripley	93	`#ifdef _AIX`
		94	`#define CONST`
		95	`#elif defined(sparc) && defined (__GNUC__) && __GNUC__ == 3`
		96	`#define CONST`
14632	tlumley	97	`#else`
21057	ripley	98	`#define CONST const`
14632	tlumley	99	`#endif`
580	ihaka	100
21057	ripley	101	`#ifdef WORDS_BIGENDIAN`
		102	`static CONST int hw = 0;`
		103	`static CONST int lw = 1;`
		104	`#else /* !WORDS_BIGENDIAN */`
		105	`static CONST int hw = 1;`
		106	`static CONST int lw = 0;`
		107	`#endif /* WORDS_BIGENDIAN */`
13997	duncan	108
21057	ripley	109
580	ihaka	110	`static double R_ValueOfNA(void)`
		111	`{`
64950	ripley	112	`/* The gcc shipping with Fedora 9 gets this wrong without`
23894	pd	113	`* the volatile declaration. Thanks to Marc Schwartz. */`
23886	pd	114	`volatile ieee_double x;`
1881	ihaka	115	`x.word[hw] = 0x7ff00000;`
		116	`x.word[lw] = 1954;`
		117	`return x.value;`
580	ihaka	118	`}`
		119
		120	`int R_IsNA(double x)`
		121	`{`
3888	hornik	122	`if (isnan(x)) {`
1881	ihaka	123	`ieee_double y;`
		124	`y.value = x;`
		125	`return (y.word[lw] == 1954);`
		126	`}`
		127	`return 0;`
580	ihaka	128	`}`
1096	ihaka	129
		130	`int R_IsNaN(double x)`
		131	`{`
3888	hornik	132	`if (isnan(x)) {`
1881	ihaka	133	`ieee_double y;`
		134	`y.value = x;`
		135	`return (y.word[lw] != 1954);`
		136	`}`
		137	`return 0;`
1096	ihaka	138	`}`
7655	ripley	139
37954	maechler	140	`/* ISNAN uses isnan, which is undefined by C++ headers`
32538	tlumley	141	`This workaround is called only when ISNAN() is used`
		142	`in a user code in a file with __cplusplus defined */`
		143
		144	`int R_isnancpp(double x)`
		145	`{`
		146	`return (isnan(x)!=0);`
		147	`}`
		148
		149
42424	ripley	150	`/* Mainly for use in packages */`
7655	ripley	151	`int R_finite(double x)`
		152	`{`
29441	ripley	153	`#ifdef HAVE_WORKING_ISFINITE`
		154	`return isfinite(x);`
15252	hornik	155	`#else`
7655	ripley	156	`return (!isnan(x) & (x != R_PosInf) & (x != R_NegInf));`
		157	`#endif`
		158	`}`
		159
		160
1881	ihaka	161	`/* Arithmetic Initialization */`
580	ihaka	162
37001	ripley	163	`void attribute_hidden InitArithmetic()`
2	r	164	`{`
571	ihaka	165	`R_NaInt = INT_MIN;`
64230	ripley	166	`R_NaReal = R_ValueOfNA();`
		167	`// we assume C99, so`
		168	`#ifndef OLD`
		169	`R_NaN = NAN;`
		170	`R_PosInf = INFINITY;`
		171	`R_NegInf = -INFINITY;`
		172	`#else`
571	ihaka	173	`R_NaN = 0.0/R_Zero_Hack;`
		174	`R_PosInf = 1.0/R_Zero_Hack;`
		175	`R_NegInf = -1.0/R_Zero_Hack;`
64230	ripley	176	`#endif`
2	r	177	`}`
		178
77193	maechler	179	`#if HAVE_LONG_DOUBLE && (SIZEOF_LONG_DOUBLE > SIZEOF_DOUBLE)`
		180	`static LDOUBLE q_1_eps = 1 / LDBL_EPSILON;`
		181	`#else`
		182	`static double q_1_eps = 1 / DBL_EPSILON;`
		183	`#endif`
		184
		185	`/* Keep myfmod() and myfloor() in step */`
4847	maechler	186	`static double myfmod(double x1, double x2)`
		187	`{`
32158	ripley	188	`if (x2 == 0.0) return R_NaN;`
77193	maechler	189	`if(fabs(x2) > q_1_eps && R_FINITE(x1) && fabs(x1) <= fabs(x2)) {`
		190	`return`
77199	maechler	191	`(fabs(x1) == fabs(x2)) ? 0 :`
77193	maechler	192	`((x1 < 0 && x2 > 0) \|\|`
77199	maechler	193	`(x2 < 0 && x1 > 0))`
		194	`? x1+x2 // differing signs`
		195	`: x1 ; // "same" signs (incl. 0)`
77193	maechler	196	`}`
		197	`double q = x1 / x2;`
		198	`if(R_FINITE(q) && (fabs(q) > q_1_eps))`
32860	ripley	199	`warning(_("probable complete loss of accuracy in modulus"));`
77193	maechler	200	`LDOUBLE tmp = (LDOUBLE)x1 - floor(q) * (LDOUBLE)x2;`
77199	maechler	201	`return (double) (tmp - floorl(tmp/x2) * x2);`
4847	maechler	202	`}`
2	r	203
34301	ripley	204	`static double myfloor(double x1, double x2)`
		205	`{`
77193	maechler	206	`double q = x1 / x2;`
		207	`if (x2 == 0.0 \|\| fabs(q) > q_1_eps \|\| !R_FINITE(q))`
		208	`return q;`
		209	`if(fabs(q) < 1)`
		210	`return (q < 0) ? -1`
		211	`: ((x1 < 0 && x2 > 0) \|\|`
		212	`(x1 > 0 && x2 < 0) // differing signs`
		213	`? -1 : 0);`
		214	`LDOUBLE tmp = (LDOUBLE)x1 - floor(q) * (LDOUBLE)x2;`
77199	maechler	215	`return (double) (floor(q) + floorl(tmp/x2));`
34301	ripley	216	`}`
		217
4839	maechler	218	`double R_pow(double x, double y) /* = x ^ y */`
4834	maechler	219	`{`
52937	luke	220	`/* squaring is the most common of the specially handled cases so`
		221	`check for it first. */`
		222	`if(y == 2.0)`
		223	`return x * x;`
4839	maechler	224	`if(x == 1. \|\| y == 0.)`
		225	`return(1.);`
7449	maechler	226	`if(x == 0.) {`
		227	`if(y > 0.) return(0.);`
47051	jmc	228	`else if(y < 0) return(R_PosInf);`
		229	`else return(y); /* NA or NaN, we assert */`
7449	maechler	230	`}`
49062	maechler	231	`if (R_FINITE(x) && R_FINITE(y)) {`
64671	ripley	232	`/* There was a special case for y == 0.5 here, but`
		233	`gcc 4.3.0 -g -O2 mis-compiled it. Showed up with`
		234	`100^0.5 as 3.162278, example(pbirthday) failed. */`
69416	murdoch	235	`#ifdef USE_POWL_IN_R_POW`
74314	ripley	236	`// this is used only on 64-bit Windows (so has powl).`
		237	`return powl(x, y);`
69416	murdoch	238	`#else`
74314	ripley	239	`return pow(x, y);`
69416	murdoch	240	`#endif`
49062	maechler	241	`}`
29410	ripley	242	`if (ISNAN(x) \|\| ISNAN(y))`
4839	maechler	243	`return(x + y);`
5107	maechler	244	`if(!R_FINITE(x)) {`
6098	pd	245	`if(x > 0) /* Inf ^ y */`
41631	ripley	246	`return (y < 0.)? 0. : R_PosInf;`
6098	pd	247	`else { /* (-Inf) ^ y */`
5107	maechler	248	`if(R_FINITE(y) && y == floor(y)) /* (-Inf) ^ n */`
69128	ripley	249	`return (y < 0.) ? 0. : (myfmod(y, 2.) != 0 ? x : -x);`
4839	maechler	250	`}`
		251	`}`
5107	maechler	252	`if(!R_FINITE(y)) {`
4859	maechler	253	`if(x >= 0) {`
6098	pd	254	`if(y > 0) /* y == +Inf */`
41631	ripley	255	`return (x >= 1) ? R_PosInf : 0.;`
6098	pd	256	`else /* y == -Inf */`
41631	ripley	257	`return (x < 1) ? R_PosInf : 0.;`
4859	maechler	258	`}`
4839	maechler	259	`}`
64671	ripley	260	`return R_NaN; // all other cases: (-Inf)^{+-Inf, non-int}; (neg)^{+-Inf}`
4834	maechler	261	`}`
		262
7867	ripley	263	`double R_pow_di(double x, int n)`
		264	`{`
10866	maechler	265	`double xn = 1.0;`
7867	ripley	266
		267	`if (ISNAN(x)) return x;`
		268	`if (n == NA_INTEGER) return NA_REAL;`
58014	maechler	269
7867	ripley	270	`if (n != 0) {`
52938	luke	271	`if (!R_FINITE(x)) return R_POW(x, (double)n);`
58014	maechler	272
		273	`Rboolean is_neg = (n < 0);`
		274	`if(is_neg) n = -n;`
7867	ripley	275	`for(;;) {`
10866	maechler	276	`if(n & 01) xn *= x;`
7867	ripley	277	`if(n >>= 1) x *= x; else break;`
		278	`}`
68923	ripley	279	`if(is_neg) xn = 1. / xn;`
7867	ripley	280	`}`
10866	maechler	281	`return xn;`
7867	ripley	282	`}`
		283
		284
5451	ihaka	285	`/* General Base Logarithms */`
2	r	286
10777	luke	287	`SEXP R_unary(SEXP, SEXP, SEXP);`
		288	`SEXP R_binary(SEXP, SEXP, SEXP, SEXP);`
63327	ripley	289	`static SEXP logical_unary(ARITHOP_TYPE, SEXP, SEXP);`
41713	ripley	290	`static SEXP integer_unary(ARITHOP_TYPE, SEXP, SEXP);`
13997	duncan	291	`static SEXP real_unary(ARITHOP_TYPE, SEXP, SEXP);`
10718	maechler	292	`static SEXP real_binary(ARITHOP_TYPE, SEXP, SEXP);`
16033	luke	293	`static SEXP integer_binary(ARITHOP_TYPE, SEXP, SEXP, SEXP);`
2	r	294
13997	duncan	295	`#if 0`
2	r	296	`static int naflag;`
		297	`static SEXP lcall;`
13997	duncan	298	`#endif`
2	r	299
63578	luke	300	`/* Integer arithmetic support */`
2	r	301
63578	luke	302	`/* The tests using integer comparisons are a bit faster than the tests`
		303	`using doubles, but they depend on a two's complement representation`
		304	`(but that is almost universal). The tests that compare results to`
		305	`double's depend on being able to accurately represent all int's as`
		306	`double's. Since int's are almost universally 32 bit that should be`
		307	`OK. */`
		308
		309	`#ifndef INT_32_BITS`
		310	`/* configure checks whether int is 32 bits. If not this code will`
		311	`need to be rewritten. Since 32 bit ints are pretty much universal,`
		312	`we can worry about writing alternate code when the need arises.`
		313	`To be safe, we signal a compiler error if int is not 32 bits. */`
		314	`# error code requires that int have 32 bits`
65682	ripley	315	`#endif`
63578	luke	316
		317	`#define INTEGER_OVERFLOW_WARNING _("NAs produced by integer overflow")`
		318
63582	luke	319	`#define CHECK_INTEGER_OVERFLOW(call, ans, naflag) do { \`
		320	`if (naflag) { \`
		321	`PROTECT(ans); \`
		322	`warningcall(call, INTEGER_OVERFLOW_WARNING); \`
		323	`UNPROTECT(1); \`
		324	`} \`
		325	`} while(0)`
		326
65526	maechler	327	`#define R_INT_MAX INT_MAX`
		328	`#define R_INT_MIN -INT_MAX`
		329	`// .. relying on fact that NA_INTEGER is outside of these`
		330
63578	luke	331	`static R_INLINE int R_integer_plus(int x, int y, Rboolean *pnaflag)`
		332	`{`
		333	`if (x == NA_INTEGER \|\| y == NA_INTEGER)`
		334	`return NA_INTEGER;`
65526	maechler	335
		336	`if (((y > 0) && (x > (R_INT_MAX - y))) \|\|`
		337	`((y < 0) && (x < (R_INT_MIN - y)))) {`
		338	`if (pnaflag != NULL)`
		339	`*pnaflag = TRUE;`
		340	`return NA_INTEGER;`
64171	maechler	341	`}`
65526	maechler	342	`return x + y;`
63578	luke	343	`}`
		344
		345	`static R_INLINE int R_integer_minus(int x, int y, Rboolean *pnaflag)`
		346	`{`
		347	`if (x == NA_INTEGER \|\| y == NA_INTEGER)`
		348	`return NA_INTEGER;`
65682	ripley	349
		350	`if (((y < 0) && (x > (R_INT_MAX + y))) \|\|`
		351	`((y > 0) && (x < (R_INT_MIN + y)))) {`
		352	`if (pnaflag != NULL)`
		353	`*pnaflag = TRUE;`
		354	`return NA_INTEGER;`
63578	luke	355	`}`
65682	ripley	356	`return x - y;`
63578	luke	357	`}`
		358
65682	ripley	359	`#define GOODIPROD(x, y, z) ((double) (x) * (double) (y) == (z))`
63578	luke	360	`static R_INLINE int R_integer_times(int x, int y, Rboolean *pnaflag)`
		361	`{`
		362	`if (x == NA_INTEGER \|\| y == NA_INTEGER)`
		363	`return NA_INTEGER;`
		364	`else {`
69111	ripley	365	`int z = x * y; // UBSAN will warn if this overflows (happens in bda)`
63578	luke	366	`if (GOODIPROD(x, y, z) && z != NA_INTEGER)`
		367	`return z;`
		368	`else {`
		369	`if (pnaflag != NULL)`
		370	`*pnaflag = TRUE;`
		371	`return NA_INTEGER;`
		372	`}`
		373	`}`
		374	`}`
		375
		376	`static R_INLINE double R_integer_divide(int x, int y)`
		377	`{`
		378	`if (x == NA_INTEGER \|\| y == NA_INTEGER)`
		379	`return NA_REAL;`
		380	`else`
		381	`return (double) x / (double) y;`
64171	maechler	382	`}`
63578	luke	383
63582	luke	384	`static R_INLINE SEXP ScalarValue1(SEXP x)`
		385	`{`
65020	luke	386	`if (NO_REFERENCES(x))`
63582	luke	387	`return x;`
		388	`else`
		389	`return allocVector(TYPEOF(x), 1);`
		390	`}`
		391
		392	`static R_INLINE SEXP ScalarValue2(SEXP x, SEXP y)`
		393	`{`
65020	luke	394	`if (NO_REFERENCES(x))`
63582	luke	395	`return x;`
65020	luke	396	`else if (NO_REFERENCES(y))`
63582	luke	397	`return y;`
		398	`else`
		399	`return allocVector(TYPEOF(x), 1);`
		400	`}`
		401
1881	ihaka	402	`/* Unary and Binary Operators */`
2	r	403
36990	ripley	404	`SEXP attribute_hidden do_arith(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	405	`{`
63578	luke	406	`SEXP ans, arg1, arg2;`
		407	`int argc;`
2	r	408
63578	luke	409	`if (args == R_NilValue)`
		410	`argc = 0;`
		411	`else if (CDR(args) == R_NilValue)`
		412	`argc = 1;`
		413	`else if (CDDR(args) == R_NilValue)`
		414	`argc = 2;`
		415	`else`
		416	`argc = length(args);`
		417	`arg1 = CAR(args);`
		418	`arg2 = CADR(args);`
2	r	419
63578	luke	420	`if (ATTRIB(arg1) != R_NilValue \|\| ATTRIB(arg2) != R_NilValue) {`
		421	`if (DispatchGroup("Ops", call, op, args, env, &ans))`
		422	`return ans;`
		423	`}`
63582	luke	424	`else if (argc == 2) {`
		425	`/* Handle some scaler operations immediately */`
63616	luke	426	`if (IS_SCALAR(arg1, REALSXP)) {`
73409	luke	427	`double x1 = SCALAR_DVAL(arg1);`
63616	luke	428	`if (IS_SCALAR(arg2, REALSXP)) {`
73409	luke	429	`double x2 = SCALAR_DVAL(arg2);`
63582	luke	430	`ans = ScalarValue2(arg1, arg2);`
		431	`switch (PRIMVAL(op)) {`
73409	luke	432	`case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans;`
		433	`case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans;`
		434	`case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans;`
		435	`case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans;`
63582	luke	436	`}`
		437	`}`
63616	luke	438	`else if (IS_SCALAR(arg2, INTSXP)) {`
73409	luke	439	`int i2 = SCALAR_IVAL(arg2);`
		440	`double x2 = i2 != NA_INTEGER ? (double) i2 : NA_REAL;`
63582	luke	441	`ans = ScalarValue1(arg1);`
		442	`switch (PRIMVAL(op)) {`
73409	luke	443	`case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans;`
		444	`case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans;`
		445	`case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans;`
		446	`case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans;`
63582	luke	447	`}`
		448	`}`
		449	`}`
63616	luke	450	`else if (IS_SCALAR(arg1, INTSXP)) {`
73409	luke	451	`int i1 = SCALAR_IVAL(arg1);`
63616	luke	452	`if (IS_SCALAR(arg2, REALSXP)) {`
73409	luke	453	`double x1 = i1 != NA_INTEGER ? (double) i1 : NA_REAL;`
		454	`double x2 = SCALAR_DVAL(arg2);`
63582	luke	455	`ans = ScalarValue1(arg2);`
		456	`switch (PRIMVAL(op)) {`
73409	luke	457	`case PLUSOP: SET_SCALAR_DVAL(ans, x1 + x2); return ans;`
		458	`case MINUSOP: SET_SCALAR_DVAL(ans, x1 - x2); return ans;`
		459	`case TIMESOP: SET_SCALAR_DVAL(ans, x1 * x2); return ans;`
		460	`case DIVOP: SET_SCALAR_DVAL(ans, x1 / x2); return ans;`
63582	luke	461	`}`
		462	`}`
63616	luke	463	`else if (IS_SCALAR(arg2, INTSXP)) {`
63582	luke	464	`Rboolean naflag = FALSE;`
73409	luke	465	`int i2 = SCALAR_IVAL(arg2);`
63582	luke	466	`switch (PRIMVAL(op)) {`
		467	`case PLUSOP:`
		468	`ans = ScalarValue2(arg1, arg2);`
73409	luke	469	`SET_SCALAR_IVAL(ans, R_integer_plus(i1, i2, &naflag));`
63582	luke	470	`CHECK_INTEGER_OVERFLOW(call, ans, naflag);`
		471	`return ans;`
		472	`case MINUSOP:`
		473	`ans = ScalarValue2(arg1, arg2);`
73409	luke	474	`SET_SCALAR_IVAL(ans, R_integer_minus(i1, i2, &naflag));`
63582	luke	475	`CHECK_INTEGER_OVERFLOW(call, ans, naflag);`
		476	`return ans;`
		477	`case TIMESOP:`
		478	`ans = ScalarValue2(arg1, arg2);`
73409	luke	479	`SET_SCALAR_IVAL(ans, R_integer_times(i1, i2, &naflag));`
63582	luke	480	`CHECK_INTEGER_OVERFLOW(call, ans, naflag);`
		481	`return ans;`
		482	`case DIVOP:`
73409	luke	483	`return ScalarReal(R_integer_divide(i1, i2));`
63582	luke	484	`}`
		485	`}`
		486	`}`
		487	`}`
		488	`else if (argc == 1) {`
63616	luke	489	`if (IS_SCALAR(arg1, REALSXP)) {`
63582	luke	490	`switch(PRIMVAL(op)) {`
		491	`case PLUSOP: return(arg1);`
		492	`case MINUSOP:`
		493	`ans = ScalarValue1(arg1);`
73409	luke	494	`SET_SCALAR_DVAL(ans, -SCALAR_DVAL(arg1));`
63582	luke	495	`return ans;`
		496	`}`
		497	`}`
63616	luke	498	`else if (IS_SCALAR(arg1, INTSXP)) {`
73409	luke	499	`int ival;`
63582	luke	500	`switch(PRIMVAL(op)) {`
		501	`case PLUSOP: return(arg1);`
		502	`case MINUSOP:`
73409	luke	503	`ival = SCALAR_IVAL(arg1);`
63582	luke	504	`ans = ScalarValue1(arg1);`
73409	luke	505	`SET_SCALAR_IVAL(ans, ival == NA_INTEGER ? NA_INTEGER : -ival);`
63582	luke	506	`return ans;`
		507	`}`
		508	`}`
		509	`}`
63578	luke	510
		511	`if (argc == 2)`
		512	`return R_binary(call, op, arg1, arg2);`
		513	`else if (argc == 1)`
		514	`return R_unary(call, op, arg1);`
		515	`else`
41713	ripley	516	`errorcall(call,_("operator needs one or two arguments"));`
6098	pd	517	`return ans; /* never used; to keep -Wall happy */`
2	r	518	`}`
		519
24538	luke	520	`#define COERCE_IF_NEEDED(v, tp, vpi) do { \`
		521	`if (TYPEOF(v) != (tp)) { \`
45446	ripley	522	`int __vo__ = OBJECT(v); \`
24538	luke	523	`REPROTECT(v = coerceVector(v, (tp)), vpi); \`
		524	`if (__vo__) SET_OBJECT(v, 1); \`
		525	`} \`
		526	`} while (0)`
2	r	527
24538	luke	528	`#define FIXUP_NULL_AND_CHECK_TYPES(v, vpi) do { \`
		529	`switch (TYPEOF(v)) { \`
71172	maechler	530	`case NILSXP: REPROTECT(v = allocVector(INTSXP,0), vpi); break; \`
24538	luke	531	`case CPLXSXP: case REALSXP: case INTSXP: case LGLSXP: break; \`
71206	maechler	532	`default: errorcall(call, _("non-numeric argument to binary operator")); \`
24538	luke	533	`} \`
		534	`} while (0)`
		535
38666	ripley	536	`SEXP attribute_hidden R_binary(SEXP call, SEXP op, SEXP x, SEXP y)`
2	r	537	`{`
71206	maechler	538	`Rboolean xattr, yattr, xarray, yarray, xts, yts, xS4, yS4;`
10777	luke	539	`PROTECT_INDEX xpi, ypi;`
39866	duncan	540	`ARITHOP_TYPE oper = (ARITHOP_TYPE) PRIMVAL(op);`
24538	luke	541	`int nprotect = 2; /* x and y */`
2	r	542
		543
10777	luke	544	`PROTECT_WITH_INDEX(x, &xpi);`
		545	`PROTECT_WITH_INDEX(y, &ypi);`
		546
24538	luke	547	`FIXUP_NULL_AND_CHECK_TYPES(x, xpi);`
		548	`FIXUP_NULL_AND_CHECK_TYPES(y, ypi);`
2	r	549
71206	maechler	550	`R_xlen_t`
		551	`nx = XLENGTH(x),`
		552	`ny = XLENGTH(y);`
24538	luke	553	`if (ATTRIB(x) != R_NilValue) {`
		554	`xattr = TRUE;`
		555	`xarray = isArray(x);`
		556	`xts = isTs(x);`
51702	jmc	557	`xS4 = isS4(x);`
1155	maechler	558	`}`
71206	maechler	559	`else xattr = xarray = xts = xS4 = FALSE;`
24538	luke	560	`if (ATTRIB(y) != R_NilValue) {`
		561	`yattr = TRUE;`
		562	`yarray = isArray(y);`
		563	`yts = isTs(y);`
51702	jmc	564	`yS4 = isS4(y);`
24538	luke	565	`}`
71206	maechler	566	`else yattr = yarray = yts = yS4 = FALSE;`
2	r	567
71360	maechler	568	`#define R_ARITHMETIC_ARRAY_1_SPECIAL`
		569
		570	`#ifdef R_ARITHMETIC_ARRAY_1_SPECIAL`
3865	pd	571	`/* If either x or y is a matrix with length 1 and the other is a`
71360	maechler	572	`vector of a different length, we want to coerce the matrix to be a vector.`
28537	ripley	573	`Do we want to? We don't do it! BDR 2004-03-06`
71206	maechler	574
		575	`From 3.4.0 (Sep. 2016), this signals a warning,`
		576	`and in the future we will disable these 2 clauses,`
		577	`so it will give an error.`
28537	ripley	578	`*/`
1010	ihaka	579
1881	ihaka	580	`/* FIXME: Danger Will Robinson.`
1155	maechler	581	`* ----- We might be trashing arguments here.`
		582	`*/`
1881	ihaka	583	`if (xarray != yarray) {`
71206	maechler	584	`if (xarray && nx==1 && ny!=1) {`
71222	maechler	585	`if(ny != 0)`
71360	maechler	586	`warningcall(call, _(`
		587	`"Recycling array of length 1 in array-vector arithmetic is deprecated.\n\`
71361	maechler	588	`Use c() or as.vector() instead.\n"));`
71206	maechler	589	`REPROTECT(x = duplicate(x), xpi);`
		590	`setAttrib(x, R_DimSymbol, R_NilValue);`
		591	`}`
		592	`if (yarray && ny==1 && nx!=1) {`
71222	maechler	593	`if(nx != 0)`
71360	maechler	594	`warningcall(call, _(`
		595	`"Recycling array of length 1 in vector-array arithmetic is deprecated.\n\`
71361	maechler	596	`Use c() or as.vector() instead.\n"));`
71206	maechler	597	`REPROTECT(y = duplicate(y), ypi);`
		598	`setAttrib(y, R_DimSymbol, R_NilValue);`
		599	`}`
785	pd	600	`}`
71360	maechler	601	`#endif`
785	pd	602
71206	maechler	603	`SEXP dims, xnames, ynames;`
571	ihaka	604	`if (xarray \|\| yarray) {`
71222	maechler	605	`/* if one is a length-atleast-1-array and the`
		606	`* other is a length-0 nonarray, then do not use array treatment */`
571	ihaka	607	`if (xarray && yarray) {`
		608	`if (!conformable(x, y))`
71206	maechler	609	`errorcall(call, _("non-conformable arrays"));`
71222	maechler	610	`PROTECT(dims = getAttrib(x, R_DimSymbol)); nprotect++;`
2	r	611	`}`
71222	maechler	612	`else if (xarray && (ny != 0 \|\| nx == 0)) {`
		613	`PROTECT(dims = getAttrib(x, R_DimSymbol)); nprotect++;`
2	r	614	`}`
71222	maechler	615	`else if (yarray && (nx != 0 \|\| ny == 0)) {`
		616	`PROTECT(dims = getAttrib(y, R_DimSymbol)); nprotect++;`
		617	`} else`
		618	`dims = R_NilValue;`
24538	luke	619	`if (xattr) {`
		620	`PROTECT(xnames = getAttrib(x, R_DimNamesSymbol));`
		621	`nprotect++;`
		622	`}`
		623	`else xnames = R_NilValue;`
		624	`if (yattr) {`
		625	`PROTECT(ynames = getAttrib(y, R_DimNamesSymbol));`
		626	`nprotect++;`
		627	`}`
		628	`else ynames = R_NilValue;`
571	ihaka	629	`}`
		630	`else {`
24538	luke	631	`dims = R_NilValue;`
		632	`if (xattr) {`
		633	`PROTECT(xnames = getAttrib(x, R_NamesSymbol));`
		634	`nprotect++;`
		635	`}`
		636	`else xnames = R_NilValue;`
		637	`if (yattr) {`
		638	`PROTECT(ynames = getAttrib(y, R_NamesSymbol));`
		639	`nprotect++;`
		640	`}`
		641	`else ynames = R_NilValue;`
571	ihaka	642	`}`
2	r	643
71206	maechler	644	`SEXP klass = NULL, tsp = NULL; // -Wall`
571	ihaka	645	`if (xts \|\| yts) {`
		646	`if (xts && yts) {`
74304	kalibera	647	`/* could check ts conformance here */`
571	ihaka	648	`PROTECT(tsp = getAttrib(x, R_TspSymbol));`
39866	duncan	649	`PROTECT(klass = getAttrib(x, R_ClassSymbol));`
2	r	650	`}`
571	ihaka	651	`else if (xts) {`
24538	luke	652	`if (nx < ny)`
71206	maechler	653	`ErrorMessage(call, ERROR_TSVEC_MISMATCH);`
571	ihaka	654	`PROTECT(tsp = getAttrib(x, R_TspSymbol));`
39866	duncan	655	`PROTECT(klass = getAttrib(x, R_ClassSymbol));`
2	r	656	`}`
6203	maechler	657	`else { /* (yts) */`
24538	luke	658	`if (ny < nx)`
71206	maechler	659	`ErrorMessage(call, ERROR_TSVEC_MISMATCH);`
571	ihaka	660	`PROTECT(tsp = getAttrib(y, R_TspSymbol));`
39866	duncan	661	`PROTECT(klass = getAttrib(y, R_ClassSymbol));`
2	r	662	`}`
24538	luke	663	`nprotect += 2;`
571	ihaka	664	`}`
4719	maechler	665
71206	maechler	666	`if (nx > 0 && ny > 0 &&`
		667	`((nx > ny) ? nx % ny : ny % nx) != 0) // mismatch`
		668	`warningcall(call,`
41680	ripley	669	`_("longer object length is not a multiple of shorter object length"));`
2	r	670
71206	maechler	671	`SEXP val;`
24538	luke	672	`/* need to preserve object here, as _binary copies class attributes /`
571	ihaka	673	`if (TYPEOF(x) == CPLXSXP \|\| TYPEOF(y) == CPLXSXP) {`
24538	luke	674	`COERCE_IF_NEEDED(x, CPLXSXP, xpi);`
		675	`COERCE_IF_NEEDED(y, CPLXSXP, ypi);`
39866	duncan	676	`val = complex_binary(oper, x, y);`
571	ihaka	677	`}`
24538	luke	678	`else if (TYPEOF(x) == REALSXP \|\| TYPEOF(y) == REALSXP) {`
63461	luke	679	`/* real_binary can handle REALSXP or INTSXP operand, but not LGLSXP. */`
		680	`/* Can get a LGLSXP. In base-Ex.R on 24 Oct '06, got 8 of these. */`
		681	`if (TYPEOF(x) != INTSXP) COERCE_IF_NEEDED(x, REALSXP, xpi);`
		682	`if (TYPEOF(y) != INTSXP) COERCE_IF_NEEDED(y, REALSXP, ypi);`
39866	duncan	683	`val = real_binary(oper, x, y);`
24538	luke	684	`}`
71206	maechler	685	`else val = integer_binary(oper, x, y, call);`
2	r	686
24539	luke	687	`/* quick return if there are no attributes */`
24538	luke	688	`if (! xattr && ! yattr) {`
		689	`UNPROTECT(nprotect);`
		690	`return val;`
		691	`}`
		692
		693	`PROTECT(val);`
		694	`nprotect++;`
		695
571	ihaka	696	`if (dims != R_NilValue) {`
24538	luke	697	`setAttrib(val, R_DimSymbol, dims);`
1881	ihaka	698	`if (xnames != R_NilValue)`
24538	luke	699	`setAttrib(val, R_DimNamesSymbol, xnames);`
1881	ihaka	700	`else if (ynames != R_NilValue)`
24538	luke	701	`setAttrib(val, R_DimNamesSymbol, ynames);`
571	ihaka	702	`}`
		703	`else {`
63461	luke	704	`if (XLENGTH(val) == xlength(xnames))`
24538	luke	705	`setAttrib(val, R_NamesSymbol, xnames);`
63461	luke	706	`else if (XLENGTH(val) == xlength(ynames))`
24538	luke	707	`setAttrib(val, R_NamesSymbol, ynames);`
571	ihaka	708	`}`
		709
6098	pd	710	`if (xts \|\| yts) { /* must set after dims! */`
24538	luke	711	`setAttrib(val, R_TspSymbol, tsp);`
39866	duncan	712	`setAttrib(val, R_ClassSymbol, klass);`
4562	pd	713	`}`
		714
51702	jmc	715	`if(xS4 \|\| yS4) { /* Only set the bit: no method defined! */`
68923	ripley	716	`val = asS4(val, TRUE, TRUE);`
51702	jmc	717	`}`
24538	luke	718	`UNPROTECT(nprotect);`
		719	`return val;`
2	r	720	`}`
		721
38666	ripley	722	`SEXP attribute_hidden R_unary(SEXP call, SEXP op, SEXP s1)`
2	r	723	`{`
39866	duncan	724	`ARITHOP_TYPE operation = (ARITHOP_TYPE) PRIMVAL(op);`
571	ihaka	725	`switch (TYPEOF(s1)) {`
		726	`case LGLSXP:`
63327	ripley	727	`return logical_unary(operation, s1, call);`
571	ihaka	728	`case INTSXP:`
41713	ripley	729	`return integer_unary(operation, s1, call);`
571	ihaka	730	`case REALSXP:`
39866	duncan	731	`return real_unary(operation, s1, call);`
571	ihaka	732	`case CPLXSXP:`
41713	ripley	733	`return complex_unary(operation, s1, call);`
571	ihaka	734	`default:`
33297	ripley	735	`errorcall(call, _("invalid argument to unary operator"));`
571	ihaka	736	`}`
6098	pd	737	`return s1; /* never used; to keep -Wall happy */`
2	r	738	`}`
		739
63327	ripley	740	`static SEXP logical_unary(ARITHOP_TYPE code, SEXP s1, SEXP call)`
		741	`{`
63333	ripley	742	`R_xlen_t n = XLENGTH(s1);`
		743	`SEXP ans = PROTECT(allocVector(INTSXP, n));`
		744	`SEXP names = PROTECT(getAttrib(s1, R_NamesSymbol));`
		745	`SEXP dim = PROTECT(getAttrib(s1, R_DimSymbol));`
		746	`SEXP dimnames = PROTECT(getAttrib(s1, R_DimNamesSymbol));`
		747	`if(names != R_NilValue) setAttrib(ans, R_NamesSymbol, names);`
		748	`if(dim != R_NilValue) setAttrib(ans, R_DimSymbol, dim);`
		749	`if(dimnames != R_NilValue) setAttrib(ans, R_DimNamesSymbol, dimnames);`
		750	`UNPROTECT(3);`
63327	ripley	751
73409	luke	752	`int *pa = INTEGER(ans);`
73732	luke	753	`const int *px = LOGICAL_RO(s1);`
73409	luke	754
63327	ripley	755	`switch (code) {`
		756	`case PLUSOP:`
73409	luke	757	`for (R_xlen_t i = 0; i < n; i++) pa[i] = px[i];`
63333	ripley	758	`break;`
63327	ripley	759	`case MINUSOP:`
63333	ripley	760	`for (R_xlen_t i = 0; i < n; i++) {`
73409	luke	761	`int x = px[i];`
		762	`pa[i] = (x == NA_INTEGER) ?`
63327	ripley	763	`NA_INTEGER : ((x == 0.0) ? 0 : -x);`
		764	`}`
63333	ripley	765	`break;`
63327	ripley	766	`default:`
		767	`errorcall(call, _("invalid unary operator"));`
		768	`}`
63333	ripley	769	`UNPROTECT(1);`
		770	`return ans;`
63327	ripley	771	`}`
		772
41713	ripley	773	`static SEXP integer_unary(ARITHOP_TYPE code, SEXP s1, SEXP call)`
2	r	774	`{`
59086	ripley	775	`R_xlen_t i, n;`
571	ihaka	776	`SEXP ans;`
2	r	777
1881	ihaka	778	`switch (code) {`
571	ihaka	779	`case PLUSOP:`
		780	`return s1;`
		781	`case MINUSOP:`
65020	luke	782	`ans = NO_REFERENCES(s1) ? s1 : duplicate(s1);`
73732	luke	783	`int *pa = INTEGER(ans);`
		784	`const int *px = INTEGER_RO(s1);`
59048	ripley	785	`n = XLENGTH(s1);`
571	ihaka	786	`for (i = 0; i < n; i++) {`
73732	luke	787	`int x = px[i];`
73409	luke	788	`pa[i] = (x == NA_INTEGER) ?`
571	ihaka	789	`NA_INTEGER : ((x == 0.0) ? 0 : -x);`
2	r	790	`}`
571	ihaka	791	`return ans;`
		792	`default:`
41713	ripley	793	`errorcall(call, _("invalid unary operator"));`
571	ihaka	794	`}`
6098	pd	795	`return s1; /* never used; to keep -Wall happy */`
2	r	796	`}`
		797
13997	duncan	798	`static SEXP real_unary(ARITHOP_TYPE code, SEXP s1, SEXP lcall)`
2	r	799	`{`
59048	ripley	800	`R_xlen_t i, n;`
571	ihaka	801	`SEXP ans;`
2	r	802
1881	ihaka	803	`switch (code) {`
571	ihaka	804	`case PLUSOP: return s1;`
		805	`case MINUSOP:`
65020	luke	806	`ans = NO_REFERENCES(s1) ? s1 : duplicate(s1);`
73732	luke	807	`double *pa = REAL(ans);`
		808	`const double *px = REAL_RO(s1);`
59048	ripley	809	`n = XLENGTH(s1);`
29410	ripley	810	`for (i = 0; i < n; i++)`
73409	luke	811	`pa[i] = -px[i];`
571	ihaka	812	`return ans;`
		813	`default:`
33101	ripley	814	`errorcall(lcall, _("invalid unary operator"));`
571	ihaka	815	`}`
6098	pd	816	`return s1; /* never used; to keep -Wall happy */`
2	r	817	`}`
988	maechler	818
16033	luke	819	`static SEXP integer_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2, SEXP lcall)`
2	r	820	`{`
59048	ripley	821	`R_xlen_t i, i1, i2, n, n1, n2;`
571	ihaka	822	`int x1, x2;`
		823	`SEXP ans;`
16033	luke	824	`Rboolean naflag = FALSE;`
2	r	825
59048	ripley	826	`n1 = XLENGTH(s1);`
		827	`n2 = XLENGTH(s2);`
16613	ripley	828	`/* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */`
		829	`if (n1 == 0 \|\| n2 == 0) n = 0; else n = (n1 > n2) ? n1 : n2;`
2	r	830
571	ihaka	831	`if (code == DIVOP \|\| code == POWOP)`
		832	`ans = allocVector(REALSXP, n);`
		833	`else`
63462	luke	834	`ans = R_allocOrReuseVector(s1, s2, INTSXP, n);`
63461	luke	835	`if (n == 0) return(ans);`
63290	luke	836	`PROTECT(ans);`
2	r	837
571	ihaka	838	`switch (code) {`
		839	`case PLUSOP:`
73409	luke	840	`{`
		841	`int *pa = INTEGER(ans);`
73732	luke	842	`const int *px1 = INTEGER_RO(s1);`
		843	`const int *px2 = INTEGER_RO(s2);`
73409	luke	844	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		845	`x1 = px1[i1];`
		846	`x2 = px2[i2];`
		847	`pa[i] = R_integer_plus(x1, x2, &naflag);`
		848	`});`
		849	`if (naflag)`
		850	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
		851	`}`
571	ihaka	852	`break;`
		853	`case MINUSOP:`
73409	luke	854	`{`
		855	`int *pa = INTEGER(ans);`
73732	luke	856	`const int *px1 = INTEGER_RO(s1);`
		857	`const int *px2 = INTEGER_RO(s2);`
73409	luke	858	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		859	`x1 = px1[i1];`
		860	`x2 = px2[i2];`
		861	`pa[i] = R_integer_minus(x1, x2, &naflag);`
		862	`});`
		863	`if (naflag)`
		864	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
		865	`}`
571	ihaka	866	`break;`
		867	`case TIMESOP:`
73409	luke	868	`{`
		869	`int *pa = INTEGER(ans);`
73732	luke	870	`const int *px1 = INTEGER_RO(s1);`
		871	`const int *px2 = INTEGER_RO(s2);`
73409	luke	872	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		873	`x1 = px1[i1];`
		874	`x2 = px2[i2];`
		875	`pa[i] = R_integer_times(x1, x2, &naflag);`
		876	`});`
		877	`if (naflag)`
		878	`warningcall(lcall, INTEGER_OVERFLOW_WARNING);`
		879	`}`
571	ihaka	880	`break;`
		881	`case DIVOP:`
73409	luke	882	`{`
		883	`double *pa = REAL(ans);`
73732	luke	884	`const int *px1 = INTEGER_RO(s1);`
		885	`const int *px2 = INTEGER_RO(s2);`
73409	luke	886	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		887	`x1 = px1[i1];`
		888	`x2 = px2[i2];`
		889	`pa[i] = R_integer_divide(x1, x2);`
		890	`});`
		891	`}`
571	ihaka	892	`break;`
		893	`case POWOP:`
73409	luke	894	`{`
		895	`double *pa = REAL(ans);`
73732	luke	896	`const int *px1 = INTEGER_RO(s1);`
		897	`const int *px2 = INTEGER_RO(s2);`
73409	luke	898	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		899	`if((x1 = px1[i1]) == 1 \|\| (x2 = px2[i2]) == 0)`
		900	`pa[i] = 1.;`
		901	`else if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER)`
		902	`pa[i] = NA_REAL;`
		903	`else`
		904	`pa[i] = R_POW((double) x1, (double) x2);`
		905	`});`
		906	`}`
571	ihaka	907	`break;`
		908	`case MODOP:`
73409	luke	909	`{`
		910	`int *pa = INTEGER(ans);`
73732	luke	911	`const int *px1 = INTEGER_RO(s1);`
		912	`const int *px2 = INTEGER_RO(s2);`
73409	luke	913	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		914	`x1 = px1[i1];`
		915	`x2 = px2[i2];`
		916	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER \|\| x2 == 0)`
		917	`pa[i] = NA_INTEGER;`
		918	`else {`
		919	`pa[i] = /* till 0.63.2: x1 % x2 */`
		920	`(x1 >= 0 && x2 > 0) ? x1 % x2 :`
		921	`(int)myfmod((double)x1,(double)x2);`
		922	`}`
		923	`});`
		924	`}`
571	ihaka	925	`break;`
		926	`case IDIVOP:`
73409	luke	927	`{`
		928	`int *pa = INTEGER(ans);`
73732	luke	929	`const int *px1 = INTEGER_RO(s1);`
		930	`const int *px2 = INTEGER_RO(s2);`
73409	luke	931	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2, {`
		932	`x1 = px1[i1];`
		933	`x2 = px2[i2];`
		934	`/* This had x %/% 0 == 0 prior to 2.14.1, but`
		935	`it seems conventionally to be undefined */`
		936	`if (x1 == NA_INTEGER \|\| x2 == NA_INTEGER \|\| x2 == 0)`
		937	`pa[i] = NA_INTEGER;`
		938	`else`
		939	`pa[i] = (int) floor((double)x1 / (double)x2);`
		940	`});`
		941	`}`
571	ihaka	942	`break;`
		943	`}`
63290	luke	944	`UNPROTECT(1);`
1010	ihaka	945
24539	luke	946	`/* quick return if there are no attributes */`
		947	`if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)`
		948	`return ans;`
		949
16613	ripley	950	`/* Copy attributes from longer argument. */`
1010	ihaka	951
63462	luke	952	`if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue)`
1010	ihaka	953	`copyMostAttrib(s2, ans);`
63462	luke	954	`if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue)`
63461	luke	955	`copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */`
1010	ihaka	956
571	ihaka	957	`return ans;`
2	r	958	`}`
		959
73409	luke	960	`#define R_INTEGER(x) (double) ((x) == NA_INTEGER ? NA_REAL : (x))`
39742	duncan	961
10718	maechler	962	`static SEXP real_binary(ARITHOP_TYPE code, SEXP s1, SEXP s2)`
2	r	963	`{`
59048	ripley	964	`R_xlen_t i, i1, i2, n, n1, n2;`
571	ihaka	965	`SEXP ans;`
2	r	966
1010	ihaka	967	`/* Note: "s1" and "s2" are protected above. */`
59048	ripley	968	`n1 = XLENGTH(s1);`
		969	`n2 = XLENGTH(s2);`
16613	ripley	970
		971	`/* S4-compatibility change: if n1 or n2 is 0, result is of length 0 */`
		972	`if (n1 == 0 \|\| n2 == 0) return(allocVector(REALSXP, 0));`
		973
18510	maechler	974	`n = (n1 > n2) ? n1 : n2;`
63462	luke	975	`PROTECT(ans = R_allocOrReuseVector(s1, s2, REALSXP, n));`
2	r	976
571	ihaka	977	`switch (code) {`
		978	`case PLUSOP:`
39742	duncan	979	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
63918	luke	980	`double *da = REAL(ans);`
73732	luke	981	`const double *dx = REAL_RO(s1);`
		982	`const double *dy = REAL_RO(s2);`
68923	ripley	983	`if (n2 == 1) {`
63918	luke	984	`double tmp = dy[0];`
		985	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] + tmp;);`
		986	`}`
68923	ripley	987	`else if (n1 == 1) {`
63918	luke	988	`double tmp = dx[0];`
		989	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp + dy[i];);`
		990	`}`
68923	ripley	991	`else if (n1 == n2)`
63918	luke	992	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] + dy[i];);`
68923	ripley	993	`else`
68663	luke	994	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
63918	luke	995	`da[i] = dx[i1] + dy[i2];);`
2	r	996	`}`
73409	luke	997	`else if(TYPEOF(s1) == INTSXP ) {`
		998	`double *da = REAL(ans);`
73732	luke	999	`const int *px1 = INTEGER_RO(s1);`
		1000	`const double *px2 = REAL_RO(s2);`
68663	luke	1001	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1002	`da[i] = R_INTEGER(px1[i1]) + px2[i2];);`
		1003	`}`
		1004	`else if(TYPEOF(s2) == INTSXP ) {`
		1005	`double *da = REAL(ans);`
73732	luke	1006	`const double *px1 = REAL_RO(s1);`
		1007	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1008	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1009	`da[i] = px1[i1] + R_INTEGER(px2[i2]););`
		1010	`}`
571	ihaka	1011	`break;`
		1012	`case MINUSOP:`
39742	duncan	1013	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
63918	luke	1014	`double *da = REAL(ans);`
73732	luke	1015	`const double *dx = REAL_RO(s1);`
		1016	`const double *dy = REAL_RO(s2);`
68923	ripley	1017	`if (n2 == 1) {`
		1018	`double tmp = dy[0];`
63918	luke	1019	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] - tmp;);`
68923	ripley	1020	`}`
		1021	`else if (n1 == 1) {`
		1022	`double tmp = dx[0];`
63918	luke	1023	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp - dy[i];);`
68923	ripley	1024	`}`
		1025	`else if (n1 == n2)`
63918	luke	1026	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] - dy[i];);`
68923	ripley	1027	`else`
68663	luke	1028	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
63918	luke	1029	`da[i] = dx[i1] - dy[i2];);`
571	ihaka	1030	`}`
73409	luke	1031	`else if(TYPEOF(s1) == INTSXP ) {`
		1032	`double *da = REAL(ans);`
73732	luke	1033	`const int *px1 = INTEGER_RO(s1);`
		1034	`const double *px2 = REAL_RO(s2);`
68663	luke	1035	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1036	`da[i] = R_INTEGER(px1[i1]) - px2[i2];);`
		1037	`}`
		1038	`else if(TYPEOF(s2) == INTSXP ) {`
		1039	`double *da = REAL(ans);`
73732	luke	1040	`const double *px1 = REAL_RO(s1);`
		1041	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1042	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1043	`da[i] = px1[i1] - R_INTEGER(px2[i2]););`
		1044	`}`
571	ihaka	1045	`break;`
		1046	`case TIMESOP:`
39742	duncan	1047	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
63918	luke	1048	`double *da = REAL(ans);`
73732	luke	1049	`const double *dx = REAL_RO(s1);`
		1050	`const double *dy = REAL_RO(s2);`
68923	ripley	1051	`if (n2 == 1) {`
		1052	`double tmp = dy[0];`
63918	luke	1053	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] * tmp;);`
68923	ripley	1054	`}`
		1055	`else if (n1 == 1) {`
73409	luke	1056	`double tmp = dx[0];`
63918	luke	1057	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp * dy[i];);`
68923	ripley	1058	`}`
		1059	`else if (n1 == n2)`
63918	luke	1060	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] * dy[i];);`
68923	ripley	1061	`else`
68663	luke	1062	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
63918	luke	1063	`da[i] = dx[i1] * dy[i2];);`
571	ihaka	1064	`}`
73409	luke	1065	`else if(TYPEOF(s1) == INTSXP ) {`
		1066	`double *da = REAL(ans);`
73732	luke	1067	`const int *px1 = INTEGER_RO(s1);`
		1068	`const double *px2 = REAL_RO(s2);`
68663	luke	1069	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1070	`da[i] = R_INTEGER(px1[i1]) * px2[i2];);`
		1071	`}`
		1072	`else if(TYPEOF(s2) == INTSXP ) {`
		1073	`double *da = REAL(ans);`
73732	luke	1074	`const double *px1 = REAL_RO(s1);`
		1075	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1076	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1077	`da[i] = px1[i1] * R_INTEGER(px2[i2]););`
		1078	`}`
571	ihaka	1079	`break;`
		1080	`case DIVOP:`
39742	duncan	1081	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
63918	luke	1082	`double *da = REAL(ans);`
73732	luke	1083	`const double *dx = REAL_RO(s1);`
		1084	`const double *dy = REAL_RO(s2);`
68923	ripley	1085	`if (n2 == 1) {`
		1086	`double tmp = dy[0];`
63918	luke	1087	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] / tmp;);`
68923	ripley	1088	`}`
		1089	`else if (n1 == 1) {`
		1090	`double tmp = dx[0];`
63918	luke	1091	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = tmp / dy[i];);`
68923	ripley	1092	`}`
		1093	`else if (n1 == n2)`
63918	luke	1094	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = dx[i] / dy[i];);`
68923	ripley	1095	`else`
68663	luke	1096	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
63918	luke	1097	`da[i] = dx[i1] / dy[i2];);`
571	ihaka	1098	`}`
73409	luke	1099	`else if(TYPEOF(s1) == INTSXP ) {`
		1100	`double *da = REAL(ans);`
73732	luke	1101	`const int *px1 = INTEGER_RO(s1);`
		1102	`const double *px2 = REAL_RO(s2);`
68663	luke	1103	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1104	`da[i] = R_INTEGER(px1[i1]) / px2[i2];);`
		1105	`}`
		1106	`else if(TYPEOF(s2) == INTSXP ) {`
		1107	`double *da = REAL(ans);`
73732	luke	1108	`const double *px1 = REAL_RO(s1);`
		1109	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1110	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1111	`da[i] = px1[i1] / R_INTEGER(px2[i2]););`
		1112	`}`
571	ihaka	1113	`break;`
		1114	`case POWOP:`
39742	duncan	1115	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
63918	luke	1116	`double *da = REAL(ans);`
73732	luke	1117	`const double *dx = REAL_RO(s1);`
		1118	`const double *dy = REAL_RO(s2);`
68923	ripley	1119	`if (n2 == 1) {`
		1120	`double tmp = dy[0];`
63918	luke	1121	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(dx[i], tmp););`
68923	ripley	1122	`}`
		1123	`else if (n1 == 1) {`
		1124	`double tmp = dx[0];`
63918	luke	1125	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(tmp, dy[i]););`
68923	ripley	1126	`}`
		1127	`else if (n1 == n2)`
63918	luke	1128	`R_ITERATE_CHECK(NINTERRUPT, n, i, da[i] = R_POW(dx[i], dy[i]););`
68923	ripley	1129	`else`
68663	luke	1130	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
63918	luke	1131	`da[i] = R_POW(dx[i1], dy[i2]););`
571	ihaka	1132	`}`
73409	luke	1133	`else if(TYPEOF(s1) == INTSXP ) {`
		1134	`double *da = REAL(ans);`
73732	luke	1135	`const int *px1 = INTEGER_RO(s1);`
		1136	`const double *px2 = REAL_RO(s2);`
68663	luke	1137	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1138	`da[i] = R_POW( R_INTEGER(px1[i1]), px2[i2]););`
		1139	`}`
		1140	`else if(TYPEOF(s2) == INTSXP ) {`
		1141	`double *da = REAL(ans);`
73732	luke	1142	`const double *px1 = REAL_RO(s1);`
		1143	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1144	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1145	`da[i] = R_POW(px1[i1], R_INTEGER(px2[i2])););`
		1146	`}`
571	ihaka	1147	`break;`
		1148	`case MODOP:`
73409	luke	1149	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
		1150	`double *da = REAL(ans);`
73732	luke	1151	`const double *px1 = REAL_RO(s1);`
		1152	`const double *px2 = REAL_RO(s2);`
68663	luke	1153	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1154	`da[i] = myfmod(px1[i1], px2[i2]););`
		1155	`}`
		1156	`else if(TYPEOF(s1) == INTSXP ) {`
		1157	`double *da = REAL(ans);`
73732	luke	1158	`const int *px1 = INTEGER_RO(s1);`
		1159	`const double *px2 = REAL_RO(s2);`
68663	luke	1160	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1161	`da[i] = myfmod(R_INTEGER(px1[i1]), px2[i2]););`
		1162	`}`
		1163	`else if(TYPEOF(s2) == INTSXP ) {`
		1164	`double *da = REAL(ans);`
73732	luke	1165	`const double *px1 = REAL_RO(s1);`
		1166	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1167	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1168	`da[i] = myfmod(px1[i1], R_INTEGER(px2[i2])););`
		1169	`}`
571	ihaka	1170	`break;`
		1171	`case IDIVOP:`
73409	luke	1172	`if(TYPEOF(s1) == REALSXP && TYPEOF(s2) == REALSXP) {`
		1173	`double *da = REAL(ans);`
73732	luke	1174	`const double *px1 = REAL_RO(s1);`
		1175	`const double *px2 = REAL_RO(s2);`
68663	luke	1176	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1177	`da[i] = myfloor(px1[i1], px2[i2]););`
		1178	`}`
		1179	`else if(TYPEOF(s1) == INTSXP ) {`
		1180	`double *da = REAL(ans);`
73732	luke	1181	`const int *px1 = INTEGER_RO(s1);`
		1182	`const double *px2 = REAL_RO(s2);`
68663	luke	1183	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1184	`da[i] = myfloor(R_INTEGER(px1[i1]), px2[i2]););`
		1185	`}`
		1186	`else if(TYPEOF(s2) == INTSXP ) {`
		1187	`double *da = REAL(ans);`
73732	luke	1188	`const double *px1 = REAL_RO(s1);`
		1189	`const int *px2 = INTEGER_RO(s2);`
68663	luke	1190	`MOD_ITERATE2_CHECK(NINTERRUPT, n, n1, n2, i, i1, i2,`
73409	luke	1191	`da[i] = myfloor(px1[i1], R_INTEGER(px2[i2])););`
		1192	`}`
571	ihaka	1193	`break;`
		1194	`}`
63461	luke	1195	`UNPROTECT(1);`
1010	ihaka	1196
24539	luke	1197	`/* quick return if there are no attributes */`
63461	luke	1198	`if (ATTRIB(s1) == R_NilValue && ATTRIB(s2) == R_NilValue)`
24539	luke	1199	`return ans;`
37954	maechler	1200
16613	ripley	1201	`/* Copy attributes from longer argument. */`
1010	ihaka	1202
63462	luke	1203	`if (ans != s2 && n == n2 && ATTRIB(s2) != R_NilValue)`
1010	ihaka	1204	`copyMostAttrib(s2, ans);`
63462	luke	1205	`if (ans != s1 && n == n1 && ATTRIB(s1) != R_NilValue)`
63461	luke	1206	`copyMostAttrib(s1, ans); /* Done 2nd so s1's attrs overwrite s2's */`
1010	ihaka	1207
571	ihaka	1208	`return ans;`
2	r	1209	`}`
		1210
7615	maechler	1211
1881	ihaka	1212	`/* Mathematical Functions of One Argument */`
2	r	1213
39866	duncan	1214	`static SEXP math1(SEXP sa, double(*f)(double), SEXP lcall)`
2	r	1215	`{`
571	ihaka	1216	`SEXP sy;`
59048	ripley	1217	`R_xlen_t i, n;`
13997	duncan	1218	`int naflag;`
2	r	1219
7615	maechler	1220	`if (!isNumeric(sa))`
9147	maechler	1221	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	1222
63461	luke	1223	`n = XLENGTH(sa);`
21984	ripley	1224	`/* coercion can lose the object bit */`
7615	maechler	1225	`PROTECT(sa = coerceVector(sa, REALSXP));`
65020	luke	1226	`PROTECT(sy = NO_REFERENCES(sa) ? sa : allocVector(REALSXP, n));`
73732	luke	1227	`const double *a = REAL_RO(sa);`
		1228	`double *y = REAL(sy);`
7615	maechler	1229	`naflag = 0;`
		1230	`for (i = 0; i < n; i++) {`
69955	luke	1231	`double x = a[i]; /* in case y == a */`
		1232	`/* This code assumes that ISNAN(x) implies ISNAN(f(x)), so we`
		1233	`only need to check ISNAN(x) if ISNAN(f(x)) is true. */`
		1234	`y[i] = f(x);`
		1235	`if (ISNAN(y[i])) {`
		1236	`if (ISNAN(x))`
		1237	`y[i] = x; /* make sure the incoming NaN is preserved */`
		1238	`else`
		1239	`naflag = 1;`
2	r	1240	`}`
571	ihaka	1241	`}`
59447	ripley	1242	`/* These are primitives, so need to use the call */`
		1243	`if(naflag) warningcall(lcall, R_MSG_NA);`
7615	maechler	1244
63462	luke	1245	`if (sa != sy && ATTRIB(sa) != R_NilValue)`
69108	luke	1246	`SHALLOW_DUPLICATE_ATTRIB(sy, sa);`
7615	maechler	1247	`UNPROTECT(2);`
		1248	`return sy;`
2	r	1249	`}`
		1250
66651	ripley	1251
36990	ripley	1252	`SEXP attribute_hidden do_math1(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1253	`{`
571	ihaka	1254	`SEXP s;`
2	r	1255
41221	ripley	1256	`checkArity(op, args);`
51267	ripley	1257	`check1arg(args, call, "x");`
41221	ripley	1258
44894	ripley	1259	`if (DispatchGroup("Math", call, op, args, env, &s))`
571	ihaka	1260	`return s;`
2	r	1261
571	ihaka	1262	`if (isComplex(CAR(args)))`
		1263	`return complex_math1(call, op, args, env);`
2	r	1264
13997	duncan	1265	`#define MATH1(x) math1(CAR(args), x, call);`
571	ihaka	1266	`switch (PRIMVAL(op)) {`
13997	duncan	1267	`case 1: return MATH1(floor);`
		1268	`case 2: return MATH1(ceil);`
		1269	`case 3: return MATH1(sqrt);`
		1270	`case 4: return MATH1(sign);`
41221	ripley	1271	`/* case 5: return MATH1(trunc); separate from 2.6.0 */`
2	r	1272
13997	duncan	1273	`case 10: return MATH1(exp);`
18510	maechler	1274	`case 11: return MATH1(expm1);`
13997	duncan	1275	`case 12: return MATH1(log1p);`
77193	maechler	1276
13997	duncan	1277	`case 20: return MATH1(cos);`
		1278	`case 21: return MATH1(sin);`
		1279	`case 22: return MATH1(tan);`
		1280	`case 23: return MATH1(acos);`
		1281	`case 24: return MATH1(asin);`
51270	ripley	1282	`case 25: return MATH1(atan);`
14701	ripley	1283
13997	duncan	1284	`case 30: return MATH1(cosh);`
		1285	`case 31: return MATH1(sinh);`
		1286	`case 32: return MATH1(tanh);`
		1287	`case 33: return MATH1(acosh);`
		1288	`case 34: return MATH1(asinh);`
		1289	`case 35: return MATH1(atanh);`
14701	ripley	1290
13997	duncan	1291	`case 40: return MATH1(lgammafn);`
		1292	`case 41: return MATH1(gammafn);`
14701	ripley	1293
13997	duncan	1294	`case 42: return MATH1(digamma);`
		1295	`case 43: return MATH1(trigamma);`
33174	ripley	1296	`/* case 44: return MATH1(tetragamma);`
37954	maechler	1297	`case 45: return MATH1(pentagamma);`
77193	maechler	1298	`removed in 2.0.0 -- rather use Math2's psigamma()`
64171	maechler	1299
		1300	`case 46: return MATH1(Rf_gamma_cody); removed in 2.8.0`
33174	ripley	1301	`*/`
64202	maechler	1302	`case 47: return MATH1(cospi);`
		1303	`case 48: return MATH1(sinpi);`
70093	ripley	1304	`#if defined(HAVE_TANPI) \|\| defined(HAVE___TANPI)`
		1305	`case 49: return MATH1(Rtanpi);`
		1306	`#else`
64171	maechler	1307	`case 49: return MATH1(tanpi);`
64645	ripley	1308	`#endif`
14701	ripley	1309
571	ihaka	1310	`default:`
33297	ripley	1311	`errorcall(call, _("unimplemented real function of 1 argument"));`
571	ihaka	1312	`}`
40141	maechler	1313	`return s; /* never used; to keep -Wall happy */`
2	r	1314	`}`
		1315
51270	ripley	1316	`/* methods are allowed to have more than one arg */`
41221	ripley	1317	`SEXP attribute_hidden do_trunc(SEXP call, SEXP op, SEXP args, SEXP env)`
		1318	`{`
		1319	`SEXP s;`
		1320	`if (DispatchGroup("Math", call, op, args, env, &s))`
		1321	`return s;`
74907	maechler	1322	`// checkArity(op, args); /* is -1 in names.c */`
51316	ripley	1323	`check1arg(args, call, "x");`
41221	ripley	1324	`if (isComplex(CAR(args)))`
		1325	`errorcall(call, _("unimplemented complex function"));`
		1326	`return math1(CAR(args), trunc, call);`
		1327	`}`
		1328
51267	ripley	1329	`/*`
58014	maechler	1330	`Note that this is slightly different from the do_math1 set,`
51267	ripley	1331	`both for integer/logical inputs and what it dispatches to for complex ones.`
		1332	`*/`
		1333
36990	ripley	1334	`SEXP attribute_hidden do_abs(SEXP call, SEXP op, SEXP args, SEXP env)`
18454	ripley	1335	`{`
51267	ripley	1336	`SEXP x, s = R_NilValue /* -Wall */;`
		1337
		1338	`checkArity(op, args);`
		1339	`check1arg(args, call, "x");`
		1340	`x = CAR(args);`
		1341
18454	ripley	1342	`if (DispatchGroup("Math", call, op, args, env, &s))`
		1343	`return s;`
51267	ripley	1344
		1345	`if (isInteger(x) \|\| isLogical(x)) {`
		1346	`/* integer or logical ==> return integer,`
		1347	`factor was covered by Math.factor. */`
63461	luke	1348	`R_xlen_t i, n = XLENGTH(x);`
65020	luke	1349	`s = (NO_REFERENCES(x) && TYPEOF(x) == INTSXP) ?`
		1350	`x : allocVector(INTSXP, n);`
63462	luke	1351	`PROTECT(s);`
40141	maechler	1352	`/* Note: relying on INTEGER(.) === LOGICAL(.) : */`
73409	luke	1353	`int *pa = INTEGER(s);`
73732	luke	1354	`const int *px = INTEGER_RO(x);`
65563	ripley	1355	`for(i = 0 ; i < n ; i++) {`
73409	luke	1356	`int xi = px[i];`
		1357	`pa[i] = (xi == NA_INTEGER) ? xi : abs(xi);`
68923	ripley	1358	`}`
51267	ripley	1359	`} else if (TYPEOF(x) == REALSXP) {`
63461	luke	1360	`R_xlen_t i, n = XLENGTH(x);`
65020	luke	1361	`PROTECT(s = NO_REFERENCES(x) ? x : allocVector(REALSXP, n));`
73409	luke	1362	`double *pa = REAL(s);`
73732	luke	1363	`const double *px = REAL_RO(x);`
51267	ripley	1364	`for(i = 0 ; i < n ; i++)`
73409	luke	1365	`pa[i] = fabs(px[i]);`
51267	ripley	1366	`} else if (isComplex(x)) {`
68923	ripley	1367	`SET_TAG(args, R_NilValue); /* cmathfuns want "z"; we might have "x" PR#16047 */`
51267	ripley	1368	`return do_cmathfuns(call, op, args, env);`
		1369	`} else`
		1370	`errorcall(call, R_MSG_NONNUM_MATH);`
63462	luke	1371
		1372	`if (x != s && ATTRIB(x) != R_NilValue)`
69108	luke	1373	`SHALLOW_DUPLICATE_ATTRIB(s, x);`
63462	luke	1374	`UNPROTECT(1);`
51267	ripley	1375	`return s;`
18454	ripley	1376	`}`
		1377
7615	maechler	1378	`/* Mathematical Functions of Two Numeric Arguments (plus 1 int) */`
		1379
59448	ripley	1380	`/* math2_1 and math2_2 and related can be removed once the byte`
		1381	`compiler knows how to optimize to .External rather than`
		1382	`.Internal */`
		1383
7615	maechler	1384	`#define if_NA_Math2_set(y,a,b) \`
		1385	`if (ISNA (a) \|\| ISNA (b)) y = NA_REAL; \`
		1386	`else if (ISNAN(a) \|\| ISNAN(b)) y = R_NaN;`
		1387
39866	duncan	1388	`static SEXP math2(SEXP sa, SEXP sb, double (*f)(double, double),`
		1389	`SEXP lcall)`
2	r	1390	`{`
571	ihaka	1391	`SEXP sy;`
59048	ripley	1392	`R_xlen_t i, ia, ib, n, na, nb;`
73732	luke	1393	`double ai, bi, *y;`
		1394	`const double a, b;`
13997	duncan	1395	`int naflag;`
2	r	1396
571	ihaka	1397	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1398	`errorcall(lcall, R_MSG_NONNUM_MATH);`
2	r	1399
27167	maechler	1400	`/* for 0-length a we want the attributes of a, not those of b`
27095	ripley	1401	`as no recycling will occur */`
69108	luke	1402	`#define SETUP_Math2 \`
		1403	`na = XLENGTH(sa); \`
		1404	`nb = XLENGTH(sb); \`
		1405	`if ((na == 0) \|\| (nb == 0)) { \`
		1406	`PROTECT(sy = allocVector(REALSXP, 0)); \`
		1407	`if (na == 0) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \`
		1408	`UNPROTECT(1); \`
		1409	`return(sy); \`
		1410	`} \`
		1411	`n = (na < nb) ? nb : na; \`
		1412	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1413	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1414	`PROTECT(sy = allocVector(REALSXP, n)); \`
73732	luke	1415	`a = REAL_RO(sa); \`
		1416	`b = REAL_RO(sb); \`
69108	luke	1417	`y = REAL(sy); \`
7615	maechler	1418	`naflag = 0`
		1419
		1420	`SETUP_Math2;`
		1421
68663	luke	1422	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
60336	ripley	1423	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
6098	pd	1424	`ai = a[ia];`
		1425	`bi = b[ib];`
7615	maechler	1426	`if_NA_Math2_set(y[i], ai, bi)`
6098	pd	1427	`else {`
29412	ripley	1428	`y[i] = f(ai, bi);`
7615	maechler	1429	`if (ISNAN(y[i])) naflag = 1;`
2	r	1430	`}`
68663	luke	1431	`});`
7615	maechler	1432
69108	luke	1433	`#define FINISH_Math2 \`
		1434	`if(naflag) warning(R_MSG_NA); \`
		1435	`if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \`
		1436	`else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \`
7615	maechler	1437	`UNPROTECT(3)`
		1438
		1439	`FINISH_Math2;`
		1440
		1441	`return sy;`
		1442	`} /* math2() */`
		1443
40141	maechler	1444	`static SEXP math2_1(SEXP sa, SEXP sb, SEXP sI,`
39866	duncan	1445	`double (*f)(double, double, int), SEXP lcall)`
7615	maechler	1446	`{`
		1447	`SEXP sy;`
59090	ripley	1448	`R_xlen_t i, ia, ib, n, na, nb;`
73732	luke	1449	`double ai, bi, *y;`
		1450	`const double a, b;`
7615	maechler	1451	`int m_opt;`
13997	duncan	1452	`int naflag;`
7615	maechler	1453
		1454	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1455	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	1456
		1457	`SETUP_Math2;`
		1458	`m_opt = asInteger(sI);`
		1459
68663	luke	1460	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
60336	ripley	1461	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	1462	`ai = a[ia];`
		1463	`bi = b[ib];`
		1464	`if_NA_Math2_set(y[i], ai, bi)`
		1465	`else {`
29412	ripley	1466	`y[i] = f(ai, bi, m_opt);`
7615	maechler	1467	`if (ISNAN(y[i])) naflag = 1;`
		1468	`}`
68663	luke	1469	`});`
7615	maechler	1470	`FINISH_Math2;`
		1471	`return sy;`
		1472	`} /* math2_1() */`
		1473
40141	maechler	1474	`static SEXP math2_2(SEXP sa, SEXP sb, SEXP sI1, SEXP sI2,`
39866	duncan	1475	`double (*f)(double, double, int, int), SEXP lcall)`
7615	maechler	1476	`{`
		1477	`SEXP sy;`
59090	ripley	1478	`R_xlen_t i, ia, ib, n, na, nb;`
73732	luke	1479	`double ai, bi, *y;`
		1480	`const double a, b;`
7615	maechler	1481	`int i_1, i_2;`
13997	duncan	1482	`int naflag;`
7615	maechler	1483	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
9147	maechler	1484	`errorcall(lcall, R_MSG_NONNUM_MATH);`
7615	maechler	1485
		1486	`SETUP_Math2;`
		1487	`i_1 = asInteger(sI1);`
		1488	`i_2 = asInteger(sI2);`
		1489
68663	luke	1490	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
60336	ripley	1491	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	1492	`ai = a[ia];`
		1493	`bi = b[ib];`
		1494	`if_NA_Math2_set(y[i], ai, bi)`
		1495	`else {`
29412	ripley	1496	`y[i] = f(ai, bi, i_1, i_2);`
7615	maechler	1497	`if (ISNAN(y[i])) naflag = 1;`
		1498	`}`
68663	luke	1499	`});`
7615	maechler	1500	`FINISH_Math2;`
571	ihaka	1501	`return sy;`
7615	maechler	1502	`} /* math2_2() */`
2	r	1503
63171	ripley	1504	`/* This is only used directly by .Internal for Bessel functions,`
		1505	`so managing R_alloc stack is only prudence */`
51657	luke	1506	`static SEXP math2B(SEXP sa, SEXP sb, double (f)(double, double, double ),`
		1507	`SEXP lcall)`
		1508	`{`
		1509	`SEXP sy;`
59090	ripley	1510	`R_xlen_t i, ia, ib, n, na, nb;`
73732	luke	1511	`double ai, bi, *y;`
		1512	`const double a, b;`
51657	luke	1513	`int naflag;`
		1514	`double amax, *work;`
64677	ripley	1515	`size_t nw;`
51657	luke	1516
67884	maechler	1517	`#define besselJY_max_nu 1e7`
		1518
51657	luke	1519	`if (!isNumeric(sa) \|\| !isNumeric(sb))`
		1520	`errorcall(lcall, R_MSG_NONNUM_MATH);`
		1521
		1522	`/* for 0-length a we want the attributes of a, not those of b`
		1523	`as no recycling will occur */`
		1524	`SETUP_Math2;`
		1525
		1526	`/* allocate work array for BesselJ, BesselY large enough for all`
		1527	`arguments */`
		1528	`amax = 0.0;`
		1529	`for (i = 0; i < nb; i++) {`
		1530	`double av = b[i] < 0 ? -b[i] : b[i];`
67884	maechler	1531	`if (amax < av)`
		1532	`amax = av;`
51657	luke	1533	`}`
67884	maechler	1534	`if (amax > besselJY_max_nu)`
		1535	`amax = besselJY_max_nu; // and warning will happen in ../nmath/bessel_[jy].c`
63171	ripley	1536	`const void *vmax = vmaxget();`
64677	ripley	1537	`nw = 1 + (size_t)floor(amax);`
		1538	`work = (double *) R_alloc(nw, sizeof(double));`
51657	luke	1539
68663	luke	1540	`MOD_ITERATE2(n, na, nb, i, ia, ib, {`
60336	ripley	1541	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
51657	luke	1542	`ai = a[ia];`
		1543	`bi = b[ib];`
		1544	`if_NA_Math2_set(y[i], ai, bi)`
		1545	`else {`
		1546	`y[i] = f(ai, bi, work);`
		1547	`if (ISNAN(y[i])) naflag = 1;`
		1548	`}`
68663	luke	1549	`});`
51657	luke	1550
63171	ripley	1551	`vmaxset(vmax);`
51657	luke	1552	`FINISH_Math2;`
		1553
		1554	`return sy;`
		1555	`} /* math2B() */`
		1556
13997	duncan	1557	`#define Math2(A, FUN) math2(CAR(A), CADR(A), FUN, call);`
		1558	`#define Math2_1(A, FUN) math2_1(CAR(A), CADR(A), CADDR(A), FUN, call);`
		1559	`#define Math2_2(A, FUN) math2_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call)`
51657	luke	1560	`#define Math2B(A, FUN) math2B(CAR(A), CADR(A), FUN, call);`
2	r	1561
36990	ripley	1562	`SEXP attribute_hidden do_math2(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1563	`{`
71820	luke	1564	`/* For .Internals, fix up call so errorcall() behaves like error(). */`
		1565	`if (TYPEOF(CAR(call)) == SYMSXP && INTERNAL(CAR(call)) == op)`
		1566	`call = R_CurrentExpression;`
		1567
571	ihaka	1568	`checkArity(op, args);`
2	r	1569
37600	ripley	1570	`if (isComplex(CAR(args)) \|\|`
		1571	`(PRIMVAL(op) == 0 && isComplex(CADR(args))))`
571	ihaka	1572	`return complex_math2(call, op, args, env);`
2	r	1573
14701	ripley	1574
571	ihaka	1575	`switch (PRIMVAL(op)) {`
13997	duncan	1576
7615	maechler	1577	`case 0: return Math2(args, atan2);`
67524	maechler	1578	`case 10001: return Math2(args, fround);// round(), ../nmath/fround.c`
		1579	`case 10004: return Math2(args, fprec); // signif(), ../nmath/fprec.c`
2	r	1580
7615	maechler	1581	`case 2: return Math2(args, lbeta);`
		1582	`case 3: return Math2(args, beta);`
		1583	`case 4: return Math2(args, lchoose);`
		1584	`case 5: return Math2(args, choose);`
2	r	1585
7671	maechler	1586	`case 6: return Math2_1(args, dchisq);`
		1587	`case 7: return Math2_2(args, pchisq);`
		1588	`case 8: return Math2_2(args, qchisq);`
14701	ripley	1589
7671	maechler	1590	`case 9: return Math2_1(args, dexp);`
		1591	`case 10: return Math2_2(args, pexp);`
		1592	`case 11: return Math2_2(args, qexp);`
14701	ripley	1593
7671	maechler	1594	`case 12: return Math2_1(args, dgeom);`
		1595	`case 13: return Math2_2(args, pgeom);`
		1596	`case 14: return Math2_2(args, qgeom);`
14701	ripley	1597
7615	maechler	1598	`case 15: return Math2_1(args, dpois);`
		1599	`case 16: return Math2_2(args, ppois);`
		1600	`case 17: return Math2_2(args, qpois);`
14701	ripley	1601
7671	maechler	1602	`case 18: return Math2_1(args, dt);`
		1603	`case 19: return Math2_2(args, pt);`
		1604	`case 20: return Math2_2(args, qt);`
14701	ripley	1605
7671	maechler	1606	`case 21: return Math2_1(args, dsignrank);`
		1607	`case 22: return Math2_2(args, psignrank);`
		1608	`case 23: return Math2_2(args, qsignrank);`
2239	hornik	1609
51657	luke	1610	`case 24: return Math2B(args, bessel_j_ex);`
		1611	`case 25: return Math2B(args, bessel_y_ex);`
28529	maechler	1612	`case 26: return Math2(args, psigamma);`
2632	maechler	1613
571	ihaka	1614	`default:`
70830	luke	1615	`error(_("unimplemented real function of %d numeric arguments"), 2);`
571	ihaka	1616	`}`
6098	pd	1617	`return op; /* never used; to keep -Wall happy */`
2	r	1618	`}`
		1619
13997	duncan	1620
22644	ripley	1621	`/* The S4 Math2 group, round and signif */`
51245	ripley	1622	`/* This is a primitive SPECIALSXP with internal argument matching */`
36990	ripley	1623	`SEXP attribute_hidden do_Math2(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1624	`{`
66326	luke	1625	`SEXP res, call2;`
46170	ripley	1626	`int n, nprotect = 2;`
66326	luke	1627	`static SEXP do_Math2_formals = NULL;`
13997	duncan	1628
45446	ripley	1629	`if (length(args) >= 2 &&`
		1630	`isSymbol(CADR(args)) && R_isMissing(CADR(args), env)) {`
		1631	`double digits = 0;`
62648	maechler	1632	`if(PRIMVAL(op) == 10004) digits = 6.0; // for signif()`
43757	ripley	1633	`PROTECT(args = list2(CAR(args), ScalarReal(digits))); nprotect++;`
42643	ripley	1634	`}`
		1635
		1636	`PROTECT(args = evalListKeepMissing(args, env));`
46170	ripley	1637	`PROTECT(call2 = lang2(CAR(call), R_NilValue));`
		1638	`SETCDR(call2, args);`
42643	ripley	1639
		1640	`n = length(args);`
63461	luke	1641	`if (n != 1 && n != 2)`
68923	ripley	1642	`error(ngettext("%d argument passed to '%s' which requires 1 or 2 arguments",`
		1643	`"%d arguments passed to '%s'which requires 1 or 2 arguments", n),`
		1644	`n, PRIMNAME(op));`
45446	ripley	1645
46170	ripley	1646	`if (! DispatchGroup("Math", call2, op, args, env, &res)) {`
42643	ripley	1647	`if(n == 1) {`
		1648	`double digits = 0.0;`
		1649	`if(PRIMVAL(op) == 10004) digits = 6.0;`
		1650	`SETCDR(args, CONS(ScalarReal(digits), R_NilValue));`
44313	ripley	1651	`} else {`
		1652	`/* If named, do argument matching by name */`
		1653	`if (TAG(args) != R_NilValue \|\| TAG(CDR(args)) != R_NilValue) {`
68923	ripley	1654	`if (do_Math2_formals == NULL)`
		1655	`do_Math2_formals = allocFormalsList2(install("x"),`
66326	luke	1656	`install("digits"));`
77452	luke	1657	`PROTECT(args = matchArgs_NR(do_Math2_formals, args, call));`
66326	luke	1658	`nprotect++;`
44313	ripley	1659	`}`
		1660	`if (length(CADR(args)) == 0)`
		1661	`errorcall(call, _("invalid second argument of length 0"));`
		1662	`}`
42095	ripley	1663	`res = do_math2(call, op, args, env);`
		1664	`}`
43757	ripley	1665	`UNPROTECT(nprotect);`
42095	ripley	1666	`return res;`
		1667	`}`
		1668
74907	maechler	1669	`/* log{2,10}() builtins : */`
42095	ripley	1670	`SEXP attribute_hidden do_log1arg(SEXP call, SEXP op, SEXP args, SEXP env)`
		1671	`{`
51270	ripley	1672	`SEXP res, call2, args2, tmp = R_NilValue /* -Wall */;`
42095	ripley	1673
22644	ripley	1674	`checkArity(op, args);`
51267	ripley	1675	`check1arg(args, call, "x");`
42095	ripley	1676
		1677	`if (DispatchGroup("Math", call, op, args, env, &res)) return res;`
		1678
68092	luke	1679	`SEXP sLog = install("log");`
42095	ripley	1680	`if(PRIMVAL(op) == 10) tmp = ScalarReal(10.0);`
		1681	`if(PRIMVAL(op) == 2) tmp = ScalarReal(2.0);`
		1682
68092	luke	1683	`PROTECT(call2 = lang3(sLog, CAR(args), tmp));`
51270	ripley	1684	`PROTECT(args2 = lang2(CAR(args), tmp));`
		1685	`if (! DispatchGroup("Math", call2, op, args2, env, &res)) {`
		1686	`if (isComplex(CAR(args)))`
		1687	`res = complex_math2(call2, op, args2, env);`
		1688	`else`
		1689	`res = math2(CAR(args), tmp, logbase, call);`
		1690	`}`
		1691	`UNPROTECT(2);`
42095	ripley	1692	`return res;`
2	r	1693	`}`
		1694
66434	luke	1695	`#ifdef M_E`
		1696	`# define DFLT_LOG_BASE M_E`
		1697	`#else`
		1698	`# define DFLT_LOG_BASE exp(1.)`
		1699	`#endif`
51270	ripley	1700
66605	luke	1701	`/* do_log is a primitive SPECIALSXP with internal argument`
		1702	`matching. do_log_builtin is the BUILTIN version that expects`
		1703	`evaluated arguments to be passed as 'args', expect that these may`
		1704	`contain missing arguments. */`
36990	ripley	1705	`SEXP attribute_hidden do_log(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1706	`{`
66605	luke	1707	`args = evalListKeepMissing(args, env);`
		1708	`return do_log_builtin(call, op, args, env);`
67524	maechler	1709	`}`
66605	luke	1710
		1711	`SEXP attribute_hidden do_log_builtin(SEXP call, SEXP op, SEXP args, SEXP env)`
		1712	`{`
		1713	`PROTECT(args);`
66156	luke	1714	`int n = length(args);`
66431	luke	1715	`SEXP res;`
66156	luke	1716
		1717	`if (n == 1 && TAG(args) == R_NilValue) {`
66431	luke	1718	`/* log(x) is handled here */`
66156	luke	1719	`SEXP x = CAR(args);`
66431	luke	1720	`if (x != R_MissingArg && ! OBJECT(x)) {`
66156	luke	1721	`if (isComplex(x))`
		1722	`res = complex_math1(call, op, args, env);`
		1723	`else`
		1724	`res = math1(x, R_log, call);`
		1725	`UNPROTECT(1);`
		1726	`return res;`
		1727	`}`
		1728	`}`
		1729	`else if (n == 2 &&`
		1730	`TAG(args) == R_NilValue &&`
67524	maechler	1731	`(TAG(CDR(args)) == R_NilValue \|\| TAG(CDR(args)) == R_BaseSymbol)) {`
66156	luke	1732	`/* log(x, y) or log(x, base = y) are handled here */`
		1733	`SEXP x = CAR(args);`
		1734	`SEXP y = CADR(args);`
66431	luke	1735	`if (x != R_MissingArg && y != R_MissingArg &&`
		1736	`! OBJECT(x) && ! OBJECT(y)) {`
66156	luke	1737	`if (isComplex(x) \|\| isComplex(y))`
		1738	`res = complex_math2(call, op, args, env);`
		1739	`else`
		1740	`res = math2(x, y, logbase, call);`
		1741	`UNPROTECT(1);`
		1742	`return res;`
		1743	`}`
		1744	`}`
		1745
66431	luke	1746	`static SEXP do_log_formals = NULL;`
		1747	`static SEXP R_x_Symbol = NULL;`
		1748	`if (do_log_formals == NULL) {`
		1749	`R_x_Symbol = install("x");`
67524	maechler	1750	`do_log_formals = allocFormalsList2(R_x_Symbol, R_BaseSymbol);`
66431	luke	1751	`}`
		1752
66434	luke	1753	`if (n == 1) {`
		1754	`if (CAR(args) == R_MissingArg \|\|`
		1755	`(TAG(args) != R_NilValue && TAG(args) != R_x_Symbol))`
		1756	`error(_("argument \"%s\" is missing, with no default"), "x");`
66432	luke	1757
66435	luke	1758	`if (! DispatchGroup("Math", call, op, args, env, &res)) {`
42095	ripley	1759	`if (isComplex(CAR(args)))`
		1760	`res = complex_math1(call, op, args, env);`
		1761	`else`
		1762	`res = math1(CAR(args), R_log, call);`
66434	luke	1763	`}`
66435	luke	1764	`UNPROTECT(1);`
66434	luke	1765	`return res;`
		1766	`}`
		1767	`else {`
		1768	`/* match argument names if supplied */`
		1769	`/* will signal an error unless there are one or two arguments */`
		1770	`/* after the match, length(args) will be 2 */`
77452	luke	1771	`PROTECT(args = matchArgs_NR(do_log_formals, args, call));`
66434	luke	1772
		1773	`if(CAR(args) == R_MissingArg)`
		1774	`error(_("argument \"%s\" is missing, with no default"), "x");`
		1775	`if (CADR(args) == R_MissingArg)`
		1776	`SETCADR(args, ScalarReal(DFLT_LOG_BASE));`
		1777
66435	luke	1778	`if (! DispatchGroup("Math", call, op, args, env, &res)) {`
42095	ripley	1779	`if (length(CADR(args)) == 0)`
		1780	`errorcall(call, _("invalid argument 'base' of length 0"));`
42101	ripley	1781	`if (isComplex(CAR(args)) \|\| isComplex(CADR(args)))`
42095	ripley	1782	`res = complex_math2(call, op, args, env);`
		1783	`else`
		1784	`res = math2(CAR(args), CADR(args), logbase, call);`
		1785	`}`
66435	luke	1786	`UNPROTECT(2);`
66434	luke	1787	`return res;`
42093	ripley	1788	`}`
2	r	1789	`}`
		1790
13997	duncan	1791
7615	maechler	1792	`/* Mathematical Functions of Three (Real) Arguments */`
		1793
59448	ripley	1794	`/* math3_1 and math3_2 and related can be removed once the byte`
		1795	`compiler knows how to optimize to .External rather than`
		1796	`.Internal */`
		1797
		1798
7615	maechler	1799	`#define if_NA_Math3_set(y,a,b,c) \`
		1800	`if (ISNA (a) \|\| ISNA (b)\|\| ISNA (c)) y = NA_REAL; \`
		1801	`else if (ISNAN(a) \|\| ISNAN(b)\|\| ISNAN(c)) y = R_NaN;`
		1802
		1803	`#define SETUP_Math3 \`
		1804	`if (!isNumeric(sa) \|\| !isNumeric(sb) \|\| !isNumeric(sc)) \`
70830	luke	1805	`error(R_MSG_NONNUM_MATH); \`
7615	maechler	1806	`\`
59048	ripley	1807	`na = XLENGTH(sa); \`
		1808	`nb = XLENGTH(sb); \`
		1809	`nc = XLENGTH(sc); \`
7615	maechler	1810	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0)) \`
		1811	`return(allocVector(REALSXP, 0)); \`
		1812	`n = na; \`
		1813	`if (n < nb) n = nb; \`
		1814	`if (n < nc) n = nc; \`
		1815	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		1816	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		1817	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		1818	`PROTECT(sy = allocVector(REALSXP, n)); \`
73732	luke	1819	`a = REAL_RO(sa); \`
		1820	`b = REAL_RO(sb); \`
		1821	`c = REAL_RO(sc); \`
7615	maechler	1822	`y = REAL(sy); \`
		1823	`naflag = 0`
834	maechler	1824
69108	luke	1825	`#define FINISH_Math3 \`
		1826	`if(naflag) warning(R_MSG_NA); \`
69175	ripley	1827	`\`
69108	luke	1828	`if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \`
		1829	`else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \`
		1830	`else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \`
51790	ripley	1831	`UNPROTECT(4)`
7615	maechler	1832
40141	maechler	1833	`static SEXP math3_1(SEXP sa, SEXP sb, SEXP sc, SEXP sI,`
39866	duncan	1834	`double (*f)(double, double, double, int), SEXP lcall)`
7615	maechler	1835	`{`
		1836	`SEXP sy;`
59048	ripley	1837	`R_xlen_t i, ia, ib, ic, n, na, nb, nc;`
73732	luke	1838	`double ai, bi, ci, *y;`
		1839	`const double a, b, *c;`
7615	maechler	1840	`int i_1;`
13997	duncan	1841	`int naflag;`
7615	maechler	1842
		1843	`SETUP_Math3;`
		1844	`i_1 = asInteger(sI);`
		1845
68663	luke	1846	`MOD_ITERATE3(n, na, nb, nc, i, ia, ib, ic, {`
60336	ripley	1847	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	1848	`ai = a[ia];`
		1849	`bi = b[ib];`
		1850	`ci = c[ic];`
		1851	`if_NA_Math3_set(y[i], ai,bi,ci)`
		1852	`else {`
29412	ripley	1853	`y[i] = f(ai, bi, ci, i_1);`
7615	maechler	1854	`if (ISNAN(y[i])) naflag = 1;`
		1855	`}`
68663	luke	1856	`});`
7615	maechler	1857
		1858	`FINISH_Math3;`
		1859	`return sy;`
		1860	`} /* math3_1 */`
		1861
40141	maechler	1862	`static SEXP math3_2(SEXP sa, SEXP sb, SEXP sc, SEXP sI, SEXP sJ,`
39866	duncan	1863	`double (*f)(double, double, double, int, int), SEXP lcall)`
7615	maechler	1864	`{`
		1865	`SEXP sy;`
59048	ripley	1866	`R_xlen_t i, ia, ib, ic, n, na, nb, nc;`
73732	luke	1867	`double ai, bi, ci, *y;`
		1868	`const double a, b, *c;`
7615	maechler	1869	`int i_1,i_2;`
13997	duncan	1870	`int naflag;`
7615	maechler	1871
		1872	`SETUP_Math3;`
		1873	`i_1 = asInteger(sI);`
		1874	`i_2 = asInteger(sJ);`
		1875
68663	luke	1876	`MOD_ITERATE3 (n, na, nb, nc, i, ia, ib, ic, {`
60336	ripley	1877	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	1878	`ai = a[ia];`
		1879	`bi = b[ib];`
		1880	`ci = c[ic];`
		1881	`if_NA_Math3_set(y[i], ai,bi,ci)`
		1882	`else {`
29412	ripley	1883	`y[i] = f(ai, bi, ci, i_1, i_2);`
7615	maechler	1884	`if (ISNAN(y[i])) naflag = 1;`
		1885	`}`
68663	luke	1886	`});`
7615	maechler	1887
		1888	`FINISH_Math3;`
571	ihaka	1889	`return sy;`
7615	maechler	1890	`} /* math3_2 */`
2	r	1891
63171	ripley	1892	`/* This is only used directly by .Internal for Bessel functions,`
		1893	`so managing R_alloc stack is only prudence */`
51657	luke	1894	`static SEXP math3B(SEXP sa, SEXP sb, SEXP sc,`
		1895	`double (f)(double, double, double, double ), SEXP lcall)`
		1896	`{`
		1897	`SEXP sy;`
59048	ripley	1898	`R_xlen_t i, ia, ib, ic, n, na, nb, nc;`
73732	luke	1899	`double ai, bi, ci, *y;`
		1900	`const double a, b, *c;`
51657	luke	1901	`int naflag;`
		1902	`double amax, *work;`
64677	ripley	1903	`size_t nw;`
51657	luke	1904
		1905	`SETUP_Math3;`
		1906
		1907	`/* allocate work array for BesselI, BesselK large enough for all`
		1908	`arguments */`
		1909	`amax = 0.0;`
		1910	`for (i = 0; i < nb; i++) {`
		1911	`double av = b[i] < 0 ? -b[i] : b[i];`
		1912	`if (av > amax) amax = av;`
		1913	`}`
63171	ripley	1914	`const void *vmax = vmaxget();`
64677	ripley	1915	`nw = 1 + (size_t)floor(amax);`
		1916	`work = (double *) R_alloc(nw, sizeof(double));`
51657	luke	1917
68663	luke	1918	`MOD_ITERATE3 (n, na, nb, nc, i, ia, ib, ic, {`
60336	ripley	1919	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
51657	luke	1920	`ai = a[ia];`
		1921	`bi = b[ib];`
		1922	`ci = c[ic];`
		1923	`if_NA_Math3_set(y[i], ai,bi,ci)`
		1924	`else {`
		1925	`y[i] = f(ai, bi, ci, work);`
		1926	`if (ISNAN(y[i])) naflag = 1;`
		1927	`}`
68663	luke	1928	`});`
51657	luke	1929
		1930	`FINISH_Math3;`
63171	ripley	1931	`vmaxset(vmax);`
51657	luke	1932
		1933	`return sy;`
		1934	`} /* math3B */`
		1935
13997	duncan	1936	`#define Math3_1(A, FUN) math3_1(CAR(A), CADR(A), CADDR(A), CADDDR(A), FUN, call);`
		1937	`#define Math3_2(A, FUN) math3_2(CAR(A), CADR(A), CADDR(A), CADDDR(A), CAD4R(A), FUN, call)`
51657	luke	1938	`#define Math3B(A, FUN) math3B (CAR(A), CADR(A), CADDR(A), FUN, call);`
2	r	1939
36990	ripley	1940	`SEXP attribute_hidden do_math3(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	1941	`{`
571	ihaka	1942	`checkArity(op, args);`
2	r	1943
571	ihaka	1944	`switch (PRIMVAL(op)) {`
2	r	1945
7671	maechler	1946	`case 1: return Math3_1(args, dbeta);`
		1947	`case 2: return Math3_2(args, pbeta);`
		1948	`case 3: return Math3_2(args, qbeta);`
2	r	1949
7671	maechler	1950	`case 4: return Math3_1(args, dbinom);`
		1951	`case 5: return Math3_2(args, pbinom);`
		1952	`case 6: return Math3_2(args, qbinom);`
2	r	1953
7671	maechler	1954	`case 7: return Math3_1(args, dcauchy);`
		1955	`case 8: return Math3_2(args, pcauchy);`
		1956	`case 9: return Math3_2(args, qcauchy);`
2	r	1957
7671	maechler	1958	`case 10: return Math3_1(args, df);`
		1959	`case 11: return Math3_2(args, pf);`
		1960	`case 12: return Math3_2(args, qf);`
2	r	1961
7671	maechler	1962	`case 13: return Math3_1(args, dgamma);`
		1963	`case 14: return Math3_2(args, pgamma);`
		1964	`case 15: return Math3_2(args, qgamma);`
2	r	1965
7639	maechler	1966	`case 16: return Math3_1(args, dlnorm);`
		1967	`case 17: return Math3_2(args, plnorm);`
		1968	`case 18: return Math3_2(args, qlnorm);`
2	r	1969
7671	maechler	1970	`case 19: return Math3_1(args, dlogis);`
		1971	`case 20: return Math3_2(args, plogis);`
		1972	`case 21: return Math3_2(args, qlogis);`
2	r	1973
7671	maechler	1974	`case 22: return Math3_1(args, dnbinom);`
		1975	`case 23: return Math3_2(args, pnbinom);`
		1976	`case 24: return Math3_2(args, qnbinom);`
2	r	1977
7639	maechler	1978	`case 25: return Math3_1(args, dnorm);`
		1979	`case 26: return Math3_2(args, pnorm);`
		1980	`case 27: return Math3_2(args, qnorm);`
2	r	1981
7671	maechler	1982	`case 28: return Math3_1(args, dunif);`
		1983	`case 29: return Math3_2(args, punif);`
		1984	`case 30: return Math3_2(args, qunif);`
2	r	1985
7671	maechler	1986	`case 31: return Math3_1(args, dweibull);`
		1987	`case 32: return Math3_2(args, pweibull);`
		1988	`case 33: return Math3_2(args, qweibull);`
2	r	1989
7671	maechler	1990	`case 34: return Math3_1(args, dnchisq);`
		1991	`case 35: return Math3_2(args, pnchisq);`
		1992	`case 36: return Math3_2(args, qnchisq);`
2	r	1993
7671	maechler	1994	`case 37: return Math3_1(args, dnt);`
		1995	`case 38: return Math3_2(args, pnt);`
		1996	`case 39: return Math3_2(args, qnt);`
602	ihaka	1997
7671	maechler	1998	`case 40: return Math3_1(args, dwilcox);`
		1999	`case 41: return Math3_2(args, pwilcox);`
		2000	`case 42: return Math3_2(args, qwilcox);`
1276	hornik	2001
51657	luke	2002	`case 43: return Math3B(args, bessel_i_ex);`
		2003	`case 44: return Math3B(args, bessel_k_ex);`
2632	maechler	2004
46039	maechler	2005	`case 45: return Math3_1(args, dnbinom_mu);`
		2006	`case 46: return Math3_2(args, pnbinom_mu);`
		2007	`case 47: return Math3_2(args, qnbinom_mu);`
2632	maechler	2008
571	ihaka	2009	`default:`
70830	luke	2010	`error(_("unimplemented real function of %d numeric arguments"), 3);`
571	ihaka	2011	`}`
6098	pd	2012	`return op; /* never used; to keep -Wall happy */`
4719	maechler	2013	`} /* do_math3() */`
988	maechler	2014
7615	maechler	2015	`/* Mathematical Functions of Four (Real) Arguments */`
		2016
59448	ripley	2017	`/* This can be removed completely once the byte compiler knows how to`
		2018	`optimize to .External rather than .Internal */`
		2019
45446	ripley	2020	`#define if_NA_Math4_set(y,a,b,c,d) \`
7615	maechler	2021	`if (ISNA (a)\|\| ISNA (b)\|\| ISNA (c)\|\| ISNA (d)) y = NA_REAL;\`
		2022	`else if (ISNAN(a)\|\| ISNAN(b)\|\| ISNAN(c)\|\| ISNAN(d)) y = R_NaN;`
		2023
40141	maechler	2024	`static SEXP math4(SEXP sa, SEXP sb, SEXP sc, SEXP sd,`
39866	duncan	2025	`double (*f)(double, double, double, double), SEXP lcall)`
2	r	2026	`{`
571	ihaka	2027	`SEXP sy;`
59048	ripley	2028	`R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd;`
73732	luke	2029	`double ai, bi, ci, di, *y;`
		2030	`const double a, b, c, d;`
13997	duncan	2031	`int naflag;`
2	r	2032
7615	maechler	2033	`#define SETUP_Math4 \`
		2034	`if(!isNumeric(sa)\|\| !isNumeric(sb)\|\| !isNumeric(sc)\|\| !isNumeric(sd))\`
70830	luke	2035	`error(R_MSG_NONNUM_MATH); \`
7615	maechler	2036	`\`
59048	ripley	2037	`na = XLENGTH(sa); \`
		2038	`nb = XLENGTH(sb); \`
		2039	`nc = XLENGTH(sc); \`
		2040	`nd = XLENGTH(sd); \`
7615	maechler	2041	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0) \|\| (nd == 0)) \`
		2042	`return(allocVector(REALSXP, 0)); \`
		2043	`n = na; \`
		2044	`if (n < nb) n = nb; \`
		2045	`if (n < nc) n = nc; \`
		2046	`if (n < nd) n = nd; \`
		2047	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		2048	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		2049	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		2050	`PROTECT(sd = coerceVector(sd, REALSXP)); \`
		2051	`PROTECT(sy = allocVector(REALSXP, n)); \`
73732	luke	2052	`a = REAL_RO(sa); \`
		2053	`b = REAL_RO(sb); \`
		2054	`c = REAL_RO(sc); \`
		2055	`d = REAL_RO(sd); \`
7615	maechler	2056	`y = REAL(sy); \`
		2057	`naflag = 0`
834	maechler	2058
7615	maechler	2059	`SETUP_Math4;`
		2060
68663	luke	2061	`MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, {`
60336	ripley	2062	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
6098	pd	2063	`ai = a[ia];`
		2064	`bi = b[ib];`
		2065	`ci = c[ic];`
		2066	`di = d[id];`
7615	maechler	2067	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
6098	pd	2068	`else {`
29412	ripley	2069	`y[i] = f(ai, bi, ci, di);`
7615	maechler	2070	`if (ISNAN(y[i])) naflag = 1;`
2	r	2071	`}`
68663	luke	2072	`});`
7615	maechler	2073
69108	luke	2074	`#define FINISH_Math4 \`
		2075	`if(naflag) warning(R_MSG_NA); \`
69175	ripley	2076	`\`
69108	luke	2077	`if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \`
		2078	`else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \`
		2079	`else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \`
		2080	`else if (n == nd) SHALLOW_DUPLICATE_ATTRIB(sy, sd); \`
7615	maechler	2081	`UNPROTECT(5)`
		2082
		2083	`FINISH_Math4;`
		2084
		2085	`return sy;`
		2086	`} /* math4() */`
		2087
39866	duncan	2088	`static SEXP math4_1(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, double (*f)(double, double, double, double, int), SEXP lcall)`
7615	maechler	2089	`{`
		2090	`SEXP sy;`
59048	ripley	2091	`R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd;`
73732	luke	2092	`double ai, bi, ci, di, *y;`
		2093	`const double a, b, c, d;`
7615	maechler	2094	`int i_1;`
13997	duncan	2095	`int naflag;`
7615	maechler	2096
		2097	`SETUP_Math4;`
		2098	`i_1 = asInteger(sI);`
		2099
68663	luke	2100	`MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, {`
60336	ripley	2101	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	2102	`ai = a[ia];`
		2103	`bi = b[ib];`
		2104	`ci = c[ic];`
		2105	`di = d[id];`
		2106	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
		2107	`else {`
29412	ripley	2108	`y[i] = f(ai, bi, ci, di, i_1);`
7615	maechler	2109	`if (ISNAN(y[i])) naflag = 1;`
		2110	`}`
68663	luke	2111	`});`
7615	maechler	2112	`FINISH_Math4;`
		2113	`return sy;`
		2114	`} /* math4_1() */`
		2115
		2116	`static SEXP math4_2(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP sI, SEXP sJ,`
39866	duncan	2117	`double (*f)(double, double, double, double, int, int), SEXP lcall)`
7615	maechler	2118	`{`
		2119	`SEXP sy;`
59048	ripley	2120	`R_xlen_t i, ia, ib, ic, id, n, na, nb, nc, nd;`
73732	luke	2121	`double ai, bi, ci, di, *y;`
		2122	`const double a, b, c, d;`
7615	maechler	2123	`int i_1, i_2;`
13997	duncan	2124	`int naflag;`
7615	maechler	2125
		2126	`SETUP_Math4;`
		2127	`i_1 = asInteger(sI);`
		2128	`i_2 = asInteger(sJ);`
		2129
68663	luke	2130	`MOD_ITERATE4 (n, na, nb, nc, nd, i, ia, ib, ic, id, {`
60336	ripley	2131	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
7615	maechler	2132	`ai = a[ia];`
		2133	`bi = b[ib];`
		2134	`ci = c[ic];`
		2135	`di = d[id];`
		2136	`if_NA_Math4_set(y[i], ai,bi,ci,di)`
		2137	`else {`
29412	ripley	2138	`y[i] = f(ai, bi, ci, di, i_1, i_2);`
7615	maechler	2139	`if (ISNAN(y[i])) naflag = 1;`
		2140	`}`
68663	luke	2141	`});`
7615	maechler	2142	`FINISH_Math4;`
571	ihaka	2143	`return sy;`
7615	maechler	2144	`} /* math4_2() */`
2	r	2145
		2146
7615	maechler	2147	`#define CAD3R CADDDR`
		2148	`/* This is not (yet) in Rinternals.h : */`
7671	maechler	2149	`#define CAD5R(e) CAR(CDR(CDR(CDR(CDR(CDR(e))))))`
4719	maechler	2150
13997	duncan	2151	`#define Math4(A, FUN) math4 (CAR(A), CADR(A), CADDR(A), CAD3R(A), FUN, call)`
7671	maechler	2152	`#define Math4_1(A, FUN) math4_1(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \`
13997	duncan	2153	`FUN, call)`
7671	maechler	2154	`#define Math4_2(A, FUN) math4_2(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), \`
13997	duncan	2155	`CAD5R(A), FUN, call)`
7615	maechler	2156
7671	maechler	2157
36990	ripley	2158	`SEXP attribute_hidden do_math4(SEXP call, SEXP op, SEXP args, SEXP env)`
2	r	2159	`{`
571	ihaka	2160	`checkArity(op, args);`
2	r	2161
13997	duncan	2162
571	ihaka	2163	`switch (PRIMVAL(op)) {`
7671	maechler	2164
		2165	`/* Completely dummy for -Wall -- math4() at all! : */`
39866	duncan	2166	`case -99: return Math4(args, (double (*)(double, double, double, double))NULL);`
7671	maechler	2167
		2168	`case 1: return Math4_1(args, dhyper);`
		2169	`case 2: return Math4_2(args, phyper);`
		2170	`case 3: return Math4_2(args, qhyper);`
		2171
		2172	`case 4: return Math4_1(args, dnbeta);`
		2173	`case 5: return Math4_2(args, pnbeta);`
		2174	`case 6: return Math4_2(args, qnbeta);`
		2175	`case 7: return Math4_1(args, dnf);`
		2176	`case 8: return Math4_2(args, pnf);`
		2177	`case 9: return Math4_2(args, qnf);`
623	ihaka	2178	`#ifdef UNIMP`
7671	maechler	2179	`case 10: return Math4_1(args, dtukey);`
1350	ihaka	2180	`#endif`
7671	maechler	2181	`case 11: return Math4_2(args, ptukey);`
		2182	`case 12: return Math4_2(args, qtukey);`
571	ihaka	2183	`default:`
70830	luke	2184	`error(_("unimplemented real function of %d numeric arguments"), 4);`
571	ihaka	2185	`}`
6098	pd	2186	`return op; /* never used; to keep -Wall happy */`
2	r	2187	`}`
4719	maechler	2188
		2189
10866	maechler	2190	`#ifdef WHEN_MATH5_IS_THERE/* as in ./arithmetic.h */`
4719	maechler	2191
7615	maechler	2192	`/* Mathematical Functions of Five (Real) Arguments */`
		2193
		2194	`#define if_NA_Math5_set(y,a,b,c,d,e) \`
45446	ripley	2195	`if (ISNA (a)\|\| ISNA (b)\|\| ISNA (c)\|\| ISNA (d)\|\| ISNA (e)) \`
		2196	`y = NA_REAL; \`
7615	maechler	2197	`else if(ISNAN(a)\|\| ISNAN(b)\|\| ISNAN(c)\|\| ISNAN(d)\|\| ISNAN(e)) \`
		2198	`y = R_NaN;`
		2199
		2200	`static SEXP math5(SEXP sa, SEXP sb, SEXP sc, SEXP sd, SEXP se, double (*f)())`
4719	maechler	2201	`{`
		2202	`SEXP sy;`
59048	ripley	2203	`R_xlen_t i, ia, ib, ic, id, ie, n, na, nb, nc, nd, ne;`
73732	luke	2204	`double ai, bi, ci, di, ei, *y;`
		2205	`const double a, b, c, d, *e;`
4719	maechler	2206
7615	maechler	2207	`#define SETUP_Math5 \`
		2208	`if (!isNumeric(sa) \|\| !isNumeric(sb) \|\| !isNumeric(sc) \|\| \`
		2209	`!isNumeric(sd) \|\| !isNumeric(se)) \`
70830	luke	2210	`error(R_MSG_NONNUM_MATH); \`
7615	maechler	2211	`\`
59048	ripley	2212	`na = XLENGTH(sa); \`
		2213	`nb = XLENGTH(sb); \`
		2214	`nc = XLENGTH(sc); \`
		2215	`nd = XLENGTH(sd); \`
		2216	`ne = XLENGTH(se); \`
7615	maechler	2217	`if ((na == 0) \|\| (nb == 0) \|\| (nc == 0) \|\| (nd == 0) \|\| (ne == 0)) \`
		2218	`return(allocVector(REALSXP, 0)); \`
		2219	`n = na; \`
		2220	`if (n < nb) n = nb; \`
		2221	`if (n < nc) n = nc; \`
		2222	`if (n < nd) n = nd; \`
		2223	`if (n < ne) n = ne; /* n = max(na,nb,nc,nd,ne) */ \`
		2224	`PROTECT(sa = coerceVector(sa, REALSXP)); \`
		2225	`PROTECT(sb = coerceVector(sb, REALSXP)); \`
		2226	`PROTECT(sc = coerceVector(sc, REALSXP)); \`
		2227	`PROTECT(sd = coerceVector(sd, REALSXP)); \`
		2228	`PROTECT(se = coerceVector(se, REALSXP)); \`
		2229	`PROTECT(sy = allocVector(REALSXP, n)); \`
73732	luke	2230	`a = REAL_RO(sa); \`
		2231	`b = REAL_RO(sb); \`
		2232	`c = REAL_RO(sc); \`
		2233	`d = REAL_RO(sd); \`
		2234	`e = REAL_RO(se); \`
7615	maechler	2235	`y = REAL(sy); \`
		2236	`naflag = 0`
4719	maechler	2237
7615	maechler	2238	`SETUP_Math5;`
		2239
68663	luke	2240	`MOD_ITERATE5 (n, na, nb, nc, nd, ne,`
		2241	`i, ia, ib, ic, id, ie, {`
60336	ripley	2242	`// if ((i+1) % NINTERRUPT == 0) R_CheckUserInterrupt();`
6098	pd	2243	`ai = a[ia];`
		2244	`bi = b[ib];`
		2245	`ci = c[ic];`
		2246	`di = d[id];`
		2247	`ei = e[ie];`
7615	maechler	2248	`if_NA_Math5_set(y[i], ai,bi,ci,di,ei)`
6098	pd	2249	`else {`
29412	ripley	2250	`y[i] = f(ai, bi, ci, di, ei);`
7615	maechler	2251	`if (ISNAN(y[i])) naflag = 1;`
4719	maechler	2252	`}`
68663	luke	2253	`});`
7615	maechler	2254
69108	luke	2255	`#define FINISH_Math5 \`
		2256	`if(naflag) warning(R_MSG_NA); \`
69175	ripley	2257	`\`
69108	luke	2258	`if (n == na) SHALLOW_DUPLICATE_ATTRIB(sy, sa); \`
		2259	`else if (n == nb) SHALLOW_DUPLICATE_ATTRIB(sy, sb); \`
		2260	`else if (n == nc) SHALLOW_DUPLICATE_ATTRIB(sy, sc); \`
		2261	`else if (n == nd) SHALLOW_DUPLICATE_ATTRIB(sy, sd); \`
		2262	`else if (n == ne) SHALLOW_DUPLICATE_ATTRIB(sy, se); \`
7615	maechler	2263	`UNPROTECT(6)`
		2264
		2265	`FINISH_Math5;`
		2266
4719	maechler	2267	`return sy;`
		2268	`} /* math5() */`
		2269
7615	maechler	2270	`#define Math5(A, FUN) \`
		2271	`math5(CAR(A), CADR(A), CADDR(A), CAD3R(A), CAD4R(A), FUN);`
4719	maechler	2272
36990	ripley	2273	`SEXP attribute_hidden do_math5(SEXP call, SEXP op, SEXP args, SEXP env)`
4719	maechler	2274	`{`
		2275	`checkArity(op, args);`
4859	maechler	2276	`lcall = call;`
4719	maechler	2277
		2278	`switch (PRIMVAL(op)) {`
7671	maechler	2279
		2280	`/* Completely dummy for -Wall -- use math5() at all! : */`
		2281	`case -99: return Math5(args, dhyper);`
4719	maechler	2282	`#ifdef UNIMP`
7671	maechler	2283	`case 2: return Math5(args, p...);`
		2284	`case 3: return Math5(args, q...);`
4719	maechler	2285	`#endif`
		2286	`default:`
70830	luke	2287	`error(_("unimplemented real function of %d numeric arguments"), 5);`
4719	maechler	2288	`}`
6098	pd	2289	`return op; /* never used; to keep -Wall happy */`
4719	maechler	2290	`} /* do_math5() */`
		2291
6098	pd	2292	`#endif /* Math5 is there */`
42439	luke	2293
		2294	`/* This is used for experimenting with parallelized nmath functions -- LT */`
		2295	`CCODE R_get_arith_function(int which)`
		2296	`{`
		2297	`switch (which) {`
		2298	`case 1: return do_math1;`
		2299	`case 2: return do_math2;`
		2300	`case 3: return do_math3;`
		2301	`case 4: return do_math4;`
		2302	`case 11: return complex_math1;`
		2303	`case 12: return complex_math2;`
		2304	`default: error("bad arith function index"); return NULL;`
		2305	`}`
		2306	`}`

The R Project SVN R

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