The R Project SVN R

Rev

Details | Last modification | View Log | RSS feed

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