lanczos_generator.cpp

Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work

//  (C) Copyright John Maddock 2006.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#include <boost/math/bindings/rr.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/math/tools/polynomial.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/numeric/ublas/operation.hpp>
#include <boost/cstdint.hpp>
#include <boost/bind.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/array.hpp>
#include <boost/test/included/test_exec_monitor.hpp>

//
// this is a sort of recursive include, since this file
// is used to create the contents of gamma.hpp.
// However, we will be using the "generic" version of gamma
// which doesn't use a lanczos approximation.
//
#include <boost/math/special_functions/gamma.hpp>

#include <iostream>
#include <iomanip>
#include <vector>

using boost::numeric::ublas::matrix;
using boost::numeric::ublas::vector;

struct lanczos_spot_data
{
   int N;
   double g;
   double err;
};

lanczos_spot_data sweet_spots[] = {
// a few selected spots from Pugh's thesis:
3, 2.03209, 6.4e-7,
4, 3.655180, 8.5e-8,
9, 8.406094, 6.9e-15,
10, 9.656578, 2.1e-16,
11, 10.900511, 6.1e-18,
12, 12.066012, 1.1e-19,
13, 13.144565, 5.2e-21,
14, 13.726821, 4.0e-22,
15, 14.977863, 1.2e-23,
20, 20.298892, 7.8e-31,
21, 21.508926, 2.1e-32,
22, 22.618910, 1.8e-34,
23, 23.118012, 5.2e-35,

// some more we've found, these are all the points where the first
// negleted term from the Lanczos series changes sign, there is one 
// point just above that point, and one just below:

3, 0.58613894134759903, 0.00036580426080686315,
3, 0.58613894879817963, 0.00036580416320668711,
3, 1.249951496720314, 4.2305536798953929e-005,
3, 1.2499515041708946, 4.2305534015367749e-005,
3, 1.8384982049465179, 1.164778449938813e-005,
3, 1.8384982123970985, 1.1647786529958048e-005,
3, 2.3276706635951996, 5.3732582230468998e-006,
3, 2.3276706710457802, 5.3732600294516588e-006,
3, 2.6038404405117035, 7.4060771171390379e-007,
3, 2.6038404479622841, 7.406058543248471e-007,
4, 0.62704569846391678, 0.00012800533557793425,
4, 0.62704570591449738, 0.00012800531551557324,
4, 1.3344274088740349, 9.0166140418663663e-006,
4, 1.3344274163246155, 9.0166116338646694e-006,
4, 1.9880444705486298, 1.4955557778725745e-006,
4, 1.9880444779992104, 1.4955557079994109e-006,
4, 2.5835030898451805, 4.5877676859600527e-007,
4, 2.5835030972957611, 4.5877667876325684e-007,
4, 3.1086832508444786, 2.1441725587707906e-007,
4, 3.1086832582950592, 2.1441723688775765e-007,
4, 3.5476611852645874, 1.0813760990515544e-007,
4, 3.547661192715168, 1.0813766234070327e-007,
4, 3.6492579728364944, 1.047642251124206e-007,
4, 3.649257980287075, 1.0476415796074403e-007,
5, 0.65435918420553207, 5.5274429026631278e-005,
5, 0.65435919165611267, 5.527441185903385e-005,
5, 1.389187254011631, 2.6307833416554052e-006,
5, 1.3891872614622116, 2.6307826816540255e-006,
5, 2.0801975876092911, 2.9014694764756466e-007,
5, 2.0801975950598717, 2.9014696142909125e-007,
5, 2.7272208333015442, 5.9176544536206063e-008,
5, 2.7272208407521248, 5.9176528019009275e-008,
5, 3.3269794136285782, 1.9379184595472579e-008,
5, 3.3269794210791588, 1.9379189057179851e-008,
5, 3.8722873777151108, 9.0123694745032944e-009,
5, 3.8722873851656914, 9.0123661007867372e-009,
5, 4.3516943231225014, 4.9816097608842378e-009,
5, 4.351694330573082, 4.9816087871575021e-009,
6, 0.67425807565450668, 2.7371903086326379e-005,
6, 0.67425808310508728, 2.7371893712684938e-005,
6, 1.4284561350941658, 9.4138428977455083e-007,
6, 1.4284561425447464, 9.4138417896969612e-007,
6, 2.1447814926505089, 7.4109049995920916e-008,
6, 2.1447815001010895, 7.4109044064461768e-008,
6, 2.8244904428720474, 1.0708333024058614e-008,
6, 2.824490450322628, 1.0708333661155651e-008,
6, 3.4669468700885773, 2.4940116660002578e-009,
6, 3.4669468775391579, 2.4940118560592828e-009,
6, 4.0696377158164978, 8.5267622527507512e-010,
6, 4.0696377232670784, 8.5267624048985822e-010,
6, 4.6278460025787354, 3.9323405231792876e-010,
6, 4.6278460100293159, 3.9323393550021793e-010,
6, 5.1345816180109978, 2.1943462473419721e-010,
6, 5.1345816254615784, 2.1943458105218707e-010,
6, 5.5917193591594696, 1.242916896389048e-010,
6, 5.5917193666100502, 1.2429175679293363e-010,
7, 0.68959406018257141, 1.4934177419305624e-005,
7, 0.68959406763315201, 1.4934171696393048e-005,
7, 1.4584139510989189, 3.8891482410396904e-007,
7, 1.4584139585494995, 3.8891475436265016e-007,
7, 2.1934122517704964, 2.2977989848570384e-008,
7, 2.193412259221077, 2.2977982819591052e-008,
7, 2.8964230641722679, 2.4742239281355561e-009,
7, 2.8964230716228485, 2.4742234478616864e-009,
7, 3.5678729563951492, 4.2816138965285379e-010,
7, 3.5678729638457298, 4.2816128067858218e-010,
7, 4.2069913893938065, 1.0931003472338519e-010,
7, 4.2069913968443871, 1.0931004659777642e-010,
7, 4.8118086308240891, 3.8518213120095777e-011,
7, 4.8118086382746696, 3.8518212172324096e-011,
7, 5.3790060579776764, 1.7613602961109952e-011,
7, 5.379006065428257, 1.7613603275796266e-011,
7, 5.9039439857006073, 9.7098281187573355e-012,
7, 5.9039439931511879, 9.7098271522398387e-012,
7, 6.3838147521018982, 5.7878801570305524e-012,
7, 6.3838147595524788, 5.7878771973838702e-012,
7, 6.7833591029047966, 2.829092368392454e-012,
7, 6.7833591103553772, 2.8290901754246269e-012,
8, 0.70188850909471512, 8.7610752024780369e-006,
8, 0.70188851654529572, 8.7610704217796693e-006,
8, 1.4822541922330856, 1.7881701068328982e-007,
8, 1.4822541996836662, 1.7881694190606741e-007,
8, 2.231784276664257, 8.2272386690167009e-009,
8, 2.2317842841148376, 8.2272367814664617e-009,
8, 2.9525627419352531, 6.8591536474980462e-010,
8, 2.9525627493858337, 6.8591515664364383e-010,
8, 3.6455128863453865, 9.1549235855884194e-011,
8, 3.6455128937959671, 9.1549217088498416e-011,
8, 4.3106606528162956, 1.8015466799942756e-011,
8, 4.3106606602668762, 1.8015467001952356e-011,
8, 4.9472371935844421, 4.918149362493385e-012,
8, 4.9472372010350227, 4.9181495178774722e-012,
8, 5.553666315972805, 1.7720036810370511e-012,
8, 5.5536663234233856, 1.7720035779693502e-012,
8, 6.1274929493665695, 8.0403304333886761e-013,
8, 6.1274929568171501, 8.0403274526493467e-013,
8, 6.6654317826032639, 4.3551924163930176e-013,
8, 6.6654317900538445, 4.3551918580885179e-013,
8, 7.1646013185381889, 2.6142608803934048e-013,
8, 7.1646013259887695, 2.6142610866873453e-013,
8, 7.6312204822897911, 1.5767278184484503e-013,
8, 7.6312204897403717, 1.576728292094547e-013,
8, 7.8779716044664383, 6.0579161144149095e-014,
8, 7.8779716119170189, 6.0579145495868576e-014,
9, 0.7120363786816597, 5.4383680639900953e-006,
9, 0.7120363861322403, 5.4383647430909934e-006,
9, 1.5018193200230598, 8.9331552269793478e-008,
9, 1.5018193274736404, 8.9331529217486569e-008,
9, 2.2630845010280609, 3.2922940247769628e-009,
9, 2.2630845084786415, 3.2922936681951672e-009,
9, 2.9980199560523033, 2.1887825488935961e-010,
9, 2.9980199635028839, 2.1887818854849753e-010,
9, 3.7078000977635384, 2.3213138957857667e-011,
9, 3.707800105214119, 2.3213132917173771e-011,
9, 4.3928612694144249, 3.6226509291513855e-012,
9, 4.3928612768650055, 3.6226507116549233e-012,
9, 5.0530400574207306, 7.8472620106742312e-013,
9, 5.0530400648713112, 7.8472599177745781e-013,
9, 5.6876183524727821, 2.2542091218064693e-013,
9, 5.6876183599233627, 2.2542085472082343e-013,
9, 6.2953110262751579, 8.261026286619692e-014,
9, 6.2953110337257385, 8.2610267677539022e-014,
9, 6.8742315545678139, 3.7226229880661886e-014,
9, 6.8742315620183945, 3.722621661830273e-014,
9, 7.4219388365745544, 1.9804198141916339e-014,
9, 7.421938844025135, 1.9804197841080089e-014,
9, 7.936030775308609, 1.1800366216622659e-014,
9, 7.9360307827591896, 1.1800362381207757e-014,
9, 8.4175376817584038, 7.341772297716249e-015,
9, 8.4175376892089844, 7.3417744970482286e-015,
10, 0.72060277312994003, 3.5324438420494327e-006,
10, 0.72060278058052063, 3.5324425213780505e-006,
10, 1.5182576552033424, 4.7690979892410133e-008,
10, 1.518257662653923, 4.7690968887170632e-008,
10, 2.2892609387636185, 1.4396365431590603e-009,
10, 2.2892609462141991, 1.4396363283470758e-009,
10, 3.0358340814709663, 7.8118059457005561e-011,
10, 3.0358340889215469, 7.8118048721726738e-011,
10, 3.7592824399471283, 6.7417155809106958e-012,
10, 3.7592824473977089, 6.7417147580480395e-012,
10, 4.4602697342634201, 8.5434024476488055e-013,
10, 4.4602697417140007, 8.5434000006256964e-013,
10, 5.1389683187007904, 1.5015237157428729e-013,
10, 5.138968326151371, 1.5015232503714518e-013,
10, 5.795127771794796, 3.5035278065891386e-014,
10, 5.7951277792453766, 3.5035266685032581e-014,
10, 6.4280949085950851, 1.0474518584578706e-014,
10, 6.4280949160456657, 1.047451900091486e-014,
10, 7.0368044599890709, 3.8902615249132166e-015,
10, 7.0368044674396515, 3.8902605386615012e-015,
10, 7.6197631657123566, 1.742442940023884e-015,
10, 7.6197631731629372, 1.7424426824540998e-015,
10, 8.1750892177224159, 9.1144555287546294e-016,
10, 8.1750892251729965, 9.1144560389833963e-016,
10, 8.7008367404341698, 5.3541847847332351e-016,
10, 8.7008367478847504, 5.3541844723127514e-016,
10, 9.1965094059705734, 3.3560389199894529e-016,
10, 9.196509413421154, 3.3560370421140035e-016,
10, 9.6678015515208244, 2.1059571242602905e-016,
10, 9.667801558971405, 2.1059579862533338e-016,
11, 0.72796426713466644, 2.3816013324880167e-006,
11, 0.72796427458524704, 2.3816003939892356e-006,
11, 1.5323273167014122, 2.6884539057733283e-008,
11, 1.5323273241519928, 2.6884532967611841e-008,
11, 2.3115825653076172, 6.7682186491374862e-010,
11, 2.3115825727581978, 6.768217495558193e-010,
11, 3.0679489448666573, 3.0544688315708309e-011,
11, 3.0679489523172379, 3.054468452617125e-011,
11, 3.8027990683913231, 2.186949692679853e-012,
11, 3.8027990758419037, 2.1869496552307528e-012,
11, 4.5169280469417572, 2.2947351040888336e-013,
11, 4.5169280543923378, 2.2947348773371414e-013,
11, 5.2107065990567207, 3.3353640359739479e-014,
11, 5.2107066065073013, 3.335363154144172e-014,
11, 5.8841593414545059, 6.4351159890193404e-015,
11, 5.8841593489050865, 6.435115249418266e-015,
11, 6.5370001494884491, 1.5930548620019077e-015,
11, 6.5370001569390297, 1.59305482454352e-015,
11, 7.1686412021517754, 4.9183710858836462e-016,
11, 7.168641209602356, 4.9183713116224844e-016,
11, 7.7781866490840912, 1.8464804238726119e-016,
11, 7.7781866565346718, 1.8464801709489452e-016,
11, 8.3644260391592979, 8.2263320775918765e-017,
11, 8.3644260466098785, 8.2263297787354186e-017,
11, 8.9258662387728691, 4.2375346461043962e-017,
11, 8.9258662462234497, 4.2375329267920249e-017,
11, 9.4609249979257584, 2.4482616150935944e-017,
11, 9.460925005376339, 2.4482619419420328e-017,
11, 9.968697115778923, 1.5254006021629585e-017,
11, 9.9686971232295036, 1.5254009931381535e-017,
11, 10.451380550861359, 9.7234265908062383e-018,
11, 10.451380558311939, 9.7234241296957717e-018,
11, 10.909898526966572, 6.3147967493977272e-018,
11, 10.909898534417152, 6.3147961332742783e-018,
12, 0.73438271880149841, 1.6566052109275983e-006,
12, 0.73438272625207901, 1.6566045155548952e-006,
12, 1.5445517003536224, 1.5860666784152459e-008,
12, 1.544551707804203, 1.5860658867525684e-008,
12, 2.33091701567173, 3.380059232629492e-010,
12, 2.3309170231223106, 3.3800585031437757e-010,
12, 3.095676377415657, 1.2884951580384992e-011,
12, 3.0956763848662376, 1.2884947453406877e-011,
12, 3.8402339816093445, 7.7768150022424858e-013,
12, 3.8402339890599251, 7.776814490345326e-013,
12, 4.5654622465372086, 6.8670097741098332e-014,
12, 4.5654622539877892, 6.8670083429297199e-014,
12, 5.2718541026115417, 8.38897478481124e-015,
12, 5.2718541100621223, 8.3889742372884196e-015,
12, 5.9596048966050148, 1.3595098759195538e-015,
12, 5.9596049040555954, 1.3595098086325987e-015,
12, 6.6286550611257553, 2.827608125932321e-016,
12, 6.6286550685763359, 2.8276074709416868e-016,
12, 7.2787091657519341, 7.3455425966950088e-017,