Subversion Repositories gng

/*
 * Piecewise polynomial approximation of
 * inverse of the normal cumulative distribution function (CDF)
 */
 
/*
 * Copyright (C) 2014, Guangxi Liu <guangxi.liu@opencores.org>
 *
 * This source file may be used and distributed without restriction provided
 * that this copyright statement is not removed from the file and that any
 * derivative work contains the original copyright notice and the associated
 * disclaimer.
 *
 * This source file 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 2.1 of the License,
 * or (at your option) any later version.
 *
 * This source 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 Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this source; if not, download it from
 * http://www.opencores.org/lgpl.shtml
 */
 
 
#include "icdf.h"
 
 
/* Coefficients table */
static const long C0[248] = {
    11049, 8007, 5220, 2577, 18846, 16547, 14535, 12721, 25134, 23229,
    21594, 20150, 30518, 28863, 27458, 26231, 35288, 33808, 32562, 31479,
    39608, 38260, 37130, 36152, 43582, 42336, 41296, 40399, 47277, 46115,
    45147, 44314, 50745, 49651, 48742, 47963, 54020, 52985, 52126, 51390,
    57132, 56147, 55331, 54632, 60101, 59160, 58381, 57715, 62945, 62043,
    61296, 60659, 65679, 64811, 64093, 63481, 68314, 67476, 66784, 66194,
    70859, 70049, 69381, 68811, 73323, 72538, 71891, 71340, 75714, 74952,
    74324, 73790, 78036, 77296, 76686, 76167, 80297, 79576, 78982, 78477,
    82500, 81797, 81218, 80726, 84649, 83963, 83398, 82918, 86748, 86078,
    85526, 85057, 88800, 88145, 87606, 87147, 90809, 90167, 89640, 89191,
    92777, 92148, 91631, 91192, 94706, 94090, 93583, 93152, 96599, 95994,
    95496, 95074, 98457, 97863, 97374, 96960, 100282, 99698, 99219, 98811,
    102076, 101502, 101031, 100630, 103841, 103276, 102812, 102419, 105577, 105021,
    104565, 104178, 107286, 106739, 106290, 105909, 108970, 108431, 107989, 107613,
    110629, 110098, 109662, 109292, 112265, 111741, 111311, 110946, 113878, 113361,
    112937, 112578, 115469, 114959, 114541, 114186, 117040, 116536, 116124, 115774,
    118590, 118093, 117686, 117340, 120121, 119630, 119228, 118887, 121633, 121149,
    120751, 120414, 123128, 122649, 122256, 121923, 124605, 124132, 123743, 123414,
    126065, 125597, 125213, 124888, 127510, 127047, 126667, 126345, 128938, 128480,
    128105, 127786, 130351, 129898, 129527, 129212, 131750, 131301, 130934, 130622,
    133134, 132690, 132326, 132018, 134505, 134065, 133705, 133400, 135862, 135426,
    135070, 134767, 137206, 136775, 136421, 136122, 138537, 138110, 137760, 137463,
    139856, 139433, 139086, 138792, 141163, 140743, 140400, 140109, 142458, 142042,
    141702, 141413, 143742, 143330, 142992, 142706, 145015, 144606, 144272, 143988,
    146277, 145872, 145540, 145259, 150000, 148769, 147528, 146797
};
 
static const long C1[248] = {
    -102514, -92199, -86162, -82951, -79070, -68105, -60693, -55396, -66134, -55839,
    -48783, -43639, -57755, -48215, -41670, -36890, -51803, -42932, -36856, -32421,
    -47313, -39012, -33338, -29203, -43778, -35964, -30631, -26751, -40907, -33511,
    -28470, -24807, -38519, -31483, -26695, -23220, -36492, -29774, -25206, -21893,
    -34747, -28307, -23933, -20765, -33223, -27032, -22831, -19789, -31878, -25911,
    -21864, -18936, -30681, -24915, -21007, -18182, -29606, -24023, -20242, -17509,
    -28634, -23218, -19552, -16903, -27750, -22488, -18927, -16356, -26942, -21821,
    -18357, -15857, -26198, -21208, -17835, -15400, -25512, -20644, -17354, -14980,
    -24876, -20121, -16909, -14591, -24285, -19636, -16496, -14231, -23733, -19184,
    -16111, -13896, -23217, -18760, -15752, -13583, -22732, -18364, -15415, -13290,
    -22276, -17991, -15099, -13015, -21846, -17639, -14801, -12756, -21440, -17307,
    -14520, -12512, -21055, -16993, -14253, -12281, -20690, -16695, -14001, -12062,
    -20342, -16412, -13762, -11854, -20012, -16143, -13534, -11656, -19697, -15886,
    -13317, -11468, -19396, -15641, -13111, -11289, -19109, -15407, -12913, -11118,
    -18834, -15183, -12724, -10954, -18570, -14969, -12543, -10797, -18317, -14763,
    -12369, -10646, -18073, -14565, -12202, -10502, -17839, -14375, -12042, -10363,
    -17614, -14193, -11888, -10230, -17397, -14016, -11739, -10102, -17188, -13847,
    -11596, -9978, -16986, -13683, -11458, -9858, -16791, -13525, -11325, -9743,
    -16603, -13372, -11196, -9632, -16421, -13224, -11072, -9524, -16244, -13081,
    -10951, -9420, -16073, -12942, -10835, -9319, -15907, -12808, -10722, -9222,
    -15747, -12678, -10612, -9127, -15591, -12552, -10506, -9035, -15440, -12429,
    -10403, -8946, -15294, -12311, -10303, -8860, -15154, -12196, -10206, -8776,
    -15020, -12086, -10113, -8695, -14894, -11982, -10024, -8618, -14779, -11885,
    -9940, -8544, -14681, -11797, -9863, -8475, -13792, -11190, -9418, -8133,
    0, 0, 0, 0, 0, 0, 0, 0
};
 
static const long C2[248] = {
    83840, 48883, 25920, 8147, 88988, 59916, 42707, 31577, 83526, 57019,
    41464, 31523, 77411, 52885, 38533, 29386, 71990, 49106, 35743, 27241,
    67363, 45865, 33332, 25372, 63422, 43107, 31282, 23780, 60040, 40746,
    29530, 22422, 57109, 38706, 28020, 21254, 54545, 36926, 26705, 20239,
    52279, 35358, 25549, 19349, 50262, 33965, 24525, 18561, 48452, 32718,
    23609, 17858, 46818, 31593, 22785, 17226, 45333, 30573, 22038, 16654,
    43976, 29643, 21358, 16134, 42731, 28790, 20735, 15658, 41583, 28005,
    20163, 15221, 40521, 27279, 19634, 14817, 39534, 26606, 19143, 14444,
    38614, 25978, 18687, 14096, 37754, 25392, 18261, 13772, 36948, 24844,
    17862, 13468, 36190, 24328, 17488, 13184, 35477, 23843, 17136, 12916,
    34803, 23385, 16804, 12664, 34165, 22952, 16490, 12426, 33561, 22542,
    16193, 12200, 32987, 22153, 15911, 11986, 32441, 21783, 15643, 11783,
    31921, 21431, 15388, 11590, 31425, 21094, 15145, 11406, 30951, 20773,
    14913, 11230, 30497, 20466, 14691, 11061, 30063, 20173, 14479, 10901,
    29646, 19891, 14275, 10746, 29246, 19620, 14080, 10599, 28861, 19360,
    13892, 10456, 28492, 19110, 13711, 10320, 28135, 18870, 13538, 10189,
    27792, 18638, 13370, 10062, 27461, 18414, 13209, 9940, 27141, 18199,
    13053, 9822, 26832, 17990, 12903, 9709, 26534, 17788, 12758, 9599,
    26245, 17593, 12617, 9492, 25965, 17405, 12481, 9390, 25694, 17222,
    12349, 9290, 25432, 17045, 12222, 9194, 25178, 16874, 12098, 9100,
    24933, 16708, 11979, 9010, 24697, 16549, 11864, 8923, 24474, 16397,
    11754, 8840, 24268, 16255, 11650, 8761, 24088, 16129, 11557, 8689,
    23954, 16029, 11480, 8628, 23908, 15977, 11432, 8586, 24029, 16017,
    11440, 8581, 24481, 16238, 11559, 8648, 0, 0, 0, 0,
    0, 0, 0, 0, 0, 0, 0, 0
};
 
 
/* Calculate Inverse of the normal CDF */
long icdf(unsigned long long n)
{
    unsigned long long t;
    int i;
    unsigned num_lzd;
    unsigned addr;
    long c0, c1, c2;
    long x;
    long long y;
    int carry;
 
    t = n;
    num_lzd = 0;
    for (i = 0; i < 61; i++) {
        if (t & 0x8000000000000000ULL)
            break;
        else {
            ++num_lzd;
            t <<= 1;
        }
    }
    addr = num_lzd << 2;
    if (n & 2ULL)
        addr += 2;
    if (n & 4ULL)
        ++addr;
 
    c0 = C0[addr];
    c1 = C1[addr];
    c2 = C2[addr];
 
    if (num_lzd <= 60)
        n &= ~(1ULL << (63 - num_lzd));
 
    t = n >> 3;
    x = 0;
    for (i = 0; i < 15; i++) {
        x <<= 1;
        if (t & 1ULL)
            ++x;
        t >>= 1;
    }
 
    y = (long long)c2 * (long long)x;
    y += ((long long)c1 << 19);
    y >>= 20;
    y *= (long long)x;
    y >>= 19;
    y += (long long)c0;
    carry = (y & 4ULL) ? 1 : 0;
    y >>= 3;
    y += carry;
    if (n & 1ULL)
        y = -y;
 
    return (long)y;
}