rmathd.exponential source code

1 module rmathd.exponential;
2 
3 public import rmathd.common;
4 
5 /*
6 ** dmd exponential.d common.d && ./exponential
7 */
8 
9 T dexp(T: double)(T x, T scale, int give_log)
10 {
11     /* NaNs propagated correctly */
12     if (isNaN(x) || isNaN(scale))
13     	return x + scale;
14     if (scale <= 0.0)
15     	return T.nan;
16 
17     mixin R_D__0!(give_log);
18     if (x < 0.)
19 	    return R_D__0;
20     return (give_log ? (-x / scale) - log(scale) : exp(-x / scale) / scale);
21 }
22 
23 
24 T pexp(T: double)(T x, T scale, int lower_tail, int log_p)
25 {
26     if (isNaN(x) || isNaN(scale))
27         return x + scale;
28     if (scale < 0)
29         return T.nan;
30 
31     if (x <= 0.)
32         return R_DT_0!T(lower_tail, log_p);
33     /* same as weibull( shape = 1): */
34     x = -(x / scale);
35     return lower_tail ? (log_p ? R_Log1_Exp!T(x) : -expm1(x)) : R_D_exp!T(x, log_p);
36 }
37 
38 
39 T qexp(T: double)(T p, T scale, int lower_tail, int log_p)
40 {
41     if (isNaN(p) || isNaN(scale))
42         return p + scale;
43 
44     if (scale < 0)
45         return T.nan;
46 
47     mixin (R_Q_P01_check!p());
48     if (p == R_DT_0!T(lower_tail, log_p))
49         return 0;
50 
51     return - scale * R_DT_Clog!T(p, lower_tail, log_p);
52 }
53 
54 
55 T rexp(T: double)(T scale)
56 {
57     if (!isFinite(scale) || scale <= 0.0) {
58         if(scale == 0.)
59             return 0.;
60         /* else */
61         return T.nan;
62     }
63     return scale * exp_rand!T(); // --> in ./sexp.c
64 }
65 
66 
67 
68 void test_exp()
69 {
70 	import std.stdio: write, writeln;
71 	writeln("dexp: ", dexp(2., 1., 0));
72     writeln("pexp: ", pexp(2., 1., 1, 0));
73     writeln("qexp: ", qexp(.8, 1., 1, 0));
74     writeln("rexp: ", rexp(1.), ", rexp: ", rexp(1.),
75             ", rexp: ", rexp(1.), ", rexp: ", rexp(1.));
76 }