A library for working with phylogenetic data.
v0.25.0
common.cpp
Go to the documentation of this file.
1 /*
2  Genesis - A toolkit for working with phylogenetic data.
3  Copyright (C) 2014-2020 Lucas Czech
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 3 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, see <http://www.gnu.org/licenses/>.
17 
18  Contact:
19  Lucas Czech <lucas.czech@h-its.org>
20  Exelixis Lab, Heidelberg Institute for Theoretical Studies
21  Schloss-Wolfsbrunnenweg 35, D-69118 Heidelberg, Germany
22 */
23 
32 
33 #include <array>
34 #include <cassert>
35 #include <iostream>
36 #include <limits>
37 #include <sstream>
38 #include <stdexcept>
39 #include <string>
40 
41 namespace genesis {
42 namespace utils {
43 
44 // =============================================================================
45 // Binomial Coefficient
46 // =============================================================================
47 
53 static const std::array<double, 1024> log_factorial_lookup_ = {{
54  0.0000000000000000, 0.0000000000000000, 0.6931471805599452, 1.7917594692280550, 3.1780538303479453,
55  4.7874917427820458, 6.5792512120101012, 8.5251613610654147, 10.604602902745251, 12.801827480081471,
56  15.104412573075518, 17.502307845873887, 19.987214495661888, 22.552163853123425, 25.191221182738683,
57  27.899271383840894, 30.671860106080675, 33.505073450136891, 36.395445208033053, 39.339884187199495,
58  42.335616460753485, 45.380138898476908, 48.471181351835227, 51.606675567764377, 54.784729398112326,
59  58.003605222980525, 61.261701761002008, 64.557538627006338, 67.889743137181540, 71.257038967168015,
60  74.658236348830172, 78.092223553315321, 81.557959456115043, 85.054467017581530, 88.580827542197696,
61  92.136175603687107, 95.719694542143216, 99.330612454787442, 102.96819861451382, 106.63176026064347,
62  110.32063971475741, 114.03421178146171, 117.77188139974507, 121.53308151543864, 125.31727114935690,
63  129.12393363912722, 132.95257503561632, 136.80272263732638, 140.67392364823428, 144.56574394634490,
64  148.47776695177305, 152.40959258449737, 156.36083630307880, 160.33112821663093, 164.32011226319520,
65  168.32744544842768, 172.35279713916282, 176.39584840699737, 180.45629141754378, 184.53382886144951,
66  188.62817342367160, 192.73904728784490, 196.86618167288998, 201.00931639928152, 205.16819948264120,
67  209.34258675253685, 213.53224149456327, 217.73693411395422, 221.95644181913033, 226.19054832372760,
68  230.43904356577696, 234.70172344281826, 238.97838956183432, 243.26884900298271, 247.57291409618688,
69  251.89040220972319, 256.22113555000954, 260.56494097186322, 264.92164979855283, 269.29109765101987,
70  273.67312428569375, 278.06757344036617, 282.47429268763045, 286.89313329542705, 291.32395009427034,
71  295.76660135076065, 300.22094864701415, 304.68685676566872, 309.16419358014696, 313.65282994987911,
72  318.15263962020936, 322.66349912672621, 327.18528770377526, 331.71788719692853, 336.26118197919851,
73  340.81505887079908, 345.37940706226692, 349.95411804077031, 354.53908551944090, 359.13420536957551,
74  363.73937555556358, 368.35449607240486, 372.97946888568913, 377.61419787391878, 382.25858877306018,
75  386.91254912321773, 391.57598821732978, 396.24881705179166, 400.93094827891588, 405.62229616114502,
76  410.32277652693745, 415.03230672824981, 419.75080559954489, 424.47819341825721, 429.21439186665168,
77  433.95932399501493, 438.71291418612128, 443.47508812091905, 448.24577274538473, 453.02489623849624,
78  457.81238798127828, 462.60817852687501, 467.41219957160826, 472.22438392698069, 477.04466549258575,
79  481.87297922988807, 486.70926113683953, 491.55344822329812, 496.40547848721775, 501.26529089157941,
80  506.13282534203501, 511.00802266523618, 515.89082458782252, 520.78117371604424, 525.67901351599517,
81  530.58428829443358, 535.49694318016964, 540.41692410599774, 545.34417779115495, 550.27865172428562,
82  555.22029414689496, 560.16905403727310, 565.12488109487435, 570.08772572513431, 575.05753902471031,
83  580.03427276713091, 585.01787938883922, 590.00831197561797, 595.00552424938212, 600.00947055532754,
84  605.02010584942377, 610.03738568623874, 615.06126620708505, 620.09170412847755, 625.12865673089118,
85  630.17208184781043, 635.22193785505999, 640.27818366040833, 645.34077869343525, 650.40968289565546,
86  655.48485671088929, 660.56626107587374, 665.65385741110617, 670.74760761191294, 675.84747403973711,
87  680.95341951363764, 686.06540730199413, 691.18340111441091, 696.30736509381416, 701.43726380873727,
88  706.57306224578758, 711.71472580229022, 716.86222027910367, 722.01551187360144, 727.17456717281596,
89  732.33935314673943, 737.50983714177755, 742.68598687435133, 747.86777042464337, 753.05515623048416,
90  758.24811308137441, 763.44661011264020, 768.65061679971700, 773.86010295255846, 779.07503871016741,
91  784.29539453524569, 789.52114120895885, 794.75224982581346, 799.98869178864345, 805.23043880370312,
92  810.47746287586358, 815.72973630391016, 820.98723167593789, 826.24992186484280, 831.51778002390608,
93  836.79077958246978, 842.06889424170026, 847.35209797043831, 852.64036500113286, 857.93366982585735,
94  863.23198719240543, 868.53529210046452, 873.84355979786574, 879.15676577690749, 884.47488577075171,
95  889.79789574989013, 895.12577191867967, 900.45849071194505, 905.79602879164634, 911.13836304361109,
96  916.48547057432859, 921.83732870780466, 927.19391498247671, 932.55520714818613, 937.92118316320796,
97  943.29182119133566, 948.66709959901982, 954.04699695256033, 959.43149201534948, 964.82056374516594,
98  970.21419129151832, 975.61235399303609, 981.01503137490840, 986.42220314636847, 991.83384919822356,
99  997.24994960042795, 1002.6704845997002, 1008.0954346171816, 1013.5247802461361, 1018.9585022496904,
100  1024.3965815586137, 1029.8389992691355, 1035.2857366408018, 1040.7367750943674, 1046.1920962097252,
101  1051.6516817238694, 1057.1155135288950, 1062.5835736700301, 1068.0558443437017, 1073.5323078956333,
102  1079.0129468189753, 1084.4977437524658, 1089.9866814786226, 1095.4797429219632, 1100.9769111472565,
103  1106.4781693578011, 1111.9835008937334, 1117.4928892303615, 1123.0063179765264, 1128.5237708729910,
104  1134.0452317908532, 1139.5706847299850, 1145.1001138174965, 1150.6335033062242, 1156.1708375732428,
105  1161.7121011184013, 1167.2572785628809, 1172.8063546477761, 1178.3593142326977, 1183.9161422943971,
106  1189.4768239254126, 1195.0413443327354, 1200.6096888364966, 1206.1818428686745, 1211.7577919718208,
107  1217.3375217978071, 1222.9210181065887, 1228.5082667649890, 1234.0992537454999, 1239.6939651251018,
108  1245.2923870841003, 1250.8945059049800, 1256.5003079712760, 1262.1097797664611, 1267.7229078728492,
109  1273.3396789705157, 1278.9600798362328, 1284.5840973424201, 1290.2117184561107, 1295.8429302379320,
110  1301.4777198411014, 1307.1160745104351, 1312.7579815813731, 1318.4034284790164, 1324.0524027171775,
111  1329.7048918974463, 1335.3608837082661, 1341.0203659240258, 1346.6833264041618, 1352.3497530922741,
112  1358.0196340152547, 1363.6929572824263, 1369.3697110846945, 1375.0498836937115, 1380.7334634610502,
113  1386.4204388173900, 1392.1107982717142, 1397.8045304105169, 1403.5016238970222, 1409.2020674704129,
114  1414.9058499450691, 1420.6129602098181, 1426.3233872271931, 1432.0371200327024, 1437.7541477341088,
115  1443.4744595107161, 1449.1980446126686, 1454.9248923602559, 1460.6549921432295, 1466.3883334201273,
116  1472.1249057176065, 1477.8646986297856, 1483.6077018175952, 1489.3539050081354, 1495.1032979940437,
117  1500.8558706328693, 1506.6116128464562, 1512.3705146203336, 1518.1325660031137, 1523.8977571058986,
118  1529.6660781016924, 1535.4375192248224, 1541.2120707703668, 1546.9897230935894, 1552.7704666093819,
119  1558.5542917917116, 1564.3411891730784, 1570.1311493439757, 1575.9241629523599, 1581.7202207031253,
120  1587.5193133575858, 1593.3214317329630, 1599.1265667018795, 1604.9347091918598, 1610.7458501848366,
121  1616.5599807166616, 1622.3770918766247, 1628.1971748069770, 1634.0202207024602, 1639.8462208098406,
122  1645.6751664274509, 1651.5070489047343, 1657.3418596417969, 1663.1795900889629, 1669.0202317463363,
123  1674.8637761633677, 1680.7102149384255, 1686.5595397183722, 1692.4117421981466, 1698.2668141203490,
124  1704.1247472748323, 1709.9855334982981, 1715.8491646738962, 1721.7156327308296, 1727.5849296439633,
125  1733.4570474334387, 1739.3319781642906, 1745.2097139460702, 1751.0902469324710, 1756.9735693209593,
126  1762.8596733524096, 1768.7485513107424, 1774.6401955225681, 1780.5345983568329, 1786.4317522244696,
127  1792.3316495780521, 1798.2342829114534, 1804.1396447595080, 1810.0477276976769, 1815.9585243417175,
128  1821.8720273473557, 1827.7882294099632, 1833.7071232642363, 1839.6287016838801, 1845.5529574812947,
129  1851.4798835072652, 1857.4094726506551, 1863.3417178381030, 1869.2766120337226, 1875.2141482388049,
130  1881.1543194915253, 1887.0971188666520, 1893.0425394752585, 1898.9905744644391, 1904.9412170170267,
131  1910.8944603513146, 1916.8502977207795, 1922.8087224138094, 1928.7697277534326, 1934.7333070970510,
132  1940.6994538361746, 1946.6681613961600, 1952.6394232359505, 1958.6132328478197, 1964.5895837571177,
133  1970.5684695220189, 1976.5498837332734, 1982.5338200139606, 1988.5202720192451, 1994.5092334361350,
134  2000.5006979832431, 2006.4946594105497, 2012.4911114991687, 2018.4900480611154, 2024.4914629390767,
135  2030.4953500061831, 2036.5017031657849, 2042.5105163512276, 2048.5217835256317, 2054.5354986816747,
136  2060.5516558413728, 2066.5702490558692, 2072.5912724052187, 2078.6147199981797, 2084.6405859720048,
137  2090.6688644922356, 2096.6995497524968, 2102.7326359742956, 2108.7681174068202, 2114.8059883267424,
138  2120.8462430380200, 2126.8888758717026, 2132.9338811857388, 2138.9812533647851, 2145.0309868200170,
139  2151.0830759889413, 2157.1375153352105, 2163.1942993484390, 2169.2534225440209, 2175.3148794629487,
140  2181.3786646716362, 2187.4447727617398, 2193.5131983499841, 2199.5839360779864, 2205.6569806120870,
141  2211.7323266431758, 2217.8099688865250, 2223.8899020816207, 2229.9721209919971, 2236.0566204050724,
142  2242.1433951319846, 2248.2324400074313, 2254.3237498895091, 2260.4173196595543, 2266.5131442219867,
143  2272.6112185041529, 2278.7115374561731, 2284.8140960507867, 2290.9188892832017, 2297.0259121709441,
144  2303.1351597537086, 2309.2466270932114, 2315.3603092730436, 2321.4762013985269, 2327.5942985965680,
145  2333.7145960155190, 2339.8370888250333, 2345.9617722159273, 2352.0886414000415, 2358.2176916101021,
146  2364.3489180995853, 2370.4823161425820, 2376.6178810336637, 2382.7556080877498, 2388.8954926399761,
147  2395.0375300455635, 2401.1817156796892, 2407.3280449373583, 2413.4765132332759, 2419.6271160017222,
148  2425.7798486964261, 2431.9347067904428, 2438.0916857760285, 2444.2507811645205, 2450.4119884862157,
149  2456.5753032902503, 2462.7407211444815, 2468.9082376353699, 2475.0778483678614, 2481.2495489652724,
150  2487.4233350691743, 2493.5992023392801, 2499.7771464533307, 2505.9571631069834, 2512.1392480137001,
151  2518.3233969046378, 2524.5096055285385, 2530.6978696516212, 2536.8881850574744, 2543.0805475469492,
152  2549.2749529380540, 2555.4713970658486, 2561.6698757823410, 2567.8703849563835, 2574.0729204735712,
153  2580.2774782361398, 2586.4840541628646, 2592.6926441889614, 2598.9032442659859, 2605.1158503617376,
154  2611.3304584601597, 2617.5470645612445, 2623.7656646809364, 2629.9862548510359, 2636.2088311191073,
155  2642.4333895483828, 2648.6599262176701, 2654.8884372212615, 2661.1189186688398, 2667.3513666853905,
156  2673.5857774111091, 2679.8221470013127, 2686.0604716263524, 2692.3007474715232, 2698.5429707369785,
157  2704.7871376376424, 2711.0332444031242, 2717.2812872776326, 2723.5312625198922, 2729.7831664030582,
158  2736.0369952146339, 2742.2927452563872, 2748.5504128442699, 2754.8099943083348, 2761.0714859926557,
159  2767.3348842552473, 2773.6001854679848, 2779.8673860165263, 2786.1364823002327, 2792.4074707320910,
160  2798.6803477386370, 2804.9551097598787, 2811.2317532492202, 2817.5102746733860, 2823.7906705123460,
161  2830.0729372592418, 2836.3570714203124, 2842.6430695148215, 2848.9309280749831, 2855.2206436458919,
162  2861.5122127854502, 2867.8056320642968, 2874.1008980657366, 2880.3980073856706, 2886.6969566325265,
163  2892.9977424271897, 2899.3003614029344, 2905.6048102053564, 2911.9110854923047, 2918.2191839338143,
164  2924.5291022120409, 2930.8408370211937, 2937.1543850674707, 2943.4697430689930, 2949.7869077557402,
165  2956.1058758694867, 2962.4266441637374, 2968.7492094036647, 2975.0735683660459, 2981.3997178392010,
166  2987.7276546229300, 2994.0573755284527, 3000.3888773783465, 3006.7221570064862, 3013.0572112579844,
167  3019.3940369891307, 3025.7326310673338, 3032.0729903710617, 3038.4151117897827, 3044.7589922239090,
168  3051.1046285847374, 3057.4520177943932, 3063.8011567857729, 3070.1520425024878, 3076.5046718988074,
169  3082.8590419396046, 3089.2151496003003, 3095.5729918668085, 3101.9325657354807, 3108.2938682130539,
170  3114.6568963165942, 3121.0216470734463, 3127.3881175211777, 3133.7563047075282, 3140.1262056903565,
171  3146.4978175375882, 3152.8711373271653, 3159.2461621469934, 3165.6228890948919, 3172.0013152785436,
172  3178.3814378154434, 3184.7632538328494, 3191.1467604677332, 3197.5319548667308, 3203.9188341860936,
173  3210.3073955916393, 3216.6976362587047, 3223.0895533720973, 3229.4831441260480, 3235.8784057241637,
174  3242.2753353793801, 3248.6739303139152, 3255.0741877592241, 3261.4761049559511, 3267.8796791538857,
175  3274.2849076119164, 3280.6917875979857, 3287.1003163890450, 3293.5104912710112, 3299.9223095387210,
176  3306.3357684958883, 3312.7508654550597, 3319.1675977375721, 3325.5859626735082, 3332.0059576016552,
177  3338.4275798694616, 3344.8508268329952, 3351.2756958569007, 3357.7021843143584, 3364.1302895870431,
178  3370.5600090650823, 3376.9913401470158, 3383.4242802397548, 3389.8588267585424, 3396.2949771269118,
179  3402.7327287766484, 3409.1720791477487, 3415.6130256883816, 3422.0555658548496, 3428.4996971115502,
180  3434.9454169309356, 3441.3927227934769, 3447.8416121876239, 3454.2920826097679, 3460.7441315642050,
181  3467.1977565630978, 3473.6529551264380, 3480.1097247820103, 3486.5680630653551, 3493.0279675197326,
182  3499.4894356960863, 3505.9524651530069, 3512.4170534566970, 3518.8831981809344, 3525.3508969070390,
183  3531.8201472238347, 3538.2909467276172, 3544.7632930221180, 3551.2371837184701, 3557.7126164351744,
184  3564.1895887980641, 3570.6680984402728, 3577.1481430021995, 3583.6297201314760, 3590.1128274829334,
185  3596.5974627185687, 3603.0836235075130, 3609.5713075259978, 3616.0605124573231, 3622.5512359918257,
186  3629.0434758268461, 3635.5372296666978, 3642.0324952226347, 3648.5292702128204, 3655.0275523622968,
187  3661.5273394029527, 3668.0286290734930, 3674.5314191194088, 3681.0357072929455, 3687.5414913530735,
188  3694.0487690654586, 3700.5575382024304, 3707.0677965429536, 3713.5795418725984, 3720.0927719835108,
189  3726.6074846743832, 3733.1236777504259, 3739.6413490233381, 3746.1604963112786, 3752.6811174388372,
190  3759.2032102370072, 3765.7267725431566, 3772.2518022009999, 3778.7782970605708, 3785.3062549781935,
191  3791.8356738164557, 3798.3665514441814, 3804.8988857364038, 3811.4326745743374, 3817.9679158453509,
192  3824.5046074429424, 3831.0427472667102, 3837.5823332223281, 3844.1233632215180, 3850.6658351820247,
193  3857.2097470275894, 3863.7550966879239, 3870.3018820986845, 3876.8501012014467, 3883.3997519436807,
194  3889.9508322787242, 3896.5033401657588, 3903.0572735697847, 3909.6126304615955, 3916.1694088177537,
195  3922.7276066205659, 3929.2872218580592, 3935.8482525239556, 3942.4106966176491, 3948.9745521441814,
196  3955.5398171142169, 3962.1064895440200, 3968.6745674554318, 3975.2440488758461, 3981.8149318381857,
197  3988.3872143808799, 3994.9608945478403, 4001.5359703884401, 4008.1124399574883, 4014.6903013152091,
198  4021.2695525272193, 4027.8501916645041, 4034.4322168033968, 4041.0156260255558, 4047.6004174179416,
199  4054.1865890727963, 4060.7741390876213, 4067.3630655651550, 4073.9533666133516, 4080.5450403453601,
200  4087.1380848795025, 4093.7324983392523, 4100.3282788532133, 4106.9254245551001, 4113.5239335837150,
201  4120.1238040829276, 4126.7250342016569, 4133.3276220938460, 4139.9315659184467, 4146.5368638393948,
202  4153.1435140255926, 4159.7515146508886, 4166.3608638940559, 4172.9715599387737, 4179.5836009736067,
203  4186.1969851919866, 4192.8117107921908, 4199.4277759773240, 4206.0451789552981, 4212.6639179388148,
204  4219.2839911453448, 4225.9053967971086, 4232.5281331210581, 4239.1521983488583, 4245.7775907168661,
205  4252.4043084661153, 4259.0323498422949, 4265.6617130957320, 4272.2923964813745, 4278.9243982587705,
206  4285.5577166920511, 4292.1923500499124, 4298.8282966055995, 4305.4655546368840, 4312.1041224260507,
207  4318.7439982598771, 4325.3851804296173, 4332.0276672309847, 4338.6714569641326, 4345.3165479336385,
208  4351.9629384484861, 4358.6106268220492, 4365.2596113720738, 4371.9098904206612, 4378.5614622942512,
209  4385.2143253236045, 4391.8684778437873, 4398.5239181941552, 4405.1806447183335, 4411.8386557642043,
210  4418.4979496838878, 4425.1585248337278, 4431.8203795742729, 4438.4835122702634, 4445.1479212906142,
211  4451.8136050083967, 4458.4805618008259, 4465.1487900492430, 4471.8182881391012, 4478.4890544599475,
212  4485.1610874054086, 4491.8343853731758, 4498.5089467649905, 4505.1847699866257, 4511.8618534478728,
213  4518.5401955625275, 4525.2197947483719, 4531.9006494271625, 4538.5827580246123, 4545.2661189703786,
214  4551.9507306980468, 4558.6365916451150, 4565.3237002529813, 4572.0120549669282, 4578.7016542361071,
215  4585.3924965135257, 4592.0845802560325, 4598.7779039243023, 4605.4724659828235, 4612.1682648998822,
216  4618.8652991475483, 4625.5635672016633, 4632.2630675418250, 4638.9637986513726, 4645.6657590173754,
217  4652.3689471306161, 4659.0733614855799, 4665.7790005804400, 4672.4858629170431, 4679.1939470008965,
218  4685.9032513411548, 4692.6137744506077, 4699.3255148456637, 4706.0384710463404, 4712.7526415762495,
219  4719.4680249625844, 4726.1846197361056, 4732.9024244311295, 4739.6214375855143, 4746.3416577406497,
220  4753.0630834414405, 4759.7857132362960, 4766.5095456771169, 4773.2345793192835, 4779.9608127216425,
221  4786.6882444464936, 4793.4168730595784, 4800.1466971300679, 4806.8777152305502, 4813.6099259370176,
222  4820.3433278288549, 4827.0779194888282, 4833.8136995030709, 4840.5506664610730, 4847.2888189556688,
223  4854.0281555830261, 4860.7686749426321, 4867.5103756372837, 4874.2532562730757, 4880.9973154593872,
224  4887.7425518088712, 4894.4889639374442, 4901.2365504642739, 4907.9853100117653, 4914.7352412055534,
225  4921.4863426744905, 4928.2386130506320, 4934.9920509692301, 4941.7466550687177, 4948.5024239907016,
226  4955.2593563799492, 4962.0174508843766, 4968.7767061550403, 4975.5371208461238, 4982.2986936149282,
227  4989.0614231218597, 4995.8253080304221, 5002.5903470072026, 5009.3565387218632, 5016.1238818471284,
228  5022.8923750587774, 5029.6620170356300, 5036.4328064595393, 5043.2047420153785, 5049.9778223910344,
229  5056.7520462773919, 5063.5274123683284, 5070.3039193607010, 5077.0815659543359, 5083.8603508520209,
230  5090.6402727594932, 5097.4213303854294, 5104.2035224414358, 5110.9868476420397, 5117.7713047046773,
231  5124.5568923496849, 5131.3436093002902, 5138.1314542825994, 5144.9204260255920, 5151.7105232611057,
232  5158.5017447238315, 5165.2940891513026, 5172.0875552838825, 5178.8821418647585, 5185.6778476399322,
233  5192.4746713582072, 5199.2726117711818, 5206.0716676332404, 5212.8718377015430, 5219.6731207360144,
234  5226.4755154993391, 5233.2790207569478, 5240.0836352770102, 5246.8893578304269, 5253.6961871908188,
235  5260.5041221345191, 5267.3131614405620, 5274.1233038906776, 5280.9345482692788, 5287.7468933634564,
236  5294.5603379629674, 5301.3748808602277, 5308.1905208503022, 5315.0072567308971, 5321.8250873023517,
237  5328.6440113676272, 5335.4640277323015, 5342.2851352045582, 5349.1073325951784, 5355.9306187175343,
238  5362.7549923875777, 5369.5804524238329, 5376.4069976473893, 5383.2346268818919, 5390.0633389535333,
239  5396.8931326910460, 5403.7240069256923, 5410.5559604912578, 5417.3889922240442, 5424.2231009628576,
240  5431.0582855490047, 5437.8945448262821, 5444.7318776409675, 5451.5702828418152, 5458.4097592800445,
241  5465.2503058093334, 5472.0919212858107, 5478.9346045680495, 5485.7783545170560, 5492.6231699962646,
242  5499.4690498715290, 5506.3159930111142, 5513.1639982856905, 5520.0130645683239, 5526.8631907344698,
243  5533.7143756619635, 5540.5666182310151, 5547.4199173242014, 5554.2742718264562, 5561.1296806250666,
244  5567.9861426096613, 5574.8436566722066, 5581.7022217069980, 5588.5618366106519, 5595.4225002821004,
245  5602.2842116225811, 5609.1469695356327, 5616.0107729270858, 5622.8756207050565, 5629.7415117799401,
246  5636.6084450644021, 5643.4764194733725, 5650.3454339240379, 5657.2154873358359, 5664.0865786304466,
247  5670.9587067317852, 5677.8318705659976, 5684.7060690614508, 5691.5813011487271, 5698.4575657606183,
248  5705.3348618321161, 5712.2131883004076, 5719.0925441048685, 5725.9729281870541, 5732.8543394906965,
249  5739.7367769616940, 5746.6202395481068, 5753.5047262001499, 5760.3902358701844, 5767.2767675127152,
250  5774.1643200843801, 5781.0528925439457, 5787.9424838523000, 5794.8330929724470, 5801.7247188694992,
251  5808.6173605106715, 5815.5110168652745, 5822.4056869047081, 5829.3013696024564, 5836.1980639340791,
252  5843.0957688772078, 5849.9944834115377, 5856.8942065188221, 5863.7949371828672, 5870.6966743895237,
253  5877.5994171266821, 5884.5031643842667, 5891.4079151542282, 5898.3136684305400, 5905.2204232091890,
254  5912.1281784881712, 5919.0369332674863, 5925.9466865491313, 5932.8574373370930, 5939.7691846373446,
255  5946.6819274578374, 5953.5956648084975, 5960.5103957012161, 5967.4261191498472, 5974.3428341702011,
256  5981.2605397800362, 5988.1792349990565, 5995.0989188489038, 6002.0195903531521, 6008.9412485373032,
257  6015.8638924287789, 6022.7875210569173, 6029.7121334529656, 6036.6377286500765, 6043.5643056832996,
258  6050.4918635895783, 6057.4204014077432, 6064.3499181785064, 6071.2804129444585
259 }};
260 
261 // constexpr version to generate the above table. Does only work in C++ >14, but we use C++11 here.
262 // We could probably rewrite this to use some template loop unrolling at compile time via
263 // variadic template argument lists. But just not worth the effort for now. Also, we might run into
264 // issues with the maximum template instantiation depth...
265 
266 // template<size_t n>
267 // struct LogFactorial{
268 // std::array<double, n> data;
269 //
270 // constexpr LogFactorial()
271 // : data( 0.0 )
272 // {
273 // // That loop needs C++14 or above...
274 // for( size_t i = 1; i < n; ++i ) {
275 // data[i] = data[i-1] + constexpr_log(i);
276 // }
277 // }
278 //
279 // constexpr double constexpr_log( double x ) const
280 // {
281 // // We need to wrap this, as std::log is not marked constexpr, but is.
282 // return std::log(x);
283 // }
284 //
285 // double operator[]( size_t i ) const
286 // {
287 // // Only called locally from within this source file, so we assume that the range is valid.
288 // // Adapt to proper boundary check if this is also exposed to the outside at some point.
289 // assert( i < n );
290 // return data[i];
291 // }
292 //
293 // size_t size() const
294 // {
295 // return n;
296 // }
297 // };
298 //
299 // Create a lookup table that behaves the same as the static std::array, but is filled at compile time.
300 // constexpr auto const log_factorial_lookup_ = LogFactorial<1024>();
301 
310 std::vector<double> log_factorial_lookup_generator_()
311 {
312  // It's stupidly hard to get double to string with a fixed number of digits (so that our
313  // printed out table is neatly aligned), while also ensuring correct rounding...
314  auto double_to_string = []( double val ){
315  std::ostringstream s;
316  s.precision( 17 );
317  s << val;
318  if( val == 0.0 ) {
319  return std::string( "0.0000000000000000" );
320  } else if( s.str().size() > 18 ) {
321  return s.str().substr( 0, 18 );
322  } else if( s.str().size() < 18 ) {
323  return s.str() + std::string( 18 - s.str().size(), '0' );
324  } else {
325  return s.str();
326  }
327  };
328 
329  // Prepare result.
330  size_t m = 1024;
331  auto logf = std::vector<double>(m);
332  logf[0] = 0.0;
333 
334  // Loop and add up log(x!) = log((x-1)!) + log(x)
335  std::cout << " " << double_to_string( logf[0] ) << ", ";
336  for( size_t i = 1; i < m; ++i) {
337  logf[i] = logf[i-1] + std::log(i);
338 
339  std::cout << double_to_string( logf[i] ) << ", ";
340  if( (i+1) % 5 == 0 ) {
341  std::cout << "\n ";
342  }
343  }
344  std::cout << "\n";
345  return logf;
346 }
347 
348 double log_factorial( size_t n )
349 {
350  // For small numbers, use or precise lookup table.
351  if( n < log_factorial_lookup_.size() ) {
352  return log_factorial_lookup_[ n ];
353  }
354 
355  // Stirling's approximation for larger values. Highly accurate, although it doesn't look like it.
356  // See https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
357  double const x = static_cast<double>( n + 1 );
358  return (x - 0.5) * std::log(x) - x + 0.5 * std::log(2 * PI) + 1.0 / (12.0 * x);
359 }
360 
361 size_t binomial_coefficient( size_t n, size_t k )
362 {
363  // Error cases.
364  if( n == 0 ) {
365  throw std::invalid_argument( "Cannot compute binomial coefficient with n == 0" );
366  }
367  if( k == 0 ) {
368  throw std::invalid_argument( "Cannot compute binomial coefficient with k == 0" );
369  }
370  if( k > n ) {
371  throw std::invalid_argument( "Cannot compute binomial coefficient with k > n" );
372  }
373 
374  // Since C(n, k) = C(n, n-k), we can shortcut.
375  if( k > n - k ) {
376  k = n - k;
377  }
378  assert( k <= n/2 );
379 
380  // We follow http://csharphelper.com/blog/2014/08/calculate-the-binomial-coefficient-n-choose-k-efficiently-in-c/
381  // In short, C(n, k) = n/k * C(n-1, k-1), which we can continue until C(n-(k-1), 1) to get
382  // a long chain product. Reading the factors from right to left, we always have the binomial
383  // coefficients of smaller n,k, and use these as our intermediates to ensure that we don't
384  // get an overflow here while multiplying them up.
385 
386  size_t result = 1;
387  for( size_t i = 1; i <= k; i++)
388  {
389  // Multiply, with overflow check.
390  assert( n > k - i );
391  if( result >= std::numeric_limits<size_t>::max() / (n - (k - i)) ) {
392  throw std::runtime_error(
393  "Cannot compute binomial coefficient with n == " + std::to_string(n) +
394  " and k == " + std::to_string(k) + " due to numerical overflow. " +
395  "Use binomial_coefficient_approx() instead for large values."
396  );
397  }
398  result *= n - (k - i);
399 
400  // As we are computing the product "backwards", we always have a valid intermediate value,
401  // that is, the division by i is always an integer division (which we assert).
402  assert( result % i == 0 );
403  result /= i;
404  }
405  return result;
406 }
407 
408 double binomial_coefficient_approx( size_t n, size_t k )
409 {
410  // Shorthand.
411  auto const& logf = log_factorial_lookup_;
412 
413  // Error cases.
414  if( n == 0 ) {
415  throw std::invalid_argument( "Cannot compute binomial coefficient with n == 0" );
416  }
417  if( k == 0 ) {
418  throw std::invalid_argument( "Cannot compute binomial coefficient with k == 0" );
419  }
420  if( k > n ) {
421  throw std::invalid_argument( "Cannot compute binomial coefficient with k > n" );
422  }
423  if( n >= logf.size() ) {
424  throw std::invalid_argument(
425  "Cannot compute binomial coefficient with n >= " + std::to_string( logf.size() ) +
426  " due to numerical overflow"
427  );
428  }
429  assert( n < logf.size() );
430  assert( k < logf.size() );
431  assert( k <= n );
432 
433  // We use a log factorial lookup table, see https://stackoverflow.com/a/37715980/4184258
434  // Basically, what we are computing is C(n, k) = n!/((n-k)! * k!), which in log form is
435  // log(n!) - log((n-k)!) - log(k!). We use a lookup table for these log factorial values,
436  // and then finally only need to exponentiate them back to get the result.
437  return std::exp( logf[n] - logf[n-k] - logf[k] );
438 }
439 
440 } // namespace utils
441 } // namespace genesis
genesis::utils::log_factorial
double log_factorial(size_t n)
Return the logarithm (base e) of the factorial of n, that is log(n!).
Definition: common.cpp:348
common.hpp
genesis::utils::log_factorial_lookup_generator_
std::vector< double > log_factorial_lookup_generator_()
Helper function to generate above lookup table.
Definition: common.cpp:310
genesis::utils::binomial_coefficient
size_t binomial_coefficient(size_t n, size_t k)
Compute the binomial coefficient, that is n choose k, for two integer numbers.
Definition: common.cpp:361
genesis::utils::PI
constexpr double PI
Make the world go round.
Definition: common.hpp:55
genesis
Container namespace for all symbols of genesis in order to keep them separate when used as a library.
Definition: placement/formats/edge_color.cpp:42
genesis::population::to_string
std::string to_string(GenomeRegion const &region)
Definition: functions/genome_region.cpp:55
genesis::utils::binomial_coefficient_approx
double binomial_coefficient_approx(size_t n, size_t k)
Compute the binomial coefficient, that is n choose k, for two integer numbers, for large numbers.
Definition: common.cpp:408
genesis::utils::log_factorial_lookup_
static const std::array< double, 1024 > log_factorial_lookup_
Log-factorial log(x!) lookup table for fast computation of binomial coefficients. See binomial_coeffi...
Definition: common.cpp:53