#include <boost/cstdint.hpp>
#include <boost/bind.hpp>
#include <boost/math/special_functions/log1p.hpp>
+#include <boost/math/special_functions/relative_difference.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/array.hpp>
+#include <boost/random.hpp>
+#include <boost/chrono.hpp>
#include "mp_t.hpp"
//
+// Define MP_TYPE
+// to a specific precision multiprecision type if you want to
+// create a lanczos approximation for that type.
+//
+
+#define MP_TYPE boost::multiprecision::number<boost::multiprecision::cpp_bin_float<100> >
+
+//
// 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
48, 36.410855412483215, 1.5004189889736636e-65,
48, 37.011237367987633, 6.8526723104762059e-66,
48, 37.011237375438213, 6.8526715926507402e-66,
+49, 32.804697265625, 1.43398975614025e-77,
+49, 32.80474609375, 1.44673437055229e-77,
+49, 33.440634765625, 4.27570575602756e-78,
+49, 33.44068359375, 4.31374583006685e-78,
+49, 34.071689453125, 1.3311629562209e-78,
+49, 34.07173828125, 1.34302740537616e-78,
+49, 34.6978125, 4.3180939436597e-79,
+49, 34.697861328125, 4.35671988075901e-79,
+49, 35.31900390625, 1.47025173244693e-79,
+49, 35.319052734375, 1.48336146723368e-79,
+49, 35.9351171875, 5.15604524622028e-80,
+49, 35.935166015625, 5.20239853250623e-80,
+49, 36.546201171875, 1.89501379810394e-80,
+49, 36.54625, 1.9120635932297e-80,
+49, 37.152109375, 7.16858816674083e-81,
+49, 37.152158203125, 7.23377831202289e-81,
+49, 37.752890625, 2.84316453204977e-81,
+49, 37.752939453125, 2.86903747838692e-81,
+49, 38.3483984375, 1.16185980624431e-81,
+49, 38.348447265625, 1.17251002885751e-81,
+49, 38.938583984375, 4.89658368061457e-82,
+49, 38.9386328125, 4.94200005973042e-82,
+49, 39.5233984375, 2.13065629695396e-82,
+49, 39.523447265625, 2.15069572328143e-82,
+49, 40.102841796875, 9.67604255181201e-83,
+49, 40.102890625, 9.76740702977543e-83,
+49, 40.676767578125, 4.50630247183498e-83,
+49, 40.67681640625, 4.54930080744076e-83,
+49, 41.245126953125, 2.15387681720274e-83,
+49, 41.24517578125, 2.17473706057203e-83,
+49, 41.80787109375, 1.05682121783075e-83,
+49, 41.807919921875, 1.06723874854161e-83,
+49, 42.365, 5.3732505406968e-84,
+49, 42.365048828125, 5.42671775522331e-84,
+49, 42.916416015625, 2.79825365457182e-84,
+49, 42.91646484375, 2.82641567518989e-84,
+49, 43.4620703125, 1.48702015592533e-84,
+49, 43.462119140625, 1.50221776839907e-84,
+49, 44.0019140625, 8.01177208266458e-85,
+49, 44.001962890625, 8.09565323702789e-85,
+49, 44.53599609375, 4.42566611189109e-85,
+49, 44.536044921875, 4.47291966331024e-85,
+49, 45.06431640625, 2.49957129718041e-85,
+49, 45.064365234375, 2.52668601869568e-85,
+49, 45.586875, 1.43196123451618e-85,
+49, 45.586923828125, 1.4477759848853e-85,
+49, 46.10376953125, 8.38746806134962e-86,
+49, 46.103818359375, 8.48100243483236e-86,
+49, 46.615, 4.91772041740201e-86,
+49, 46.615048828125, 4.97369333258043e-86,
+49, 47.120712890625, 2.8986790595583e-86,
+49, 47.12076171875, 2.93248475871069e-86,
+49, 47.621103515625, 1.73183181437907e-86,
+49, 47.62115234375, 1.75238561661062e-86,
+49, 48.116318359375, 1.03658031298704e-86,
+49, 48.1163671875, 1.04913081745512e-86,
+49, 48.6066015625, 6.24079697368347e-87,
+49, 48.606650390625, 6.31757720373193e-87,
+50, 33.54638671875, 1.92761599350307e-79,
+50, 33.546435546875, 1.94516650364829e-79,
+50, 34.181982421875, 5.76192865620046e-80,
+50, 34.18203125, 5.81436002400983e-80,
+50, 34.81279296875, 1.79764857133798e-80,
+50, 34.812841796875, 1.81400022919546e-80,
+50, 35.438720703125, 5.78653632499813e-81,
+50, 35.43876953125, 5.83972710207297e-81,
+50, 36.05986328125, 1.97160801805539e-81,
+50, 36.059912109375, 1.98962723354967e-81,
+50, 36.676025390625, 6.90862452044191e-82,
+50, 36.67607421875, 6.97216511712785e-82,
+50, 37.28720703125, 2.51020639219632e-82,
+50, 37.287255859375, 2.533501052945e-82,
+50, 37.893408203125, 9.54172211017555e-83,
+50, 37.89345703125, 9.63040270153392e-83,
+50, 38.494482421875, 3.7244684010854e-83,
+50, 38.49453125, 3.75949927313683e-83,
+50, 39.090478515625, 1.52160071659105e-83,
+50, 39.09052734375, 1.53593965559733e-83,
+50, 39.68125, 6.38933786856091e-84,
+50, 39.681298828125, 6.45010684683321e-84,
+50, 40.266748046875, 2.75949080957442e-84,
+50, 40.266796875, 2.78612547118161e-84,
+50, 40.84697265625, 1.23878873498316e-84,
+50, 40.847021484375, 1.25084568494597e-84,
+50, 41.421826171875, 5.73021424860067e-85,
+50, 41.421875, 5.78652397826679e-85,
+50, 41.991259765625, 2.7328651665849e-85,
+50, 41.99130859375, 2.75996267833758e-85,
+50, 42.55517578125, 1.3306120670112e-85,
+50, 42.555224609375, 1.34403199525941e-85,
+50, 43.11357421875, 6.67538976472161e-86,
+50, 43.113623046875, 6.74368293994306e-86,
+50, 43.66640625, 3.44457915579056e-86,
+50, 43.666455078125, 3.48023818015529e-86,
+50, 44.213623046875, 1.82191303634907e-86,
+50, 44.213671875, 1.84098742019443e-86,
+50, 44.75517578125, 9.8204850440831e-87,
+50, 44.755224609375, 9.9248385898419e-87,
+50, 45.291015625, 5.34608946266382e-87,
+50, 45.291064453125, 5.40437863458585e-87,
+50, 45.821240234375, 2.99827983616452e-87,
+50, 45.8212890625, 3.03145075570816e-87,
+50, 46.34580078125, 1.70330388065361e-87,
+50, 46.345849609375, 1.7224993598544e-87,
+50, 46.864794921875, 9.90800075053591e-88,
+50, 46.86484375, 1.00207004805813e-87,
+50, 47.378173828125, 5.7288229698866e-88,
+50, 47.37822265625, 5.79582150586751e-88,
+50, 47.886181640625, 3.39976963258209e-88,
+50, 47.88623046875, 3.43999283382834e-88,
+50, 48.388818359375, 2.00196389521248e-88,
+50, 48.3888671875, 2.02630021564974e-88,
+50, 48.886376953125, 1.20465247273191e-88,
+50, 48.88642578125, 1.21945216999773e-88,
+50, 49.378955078125, 7.21301501171278e-89,
+50, 49.37900390625, 7.30328914193871e-89,
+50, 49.866796875, 4.30834039532944e-89,
+50, 49.866845703125, 4.36344856692381e-89,
+51, 34.28802734375, 2.50653782215963e-81,
+51, 34.288076171875, 2.53000371457217e-81,
+51, 34.92328125, 7.49158638633263e-82,
+51, 34.923330078125, 7.5617476436076e-82,
+51, 35.55384765625, 2.33786811783389e-82,
+51, 35.553896484375, 2.35974650122535e-82,
+51, 36.179677734375, 7.59715649735025e-83,
+51, 36.1797265625, 7.66823998247397e-83,
+51, 36.800673828125, 2.54120219840682e-83,
+51, 36.80072265625, 2.5652464992018e-83,
+51, 37.416884765625, 8.90060057210517e-84,
+51, 37.41693359375, 8.98516095098549e-84,
+51, 38.02826171875, 3.26067023053424e-84,
+51, 38.028310546875, 3.29155983535777e-84,
+51, 38.634658203125, 1.22494879115424e-84,
+51, 38.63470703125, 1.23666166304672e-84,
+51, 39.23607421875, 4.76183281948702e-85,
+51, 39.236123046875, 4.80788045863205e-85,
+51, 39.832509765625, 1.93410326826023e-85,
+51, 39.83255859375, 1.95285023241101e-85,
+51, 40.423818359375, 8.05264364304082e-86,
+51, 40.4238671875, 8.13162262140464e-86,
+51, 41.01, 3.47176148694309e-86,
+51, 41.010048828125, 3.50614975628603e-86,
+51, 41.591005859375, 1.55093045558687e-86,
+51, 41.5910546875, 1.56638804232732e-86,
+51, 42.16673828125, 7.11162716693627e-87,
+51, 42.166787109375, 7.18328462123865e-87,
+51, 42.7371484375, 3.34775725369505e-87,
+51, 42.737197265625, 3.3819741026035e-87,
+51, 43.302236328125, 1.63427281479688e-87,
+51, 43.30228515625, 1.65107897871875e-87,
+51, 43.86185546875, 8.09778432115049e-88,
+51, 43.861904296875, 8.1826011369978e-88,
+51, 44.4160546875, 4.14978056128476e-88,
+51, 44.416103515625, 4.19369025316849e-88,
+51, 44.964736328125, 2.16908484161885e-88,
+51, 44.96478515625, 2.19237214768002e-88,
+51, 45.507900390625, 1.16244501700139e-88,
+51, 45.50794921875, 1.17507586155075e-88,
+51, 46.045498046875, 6.33568823869884e-89,
+51, 46.045546875, 6.4056389506585e-89,
+51, 46.577529296875, 3.50966350105409e-89,
+51, 46.577578125, 3.54914608970848e-89,
+51, 47.103994140625, 1.9666300290808e-89,
+51, 47.10404296875, 1.98929989871267e-89,
+51, 47.624990234375, 1.12975553299459e-89,
+51, 47.6250390625, 1.14296776906837e-89,
+51, 48.140517578125, 6.55256898902797e-90,
+51, 48.14056640625, 6.630573396527e-90,
+51, 48.650673828125, 3.84161435604585e-90,
+51, 48.65072265625, 3.88816278976081e-90,
+51, 49.155556640625, 2.26271437569233e-90,
+51, 49.15560546875, 2.29072792758847e-90,
+51, 49.655361328125, 1.35342362614108e-90,
+51, 49.65541015625, 1.37038547515428e-90,
+51, 50.150185546875, 8.04288108914844e-91,
+51, 50.150234375, 8.14599579489655e-91,
+51, 50.640322265625, 4.85018565373833e-91,
+51, 50.64037109375, 4.91297519323409e-91,
+52, 34.3901171875, 1.09784132477666e-82,
+52, 34.390166015625, 1.10845736512505e-82,
+52, 35.029619140625, 3.14965843034692e-83,
+52, 35.02966796875, 3.18011008403191e-83,
+52, 35.664580078125, 9.47858460950038e-84,
+52, 35.66462890625, 9.56968525920582e-84,
+52, 36.29490234375, 2.95312719534554e-84,
+52, 36.294951171875, 2.98153078194748e-84,
+52, 36.9205859375, 9.58907405270613e-85,
+52, 36.920634765625, 9.68126325492615e-85,
+52, 37.541533203125, 3.20598766092008e-85,
+52, 37.54158203125, 3.23711321193289e-85,
+52, 38.15779296875, 1.1229532977843e-85,
+52, 38.157841796875, 1.13387086000462e-85,
+52, 38.769267578125, 4.07481357929281e-86,
+52, 38.76931640625, 4.11456729015147e-86,
+52, 39.375908203125, 1.53004878535561e-86,
+52, 39.37595703125, 1.54506157346258e-86,
+52, 39.977666015625, 5.94046664964127e-87,
+52, 39.97771484375, 5.99920997458957e-87,
+52, 40.5744921875, 2.38377058782452e-87,
+52, 40.574541015625, 2.40756323337164e-87,
+52, 41.166337890625, 9.88492129468606e-88,
+52, 41.16638671875, 9.98457131300606e-88,
+52, 41.753154296875, 4.23624932947817e-88,
+52, 41.753203125, 4.27936286666186e-88,
+52, 42.334892578125, 1.87664263717547e-88,
+52, 42.33494140625, 1.89589075628466e-88,
+52, 42.91150390625, 8.59533063061209e-89,
+52, 42.911552734375, 8.6839076793046e-89,
+52, 43.482890625, 4.02833193451732e-89,
+52, 43.482939453125, 4.07030503687417e-89,
+52, 44.04900390625, 1.93062243793654e-89,
+52, 44.049052734375, 1.9510780127415e-89,
+52, 44.60984375, 9.55143090611608e-90,
+52, 44.609892578125, 9.6538161099766e-90,
+52, 45.1653125, 4.81625778814193e-90,
+52, 45.165361328125, 4.86882655159775e-90,
+52, 45.715458984375, 2.52110578486338e-90,
+52, 45.7155078125, 2.54874785641837e-90,
+52, 46.26013671875, 1.33336774124782e-90,
+52, 46.260185546875, 1.34823517438449e-90,
+52, 46.79939453125, 7.22810225870508e-91,
+52, 46.799443359375, 7.30975111869956e-91,
+52, 47.33318359375, 3.97344561651893e-91,
+52, 47.333232421875, 4.01915570075938e-91,
+52, 47.86150390625, 2.20743711754887e-91,
+52, 47.861552734375, 2.23347698762707e-91,
+52, 48.384453125, 1.25905367991677e-91,
+52, 48.384501953125, 1.27411715811031e-91,
+52, 48.90203125, 7.27905637066785e-92,
+52, 48.902080078125, 7.36737320296161e-92,
+52, 49.414287109375, 4.23381379773091e-92,
+52, 49.4143359375, 4.28618644566405e-92,
+52, 49.9213671875, 2.50298615147254e-92,
+52, 49.921416015625, 2.53432810782328e-92,
+52, 50.4233203125, 1.47126510667252e-92,
+52, 50.423369140625, 1.49015523321499e-92,
+52, 50.920390625, 8.77981976716807e-93,
+52, 50.920439453125, 8.89421391691748e-93,
+52, 51.41267578125, 5.18853729438254e-93,
+52, 51.4127734375, 5.32746953609488e-93,
+52, 51.90046875, 3.10020589622689e-93,
+52, 51.90056640625, 3.18462954145289e-93,
+53, 35.132099609375, 1.34165815167586e-84,
+53, 35.1321484375, 1.35500182297823e-84,
+53, 35.7712109375, 3.8586346438587e-85,
+53, 35.771259765625, 3.89696659256537e-85,
+53, 36.405830078125, 1.15515485184e-85,
+53, 36.40587890625, 1.16663065451713e-85,
+53, 37.03595703125, 3.62188787444457e-86,
+53, 37.036005859375, 3.65765415366427e-86,
+53, 37.6614453125, 1.16222508687862e-86,
+53, 37.661494140625, 1.17382378602205e-86,
+53, 38.282392578125, 3.9213181641527e-87,
+53, 38.28244140625, 3.96039919460395e-87,
+53, 38.898701171875, 1.37417257661721e-87,
+53, 38.89875, 1.38784484660129e-87,
+53, 39.5102734375, 4.94269999592506e-88,
+53, 39.510322265625, 4.99232702818206e-88,
+53, 40.117109375, 1.84017054708022e-88,
+53, 40.117158203125, 1.85883994034502e-88,
+53, 40.719208984375, 7.15588091788769e-89,
+53, 40.7192578125, 7.22859337586009e-89,
+53, 41.316474609375, 2.87513573445484e-89,
+53, 41.3165234375, 2.90443251849943e-89,
+53, 41.90880859375, 1.18074666272332e-89,
+53, 41.908857421875, 1.19294813450042e-89,
+53, 42.496259765625, 5.05747920130745e-90,
+53, 42.49630859375, 5.10993883553687e-90,
+53, 43.078681640625, 2.21264357218252e-90,
+53, 43.07873046875, 2.23591087175254e-90,
+53, 43.65607421875, 9.98793944641443e-91,
+53, 43.656123046875, 1.00942684217088e-90,
+53, 44.228388671875, 4.65040411340808e-91,
+53, 44.2284375, 4.7004152251227e-91,
+53, 44.795576171875, 2.23187504561532e-91,
+53, 44.795625, 2.25605716326726e-91,
+53, 45.3575390625, 1.09071767039123e-91,
+53, 45.357587890625, 1.10272568727892e-91,
+53, 45.91427734375, 5.47462949180705e-92,
+53, 45.914326171875, 5.53577914514396e-92,
+53, 46.4657421875, 2.81267854173599e-92,
+53, 46.465791015625, 2.84457112046934e-92,
+53, 47.01193359375, 1.48845783243062e-92,
+53, 47.011982421875, 1.50546780296832e-92,
+53, 47.55275390625, 7.96299241886377e-93,
+53, 47.552802734375, 8.05564452331045e-93,
+53, 48.088251953125, 4.36254557224574e-93,
+53, 48.08830078125, 4.41399477393221e-93,
+53, 48.61837890625, 2.41436230087065e-93,
+53, 48.618427734375, 2.44344138971554e-93,
+53, 49.143232421875, 1.3742563864423e-93,
+53, 49.14328125, 1.39095171064196e-93,
+53, 49.662763671875, 7.8627249720243e-94,
+53, 49.6628125, 7.95992783735172e-94,
+53, 50.1770703125, 4.55653379867663e-94,
+53, 50.177119140625, 4.61380444528666e-94,
+53, 50.686201171875, 2.64507670539115e-94,
+53, 50.68625, 2.6791564695207e-94,
+53, 51.1903515625, 1.56964252641228e-94,
+53, 51.190400390625, 1.59007794504541e-94,
+53, 51.6895703125, 9.2823621015156e-95,
+53, 51.68966796875, 9.5288381887029e-95,
+53, 52.18400390625, 5.4283757570851e-95,
+53, 52.1841015625, 5.5775793180482e-95,
+53, 52.6739453125, 3.22377149580101e-95,
+53, 52.67404296875, 3.31422857776726e-95,
+54, 35.87408203125, 1.60299124464084e-86,
+54, 35.874130859375, 1.61925498333934e-86,
+54, 36.512802734375, 4.60744049715477e-87,
+54, 36.5128515625, 4.6542262775045e-87,
+54, 37.14712890625, 1.38230307333573e-87,
+54, 37.147177734375, 1.3963165242666e-87,
+54, 37.776962890625, 4.2676811431326e-88,
+54, 37.77701171875, 4.31135557415735e-88,
+54, 38.402353515625, 1.37886617731339e-88,
+54, 38.40240234375, 1.39301204526773e-88,
+54, 39.023251953125, 4.64806876395431e-89,
+54, 39.02330078125, 4.69564494367985e-89,
+54, 39.639560546875, 1.61376384913413e-89,
+54, 39.639609375, 1.63036800302231e-89,
+54, 40.251279296875, 5.81713197366461e-90,
+54, 40.251328125, 5.87720283929007e-90,
+54, 40.858359375, 2.17315978880267e-90,
+54, 40.858408203125, 2.19566884825944e-90,
+54, 41.460703125, 8.31413664866219e-91,
+54, 41.460751953125, 8.40143100626546e-91,
+54, 42.058359375, 3.32234094377706e-91,
+54, 42.058408203125, 3.35733728990472e-91,
+54, 42.65123046875, 1.37108630576255e-91,
+54, 42.651279296875, 1.38557871529173e-91,
+54, 43.23921875, 5.77759110628069e-92,
+54, 43.239267578125, 5.8395345444883e-92,
+54, 43.822373046875, 2.53796499025021e-92,
+54, 43.822421875, 2.56525705223881e-92,
+54, 44.400546875, 1.13685971567258e-92,
+54, 44.400595703125, 1.14924635408225e-92,
+54, 44.973740234375, 5.24494903201998e-93,
+54, 44.9737890625, 5.30278922556907e-93,
+54, 45.541904296875, 2.48974278621437e-93,
+54, 45.541953125, 2.51750108068007e-93,
+54, 46.104990234375, 1.21416817580495e-93,
+54, 46.1050390625, 1.22784384566186e-93,
+54, 46.66294921875, 6.06863156354183e-94,
+54, 46.662998046875, 6.13771586650707e-94,
+54, 47.215732421875, 3.09782695797381e-94,
+54, 47.21578125, 3.13356585781024e-94,
+54, 47.763291015625, 1.60666638154662e-94,
+54, 47.76333984375, 1.62557509885327e-94,
+54, 48.305673828125, 8.60651173329294e-95,
+54, 48.30572265625, 8.70866176467966e-95,
+54, 48.842783203125, 4.66138641992004e-95,
+54, 48.84283203125, 4.71765875807128e-95,
+54, 49.37466796875, 2.58125106144261e-95,
+54, 49.374716796875, 2.61280647106825e-95,
+54, 49.901279296875, 1.43706251427965e-95,
+54, 49.901328125, 1.4550467561572e-95,
+54, 50.42271484375, 8.17008391727145e-96,
+54, 50.422763671875, 8.27405152155282e-96,
+54, 50.9390234375, 4.73629466134273e-96,
+54, 50.939072265625, 4.79714657796633e-96,
+54, 51.45015625, 2.70637622493452e-96,
+54, 51.45025390625, 2.77837012006172e-96,
+54, 51.95640625, 1.60392443466118e-96,
+54, 51.95650390625, 1.6468673251926e-96,
+54, 52.4577734375, 9.52380966519104e-97,
+54, 52.45787109375, 9.78162481883129e-97,
+54, 52.95435546875, 5.57472185760653e-97,
+54, 52.954453125, 5.73022000928701e-97,
+54, 53.4464453125, 3.32141601130766e-97,
+54, 53.44654296875, 3.41542726435434e-97,
+54, 53.934140625, 1.95568174192375e-97,
+54, 53.93423828125, 2.01256167035286e-97,
+55, 36.615966796875, 1.8332541133448e-88,
+55, 36.616015625, 1.85247560570567e-88,
+55, 37.254345703125, 5.30481775801044e-89,
+55, 37.25439453125, 5.36017474932385e-89,
+55, 37.88837890625, 1.59132201617429e-89,
+55, 37.888427734375, 1.60790995949161e-89,
+55, 38.518017578125, 4.92735433868831e-90,
+55, 38.51806640625, 4.97903083471312e-90,
+55, 39.14326171875, 1.58522865714114e-90,
+55, 39.143310546875, 1.60194883269875e-90,
+55, 39.764111328125, 5.33698675182668e-91,
+55, 39.76416015625, 5.39311857875617e-91,
+55, 40.38046875, 1.85508041787031e-91,
+55, 40.380517578125, 1.8746204836555e-91,
+55, 40.99228515625, 6.63996228149307e-92,
+55, 40.992333984375, 6.71043727342021e-92,
+55, 41.599560546875, 2.46798186824746e-92,
+55, 41.599609375, 2.49429084756772e-92,
+55, 42.202197265625, 9.40692186686793e-93,
+55, 42.20224609375, 9.50850937940109e-93,
+55, 42.800244140625, 3.75106220500044e-93,
+55, 42.80029296875, 3.79158825843364e-93,
+55, 43.3935546875, 1.52947380015209e-93,
+55, 43.393603515625, 1.54616717565206e-93,
+55, 43.98212890625, 6.43906065579312e-94,
+55, 43.982177734375, 6.50998764609593e-94,
+55, 44.56591796875, 2.79632632865198e-94,
+55, 44.565966796875, 2.82738157386603e-94,
+55, 45.144873046875, 1.25150437492926e-94,
+55, 45.144921875, 1.26550343381352e-94,
+55, 45.7189453125, 5.76629149742495e-95,
+55, 45.718994140625, 5.83119549719855e-95,
+55, 46.288037109375, 2.70055229254301e-95,
+55, 46.2880859375, 2.73147381215483e-95,
+55, 46.852197265625, 1.31270463114748e-95,
+55, 46.85224609375, 1.32782212771267e-95,
+55, 47.411279296875, 6.45560417517413e-96,
+55, 47.411328125, 6.53138347837362e-96,
+55, 47.96533203125, 3.27538902561878e-96,
+55, 47.965380859375, 3.31427965467404e-96,
+55, 48.5142578125, 1.68599305405381e-96,
+55, 48.514306640625, 1.70640462719072e-96,
+55, 49.05810546875, 8.95772389297722e-97,
+55, 49.058154296875, 9.06711241441999e-97,
+55, 49.59677734375, 4.8089458433976e-97,
+55, 49.596826171875, 4.86873031520792e-97,
+55, 50.130322265625, 2.64138783197092e-97,
+55, 50.13037109375, 2.67465396471993e-97,
+55, 50.658740234375, 1.47989802765765e-97,
+55, 50.6587890625, 1.49871358521579e-97,
+55, 51.181982421875, 8.287936462558e-98,
+55, 51.18203125, 8.39594586787481e-98,
+55, 51.7001953125, 4.76307652889875e-98,
+55, 51.70029296875, 4.88866793750702e-98,
+55, 52.21337890625, 2.75947508065571e-98,
+55, 52.2134765625, 2.83330856381719e-98,
+55, 52.72158203125, 1.59307203690564e-98,
+55, 52.7216796875, 1.63687980385948e-98,
+55, 53.225, 9.37441146592157e-99,
+55, 53.22509765625, 9.63618543580853e-99,
+55, 53.72373046875, 5.56087800966994e-99,
+55, 53.723828125, 5.71812073882941e-99,
+55, 54.21787109375, 3.26344577772922e-99,
+55, 54.21796875, 3.3582248315885e-99,
+55, 54.70771484375, 1.95673422065812e-99,
+55, 54.7078125, 2.01394065913309e-99,
+56, 37.3578515625, 2.04647041348259e-90,
+56, 37.357900390625, 2.06847980628938e-90,
+56, 37.99583984375, 5.88243837411135e-91,
+56, 37.995888671875, 5.94590494289821e-91,
+56, 38.629580078125, 1.76021316827907e-91,
+56, 38.62962890625, 1.77923863437215e-91,
+56, 39.2590234375, 5.45779852011259e-92,
+56, 39.259072265625, 5.51704239356607e-92,
+56, 39.884169921875, 1.76524851588723e-92,
+56, 39.88421875, 1.78439348136451e-92,
+56, 40.504970703125, 5.93342330120247e-93,
+56, 40.50501953125, 5.9975792008348e-93,
+56, 41.121328125, 2.04345614663194e-93,
+56, 41.121376953125, 2.06573687666038e-93,
+56, 41.733291015625, 7.34765755919972e-94,
+56, 41.73333984375, 7.42775943660008e-94,
+56, 42.34076171875, 2.72161843366724e-94,
+56, 42.340810546875, 2.75141110970169e-94,
+56, 42.94369140625, 1.03613177248368e-94,
+56, 42.943740234375, 1.04758637844824e-94,
+56, 43.542080078125, 4.09126499345346e-95,
+56, 43.54212890625, 4.13674598350119e-95,
+56, 44.135830078125, 1.65424896423927e-95,
+56, 44.13587890625, 1.67288546036148e-95,
+56, 44.724990234375, 6.99446261996147e-96,
+56, 44.7250390625, 7.07318147969286e-96,
+56, 45.3094140625, 3.01976769844332e-96,
+56, 45.309462890625, 3.05402248905396e-96,
+56, 45.889052734375, 1.32918182143927e-96,
+56, 45.8891015625, 1.34452414236527e-96,
+56, 46.463955078125, 6.09433185040725e-97,
+56, 46.46400390625, 6.16497223474698e-97,
+56, 47.033974609375, 2.83946834327271e-97,
+56, 47.0340234375, 2.87288014917526e-97,
+56, 47.599111328125, 1.35701697542241e-97,
+56, 47.59916015625, 1.37323222143338e-97,
+56, 48.159365234375, 6.71405638326024e-98,
+56, 48.1594140625, 6.79470565448062e-98,
+56, 48.714638671875, 3.38572542994683e-98,
+56, 48.7146875, 3.42679289110478e-98,
+56, 49.2648828125, 1.73113405012133e-98,
+56, 49.264931640625, 1.75251860795891e-98,
+56, 49.81009765625, 9.02288638273931e-99,
+56, 49.810146484375, 9.13660050992263e-99,
+56, 50.35033203125, 4.87338593256003e-99,
+56, 50.350380859375, 4.93504167152289e-99,
+56, 50.885439453125, 2.62897432578778e-99,
+56, 50.88548828125, 2.66302222150824e-99,
+56, 51.41546875, 1.43071360755971e-99,
+56, 51.41556640625, 1.46894783889895e-99,
+56, 51.94056640625, 8.12763956845104e-100,
+56, 51.9406640625, 8.34551910155042e-100,
+56, 52.4605859375, 4.58251489759803e-100,
+56, 52.46068359375, 4.70835817217817e-100,
+56, 52.97572265625, 2.66457725640481e-100,
+56, 52.9758203125, 2.73809355167495e-100,
+56, 53.48587890625, 1.52026667140348e-100,
+56, 53.4859765625, 1.5636427659068e-100,
+56, 53.99134765625, 8.99584141358291e-101,
+56, 53.9914453125, 9.25372909114439e-101,
+56, 54.49212890625, 5.32825924303324e-101,
+56, 54.4922265625, 5.48250901435612e-101,
+56, 54.9883203125, 3.11278188657358e-101,
+56, 54.98841796875, 3.20543677910749e-101,
+56, 55.48021484375, 1.86269305932801e-101,
+56, 55.4803125, 1.91846988478933e-101,
+57, 38.0996875, 2.20476439543612e-92,
+57, 38.099736328125, 2.22917432379422e-92,
+57, 38.737333984375, 6.34433356523382e-93,
+57, 38.7373828125, 6.41479399056003e-93,
+57, 39.370830078125, 1.91018618887989e-93,
+57, 39.37087890625, 1.93131173925967e-93,
+57, 40.000078125, 5.92216057314857e-94,
+57, 40.000126953125, 5.98791393861005e-94,
+57, 40.625078125, 1.90309216004904e-94,
+57, 40.625126953125, 1.92431835469877e-94,
+57, 41.245830078125, 6.38398110903664e-95,
+57, 41.24587890625, 6.45498354446764e-95,
+57, 41.862236328125, 2.20268568579443e-95,
+57, 41.86228515625, 2.22728247906098e-95,
+57, 42.474248046875, 7.79068157690964e-96,
+57, 42.474296875, 7.87885966415e-96,
+57, 43.0819140625, 2.88050550788237e-96,
+57, 43.081962890625, 2.91318272595493e-96,
+57, 43.68513671875, 1.09806112159568e-96,
+57, 43.685185546875, 1.11057160166479e-96,
+57, 44.2838671875, 4.30562812876992e-97,
+57, 44.283916015625, 4.35507143677268e-97,
+57, 44.87810546875, 1.7530232640245e-97,
+57, 44.878154296875, 1.7731758707715e-97,
+57, 45.46775390625, 7.31286187931292e-98,
+57, 45.467802734375, 7.39751767287751e-98,
+57, 46.052763671875, 3.11976304411238e-98,
+57, 46.0528125, 3.15638274069021e-98,
+57, 46.633134765625, 1.37464423777275e-98,
+57, 46.63318359375, 1.39093953120416e-98,
+57, 47.208818359375, 6.24536630012012e-99,
+57, 47.2088671875, 6.31989226309693e-99,
+57, 47.779716796875, 2.88489889714976e-99,
+57, 47.779765625, 2.91990162728214e-99,
+57, 48.34587890625, 1.38459867556409e-99,
+57, 48.345927734375, 1.40146037072031e-99,
+57, 48.907158203125, 6.71827295305273e-100,
+57, 48.90720703125, 6.80152139672612e-100,
+57, 49.463603515625, 3.36417686184697e-100,
+57, 49.46365234375, 3.40624440071944e-100,
+57, 50.0151171875, 1.70838155932608e-100,
+57, 50.015166015625, 1.73011766047011e-100,
+57, 50.561748046875, 8.96377195112949e-101,
+57, 50.561796875, 9.07844496707364e-101,
+57, 51.1033984375, 4.7573624982879e-101,
+57, 51.103447265625, 4.81906515869594e-101,
+57, 51.64001953125, 2.52428955803692e-101,
+57, 51.6401171875, 2.59192149716193e-101,
+57, 52.1717578125, 1.39187870966602e-101,
+57, 52.17185546875, 1.42956730803746e-101,
+57, 52.698515625, 7.7317500805067e-102,
+57, 52.69861328125, 7.94501950652444e-102,
+57, 53.220390625, 4.4082211286455e-102,
+57, 53.22048828125, 4.53055204196072e-102,
+57, 53.73728515625, 2.47396166798242e-102,
+57, 53.7373828125, 2.54499159923448e-102,
+57, 54.2494921875, 1.45907071692507e-102,
+57, 54.24958984375, 1.50072529856044e-102,
+57, 54.75681640625, 8.32179730748228e-103,
+57, 54.7569140625, 8.56822893627666e-103,
+57, 55.2596484375, 4.99089804791844e-103,
+57, 55.25974609375, 5.1376157827554e-103,
+57, 55.757890625, 2.93846786006917e-103,
+57, 55.75798828125, 3.02626005750667e-103,
+57, 56.25173828125, 1.71972275925121e-103,
+57, 56.2518359375, 1.77242201439855e-103,
+57, 56.74138671875, 1.00746867642176e-103,
+57, 56.741484375, 1.03914712466117e-103,
+58, 38.841474609375, 2.28962685408854e-94,
+58, 38.8415234375, 2.31583571384595e-94,
+58, 39.478828125, 6.65279878387051e-95,
+58, 39.478876953125, 6.72851011961114e-95,
+58, 40.11203125, 1.98889585849889e-95,
+58, 40.112080078125, 2.0116031211386e-95,
+58, 40.741083984375, 6.15547181773793e-96,
+58, 40.7411328125, 6.22611457465467e-96,
+58, 41.365986328125, 1.98535700187755e-96,
+58, 41.36603515625, 2.00813333128919e-96,
+58, 41.986640625, 6.56800883846017e-97,
+58, 41.986689453125, 6.64407409629831e-97,
+58, 42.603095703125, 2.27151430941761e-97,
+58, 42.60314453125, 2.29779962665845e-97,
+58, 43.21525390625, 8.08887967205788e-98,
+58, 43.215302734375, 8.18281575271904e-98,
+58, 43.82306640625, 2.95592121751972e-98,
+58, 43.823115234375, 2.99061314771485e-98,
+58, 44.42658203125, 1.13134888356557e-98,
+58, 44.426630859375, 1.14457524376905e-98,
+58, 45.025654296875, 4.41845285569886e-99,
+58, 45.025703125, 4.47048673389212e-99,
+58, 45.620283203125, 1.77714624291313e-99,
+58, 45.62033203125, 1.79825016216838e-99,
+58, 46.210419921875, 7.34443777835212e-100,
+58, 46.21046875, 7.43260970647061e-100,
+58, 46.796064453125, 3.14964704226534e-100,
+58, 46.79611328125, 3.18755872885081e-100,
+58, 47.377119140625, 1.38201826867972e-100,
+58, 47.37716796875, 1.39878290995527e-100,
+58, 47.95353515625, 6.19027397147059e-101,
+58, 47.953583984375, 6.26644970781979e-101,
+58, 48.5253125, 2.85831172869441e-101,
+58, 48.525361328125, 2.89384123773075e-101,
+58, 49.09240234375, 1.35697005525987e-101,
+58, 49.092451171875, 1.37396450788758e-101,
+58, 49.65470703125, 6.51873789254857e-102,
+58, 49.654755859375, 6.60202992320111e-102,
+58, 50.212275390625, 3.23653359270438e-102,
+58, 50.21232421875, 3.2783090214757e-102,
+58, 50.76505859375, 1.65267733466975e-102,
+58, 50.765107421875, 1.67409628228237e-102,
+58, 51.31291015625, 8.39865374766313e-103,
+58, 51.3130078125, 8.62296916999248e-103,
+58, 51.8559765625, 4.43898857454993e-103,
+58, 51.85607421875, 4.55876608676698e-103,
+58, 52.39416015625, 2.38284566979309e-103,
+58, 52.3942578125, 2.44798937063255e-103,
+58, 52.9274609375, 1.29896683738751e-103,
+58, 52.92755859375, 1.33500351596253e-103,
+58, 53.45587890625, 7.16661264811987e-104,
+58, 53.4559765625, 7.36908508006679e-104,
+58, 53.9794140625, 3.97139434881222e-104,
+58, 53.97951171875, 4.08676376045491e-104,
+58, 54.49826171875, 2.31430475778334e-104,
+58, 54.498359375, 2.38084236991499e-104,
+58, 55.0122265625, 1.31587661343037e-104,
+58, 55.01232421875, 1.354675556019e-104,
+58, 55.52150390625, 7.55721877752536e-105,
+58, 55.5216015625, 7.78550839762316e-105,
+58, 56.02619140625, 4.38540257957703e-105,
+58, 56.0262890625, 4.52070576457663e-105,
+58, 56.526484375, 2.63611015718921e-105,
+58, 56.52658203125, 2.71673162725835e-105,
+58, 57.02228515625, 1.52363825877265e-105,
+58, 57.0223828125, 1.57187635540744e-105,
+58, 57.513984375, 9.12003014967339e-106,
+58, 57.51408203125, 9.40921222802051e-106,
};
//
return result;
}
//
+// Stopwatch for measuring performance:
+//
+template <class Clock>
+struct stopwatch
+{
+ typedef typename Clock::duration duration;
+ stopwatch()
+ {
+ m_start = Clock::now();
+ }
+ duration elapsed()
+ {
+ return Clock::now() - m_start;
+ }
+ void reset()
+ {
+ m_start = Clock::now();
+ }
+
+private:
+ typename Clock::time_point m_start;
+};
+//
// Returns the factorials and half factorials:
//
template <class T>
static std::vector<std::vector<T> > data;
if(data.empty())
{
- T fact = 1;
+ mp_t fact = 1;
int k = 1;
std::vector<T> item;
do
{
data.push_back(item);
data.back().push_back(k-1);
- data.back().push_back(fact);
- data.back().push_back(log(fact));
+ data.back().push_back(static_cast<T>(fact));
+ data.back().push_back(static_cast<T>(log(fact)));
fact = fact * T(k++);
}while(k < 100);
fact = 0.5;
- T srpi = sqrt(boost::math::constants::pi<T>());
- T mul = 1.5;
+ mp_t srpi = sqrt(boost::math::constants::pi<mp_t>());
+ mp_t mul = 1.5;
do{
data.push_back(item);
- data.back().push_back(mul-1);
- data.back().push_back(fact*srpi);
- data.back().push_back(log(fact*srpi));
+ data.back().push_back(static_cast<T>(mul-1));
+ data.back().push_back(static_cast<T>(fact*srpi));
+ data.back().push_back(static_cast<T>(log(fact*srpi)));
fact *= mul;
mul += 1;
}while(mul < 100);
return get_test_data_near_x(T(2));
}
+template <class T>
+std::vector<std::vector<T> > const& get_random_test_data()
+{
+ static std::vector<std::vector<T> > data;
+ if (data.empty())
+ {
+ boost::random::mt19937_64 gen;
+ boost::random::uniform_real_distribution<T> dist(1e-30, 100);
+ for (unsigned i = 0; i < 100; ++i)
+ {
+ T value = dist(gen);
+ T result = static_cast<T>(boost::math::tgamma(mp_t(value) + 1));
+ data.push_back(std::vector<T>({value, result}));
+ }
+ }
+ return data;
+}
+
//
// Converts Lanczos approximation into rational form via
// polynomial arithmetic:
// half factorials, returns the max error found:
//
template <class T, class R>
-T get_max_error(const lanczos_info<T>& dat, R)
+T get_max_error(const lanczos_info<T>& dat, R, double* time = nullptr)
{
- T max_error = 0;
- std::vector<std::vector<T> > const & tests = get_test_data<T>();
+ R max_error = 0;
+ std::vector<std::vector<R> > const & tests1 = get_test_data<R>();
+ std::vector<std::vector<R> > const & tests2 = get_random_test_data<R>();
+
+ std::vector<std::vector<R> > tests(tests1.begin(), tests1.end());
+ tests.insert(tests.end(), tests2.begin(), tests2.end());
+
lanczos_rational<R> rational(dat);
+ stopwatch<boost::chrono::high_resolution_clock> w;
+
for(unsigned i = 0; i < tests.size(); ++i)
{
- R input = boost::math::tools::real_cast<R>(tests[i][0]);
+ R input = tests[i][0];
+ R expected = tests[i][1];
if(std::numeric_limits<R>::is_specialized && (boost::math::tools::real_cast<R>(tests[i][1]) > (std::numeric_limits<R>::max)()))
continue;
R gamr = rational.factorial(input);
if(std::numeric_limits<R>::is_specialized && (gamr > (std::numeric_limits<R>::max)()))
continue;
- T gam = gamr;
- T exp = tests[i][1];
- T err = relative_error(gam, exp);
+ R err = boost::math::relative_difference(gamr, expected);
if(err > max_error)
max_error = err;
}
+ if (time)
+ * time = boost::chrono::duration_cast<boost::chrono::duration<double>>(w.elapsed()).count();
+
return max_error;
}
+
+template <class T>
+T get_max_error_sterling(double* time = nullptr)
+{
+ std::vector<std::vector<T> > const& tests1 = get_test_data<T>();
+ std::vector<std::vector<T> > const& tests2 = get_random_test_data<T>();
+
+ std::vector<std::vector<T> > tests(tests1.begin(), tests1.end());
+ tests.insert(tests.end(), tests2.begin(), tests2.end());
+
+ T max_err = 0;
+
+ stopwatch<boost::chrono::high_resolution_clock> w;
+
+ for (unsigned i = 0; i < tests.size(); ++i)
+ {
+ T loc = tests[i][0] + 1;
+ T found = boost::math::detail::gamma_imp(loc, boost::math::policies::policy<>(), boost::math::lanczos::undefined_lanczos());
+ T expected = tests[i][1];
+ T err = boost::math::relative_difference(expected, found);
+ if (err > max_err)
+ max_err = err;
+ }
+
+ if(time)
+ *time = boost::chrono::duration_cast<boost::chrono::duration<double>>(w.elapsed()).count();
+
+ return max_err;
+}
+
//
// This is what prints the "best" approximation out as code:
//
template <class T>
-void print_code(const lanczos_info<T>& l, const char* name)
+void print_code(const lanczos_info<T>& l, const char* name, int precision = std::numeric_limits<T>::max_digits10, int precision2 = std::numeric_limits<T>::digits)
{
- std::cout << std::scientific << std::setprecision(std::numeric_limits<T>::digits10 + 3);
+ std::cout << std::scientific << std::setprecision(precision);
using namespace std;
lanczos_info<T> l2(l);
T factor = exp(l.r);
std::cout << name;
std::cout << " precision arithmetic) " << l.err <<
"\n// Generated with compiler: " << BOOST_COMPILER << " on " << BOOST_PLATFORM << " at " << __DATE__ << "\n"
+ "// Type precision was " << precision2 << " bits or " << precision << " max_digits10\n"
"//\n"
- "struct lanczos" << l.n << name << "\n"
+ "struct lanczos" << l.n << name << " : public mpl::int_<" << precision2 << ">\n"
"{\n"
" template <class T>\n"
" static T lanczos_sum(const T& z)\n"
for(unsigned i = 0; i < rat.num.size(); ++i)
{
- std::cout << " static_cast<T>(" << rat.num[i] << "L)";
+ if(precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " static_cast<T>(" << rat.num[i] << "L)";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << rat.num[i] << "))";
if(i != rat.num.size() - 1)
std::cout << ",";
std::cout << "\n";
" static const " << denom_type << " denom[" << rat.denom.size() << "] = {\n";
for(unsigned i = 0; i < rat.denom.size(); ++i)
{
- std::cout << " " << cast_type << "(" << rat.denom[i] << suffix_type << ")";
+ if(precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " " << cast_type << "(" << rat.denom[i] << suffix_type << ")";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << rat.denom[i] << "))";
if(i != rat.denom.size() - 1)
std::cout << ",";
std::cout << "\n";
for(unsigned i = 0; i < rat.num.size(); ++i)
{
- std::cout << " static_cast<T>(" << (rat.num[i]/factor) << "L)";
+ if (precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " static_cast<T>(" << (rat.num[i]/factor) << "L)";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << (rat.num[i]/factor) << "))";
if(i != rat.num.size() - 1)
std::cout << ",";
std::cout << "\n";
" static const " << denom_type << " denom[" << rat.denom.size() << "] = {\n";
for(unsigned i = 0; i < rat.denom.size(); ++i)
{
- std::cout << " " << cast_type << "(" << rat.denom[i] << suffix_type << ")";
+ if (precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " " << cast_type << "(" << rat.denom[i] << suffix_type << ")";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << rat.denom[i] << "))";
if(i != rat.denom.size() - 1)
std::cout << ",";
std::cout << "\n";
for(int i = 1; i < l.n; ++i)
{
- std::cout << " static_cast<T>(" << l.c[i]*factor << "L)";
+ if (precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " static_cast<T>(" << l.c[i]*factor << "L)";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << l.c[i] * factor << "))";
if(i != l.n - 1)
std::cout << ",";
std::cout << "\n";
for(int i = 1; i < l.n; ++i)
{
- std::cout << " static_cast<T>(" << l.c[i]/factor << "L),\n";
+ if (precision <= std::numeric_limits<long double>::digits10)
+ std::cout << " static_cast<T>(" << l.c[i]/factor << "L),\n";
+ else
+ std::cout << " static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, " << precision2 << ", " << l.c[i]/factor << ")),\n";
}
std::cout <<
" };\n"
//
void print_test_values(const std::vector<std::vector<mp_t> >& v, const char* name, int offset = 1)
{
- std::cout << "#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))\n";
+ std::cout << std::setprecision(110);
std::cout <<
" static const boost::array<boost::array<T, 3>, " << v.size() << "> " << name << " = {\n";
for(unsigned i = 0; i < v.size(); ++i)
{
std::cout << " SC_(" << (v[i][0] + offset) << "), SC_(" << v[i][1] << "), SC_(" << v[i][2] << "),\n";
}
- std::cout << " };\n#undef SC_\n\n";
+ std::cout << " };\n\n";
}
//
// Get the error for a specific approximation, and print out it's code:
{
using namespace std;
- lanczos_info<mp_t> best;
- best.err = 100; // best case had better be better than this!
- for(int i = 0; i < sizeof(sweet_spots)/sizeof(sweet_spots[0]); ++i)
+ double exec_time;
+
+ std::cout << "Enter value for N, or 0 to scan for best approximation: ";
+ int N;
+ std::cin >> N;
+
+ if (N == 0)
{
- if((sweet_spots[i].err < eps*10) && (sweet_spots[i].N < max_scan))
+ T sterling_err = get_max_error_sterling<T>(&exec_time);
+
+ std::cout << "Max error from generic Sterling approximation was: " << static_cast<int>(sterling_err / eps) << "eps (in " << (int)(exec_time * 1000) << "ms)" << std::endl;
+
+ lanczos_info<mp_t> best;
+ best.err = 100; // best case had better be better than this!
+
+ std::cout << std::setw(20) << std::right << "N" << std::setw(20) << std::right << "g" << std::setw(20) << std::right << "eps" << std::setw(20) << std::right << "time (ms)\n";
+
+ for (int i = 0; i < sizeof(sweet_spots) / sizeof(sweet_spots[0]); ++i)
{
- lanczos_info<mp_t> info = generate_lanczos(sweet_spots[i].N, mp_t(sweet_spots[i].g));
- mp_t err = get_max_error(info, eps);
- if(err/eps < 1000)
- {
- std::cout << sweet_spots[i].N << " " << sweet_spots[i].g << " " << err/eps << std::endl;
- }
- if(err < best.err)
+ if ((sweet_spots[i].err < eps * 10) && (sweet_spots[i].N < max_scan))
{
- best = info;
- best.err = err;
+ lanczos_info<mp_t> info = generate_lanczos(sweet_spots[i].N, mp_t(sweet_spots[i].g));
+ mp_t err = get_max_error(info, eps, &exec_time);
+ if (err / eps < 1000)
+ {
+ std::cout << std::setprecision(14) << std::fixed << std::setw(20) << std::right << sweet_spots[i].N
+ << std::setw(20) << std::right << sweet_spots[i].g << std::setw(20) << std::right << static_cast<int>(err / eps)
+ << std::setw(20) << std::right << (int)(exec_time * 1000) << std::endl;
+ }
+ if (err < best.err)
+ {
+ best = info;
+ best.err = err;
+ }
}
}
+ std::cout << std::endl;
+
+ if (best.err < 100)
+ print_code(best, name, std::numeric_limits<T>::max_digits10, std::numeric_limits<T>::digits);
+ else
+ std::cout << "Sorry, no viable approximation was found!!" << std::endl;
}
- std::cout << std::endl;
+ else
+ {
+ //
+ // Test a specific N and g:
+ //
+ std::cout << "Enter a value for g: ";
+ double g;
+ std::cin >> g;
+ lanczos_info<mp_t> info = generate_lanczos(N, mp_t(g));
+ mp_t err = get_max_error(info, eps, &exec_time);
- print_code(best, name);
+ std::cout << "N = " << N << std::endl;
+ std::cout << "g = " << g << std::endl;
+ std::cout << "eps error = " << static_cast<int>(err / eps) << std::endl;
+ std::cout << "Test time (ms) = " << static_cast<int>(exec_time * 1000) << std::endl;
+ print_code(info, name, std::numeric_limits<T>::max_digits10, std::numeric_limits<T>::digits);
+ }
}
+
+void scan_for_sweet_spots(int N)
+{
+ mp_t r = N * 0.66;
+ lanczos_info<mp_t> lower, upper(generate_lanczos(N + 1, r));
+ r += 0.1;
+ while (r < N)
+ {
+ lower = upper;
+ upper = generate_lanczos(N + 1, r);
+
+ if (lower.c.back() * upper.c.back() < 0)
+ {
+ std::pair<mp_t, mp_t> location = boost::math::tools::bisect([N](const mp_t& pos)
+ {
+ return generate_lanczos(N + 1, pos).c.back();
+ }, r - 0.1, r, boost::math::tools::eps_tolerance<mp_t>(20));
+ mp_t err = get_max_error(generate_lanczos(N, location.first), std::numeric_limits<mp_t>::epsilon());
+ std::cout << std::setprecision(15) << N << ", " << location.first << ", " << err << std::endl;
+ err = get_max_error(generate_lanczos(N, location.second), std::numeric_limits<mp_t>::epsilon());
+ std::cout << std::setprecision(15) << N << ", " << location.second << ", " << err << std::endl;
+ //std::cout << std::setprecision(15) << N << ", " << location.first << ", " << location.second << std::endl;
+ }
+
+ r += 0.1;
+ }
+}
+
int main(int argc, char*argv [])
{
- bool test_double(false), test_long(false), test_float(false), test_quad(false), spots(false), test_data(false);
+ bool test_double(false), test_long(false), test_float(false), test_quad(false), test_mp(false), spots(false), test_data(false), find_sweet(false);
if(argc < 2)
{
" -float test type float for the best approximation\n"
" -double test type double for the best approximation\n"
" -long-double test type long double for the best approximation\n"
+ " -quad test quad precision for the best approximation\n"
+ " -MP test current multiprecision type for the best approximation\n"
" -spots print out the best cases found in previous runs\n"
+ " -sweet Scan for more sweet spots for the arbitary parameter G: these "
+ " will need to cut and pasted into the big table at the start of this file"
+ " in order to search for greater precision approximations than otherwise supported."
" -data print out the test data\n" << std::flush;
return 1;
}
test_long = true;
else if(std::strcmp(argv[i], "-quad-float") == 0)
test_quad = true;
+ else if(std::strcmp(argv[i], "-MP") == 0)
+ test_mp = true;
else if(std::strcmp(argv[i], "-spots") == 0)
spots = true;
else if(std::strcmp(argv[i], "-data") == 0)
test_data = true;
+ else if(std::strcmp(argv[i], "-sweet") == 0)
+ find_sweet = true;
}
if(spots)
}
if(test_quad)
{
- //find_best_lanczos("quad_float", pow(NTL::quad_float(2), NTL::quad_float(-105)));
+ find_best_lanczos("quad_float", std::numeric_limits<boost::multiprecision::cpp_bin_float_quad>::epsilon());
}
+#ifdef MP_TYPE
+ if(test_mp)
+ {
+ find_best_lanczos("MP", std::numeric_limits<MP_TYPE>::epsilon());
+ }
+#endif
if(test_data)
{
std::cout << "Test Data follows:\n\n";
print_test_values(get_test_data_near_x<mp_t>(mp_t(-10)), "near_m10", 0);
print_test_values(get_test_data_near_x<mp_t>(mp_t(-55)), "near_m55", 0);
}
+ if (find_sweet)
+ {
+ unsigned last_index = sizeof(sweet_spots) / sizeof(sweet_spots[0]) - 1;
+ for(unsigned N = sweet_spots[last_index].N + 1; N < sweet_spots[last_index].N + 11; ++N)
+ scan_for_sweet_spots(N);
+ }
return 0;
}
projects / platform / upstream / boost.git / blobdiff