1#include "test/jemalloc_test.h"
2
3#define	MAX_REL_ERR 1.0e-9
4#define	MAX_ABS_ERR 1.0e-9
5
6#include <float.h>
7
8#ifndef INFINITY
9#define	INFINITY (DBL_MAX + DBL_MAX)
10#endif
11
12static bool
13double_eq_rel(double a, double b, double max_rel_err, double max_abs_err)
14{
15	double rel_err;
16
17	if (fabs(a - b) < max_abs_err)
18		return (true);
19	rel_err = (fabs(b) > fabs(a)) ? fabs((a-b)/b) : fabs((a-b)/a);
20	return (rel_err < max_rel_err);
21}
22
23static uint64_t
24factorial(unsigned x)
25{
26	uint64_t ret = 1;
27	unsigned i;
28
29	for (i = 2; i <= x; i++)
30		ret *= (uint64_t)i;
31
32	return (ret);
33}
34
35TEST_BEGIN(test_ln_gamma_factorial)
36{
37	unsigned x;
38
39	/* exp(ln_gamma(x)) == (x-1)! for integer x. */
40	for (x = 1; x <= 21; x++) {
41		assert_true(double_eq_rel(exp(ln_gamma(x)),
42		    (double)factorial(x-1), MAX_REL_ERR, MAX_ABS_ERR),
43		    "Incorrect factorial result for x=%u", x);
44	}
45}
46TEST_END
47
48/* Expected ln_gamma([0.0..100.0] increment=0.25). */
49static const double ln_gamma_misc_expected[] = {
50	INFINITY,
51	1.28802252469807743, 0.57236494292470008, 0.20328095143129538,
52	0.00000000000000000, -0.09827183642181320, -0.12078223763524518,
53	-0.08440112102048555, 0.00000000000000000, 0.12487171489239651,
54	0.28468287047291918, 0.47521466691493719, 0.69314718055994529,
55	0.93580193110872523, 1.20097360234707429, 1.48681557859341718,
56	1.79175946922805496, 2.11445692745037128, 2.45373657084244234,
57	2.80857141857573644, 3.17805383034794575, 3.56137591038669710,
58	3.95781396761871651, 4.36671603662228680, 4.78749174278204581,
59	5.21960398699022932, 5.66256205985714178, 6.11591589143154568,
60	6.57925121201010121, 7.05218545073853953, 7.53436423675873268,
61	8.02545839631598312, 8.52516136106541467, 9.03318691960512332,
62	9.54926725730099690, 10.07315123968123949, 10.60460290274525086,
63	11.14340011995171231, 11.68933342079726856, 12.24220494005076176,
64	12.80182748008146909, 13.36802367147604720, 13.94062521940376342,
65	14.51947222506051816, 15.10441257307551943, 15.69530137706046524,
66	16.29200047656724237, 16.89437797963419285, 17.50230784587389010,
67	18.11566950571089407, 18.73434751193644843, 19.35823122022435427,
68	19.98721449566188468, 20.62119544270163018, 21.26007615624470048,
69	21.90376249182879320, 22.55216385312342098, 23.20519299513386002,
70	23.86276584168908954, 24.52480131594137802, 25.19122118273868338,
71	25.86194990184851861, 26.53691449111561340, 27.21604439872720604,
72	27.89927138384089389, 28.58652940490193828, 29.27775451504081516,
73	29.97288476399884871, 30.67186010608067548, 31.37462231367769050,
74	32.08111489594735843, 32.79128302226991565, 33.50507345013689076,
75	34.22243445715505317, 34.94331577687681545, 35.66766853819134298,
76	36.39544520803305261, 37.12659953718355865, 37.86108650896109395,
77	38.59886229060776230, 39.33988418719949465, 40.08411059791735198,
78	40.83150097453079752, 41.58201578195490100, 42.33561646075348506,
79	43.09226539146988699, 43.85192586067515208, 44.61456202863158893,
80	45.38013889847690052, 46.14862228684032885, 46.91997879580877395,
81	47.69417578616628361, 48.47118135183522014, 49.25096429545256882,
82	50.03349410501914463, 50.81874093156324790, 51.60667556776436982,
83	52.39726942748592364, 53.19049452616926743, 53.98632346204390586,
84	54.78472939811231157, 55.58568604486942633, 56.38916764371992940,
85	57.19514895105859864, 58.00360522298051080, 58.81451220059079787,
86	59.62784609588432261, 60.44358357816834371, 61.26170176100199427,
87	62.08217818962842927, 62.90499082887649962, 63.73011805151035958,
88	64.55753862700632340, 65.38723171073768015, 66.21917683354901385,
89	67.05335389170279825, 67.88974313718154008, 68.72832516833013017,
90	69.56908092082363737, 70.41199165894616385, 71.25703896716800045,
91	72.10420474200799390, 72.95347118416940191, 73.80482079093779646,
92	74.65823634883015814, 75.51370092648485866, 76.37119786778275454,
93	77.23071078519033961, 78.09222355331530707, 78.95572030266725960,
94	79.82118541361435859, 80.68860351052903468, 81.55795945611502873,
95	82.42923834590904164, 83.30242550295004378, 84.17750647261028973,
96	85.05446701758152983, 85.93329311301090456, 86.81397094178107920,
97	87.69648688992882057, 88.58082754219766741, 89.46697967771913795,
98	90.35493026581838194, 91.24466646193963015, 92.13617560368709292,
99	93.02944520697742803, 93.92446296229978486, 94.82121673107967297,
100	95.71969454214321615, 96.61988458827809723, 97.52177522288820910,
101	98.42535495673848800, 99.33061245478741341, 100.23753653310367895,
102	101.14611615586458981, 102.05634043243354370, 102.96819861451382394,
103	103.88168009337621811, 104.79677439715833032, 105.71347118823287303,
104	106.63176026064346047, 107.55163153760463501, 108.47307506906540198,
105	109.39608102933323153, 110.32063971475740516, 111.24674154146920557,
106	112.17437704317786995, 113.10353686902013237, 114.03421178146170689,
107	114.96639265424990128, 115.90007047041454769, 116.83523632031698014,
108	117.77188139974506953, 118.70999700805310795, 119.64957454634490830,
109	120.59060551569974962, 121.53308151543865279, 122.47699424143097247,
110	123.42233548443955726, 124.36909712850338394, 125.31727114935689826,
111	126.26684961288492559, 127.21782467361175861, 128.17018857322420899,
112	129.12393363912724453, 130.07905228303084755, 131.03553699956862033,
113	131.99338036494577864, 132.95257503561629164, 133.91311374698926784,
114	134.87498931216194364, 135.83819462068046846, 136.80272263732638294,
115	137.76856640092901785, 138.73571902320256299, 139.70417368760718091,
116	140.67392364823425055, 141.64496222871400732, 142.61728282114600574,
117	143.59087888505104047, 144.56574394634486680, 145.54187159633210058,
118	146.51925549072063859, 147.49788934865566148, 148.47776695177302031,
119	149.45888214327129617, 150.44122882700193600, 151.42480096657754984,
120	152.40959258449737490, 153.39559776128982094, 154.38281063467164245,
121	155.37122539872302696, 156.36083630307879844, 157.35163765213474107,
122	158.34362380426921391, 159.33678917107920370, 160.33112821663092973,
123	161.32663545672428995, 162.32330545817117695, 163.32113283808695314,
124	164.32011226319519892, 165.32023844914485267, 166.32150615984036790,
125	167.32391020678358018, 168.32744544842768164, 169.33210678954270634,
126	170.33788918059275375, 171.34478761712384198, 172.35279713916281707,
127	173.36191283062726143, 174.37212981874515094, 175.38344327348534080,
128	176.39584840699734514, 177.40934047306160437, 178.42391476654847793,
129	179.43956662288721304, 180.45629141754378111, 181.47408456550741107,
130	182.49294152078630304, 183.51285777591152737, 184.53382886144947861,
131	185.55585034552262869, 186.57891783333786861, 187.60302696672312095,
132	188.62817342367162610, 189.65435291789341932, 190.68156119837468054,
133	191.70979404894376330, 192.73904728784492590, 193.76931676731820176,
134	194.80059837318714244, 195.83288802445184729, 196.86618167288995096,
135	197.90047530266301123, 198.93576492992946214, 199.97204660246373464,
136	201.00931639928148797, 202.04757043027063901, 203.08680483582807597,
137	204.12701578650228385, 205.16819948264117102, 206.21035215404597807,
138	207.25347005962987623, 208.29754948708190909, 209.34258675253678916,
139	210.38857820024875878, 211.43552020227099320, 212.48340915813977858,
140	213.53224149456323744, 214.58201366511514152, 215.63272214993284592,
141	216.68436345542014010, 217.73693411395422004, 218.79043068359703739,
142	219.84484974781133815, 220.90018791517996988, 221.95644181913033322,
143	223.01360811766215875, 224.07168349307951871, 225.13066465172661879,
144	226.19054832372759734, 227.25133126272962159, 228.31301024565024704,
145	229.37558207242807384, 230.43904356577689896, 231.50339157094342113,
146	232.56862295546847008, 233.63473460895144740, 234.70172344281823484,
147	235.76958639009222907, 236.83832040516844586, 237.90792246359117712,
148	238.97838956183431947, 240.04971871708477238, 241.12190696702904802,
149	242.19495136964280846, 243.26884900298270509, 244.34359696498191283,
150	245.41919237324782443, 246.49563236486270057, 247.57291409618682110,
151	248.65103474266476269, 249.72999149863338175, 250.80978157713354904,
152	251.89040220972316320, 252.97185064629374551, 254.05412415488834199,
153	255.13722002152300661, 256.22113555000953511, 257.30586806178126835,
154	258.39141489572085675, 259.47777340799029844, 260.56494097186322279,
155	261.65291497755913497, 262.74169283208021852, 263.83127195904967266,
156	264.92164979855277807, 266.01282380697938379, 267.10479145686849733,
157	268.19755023675537586, 269.29109765101975427, 270.38543121973674488,
158	271.48054847852881721, 272.57644697842033565, 273.67312428569374561,
159	274.77057798174683967, 275.86880566295326389, 276.96780494052313770,
160	278.06757344036617496, 279.16810880295668085, 280.26940868320008349,
161	281.37147075030043197, 282.47429268763045229, 283.57787219260217171,
162	284.68220697654078322, 285.78729476455760050, 286.89313329542699194,
163	287.99972032146268930, 289.10705360839756395, 290.21513093526289140,
164	291.32395009427028754, 292.43350889069523646, 293.54380514276073200,
165	294.65483668152336350, 295.76660135076059532, 296.87909700685889902,
166	297.99232151870342022, 299.10627276756946458, 300.22094864701409733,
167	301.33634706277030091, 302.45246593264130297, 303.56930318639643929,
168	304.68685676566872189, 305.80512462385280514, 306.92410472600477078,
169	308.04379504874236773, 309.16419358014690033, 310.28529831966631036,
170	311.40710727801865687, 312.52961847709792664, 313.65282994987899201,
171	314.77673974032603610, 315.90134590329950015, 317.02664650446632777,
172	318.15263962020929966, 319.27932333753892635, 320.40669575400545455,
173	321.53475497761127144, 322.66349912672620803, 323.79292633000159185,
174	324.92303472628691452, 326.05382246454587403, 327.18528770377525916,
175	328.31742861292224234, 329.45024337080525356, 330.58373016603343331,
176	331.71788719692847280, 332.85271267144611329, 333.98820480709991898,
177	335.12436183088397001, 336.26118197919845443, 337.39866349777429377,
178	338.53680464159958774, 339.67560367484657036, 340.81505887079896411,
179	341.95516851178109619, 343.09593088908627578, 344.23734430290727460,
180	345.37940706226686416, 346.52211748494903532, 347.66547389743118401,
181	348.80947463481720661, 349.95411804077025408, 351.09940246744753267,
182	352.24532627543504759, 353.39188783368263103, 354.53908551944078908,
183	355.68691771819692349, 356.83538282361303118, 357.98447923746385868,
184	359.13420536957539753
185};
186
187TEST_BEGIN(test_ln_gamma_misc)
188{
189	unsigned i;
190
191	for (i = 1; i < sizeof(ln_gamma_misc_expected)/sizeof(double); i++) {
192		double x = (double)i * 0.25;
193		assert_true(double_eq_rel(ln_gamma(x),
194		    ln_gamma_misc_expected[i], MAX_REL_ERR, MAX_ABS_ERR),
195		    "Incorrect ln_gamma result for i=%u", i);
196	}
197}
198TEST_END
199
200/* Expected pt_norm([0.01..0.99] increment=0.01). */
201static const double pt_norm_expected[] = {
202	-INFINITY,
203	-2.32634787404084076, -2.05374891063182252, -1.88079360815125085,
204	-1.75068607125216946, -1.64485362695147264, -1.55477359459685305,
205	-1.47579102817917063, -1.40507156030963221, -1.34075503369021654,
206	-1.28155156554460081, -1.22652812003661049, -1.17498679206608991,
207	-1.12639112903880045, -1.08031934081495606, -1.03643338949378938,
208	-0.99445788320975281, -0.95416525314619416, -0.91536508784281390,
209	-0.87789629505122846, -0.84162123357291418, -0.80642124701824025,
210	-0.77219321418868492, -0.73884684918521371, -0.70630256284008752,
211	-0.67448975019608171, -0.64334540539291685, -0.61281299101662701,
212	-0.58284150727121620, -0.55338471955567281, -0.52440051270804067,
213	-0.49585034734745320, -0.46769879911450812, -0.43991316567323380,
214	-0.41246312944140462, -0.38532046640756751, -0.35845879325119373,
215	-0.33185334643681652, -0.30548078809939738, -0.27931903444745404,
216	-0.25334710313579978, -0.22754497664114931, -0.20189347914185077,
217	-0.17637416478086135, -0.15096921549677725, -0.12566134685507399,
218	-0.10043372051146975, -0.07526986209982976, -0.05015358346473352,
219	-0.02506890825871106, 0.00000000000000000, 0.02506890825871106,
220	0.05015358346473366, 0.07526986209982990, 0.10043372051146990,
221	0.12566134685507413, 0.15096921549677739, 0.17637416478086146,
222	0.20189347914185105, 0.22754497664114931, 0.25334710313579978,
223	0.27931903444745404, 0.30548078809939738, 0.33185334643681652,
224	0.35845879325119373, 0.38532046640756762, 0.41246312944140484,
225	0.43991316567323391, 0.46769879911450835, 0.49585034734745348,
226	0.52440051270804111, 0.55338471955567303, 0.58284150727121620,
227	0.61281299101662701, 0.64334540539291685, 0.67448975019608171,
228	0.70630256284008752, 0.73884684918521371, 0.77219321418868492,
229	0.80642124701824036, 0.84162123357291441, 0.87789629505122879,
230	0.91536508784281423, 0.95416525314619460, 0.99445788320975348,
231	1.03643338949378938, 1.08031934081495606, 1.12639112903880045,
232	1.17498679206608991, 1.22652812003661049, 1.28155156554460081,
233	1.34075503369021654, 1.40507156030963265, 1.47579102817917085,
234	1.55477359459685394, 1.64485362695147308, 1.75068607125217102,
235	1.88079360815125041, 2.05374891063182208, 2.32634787404084076
236};
237
238TEST_BEGIN(test_pt_norm)
239{
240	unsigned i;
241
242	for (i = 1; i < sizeof(pt_norm_expected)/sizeof(double); i++) {
243		double p = (double)i * 0.01;
244		assert_true(double_eq_rel(pt_norm(p), pt_norm_expected[i],
245		    MAX_REL_ERR, MAX_ABS_ERR),
246		    "Incorrect pt_norm result for i=%u", i);
247	}
248}
249TEST_END
250
251/*
252 * Expected pt_chi2(p=[0.01..0.99] increment=0.07,
253 *                  df={0.1, 1.1, 10.1, 100.1, 1000.1}).
254 */
255static const double pt_chi2_df[] = {0.1, 1.1, 10.1, 100.1, 1000.1};
256static const double pt_chi2_expected[] = {
257	1.168926411457320e-40, 1.347680397072034e-22, 3.886980416666260e-17,
258	8.245951724356564e-14, 2.068936347497604e-11, 1.562561743309233e-09,
259	5.459543043426564e-08, 1.114775688149252e-06, 1.532101202364371e-05,
260	1.553884683726585e-04, 1.239396954915939e-03, 8.153872320255721e-03,
261	4.631183739647523e-02, 2.473187311701327e-01, 2.175254800183617e+00,
262
263	0.0003729887888876379, 0.0164409238228929513, 0.0521523015190650113,
264	0.1064701372271216612, 0.1800913735793082115, 0.2748704281195626931,
265	0.3939246282787986497, 0.5420727552260817816, 0.7267265822221973259,
266	0.9596554296000253670, 1.2607440376386165326, 1.6671185084541604304,
267	2.2604828984738705167, 3.2868613342148607082, 6.9298574921692139839,
268
269	2.606673548632508, 4.602913725294877, 5.646152813924212,
270	6.488971315540869, 7.249823275816285, 7.977314231410841,
271	8.700354939944047, 9.441728024225892, 10.224338321374127,
272	11.076435368801061, 12.039320937038386, 13.183878752697167,
273	14.657791935084575, 16.885728216339373, 23.361991680031817,
274
275	70.14844087392152, 80.92379498849355, 85.53325420085891,
276	88.94433120715347, 91.83732712857017, 94.46719943606301,
277	96.96896479994635, 99.43412843510363, 101.94074719829733,
278	104.57228644307247, 107.43900093448734, 110.71844673417287,
279	114.76616819871325, 120.57422505959563, 135.92318818757556,
280
281	899.0072447849649, 937.9271278858220, 953.8117189560207,
282	965.3079371501154, 974.8974061207954, 983.4936235182347,
283	991.5691170518946, 999.4334123954690, 1007.3391826856553,
284	1015.5445154999951, 1024.3777075619569, 1034.3538789836223,
285	1046.4872561869577, 1063.5717461999654, 1107.0741966053859
286};
287
288TEST_BEGIN(test_pt_chi2)
289{
290	unsigned i, j;
291	unsigned e = 0;
292
293	for (i = 0; i < sizeof(pt_chi2_df)/sizeof(double); i++) {
294		double df = pt_chi2_df[i];
295		double ln_gamma_df = ln_gamma(df * 0.5);
296		for (j = 1; j < 100; j += 7) {
297			double p = (double)j * 0.01;
298			assert_true(double_eq_rel(pt_chi2(p, df, ln_gamma_df),
299			    pt_chi2_expected[e], MAX_REL_ERR, MAX_ABS_ERR),
300			    "Incorrect pt_chi2 result for i=%u, j=%u", i, j);
301			e++;
302		}
303	}
304}
305TEST_END
306
307/*
308 * Expected pt_gamma(p=[0.1..0.99] increment=0.07,
309 *                   shape=[0.5..3.0] increment=0.5).
310 */
311static const double pt_gamma_shape[] = {0.5, 1.0, 1.5, 2.0, 2.5, 3.0};
312static const double pt_gamma_expected[] = {
313	7.854392895485103e-05, 5.043466107888016e-03, 1.788288957794883e-02,
314	3.900956150232906e-02, 6.913847560638034e-02, 1.093710833465766e-01,
315	1.613412523825817e-01, 2.274682115597864e-01, 3.114117323127083e-01,
316	4.189466220207417e-01, 5.598106789059246e-01, 7.521856146202706e-01,
317	1.036125427911119e+00, 1.532450860038180e+00, 3.317448300510606e+00,
318
319	0.01005033585350144, 0.08338160893905107, 0.16251892949777497,
320	0.24846135929849966, 0.34249030894677596, 0.44628710262841947,
321	0.56211891815354142, 0.69314718055994529, 0.84397007029452920,
322	1.02165124753198167, 1.23787435600161766, 1.51412773262977574,
323	1.89711998488588196, 2.52572864430825783, 4.60517018598809091,
324
325	0.05741590094955853, 0.24747378084860744, 0.39888572212236084,
326	0.54394139997444901, 0.69048812513915159, 0.84311389861296104,
327	1.00580622221479898, 1.18298694218766931, 1.38038096305861213,
328	1.60627736383027453, 1.87396970522337947, 2.20749220408081070,
329	2.65852391865854942, 3.37934630984842244, 5.67243336507218476,
330
331	0.1485547402532659, 0.4657458011640391, 0.6832386130709406,
332	0.8794297834672100, 1.0700752852474524, 1.2629614217350744,
333	1.4638400448580779, 1.6783469900166610, 1.9132338090606940,
334	2.1778589228618777, 2.4868823970010991, 2.8664695666264195,
335	3.3724415436062114, 4.1682658512758071, 6.6383520679938108,
336
337	0.2771490383641385, 0.7195001279643727, 0.9969081732265243,
338	1.2383497880608061, 1.4675206597269927, 1.6953064251816552,
339	1.9291243435606809, 2.1757300955477641, 2.4428032131216391,
340	2.7406534569230616, 3.0851445039665513, 3.5043101122033367,
341	4.0575997065264637, 4.9182956424675286, 7.5431362346944937,
342
343	0.4360451650782932, 0.9983600902486267, 1.3306365880734528,
344	1.6129750834753802, 1.8767241606994294, 2.1357032436097660,
345	2.3988853336865565, 2.6740603137235603, 2.9697561737517959,
346	3.2971457713883265, 3.6731795898504660, 4.1275751617770631,
347	4.7230515633946677, 5.6417477865306020, 8.4059469148854635
348};
349
350TEST_BEGIN(test_pt_gamma_shape)
351{
352	unsigned i, j;
353	unsigned e = 0;
354
355	for (i = 0; i < sizeof(pt_gamma_shape)/sizeof(double); i++) {
356		double shape = pt_gamma_shape[i];
357		double ln_gamma_shape = ln_gamma(shape);
358		for (j = 1; j < 100; j += 7) {
359			double p = (double)j * 0.01;
360			assert_true(double_eq_rel(pt_gamma(p, shape, 1.0,
361			    ln_gamma_shape), pt_gamma_expected[e], MAX_REL_ERR,
362			    MAX_ABS_ERR),
363			    "Incorrect pt_gamma result for i=%u, j=%u", i, j);
364			e++;
365		}
366	}
367}
368TEST_END
369
370TEST_BEGIN(test_pt_gamma_scale)
371{
372	double shape = 1.0;
373	double ln_gamma_shape = ln_gamma(shape);
374
375	assert_true(double_eq_rel(
376	    pt_gamma(0.5, shape, 1.0, ln_gamma_shape) * 10.0,
377	    pt_gamma(0.5, shape, 10.0, ln_gamma_shape), MAX_REL_ERR,
378	    MAX_ABS_ERR),
379	    "Scale should be trivially equivalent to external multiplication");
380}
381TEST_END
382
383int
384main(void)
385{
386
387	return (test(
388	    test_ln_gamma_factorial,
389	    test_ln_gamma_misc,
390	    test_pt_norm,
391	    test_pt_chi2,
392	    test_pt_gamma_shape,
393	    test_pt_gamma_scale));
394}
395