The R Project SVN R

Rev

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

Rev Author Line No. Line
3279 pd 1
/*
2
 *  Mathlib : A C Library of Special Functions
5458 ripley 3
 *  Copyright (C) 1998-1999 R Development Core Team
3279 pd 4
 *
5
 *  This program is free software; you can redistribute it and/or modify
6
 *  it under the terms of the GNU General Public License as published by
7
 *  the Free Software Foundation; either version 2 of the License, or
8
 *  (at your option) any later version.
9
 *
10
 *  This program is distributed in the hope that it will be useful,
11
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13
 *  GNU General Public License for more details.
14
 *
15
 *  You should have received a copy of the GNU General Public License
16
 *  along with this program; if not, write to the Free Software
5458 ripley 17
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
3279 pd 18
 *
19
 
20
 * Mathlib.h  should contain ALL headers from R's C code in `src/nmath'
21
   ---------  such that ``the Math library'' can be used by simply
22
 
23
   ``#include "Mathlib.h" ''
24
 
25
   and nothing else.
26
*/
571 ihaka 27
#ifndef MATHLIB_H
28
#define MATHLIB_H
2 r 29
 
3279 pd 30
/*-- Mathlib as part of R --  undefine this for standalone : */
31
#define MATHLIB_IN_R
3076 pd 32
 
7108 ripley 33
#include "R_ext/Rver.h" /* for IEEE_754 */
7003 ripley 34
#include "R_ext/Arith.h"
35
/*#include "R_ext/Random.h"*/
2 r 36
 
2629 maechler 37
#ifdef FORTRAN_H
38
#error __MUST__include "Mathlib.h"  _before_  "Fortran.h"
39
#endif
40
 
571 ihaka 41
#include <errno.h>
2737 hornik 42
#include <limits.h>
571 ihaka 43
#include <float.h>
2 r 44
#include <math.h>
45
#include <stdlib.h>
46
 
2629 maechler 47
/* TRUE and FALSE conflict with the Mac --- Fortran.h still defines them... */
48
#define LTRUE	(1)
49
#define LFALSE	(0)
50
 
3244 ihaka 51
/* 30 Decimal-place constants */
52
/* Computed with bc -l (scale=32; proper round) */
2 r 53
 
3244 ihaka 54
/* SVID & X/Open Constants */
55
/* Names from Solaris math.h */
56
 
57
#ifndef M_E
58
#define M_E		2.718281828459045235360287471353	/* e */
59
#endif
60
 
61
#ifndef M_LOG2E
62
#define M_LOG2E		1.442695040888963407359924681002	/* log2(e) */
63
#endif
64
 
65
#ifndef M_LOG10E
4835 maechler 66
#define M_LOG10E	0.434294481903251827651128918917	/* log10(e) */
3244 ihaka 67
#endif
68
 
69
#ifndef M_LN2
70
#define M_LN2		0.693147180559945309417232121458	/* ln(2) */
71
#endif
72
 
73
#ifndef M_LN10
74
#define M_LN10		2.302585092994045684017991454684	/* ln(10) */
75
#endif
76
 
77
#ifndef M_PI
4835 maechler 78
#define M_PI		3.141592653589793238462643383280	/* pi */
3244 ihaka 79
#endif
80
 
81
#ifndef M_PI_2
82
#define M_PI_2		1.570796326794896619231321691640	/* pi/2 */
83
#endif
84
 
85
#ifndef M_PI_4
86
#define M_PI_4		0.785398163397448309615660845820	/* pi/4 */
87
#endif
88
 
4835 maechler 89
#ifndef M_1_PI
3244 ihaka 90
#define M_1_PI		0.318309886183790671537767526745	/* 1/pi */
91
#endif
92
 
4835 maechler 93
#ifndef M_2_PI
3244 ihaka 94
#define M_2_PI		0.636619772367581343075535053490	/* 1/pi */
95
#endif
96
 
97
#ifndef M_2_SQRTPI
98
#define M_2_SQRTPI	1.128379167095512573896158903122	/* 1/sqrt(pi) */
99
#endif
100
 
101
#ifndef M_SQRT2
102
#define M_SQRT2		1.414213562373095048801688724210	/* sqrt(2) */
103
#endif
104
 
105
#ifndef M_SQRT1_2
106
#define M_SQRT1_2	0.707106781186547524400844362105	/* 1/sqrt(2) */
107
#endif
108
 
109
/* Other, R-Specific Constants */
110
/* Note there are some repeats of values above */
111
/* Needs a cleanup */
112
 
113
#ifndef M_1_SQRT_2
777 maechler 114
#define M_1_SQRT_2	0.707106781186547524400844362105	/* 1/sqrt(2) */
3244 ihaka 115
#endif
116
 
117
#ifndef M_SQRT_32
777 maechler 118
#define M_SQRT_32	5.656854249492380195206754896838	/* sqrt(32) */
2 r 119
#endif
120
 
3244 ihaka 121
#ifndef M_LOG10_2
122
#define M_LOG10_2	0.301029995663981195213738894724	/* log10(2) */
777 maechler 123
#endif
124
 
2 r 125
#ifndef M_PI_half
3244 ihaka 126
#define M_PI_half	1.570796326794896619231321691640	/* pi/2 */
2 r 127
#endif
128
 
129
#ifndef M_SQRT_PI
3244 ihaka 130
#define M_SQRT_PI	1.772453850905516027298167483341	/* sqrt(pi) */
2 r 131
#endif
132
 
3244 ihaka 133
#ifndef M_1_SQRT_2PI
134
#define M_1_SQRT_2PI	0.398942280401432677939946059934	/* 1/sqrt(2pi) */
135
#endif
2 r 136
 
3244 ihaka 137
#ifndef M_SQRT_2dPI
138
#define M_SQRT_2dPI	0.797884560802865355879892119869	/* sqrt(2/pi) */
139
#endif
140
 
141
 
777 maechler 142
#ifndef M_LN_SQRT_PI
3244 ihaka 143
#define M_LN_SQRT_PI	0.572364942924700087071713675677	/* log(sqrt(pi)) */
2 r 144
#endif
145
 
3244 ihaka 146
#ifndef M_LN_SQRT_2PI
147
#define M_LN_SQRT_2PI	0.918938533204672741780329736406	/* log(sqrt(2*pi)) */
148
#endif
777 maechler 149
 
3244 ihaka 150
#ifndef M_LN_SQRT_PId2
151
#define M_LN_SQRT_PId2	0.225791352644727432363097614947	/* log(sqrt(pi/2)) */
152
#endif
153
 
154
 
3076 pd 155
#ifdef MATHLIB_IN_R/* Mathlib in R */
777 maechler 156
 
7007 ripley 157
#include "R_ext/Error.h"
3076 pd 158
# define MATHLIB_ERROR(fmt,x)		error(fmt,x);
159
# define MATHLIB_WARNING(fmt,x)		warning(fmt,x)
160
# define MATHLIB_WARNING2(fmt,x,x2)	warning(fmt,x,x2)
161
# define MATHLIB_WARNING3(fmt,x,x2,x3)	warning(fmt,x,x2,x3)
162
# define MATHLIB_WARNING4(fmt,x,x2,x3,x4) warning(fmt,x,x2,x3,x4)
2 r 163
 
3076 pd 164
#else/* Mathlib standalone */
165
 
166
#include <stdio.h>
3786 pd 167
# define MATHLIB_ERROR(fmt,x)	{ printf(fmt,x); exit(1) }
3076 pd 168
# define MATHLIB_WARNING(fmt,x)		printf(fmt,x)
169
# define MATHLIB_WARNING2(fmt,x,x2)	printf(fmt,x,x2)
170
# define MATHLIB_WARNING3(fmt,x,x2,x3)	printf(fmt,x,x2,x3)
171
# define MATHLIB_WARNING4(fmt,x,x2,x3,x4) printf(fmt,x,x2,x3,x4)
172
#endif
173
 
571 ihaka 174
#define ME_NONE		0
3786 pd 175
/*	no error */
571 ihaka 176
#define ME_DOMAIN	1
3786 pd 177
/*	argument out of domain */
571 ihaka 178
#define ME_RANGE	2
3786 pd 179
/*	value out of range */
180
#define ME_NOCONV	4
181
/*	process did not converge */
182
#define ME_PRECISION	8
183
/*	does not have "full" precision */
184
#define ME_UNDERFLOW	16
185
/*	and underflow occured (important for IEEE)*/
571 ihaka 186
 
187
 
188
#ifdef IEEE_754
7108 ripley 189
#ifdef OLD
3609 pd 190
# ifdef HAVE_IEEE754_H
191
#  include <ieee754.h> /* newer Linuxen */
192
# else
193
#  ifdef HAVE_IEEEFP_H
194
#   include <ieeefp.h> /* others [Solaris 2.5.x], .. */
195
#  endif
196
# endif
7108 ripley 197
#endif
995 ihaka 198
 
571 ihaka 199
extern double m_zero;
200
extern double m_one;
3548 ihaka 201
/* extern double m_tiny; */
571 ihaka 202
#define ML_ERROR(x)	/* nothing */
203
#define ML_POSINF	(m_one / m_zero)
204
#define ML_NEGINF	((-m_one) / m_zero)
205
#define ML_NAN		(m_zero / m_zero)
3548 ihaka 206
#define ML_UNDERFLOW	(DBL_MIN * DBL_MIN)
777 maechler 207
#define ML_VALID(x)	(!isnan(x))
3609 pd 208
 
209
#else/*--- NO IEEE: No +/-Inf, NAN,... ---*/
7016 ripley 210
void ml_error(int n);
571 ihaka 211
#define ML_ERROR(x)	ml_error(x)
212
#define ML_POSINF	DBL_MAX
213
#define ML_NEGINF	(-DBL_MAX)
214
#define ML_NAN		(-DBL_MAX)
215
#define ML_UNDERFLOW	0
777 maechler 216
#define ML_VALID(x)	(errno == 0)
2 r 217
#endif
218
 
571 ihaka 219
	/* Splus Compatibility */
2 r 220
 
571 ihaka 221
#define snorm	norm_rand
222
#define sunif	unif_rand
223
#define sexp	exp_rand
2 r 224
 
1394 ihaka 225
	/* Undo SGI Madness */
226
 
227
#ifdef ftrunc
5408 hornik 228
# undef ftrunc
1394 ihaka 229
#endif
230
#ifdef qexp
5408 hornik 231
# undef qexp
1394 ihaka 232
#endif
5408 hornik 233
#ifdef qgamma
234
# undef qgamma
235
#endif
1394 ihaka 236
 
571 ihaka 237
	/* Name Hiding to Avoid Clashes with Fortran */
2 r 238
 
571 ihaka 239
#ifdef HIDE_NAMES
3076 pd 240
# define d1mach	c_d1mach
241
# define i1mach	c_i1mach
571 ihaka 242
#endif
2 r 243
 
571 ihaka 244
#define	rround	fround
245
#define	prec	fprec
3244 ihaka 246
#undef trunc
571 ihaka 247
#define	trunc	ftrunc
2 r 248
 
4835 maechler 249
 
250
	/* Utilities for `dpq' handling (density/probability/quantile) */
251
 
252
#define R_D__0 (give_log ? ML_NEGINF : 0.)
253
#define R_D__1 (give_log ? 0. : 1.)
254
#define R_DT_0 (lower_tail ? R_D__0 : R_D__1)
255
#define R_DT_1 (lower_tail ? R_D__1 : R_D__0)
256
 
257
#define R_D_val(x)   (give_log	 ? log(x) : x)	      /*  x  */
258
#define R_D_log(x)   (give_log	 ?  x	  : exp(x))   /* log(x) */
259
 
260
#define R_DT_val(x)  R_D_val(lower_tail ? x	 : 1. - x) /*  x  */
261
#define R_DT_Cval(x) R_D_val(lower_tail ? 1. - x : x)	   /*  1 - x */
262
#define R_DT_log(x)  R_D_log(lower_tail ? x	 : 1. - x) /* log(x) */
263
#define R_DT_Clog(x) R_D_log(lower_tail ? 1. - x : x)	   /* log(1 - x) */
264
 
265
#define R_D_give_log(dd)    (((int)dd) >> 1) /* Extract ``give_log'' flag */
266
#define R_D_lower_tail(dd)  (((int)dd) % 2)  /* Extract ``lower_tail'' flag */
267
 
268
	/* R's version of C functions: */
269
 
270
double R_log(double x);
271
double R_pow(double x, double y);
272
 
571 ihaka 273
	/* Machine Characteristics */
2 r 274
 
571 ihaka 275
double	d1mach(int);
276
double	d1mach_(int*);
277
int	i1mach(int);
278
int	i1mach_(int*);
2 r 279
 
571 ihaka 280
	/* General Support Functions */
2 r 281
 
571 ihaka 282
int	imax2(int, int);
283
int	imin2(int, int);
284
double	fmax2(double, double);
285
double	fmin2(double, double);
286
double	fmod(double, double);
287
double	fprec(double, double);
288
double	fround(double, double);
289
double	ftrunc(double);
3076 pd 290
double	sign(double);
571 ihaka 291
double	fsign(double, double);
292
double	fsquare(double);
293
double	fcube(double);
2 r 294
 
571 ihaka 295
	/* Random Number Generators */
2 r 296
 
571 ihaka 297
double	snorm(void);
298
double	sunif(void);
299
double	sexp(void);
2 r 300
 
571 ihaka 301
	/* Chebyshev Series */
2 r 302
 
571 ihaka 303
int	chebyshev_init(double*, int, double);
304
double	chebyshev_eval(double, double *, int);
2 r 305
 
571 ihaka 306
	/* Gamma and Related Functions */
2 r 307
 
571 ihaka 308
double	logrelerr(double);
309
void	gammalims(double*, double*);
310
double	lgammacor(double);
2278 maechler 311
double	gammafn(double);
2629 maechler 312
double	gamma_cody(double);
2278 maechler 313
double	lgammafn(double);
571 ihaka 314
void	dpsifn(double, int, int, int, double*, int*, int*);
315
double	digamma(double);
316
double	trigamma(double);
317
double	tetragamma(double);
318
double	pentagamma(double);
2 r 319
 
571 ihaka 320
double	choose(double, double);
321
double	lchoose(double, double);
322
double	fastchoose(double, double);
323
double	lfastchoose(double, double);
2 r 324
 
2629 maechler 325
	/* Bessel Functions of All Kinds */
326
 
3786 pd 327
double	bessel_i(double, double, double);
328
double	bessel_j(double, double);
329
double	bessel_k(double, double, double);
330
double	bessel_y(double, double);
331
void	I_bessel(double*, double*, long*, long*, double*, long*);
332
void	J_bessel(double*, double*, long*,	 double*, long*);
333
void	K_bessel(double*, double*, long*, long*, double*, long*);
334
void	Y_bessel(double*, double*, long*,	 double*, long*);
2629 maechler 335
 
571 ihaka 336
	/* Beta and Related Functions */
2 r 337
 
571 ihaka 338
double	beta(double, double);
339
double	lbeta(double, double);
2 r 340
 
571 ihaka 341
	/* Normal Distribution */
2 r 342
 
571 ihaka 343
double	dnorm(double, double, double);
344
double	pnorm(double, double, double);
345
double	qnorm(double, double, double);
346
double	rnorm(double, double);
2 r 347
 
571 ihaka 348
	/* Uniform Distribution */
2 r 349
 
571 ihaka 350
double	dunif(double, double, double);
351
double	punif(double, double, double);
352
double	qunif(double, double, double);
353
double	runif(double, double);
2 r 354
 
571 ihaka 355
	/* Gamma Distribution */
2 r 356
 
571 ihaka 357
double	dgamma(double, double, double);
358
double	pgamma(double, double, double);
359
double	qgamma(double, double, double);
360
double	rgamma(double, double);
2 r 361
 
571 ihaka 362
	/* Beta Distribution */
2 r 363
 
571 ihaka 364
double	dbeta(double, double, double);
365
double	pbeta(double, double, double);
366
double	pbeta_raw(double, double, double);
367
double	qbeta(double, double, double);
368
double	rbeta(double, double);
2 r 369
 
571 ihaka 370
	/* Lognormal Distribution */
2 r 371
 
571 ihaka 372
double	dlnorm(double, double, double);
373
double	plnorm(double, double, double);
374
double	qlnorm(double, double, double);
375
double	rlnorm(double, double);
2 r 376
 
571 ihaka 377
	/* Chi-squared Distribution */
2 r 378
 
571 ihaka 379
double	dchisq(double, double);
380
double	pchisq(double, double);
381
double	qchisq(double, double);
382
double	rchisq(double);
2 r 383
 
571 ihaka 384
	/* Non-central Chi-squared Distribution */
2 r 385
 
605 ihaka 386
double	dnchisq(double, double, double);
571 ihaka 387
double	pnchisq(double, double, double);
388
double	qnchisq(double, double, double);
605 ihaka 389
double	rnchisq(double, double);
571 ihaka 390
 
391
	/* F Distibution */
392
 
393
double	df(double, double, double);
394
double	pf(double, double, double);
395
double	qf(double, double, double);
396
double	rf(double, double);
397
 
398
	/* Student t Distibution */
399
 
400
double	dt(double, double);
401
double	pt(double, double);
402
double	qt(double, double);
403
double	rt(double);
404
 
405
	/* Binomial Distribution */
406
 
407
double	dbinom(double, double, double);
408
double	pbinom(double, double, double);
409
double	qbinom(double, double, double);
410
double	rbinom(double, double);
411
 
412
	/* Cauchy Distribution */
413
 
414
double	dcauchy(double, double, double);
415
double	pcauchy(double, double, double);
416
double	qcauchy(double, double, double);
417
double	rcauchy(double, double);
418
 
419
	/* Exponential Distribution */
420
 
421
double	dexp(double, double);
422
double	pexp(double, double);
423
double	qexp(double, double);
424
double	rexp(double);
425
 
426
	/* Geometric Distribution */
427
 
428
double	dgeom(double, double);
429
double	pgeom(double, double);
430
double	qgeom(double, double);
431
double	rgeom(double);
432
 
433
	/* Hypergeometric Distibution */
434
 
435
double	dhyper(double, double, double, double);
436
double	phyper(double, double, double, double);
437
double	qhyper(double, double, double, double);
438
double	rhyper(double, double, double);
439
 
440
	/* Negative Binomial Distribution */
441
 
442
double	dnbinom(double, double, double);
443
double	pnbinom(double, double, double);
444
double	qnbinom(double, double, double);
445
double	rnbinom(double, double);
446
 
447
	/* Poisson Distribution */
448
 
449
double	dpois(double, double);
450
double	ppois(double, double);
451
double	qpois(double, double);
452
double	rpois(double);
453
 
454
	/* Weibull Distribution */
455
 
456
double	dweibull(double, double, double);
457
double	pweibull(double, double, double);
458
double	qweibull(double, double, double);
459
double	rweibull(double, double);
460
 
461
	/* Logistic Distribution */
462
 
463
double	dlogis(double, double, double);
464
double	plogis(double, double, double);
465
double	qlogis(double, double, double);
466
double	rlogis(double, double);
467
 
602 ihaka 468
	/* Non-central Beta Distribution */
469
 
470
double	dnbeta(double, double, double, double);
471
double	pnbeta(double, double, double, double);
472
double	qnbeta(double, double, double, double);
473
double	rnbeta(double, double, double);
474
 
475
	/* Non-central F Distribution */
476
 
477
double	dnf(double, double, double, double);
478
double	pnf(double, double, double, double);
479
double	qnf(double, double, double, double);
480
double	rnf(double, double, double);
481
 
482
	/* Non-central Student t Distribution */
483
 
484
double	dnt(double, double, double);
485
double	pnt(double, double, double);
486
double	qnt(double, double, double);
487
double	rnt(double, double);
488
 
624 ihaka 489
	/* Studentized Range Distribution */
490
 
491
double	dtukey(double, double, double, double);
492
double	ptukey(double, double, double, double);
493
double	qtukey(double, double, double, double);
494
double	rtukey(double, double, double);
495
 
2235 hornik 496
/* Wilcoxon Rank Sum Distribution */
1276 hornik 497
 
3865 pd 498
#define WILCOX_MAX 50
1276 hornik 499
double dwilcox(double, double, double);
500
double pwilcox(double, double, double);
501
double qwilcox(double, double, double);
502
double rwilcox(double, double);
503
 
2235 hornik 504
/* Wilcoxon Signed Rank Distribution */
505
 
3865 pd 506
#define SIGNRANK_MAX 50
2235 hornik 507
double dsignrank(double, double);
508
double psignrank(double, double);
509
double qsignrank(double, double);
510
double rsignrank(double);
511
 
2 r 512
#endif