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

The R Project SVN R

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