/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2021 Willem Elbers (whe@willemelbers.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
*
******************************************************************************/
/**
* @file fermi_dirac.c
*
* @brief Fast generation of pseudo-random numbers from a Fermi-Dirac
* distribution using numerical inversion (Hormann & Leydold, 2003). The
* spline tables were generated with AnyRNG. For more details, refer
* to https://github.com/wullm/AnyRNG.
*/
/* This object's header. */
#include "fermi_dirac.h"
/* Standard headers */
#include
/* Fast optimized logarithm */
#include "../../log.h"
/* Cubic spline coefficients */
struct spline {
double a0, a1, a2, a3;
};
/**
* @brief Interpolation and search tables for the Fermi-Dirac distribution
*/
struct anyrng {
/*! Number of intervals on which the interpolation is defined */
int intervalN;
/*! Endpoints of the intervals */
double *endpoints;
/*! Cubic splines of the Fermi-Dirac quantile function */
struct spline *splines;
/*! Length of the look up tables */
int tablelen;
/*! Search table to look up enclosing intervals for small u */
int *index_table_a;
/*! Search table to look up enclosing intervals for intermediate u */
int *index_table_b;
/*! Search table to look up enclosing intervals for large u */
int *index_table_c;
};
static double endpoints[118] = {
0.000000e+00, 8.812220e-16, 3.490090e-15, 1.764874e-14, 1.080345e-13,
7.475047e-13, 5.544411e-12, 4.267185e-11, 3.346975e-10, 2.649997e-09,
2.107211e-08, 7.092068e-08, 1.677779e-07, 5.644860e-07, 1.334411e-06,
2.599613e-06, 4.480999e-06, 1.057018e-05, 2.054614e-05, 3.533463e-05,
5.584298e-05, 8.296016e-05, 1.175567e-04, 1.604850e-04, 2.125789e-04,
2.746542e-04, 4.319210e-04, 6.384483e-04, 9.001106e-04, 1.222496e-03,
1.610910e-03, 2.070375e-03, 2.605637e-03, 3.221165e-03, 3.921158e-03,
4.709548e-03, 5.590000e-03, 6.565923e-03, 8.816531e-03, 1.148360e-02,
1.458549e-02, 1.813688e-02, 2.214891e-02, 2.662937e-02, 3.158292e-02,
3.701122e-02, 4.291319e-02, 4.928517e-02, 5.612111e-02, 6.341282e-02,
7.115008e-02, 8.791173e-02, 1.062920e-01, 1.261570e-01, 1.473582e-01,
1.697373e-01, 1.931305e-01, 2.173724e-01, 2.422994e-01, 2.677516e-01,
2.935760e-01, 3.196273e-01, 3.457697e-01, 3.718776e-01, 3.978360e-01,
4.235411e-01, 4.488999e-01, 4.738305e-01, 5.221310e-01, 5.679885e-01,
6.110917e-01, 6.512521e-01, 6.883826e-01, 7.224780e-01, 7.535962e-01,
7.818432e-01, 8.073587e-01, 8.303052e-01, 8.508587e-01, 8.692016e-01,
8.855172e-01, 8.999850e-01, 9.127781e-01, 9.240609e-01, 9.339876e-01,
9.427015e-01, 9.503348e-01, 9.570084e-01, 9.628322e-01, 9.679056e-01,
9.723182e-01, 9.761502e-01, 9.794730e-01, 9.823505e-01, 9.848391e-01,
9.869886e-01, 9.888431e-01, 9.904414e-01, 9.918172e-01, 9.930005e-01,
9.948897e-01, 9.962793e-01, 9.972980e-01, 9.980425e-01, 9.985850e-01,
9.989794e-01, 9.992652e-01, 9.994720e-01, 9.996213e-01, 9.998061e-01,
9.999013e-01, 9.999500e-01, 9.999748e-01, 9.999937e-01, 9.999984e-01,
9.999996e-01, 9.999999e-01, 1.000000e+00};
static struct spline splines[118] = {
{1.000000e-05, 3.177853e-05, -3.440759e-05, 1.454999e-05},
{2.192092e-05, 1.957877e-05, -1.160959e-05, 3.951739e-06},
{3.384185e-05, 4.458279e-05, -3.298529e-05, 1.224435e-05},
{5.768370e-05, 9.796089e-05, -8.223074e-05, 3.195355e-05},
{1.053674e-04, 2.077194e-04, -1.865719e-04, 7.421998e-05},
{2.007348e-04, 4.293444e-04, -3.993851e-04, 1.607755e-04},
{3.914696e-04, 8.738366e-04, -8.274556e-04, 3.350886e-04},
{7.729391e-04, 1.763374e-03, -1.684705e-03, 6.842694e-04},
{1.535878e-03, 3.542204e-03, -3.398798e-03, 1.382473e-03},
{3.061757e-03, 7.097573e-03, -6.822804e-03, 2.776987e-03},
{6.113513e-03, 4.824439e-03, -2.643433e-03, 8.707511e-04},
{9.165270e-03, 4.177172e-03, -1.553626e-03, 4.282101e-04},
{1.221703e-02, 9.643782e-03, -5.278698e-03, 1.738429e-03},
{1.832054e-02, 8.348620e-03, -3.098591e-03, 8.534846e-04},
{2.442405e-02, 7.742955e-03, -2.147031e-03, 5.075889e-04},
{3.052757e-02, 7.392976e-03, -1.625943e-03, 3.364794e-04},
{3.663108e-02, 1.666989e-02, -6.155448e-03, 1.692581e-03},
{4.883811e-02, 1.546030e-02, -4.257164e-03, 1.003894e-03},
{6.104513e-02, 1.476135e-02, -3.217990e-03, 6.636676e-04},
{7.325216e-02, 1.430652e-02, -2.570199e-03, 4.707042e-04},
{8.545918e-02, 1.398707e-02, -2.130803e-03, 3.507572e-04},
{9.766621e-02, 1.375043e-02, -1.814549e-03, 2.711427e-04},
{1.098732e-01, 1.356811e-02, -1.576708e-03, 2.156231e-04},
{1.220803e-01, 1.342333e-02, -1.391687e-03, 1.753831e-04},
{1.342873e-01, 1.330557e-02, -1.243846e-03, 1.452992e-04},
{1.464943e-01, 2.851127e-02, -5.005621e-03, 9.084052e-04},
{1.709084e-01, 2.787360e-02, -4.132288e-03, 6.727367e-04},
{1.953224e-01, 2.740090e-02, -3.503576e-03, 5.167249e-04},
{2.197365e-01, 2.703641e-02, -3.030574e-03, 4.082211e-04},
{2.441505e-01, 2.674669e-02, -2.662433e-03, 3.297933e-04},
{2.685646e-01, 2.651082e-02, -2.368087e-03, 2.713228e-04},
{2.929786e-01, 2.631498e-02, -2.127539e-03, 2.266089e-04},
{3.173927e-01, 2.614973e-02, -1.927362e-03, 1.916808e-04},
{3.418067e-01, 2.600838e-02, -1.758224e-03, 1.639006e-04},
{3.662208e-01, 2.588603e-02, -1.613441e-03, 1.414606e-04},
{3.906348e-01, 2.577907e-02, -1.488110e-03, 1.230891e-04},
{4.150489e-01, 2.568473e-02, -1.378552e-03, 1.078702e-04},
{4.394629e-01, 5.362050e-02, -5.511167e-03, 7.187774e-04},
{4.882911e-01, 5.303602e-02, -4.778531e-03, 5.706129e-04},
{5.371192e-01, 5.255849e-02, -4.191326e-03, 4.609411e-04},
{5.859473e-01, 5.216054e-02, -3.710147e-03, 3.777099e-04},
{6.347754e-01, 5.182344e-02, -3.308556e-03, 3.132206e-04},
{6.836035e-01, 5.153391e-02, -2.968181e-03, 2.623720e-04},
{7.324316e-01, 5.128230e-02, -2.675876e-03, 2.216775e-04},
{7.812597e-01, 5.106141e-02, -2.421999e-03, 1.886903e-04},
{8.300878e-01, 5.086577e-02, -2.199319e-03, 1.616529e-04},
{8.789159e-01, 5.069114e-02, -2.002315e-03, 1.392774e-04},
{9.277440e-01, 5.053420e-02, -1.826699e-03, 1.206032e-04},
{9.765721e-01, 5.039230e-02, -1.669092e-03, 1.049023e-04},
{1.025400e+00, 5.026329e-02, -1.526795e-03, 9.161505e-05},
{1.074228e+00, 1.028683e-01, -5.815751e-03, 6.036140e-04},
{1.171885e+00, 1.020332e-01, -4.846724e-03, 4.697008e-04},
{1.269541e+00, 1.013216e-01, -4.035237e-03, 3.698061e-04},
{1.367197e+00, 1.007075e-01, -3.345636e-03, 2.943847e-04},
{1.464853e+00, 1.001716e-01, -2.752379e-03, 2.369418e-04},
{1.562509e+00, 9.969996e-02, -2.236687e-03, 1.929401e-04},
{1.660166e+00, 9.928152e-02, -1.784441e-03, 1.591361e-04},
{1.757822e+00, 9.890786e-02, -1.384813e-03, 1.331654e-04},
{1.855478e+00, 9.857228e-02, -1.029346e-03, 1.132745e-04},
{1.953134e+00, 9.826940e-02, -7.113359e-04, 9.814323e-05},
{2.050790e+00, 9.799484e-02, -4.253906e-04, 8.676449e-05},
{2.148447e+00, 9.774497e-02, -1.671167e-04, 7.836135e-05},
{2.246103e+00, 9.751678e-02, 6.710498e-05, 7.232835e-05},
{2.343759e+00, 9.730773e-02, 2.802876e-04, 6.818952e-05},
{2.441415e+00, 9.711568e-02, 4.749628e-04, 6.556786e-05},
{2.539071e+00, 9.693877e-02, 6.532755e-04, 6.416282e-05},
{2.636728e+00, 9.677542e-02, 8.170569e-04, 6.373341e-05},
{2.734384e+00, 1.910291e-01, 3.767224e-03, 5.160521e-04},
{2.929696e+00, 1.899901e-01, 4.782956e-03, 5.393740e-04},
{3.125009e+00, 1.890915e-01, 5.646906e-03, 5.740339e-04},
{3.320321e+00, 1.883086e-01, 6.388338e-03, 6.155239e-04},
{3.515634e+00, 1.876218e-01, 7.029801e-03, 6.608298e-04},
{3.710946e+00, 1.870156e-01, 7.588922e-03, 7.079357e-04},
{3.906258e+00, 1.864773e-01, 8.079643e-03, 7.555006e-04},
{4.101571e+00, 1.859967e-01, 8.513103e-03, 8.026422e-04},
{4.296883e+00, 1.855653e-01, 8.898285e-03, 8.487922e-04},
{4.492196e+00, 1.851763e-01, 9.242495e-03, 8.935968e-04},
{4.687508e+00, 1.848239e-01, 9.551712e-03, 9.368494e-04},
{4.882821e+00, 1.845031e-01, 9.830867e-03, 9.784433e-04},
{5.078133e+00, 1.842100e-01, 1.008405e-02, 1.018340e-03},
{5.273445e+00, 1.839412e-01, 1.031468e-02, 1.056544e-03},
{5.468758e+00, 1.836937e-01, 1.052561e-02, 1.093093e-03},
{5.664070e+00, 1.834651e-01, 1.071926e-02, 1.128039e-03},
{5.859383e+00, 1.832533e-01, 1.089768e-02, 1.161448e-03},
{6.054695e+00, 1.830564e-01, 1.106262e-02, 1.193391e-03},
{6.250007e+00, 1.828729e-01, 1.121556e-02, 1.223940e-03},
{6.445320e+00, 1.827015e-01, 1.135780e-02, 1.253169e-03},
{6.640632e+00, 1.825408e-01, 1.149043e-02, 1.281149e-03},
{6.835945e+00, 1.823900e-01, 1.161443e-02, 1.307949e-03},
{7.031257e+00, 1.822482e-01, 1.173063e-02, 1.333635e-03},
{7.226570e+00, 1.821144e-01, 1.183975e-02, 1.358270e-03},
{7.421882e+00, 1.819881e-01, 1.194245e-02, 1.381911e-03},
{7.617194e+00, 1.818685e-01, 1.203929e-02, 1.404614e-03},
{7.812507e+00, 1.817552e-01, 1.213077e-02, 1.426431e-03},
{8.007819e+00, 1.816477e-01, 1.221734e-02, 1.447411e-03},
{8.203132e+00, 1.815454e-01, 1.229938e-02, 1.467598e-03},
{8.398444e+00, 1.814481e-01, 1.237726e-02, 1.487035e-03},
{8.593757e+00, 1.813554e-01, 1.245129e-02, 1.505762e-03},
{8.789069e+00, 1.812668e-01, 1.252176e-02, 1.523815e-03},
{8.984381e+00, 3.367030e-01, 4.147623e-02, 1.244560e-02},
{9.375006e+00, 3.361379e-01, 4.177920e-02, 1.270774e-02},
{9.765631e+00, 3.356191e-01, 4.205304e-02, 1.295269e-02},
{1.015626e+01, 3.351411e-01, 4.230171e-02, 1.318205e-02},
{1.054688e+01, 3.346991e-01, 4.252850e-02, 1.339724e-02},
{1.093751e+01, 3.342892e-01, 4.273612e-02, 1.359950e-02},
{1.132813e+01, 3.339080e-01, 4.292687e-02, 1.378995e-02},
{1.171876e+01, 3.335526e-01, 4.310271e-02, 1.396956e-02},
{1.210938e+01, 3.332203e-01, 4.326528e-02, 1.413922e-02},
{1.250000e+01, 5.722553e-01, 9.204483e-02, 1.169495e-01},
{1.328125e+01, 5.704529e-01, 9.145775e-02, 1.193390e-01},
{1.406250e+01, 5.688543e-01, 9.089620e-02, 1.214992e-01},
{1.484375e+01, 5.674266e-01, 9.036152e-02, 1.234615e-01},
{1.562500e+01, 8.513274e-01, -3.717295e-01, 1.082901e+00},
{1.718750e+01, 8.456524e-01, -3.938066e-01, 1.110653e+00},
{1.875000e+01, 8.409733e-01, -4.130061e-01, 1.134532e+00},
{2.031250e+01, 8.371720e-01, -4.305845e-01, 1.155912e+00},
{2.187500e+01, 8.346323e-01, -4.503617e-01, 1.178229e+00},
{2.343750e+01, 8.359434e-01, -4.895204e-01, 1.216076e+00}};
static int index_table_a[165] = {
0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4,
4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6,
6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9,
10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 13, 13, 13,
13, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 18, 18, 18, 19,
19, 20, 20, 21, 21, 22, 23, 24, 24, 24, 25, 25, 26, 26, 27, 27, 28, 29, 29,
30, 31, 32, 33, 34, 35, 36, 36, 37, 38, 38, 39, 40};
static int index_table_b[165] = {
41, 42, 43, 44, 45, 46, 47, 48, 48, 49, 49, 50, 50, 50, 50, 51, 51, 51, 52,
52, 52, 52, 53, 53, 53, 53, 54, 54, 54, 54, 55, 55, 55, 55, 56, 56, 56, 56,
57, 57, 57, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 59, 60, 60, 60, 60, 61,
61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65,
65, 65, 66, 66, 66, 66, 66, 66, 66, 66, 66, 67, 67, 67, 67, 67, 67, 67, 67,
68, 68, 68, 68, 68, 68, 68, 69, 69, 69, 69, 69, 69, 69, 70, 70, 70, 70, 70,
70, 70, 71, 71, 71, 71, 71, 71, 72, 72, 72, 72, 72, 73, 73, 73, 73, 73, 74,
74, 74, 74, 75, 75, 75, 75, 76, 76, 76, 76, 77, 77, 77, 78, 78, 78, 79, 79,
80, 80, 80, 81, 81, 82, 83, 83, 84, 85, 86, 87, 88};
static int index_table_c[165] = {
89, 90, 90, 91, 92, 92, 93, 93, 94, 94, 95, 96, 96, 97, 97,
98, 98, 98, 99, 99, 99, 99, 100, 100, 100, 101, 101, 101, 101, 102,
102, 102, 102, 103, 103, 103, 104, 104, 104, 104, 105, 105, 105, 105, 106,
106, 106, 107, 107, 107, 107, 107, 107, 107, 107, 108, 108, 108, 108, 108,
108, 108, 109, 109, 109, 109, 109, 109, 109, 109, 110, 110, 110, 110, 110,
110, 110, 110, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111, 111,
111, 111, 111, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112, 112,
112, 112, 112, 112, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113, 113,
113, 113, 113, 113, 113, 114, 114, 114, 114, 114, 114, 114, 114, 114, 114,
114, 114, 114, 114, 114, 114, 115, 115, 115, 115, 115, 115, 115, 115, 115,
115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115, 115};
static struct anyrng anyrng = {118, endpoints, splines, 165,
index_table_a, index_table_b, index_table_c};
/**
* @brief Transform a 64-bit unsigned integer seed into a (dimensionless)
* Fermi-Dirac momentum (units of kb*T), using cubic spline interpolation of
* the quantile function.
*
* @param seed Random seed to be transformed
*/
double neutrino_seed_to_fermi_dirac(uint64_t seed) {
/* Scramble the bits with splitmix64 */
uint64_t A = seed;
A = A + 0x9E3779B97f4A7C15;
A = (A ^ (A >> 30)) * 0xBF58476D1CE4E5B9;
A = (A ^ (A >> 27)) * 0x94D049BB133111EB;
A = A ^ (A >> 31);
/* Map the integer to the unit open interval (0, 1) */
const double u = ((double)A + 0.5) / ((double)UINT64_MAX + 1);
/* Use the hash table to find an enclosing interval */
const int tablen = anyrng.tablelen;
int index, interval;
if (u < 0.025) {
index = (int)((optimized_log10f(u) + 14.5229) / 12.9208 * tablen);
interval = anyrng.index_table_a[index < tablen ? index : tablen - 1] + 1;
} else if (u < 0.975) {
index = (int)((u - 0.025) / 0.95 * tablen);
interval = anyrng.index_table_b[index < tablen ? index : tablen - 1] + 1;
} else {
index = (int)(-(optimized_log10f(1 - u) + 1.60206) / 6.39794 * tablen);
interval = anyrng.index_table_c[index < tablen ? index : tablen - 1] + 1;
}
/* Retrieve the endpoints and cubic spline coefficients of this interval */
const double Fl = anyrng.endpoints[interval];
const double Fr = anyrng.endpoints[interval + 1];
struct spline *iv = &anyrng.splines[interval];
/* Evaluate F^-1(u) using the Hermite approximation of F in this interval */
const double u_tilde = (u - Fl) / (Fr - Fl);
const double u_tilde2 = u_tilde * u_tilde;
const double u_tilde3 = u_tilde2 * u_tilde;
return iv->a0 + iv->a1 * u_tilde + iv->a2 * u_tilde2 + iv->a3 * u_tilde3;
}
/**
* @brief Transform a 64-bit unsigned integer seed into a point on the sphere.
*
* @param seed Random seed to be transformed
*/
void neutrino_seed_to_direction(uint64_t seed, double n[3]) {
/* Skip the first draw, which is used for the magnitude */
seed += 0x9E3779B97f4A7C15;
/* Draw three numbers for the orientation of the momentum vector */
double u[3];
for (int i = 0; i < 3; i++) {
/* Scramble the bits with splitmix64 */
uint64_t N = seed += 0x9E3779B97f4A7C15;
N = (N ^ (N >> 30)) * 0xBF58476D1CE4E5B9;
N = (N ^ (N >> 27)) * 0x94D049BB133111EB;
N = N ^ (N >> 31);
/* Map the integer to the unit open interval (0,1) */
u[i] = ((double)N + 0.5) / ((double)UINT64_MAX + 1);
}
/* Map to three independent Gaussians with Box-Mueller */
const double sqrt_log_u = sqrt(-2 * log(u[0]));
n[0] = sqrt_log_u * cos(2 * M_PI * u[1]);
n[1] = sqrt_log_u * sin(2 * M_PI * u[1]);
n[2] = sqrt_log_u * cos(2 * M_PI * u[2]);
/* Normalize the vector */
const double nmag = sqrtf(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
if (nmag > 0) {
const double inv_nmag = 1.0 / nmag;
n[0] *= inv_nmag;
n[1] *= inv_nmag;
n[2] *= inv_nmag;
}
}