precomputation.h

/*********************************************************************
*
* precomputation of 
*      [        oo                                   ]
*      |       /                                     |
*      |  1    |                 2                   |
*  Real| ----  | sinc(pi Delta t)   psi(1 +  it) dt  |
*      |  pi   |                                     |
*      |       /                                     |
*      [        -oo                                  ]
*
* for Delta = 1.0, 1.005, 1.01, ..., 5.0
*
* accurate to (hopefully) within about 10^(-5)
*
* To make this list, do something like the folowing in sage:
* (will take about an hour on a fast machine)
*
* TODO: do this computation more accurately
*----------------------------------------------------------
import scipy
from scipy.special import digamma


def sinc_square(u, Delta):
    x = u * Delta * RR(pi)

    if abs(x) < .00001:
        return 1 - 1.0/3.0 * x^2 + 2.0/45.0 * x^4 - 1.0/315.0 * x^6
    else:
        a = sin(x)/x
        return a*a;

def foo2(Delta):
    f = lambda t : sinc_square(t, Delta)
    g = lambda t : RR(f(t)) * digamma(1 + CC.0 * t).real
    
    integrals = 0
    error = 0

    for n in range(100):
        X = numerical_integral(g, n * 100, (n + 1)*100)
        integrals += X[0]
        error += X[1]

    print "Error is maybe < ", error
    return RR(2) * integrals / RR(pi)



def do_precomputation():
    L = []
    for Delta in [1.00, 1.005, .., 5.0]:
        x = foo2(Delta)
        print Delta, x
        L.append(x)

    count = 0
    for x in L:
        count = count + 1
        print "{x},".format(x=x.str()[:9]),
        if count % 7 == 0:
            print

do_precomputation()

****************************************************/

#include <iostream>
#include <cmath>

double precomputed_gamma_terms[801] = {
-0.100523, -0.100436, -0.100342, -0.100245, -0.100147, -0.100042, -0.099937,
-0.099829, -0.099723, -0.099606, -0.099490, -0.099373, -0.099258, -0.099132,
-0.099008, -0.098881, -0.098755, -0.098626, -0.098495, -0.098360, -0.098223,
-0.098092, -0.097955, -0.097816, -0.097673, -0.097535, -0.097392, -0.097248,
-0.097101, -0.096956, -0.096809, -0.096660, -0.096513, -0.096360, -0.096208,
-0.096056, -0.095907, -0.095749, -0.095594, -0.095438, -0.095285, -0.095125,
-0.094967, -0.094809, -0.094652, -0.094491, -0.094331, -0.094171, -0.094008,
-0.093848, -0.093687, -0.093524, -0.093358, -0.093199, -0.093036, -0.092872,
-0.092705, -0.092545, -0.092380, -0.092216, -0.092049, -0.091886, -0.091721,
-0.091556, -0.091391, -0.091226, -0.091060, -0.090894, -0.090730, -0.090563,
-0.090397, -0.090232, -0.090069, -0.089900, -0.089735, -0.089569, -0.089406,
-0.089238, -0.089072, -0.088907, -0.088743, -0.088576, -0.088411, -0.088246,
-0.088083, -0.087916, -0.087752, -0.087587, -0.087423, -0.087259, -0.087095,
-0.086931, -0.086766, -0.086605, -0.086442, -0.086278, -0.086114, -0.085954,
-0.085791, -0.085630, -0.085464, -0.085306, -0.085145, -0.084984, -0.084824,
-0.084663, -0.084503, -0.084344, -0.084178, -0.084025, -0.083866, -0.083707,
-0.083551, -0.083392, -0.083234, -0.083076, -0.082910, -0.082763, -0.082606,
-0.082398, -0.082303, -0.082139, -0.081984, -0.081830, -0.081674, -0.081521,
-0.081368, -0.081214, -0.081055, -0.080909, -0.080754, -0.080605, -0.080463,
-0.080302, -0.080151, -0.080001, -0.079851, -0.079701, -0.079552, -0.079403,
-0.079252, -0.079107, -0.078961, -0.078811, -0.078671, -0.078518, -0.078371,
-0.078226, -0.078083, -0.077935, -0.077790, -0.077647, -0.077497, -0.077358,
-0.077215, -0.077072, -0.076928, -0.076788, -0.076646, -0.076505, -0.076342,
-0.076223, -0.076083, -0.075943, -0.075804, -0.075665, -0.075527, -0.075388,
-0.075245, -0.075113, -0.074976, -0.074839, -0.074699, -0.074567, -0.074432,
-0.074296, -0.074165, -0.074027, -0.073893, -0.073760, -0.073628, -0.073494,
-0.073361, -0.073229, -0.073095, -0.072966, -0.072835, -0.072706, -0.072570,
-0.072445, -0.072315, -0.072186, -0.072054, -0.071929, -0.071801, -0.071673,
-0.071549, -0.071419, -0.071293, -0.071167, -0.071044, -0.070916, -0.070791,
-0.070665, -0.070542, -0.070418, -0.070295, -0.070171, -0.070047, -0.069926,
-0.069804, -0.069682, -0.069560, -0.069440, -0.069319, -0.069199, -0.069081,
-0.068959, -0.068840, -0.068721, -0.068607, -0.068485, -0.068367, -0.068249,
-0.068135, -0.068015, -0.067899, -0.067783, -0.067667, -0.067551, -0.067435,
-0.067322, -0.067206, -0.067093, -0.066979, -0.066866, -0.066751, -0.066640,
-0.066526, -0.066416, -0.066303, -0.066192, -0.066081, -0.065971, -0.065863,
-0.065750, -0.065640, -0.065531, -0.065424, -0.065313, -0.065204, -0.065096,
-0.064989, -0.064881, -0.064774, -0.064667, -0.064560, -0.064454, -0.064348,
-0.064242, -0.064132, -0.064032, -0.063927, -0.063823, -0.063719, -0.063615,
-0.063511, -0.063408, -0.063298, -0.063202, -0.063100, -0.062998, -0.062902,
-0.062795, -0.062694, -0.062593, -0.062488, -0.062392, -0.062292, -0.062193,
-0.062097, -0.061994, -0.061895, -0.061797, -0.061699, -0.061601, -0.061503,
-0.061406, -0.061304, -0.061212, -0.061115, -0.061019, -0.060924, -0.060827,
-0.060732, -0.060637, -0.060539, -0.060447, -0.060353, -0.060259, -0.060159,
-0.060071, -0.059978, -0.059885, -0.059793, -0.059700, -0.059608, -0.059516,
-0.059425, -0.059332, -0.059241, -0.059150, -0.059058, -0.058969, -0.058879,
-0.058789, -0.058702, -0.058610, -0.058521, -0.058432, -0.058346, -0.058255,
-0.058167, -0.058079, -0.057990, -0.057904, -0.057817, -0.057730, -0.057641,
-0.057557, -0.057471, -0.057384, -0.057300, -0.057214, -0.057129, -0.057044,
-0.056960, -0.056874, -0.056790, -0.056706, -0.056619, -0.056539, -0.056455,
-0.056372, -0.056288, -0.056207, -0.056124, -0.056042, -0.055963, -0.055878,
-0.055797, -0.055715, -0.055638, -0.055553, -0.055473, -0.055392, -0.055313,
-0.055232, -0.055152, -0.055073, -0.054993, -0.054914, -0.054836, -0.054757,
-0.054680, -0.054600, -0.054522, -0.054444, -0.054369, -0.054289, -0.054212,
-0.054135, -0.054059, -0.053981, -0.053905, -0.053829, -0.053751, -0.053677,
-0.053601, -0.053526, -0.053449, -0.053376, -0.053301, -0.053227, -0.053152,
-0.053078, -0.053004, -0.052930, -0.052858, -0.052783, -0.052710, -0.052637,
-0.052564, -0.052491, -0.052419, -0.052347, -0.052274, -0.052203, -0.052131,
-0.052060, -0.051987, -0.051917, -0.051846, -0.051776, -0.051706, -0.051635,
-0.051564, -0.051494, -0.051428, -0.051355, -0.051285, -0.051216, -0.051149,
-0.051078, -0.051009, -0.050941, -0.050873, -0.050804, -0.050736, -0.050668,
-0.050596, -0.050533, -0.050465, -0.050398, -0.050331, -0.050264, -0.050198,
-0.050131, -0.050061, -0.049999, -0.049937, -0.049866, -0.049805, -0.049735,
-0.049670, -0.049605, -0.049537, -0.049475, -0.049410, -0.049344, -0.049281,
-0.049217, -0.049153, -0.049089, -0.049025, -0.048962, -0.048899, -0.048835,
-0.048767, -0.048709, -0.048646, -0.048584, -0.048523, -0.048459, -0.048397,
-0.048335, -0.048272, -0.048211, -0.048150, -0.048088, -0.048026, -0.047966,
-0.047905, -0.047843, -0.047787, -0.047723, -0.047663, -0.047602, -0.047544,
-0.047485, -0.047423, -0.047363, -0.047301, -0.047244, -0.047185, -0.047126,
-0.047067, -0.047008, -0.046950, -0.046891, -0.046835, -0.046775, -0.046712,
-0.046658, -0.046599, -0.047119, -0.046486, -0.046428, -0.046372, -0.046318,
-0.046257, -0.046200, -0.046146, -0.046087, -0.046031, -0.045974, -0.045918,
-0.045862, -0.045806, -0.045750, -0.045690, -0.045639, -0.045584, -0.045529,
-0.045471, -0.045419, -0.045363, -0.045309, -0.045255, -0.045200, -0.045146,
-0.045091, -0.045036, -0.044983, -0.044929, -0.044876, -0.044820, -0.044769,
-0.044715, -0.044662, -0.044598, -0.044556, -0.044503, -0.044449, -0.044400,
-0.044345, -0.044293, -0.044241, -0.044188, -0.044137, -0.044085, -0.044033,
-0.043978, -0.043930, -0.043878, -0.043827, -0.043775, -0.043725, -0.043674,
-0.043623, -0.043573, -0.043522, -0.043471, -0.043421, -0.043371, -0.043321,
-0.043271, -0.043220, -0.043170, -0.043121, -0.043072, -0.043022, -0.042974,
-0.042923, -0.042874, -0.042825, -0.042779, -0.042729, -0.042679, -0.042630,
-0.042584, -0.042533, -0.042485, -0.042437, -0.042388, -0.042341, -0.042293,
-0.042245, -0.042197, -0.042150, -0.042103, -0.042055, -0.042005, -0.041960,
-0.041915, -0.041867, -0.041822, -0.041773, -0.041727, -0.041680, -0.041632,
-0.041588, -0.041541, -0.041495, -0.041450, -0.041403, -0.041358, -0.041312,
-0.041262, -0.041221, -0.041175, -0.041130, -0.041087, -0.041040, -0.040995,
-0.040950, -0.040906, -0.040860, -0.040816, -0.040771, -0.040726, -0.040682,
-0.040638, -0.040594, -0.040552, -0.040506, -0.040461, -0.040422, -0.040371,
-0.040331, -0.040287, -0.040244, -0.040199, -0.040157, -0.040114, -0.040071,
-0.040029, -0.039985, -0.039942, -0.039900, -0.039855, -0.039815, -0.039772,
-0.039730, -0.039688, -0.039645, -0.039603, -0.039561, -0.039521, -0.039478,
-0.039436, -0.039394, -0.039350, -0.039311, -0.039270, -0.039228, -0.039188,
-0.039146, -0.039105, -0.039064, -0.039027, -0.038983, -0.038942, -0.038901,
-0.038861, -0.038820, -0.038780, -0.038740, -0.038700, -0.038659, -0.038619,
-0.038578, -0.038541, -0.038499, -0.038460, -0.038420, -0.038380, -0.038341,
-0.038302, -0.038262, -0.038221, -0.038182, -0.038145, -0.038106, -0.038068,
-0.038028, -0.037989, -0.037951, -0.037914, -0.037873, -0.037835, -0.037797,
-0.037758, -0.037720, -0.037682, -0.037643, -0.037606, -0.037568, -0.037530,
-0.037492, -0.037455, -0.037417, -0.037379, -0.037342, -0.037305, -0.037267,
-0.037230, -0.037193, -0.037155, -0.037118, -0.037081, -0.037044, -0.037008,
-0.036971, -0.036934, -0.036898, -0.036863, -0.036826, -0.036788, -0.036751,
-0.036717, -0.036678, -0.036643, -0.036607, -0.036570, -0.036535, -0.036499,
-0.036464, -0.036427, -0.036392, -0.036357, -0.036321, -0.036286, -0.036250,
-0.036215, -0.036180, -0.036144, -0.036110, -0.036075, -0.036040, -0.036003,
-0.035970, -0.035935, -0.035901, -0.035865, -0.035831, -0.035797, -0.035762,
-0.035730, -0.035694, -0.035660, -0.035625, -0.035591, -0.035557, -0.035520,
-0.035489, -0.035456, -0.035422, -0.035388, -0.035354, -0.035317, -0.035287,
-0.035254, -0.035220, -0.035189, -0.035154, -0.035121, -0.035088, -0.035054,
-0.035021, -0.034988, -0.034957, -0.034923, -0.034890, -0.034857, -0.034824,
-0.034793, -0.034759, -0.034727, -0.034694, -0.034661, -0.034630, -0.034597,
-0.034565, -0.034536, -0.034501, -0.034469, -0.034437, -0.034406, -0.034373,
-0.034341, -0.034310, -0.034277, -0.034246, -0.034215, -0.034183, -0.034152,
-0.034121, -0.034089, -0.034058, -0.034024, -0.033996, -0.033964, -0.033933,
-0.033902, -0.033872, -0.033841, -0.033810, -0.033782, -0.033749, -0.033718,
-0.033687, -0.033657, -0.033626, -0.033596, -0.033565, -0.033535, -0.033505,
-0.033474, -0.033444, -0.033415};


double gamma_terms(double Delta) {
    double delta = Delta;
    if(Delta < 1 || Delta > 5) {
        std::cerr << "WARNING in gamma_terms(), Delta was out of bounds." << std::endl;
        
        return 100000000000;    // this value should be large enough to indicate
                                // that something has gone wrong...
                                //
                                // I would return 0.0/0.0, but I don't remember right
                                // now what that does when -ffast-math is used.
    }

    if(Delta == 1) {
        return precomputed_gamma_terms[0];
    }
    if(Delta == 5) {
        return precomputed_gamma_terms[800];
    }
    else {
        // We use linear interpolation to estimate the value.
        //
        // I don't know exactly how accurate this is. Probably about
        // 4 digits?
        
        Delta = (Delta - 1.0) * 200;
        int a = (int)(std::floor(Delta));
        int b = (int)(std::ceil(Delta));

        if(a == b) {
            return precomputed_gamma_terms[a];
        }

        double x = precomputed_gamma_terms[a] * (Delta - a) + precomputed_gamma_terms[b] * (b - Delta);
        return x;
    }
}