1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
|
Description: From upstream, avoid nested C functions
--- a/CUtils/c_sources/rhyper.c
+++ b/CUtils/c_sources/rhyper.c
@@ -52,55 +52,61 @@
* If (i > 7), use Stirling's approximation, otherwise use table lookup.
*/
-static double afc(int i)
-{
- static int computed=10;
- static double al[1756] =
- {
- 0.0,
- 0,/*ln(0!)*/
- 0,/*ln(1!)*/
- 0.693147180559945309,/*ln(2!)*/
- 1.791759469228055,/*ln(3!)*/
- 3.17805383034794562,/*ln(4!)*/
- 4.78749174278204599,/*ln(5!)*/
- 6.579251212010101,/*ln(6!)*/
- 8.5251613610654143,/*ln(7!)*/
- 10.6046029027452502,/*ln(8!)*/
- 12.8018274800814696,/*ln(9!)*/
- 15.1044125730755153,/*ln(10!)*/
- };
- double compute(int n) {
- static long double cur=3628800;
- static int i=11;
- static volatile int mutex=0;
-
- while (__sync_lock_test_and_set(&mutex, 1)) {
- /* Internal loop with only read to avoid cache line ping-pong
- on multi-processors */
- while(mutex) {
- /* spinlock */
- }
- }
+struct afc_data {
+ int computed;
+ double al[1756];
+};
+
+double compute(int n, struct afc_data * __restrict__ data) {
+ static long double cur=3628800;
+ static int i=11;
+ static volatile int mutex=0;
- for(; i<=n; i++) {
- cur*=i;
- al[i+1]=logl(cur);
+ while (__sync_lock_test_and_set(&mutex, 1)) {
+ /* Internal loop with only read to avoid cache line ping-pong
+ on multi-processors */
+ while(mutex) {
+ /* spinlock */
}
- computed=n;
- __sync_lock_release(&mutex);
- return al[i];
- };
+ }
+ for(; i<=n; i++) {
+ cur*=i;
+ data->al[i+1]=logl(cur);
+ }
+ data->computed=n;
+ __sync_lock_release(&mutex);
+ return data->al[i];
+};
+
+static double afc(int i)
+{
double di, value;
+ static struct afc_data data = {
+ .computed = 10,
+ .al = {
+ 0.0,
+ 0,/*ln(0!)*/
+ 0,/*ln(1!)*/
+ 0.693147180559945309,/*ln(2!)*/
+ 1.791759469228055,/*ln(3!)*/
+ 3.17805383034794562,/*ln(4!)*/
+ 4.78749174278204599,/*ln(5!)*/
+ 6.579251212010101,/*ln(6!)*/
+ 8.5251613610654143,/*ln(7!)*/
+ 10.6046029027452502,/*ln(8!)*/
+ 12.8018274800814696,/*ln(9!)*/
+ 15.1044125730755153,/*ln(10!)*/
+ }
+ };
if (i < 0) {
fprintf(stderr, "rhyper.c: afc(i), i=%d < 0 -- SHOULD NOT HAPPEN!\n", i);
exit(1);
- } else if (i <= computed) {
- value = al[i + 1];
+ } else if (i <= data.computed) {
+ value = data.al[i + 1];
} else if (i <= 1754) {
- value = compute(i);
+ value = compute(i, &data);
} else {
di = i;
value = (di + 0.5) * log(di) - di + 0.08333333333333 / di
|
