github.com/gopherd/gonum@v0.0.4/mathext/internal/cephes/lanczos.go (about)

     1  // Derived from SciPy's special/cephes/lanczos.c
     2  // https://github.com/scipy/scipy/blob/master/scipy/special/cephes/lanczos.c
     3  
     4  // Use of this source code is governed by a BSD-style
     5  // license that can be found in the LICENSE file.
     6  // Copyright ©2006 John Maddock
     7  // Portions Copyright ©2003 Boost
     8  // Portions Copyright ©2016 The Gonum Authors. All rights reserved.
     9  
    10  package cephes
    11  
    12  // Optimal values for G for each N are taken from
    13  // http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
    14  // as are the theoretical error bounds.
    15  
    16  // Constants calculated using the method described by Godfrey
    17  // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
    18  // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
    19  
    20  var lanczosNum = [...]float64{
    21  	2.506628274631000270164908177133837338626,
    22  	210.8242777515793458725097339207133627117,
    23  	8071.672002365816210638002902272250613822,
    24  	186056.2653952234950402949897160456992822,
    25  	2876370.628935372441225409051620849613599,
    26  	31426415.58540019438061423162831820536287,
    27  	248874557.8620541565114603864132294232163,
    28  	1439720407.311721673663223072794912393972,
    29  	6039542586.35202800506429164430729792107,
    30  	17921034426.03720969991975575445893111267,
    31  	35711959237.35566804944018545154716670596,
    32  	42919803642.64909876895789904700198885093,
    33  	23531376880.41075968857200767445163675473,
    34  }
    35  
    36  var lanczosDenom = [...]float64{
    37  	1,
    38  	66,
    39  	1925,
    40  	32670,
    41  	357423,
    42  	2637558,
    43  	13339535,
    44  	45995730,
    45  	105258076,
    46  	150917976,
    47  	120543840,
    48  	39916800,
    49  	0,
    50  }
    51  
    52  var lanczosSumExpgScaledNum = [...]float64{
    53  	0.006061842346248906525783753964555936883222,
    54  	0.5098416655656676188125178644804694509993,
    55  	19.51992788247617482847860966235652136208,
    56  	449.9445569063168119446858607650988409623,
    57  	6955.999602515376140356310115515198987526,
    58  	75999.29304014542649875303443598909137092,
    59  	601859.6171681098786670226533699352302507,
    60  	3481712.15498064590882071018964774556468,
    61  	14605578.08768506808414169982791359218571,
    62  	43338889.32467613834773723740590533316085,
    63  	86363131.28813859145546927288977868422342,
    64  	103794043.1163445451906271053616070238554,
    65  	56906521.91347156388090791033559122686859,
    66  }
    67  
    68  var lanczosSumExpgScaledDenom = [...]float64{
    69  	1,
    70  	66,
    71  	1925,
    72  	32670,
    73  	357423,
    74  	2637558,
    75  	13339535,
    76  	45995730,
    77  	105258076,
    78  	150917976,
    79  	120543840,
    80  	39916800,
    81  	0,
    82  }
    83  
    84  var lanczosSumNear1D = [...]float64{
    85  	0.3394643171893132535170101292240837927725e-9,
    86  	-0.2499505151487868335680273909354071938387e-8,
    87  	0.8690926181038057039526127422002498960172e-8,
    88  	-0.1933117898880828348692541394841204288047e-7,
    89  	0.3075580174791348492737947340039992829546e-7,
    90  	-0.2752907702903126466004207345038327818713e-7,
    91  	-0.1515973019871092388943437623825208095123e-5,
    92  	0.004785200610085071473880915854204301886437,
    93  	-0.1993758927614728757314233026257810172008,
    94  	1.483082862367253753040442933770164111678,
    95  	-3.327150580651624233553677113928873034916,
    96  	2.208709979316623790862569924861841433016,
    97  }
    98  
    99  var lanczosSumNear2D = [...]float64{
   100  	0.1009141566987569892221439918230042368112e-8,
   101  	-0.7430396708998719707642735577238449585822e-8,
   102  	0.2583592566524439230844378948704262291927e-7,
   103  	-0.5746670642147041587497159649318454348117e-7,
   104  	0.9142922068165324132060550591210267992072e-7,
   105  	-0.8183698410724358930823737982119474130069e-7,
   106  	-0.4506604409707170077136555010018549819192e-5,
   107  	0.01422519127192419234315002746252160965831,
   108  	-0.5926941084905061794445733628891024027949,
   109  	4.408830289125943377923077727900630927902,
   110  	-9.8907772644920670589288081640128194231,
   111  	6.565936202082889535528455955485877361223,
   112  }
   113  
   114  const lanczosG = 6.024680040776729583740234375
   115  
   116  func lanczosSum(x float64) float64 {
   117  	return ratevl(x,
   118  		lanczosNum[:],
   119  		len(lanczosNum)-1,
   120  		lanczosDenom[:],
   121  		len(lanczosDenom)-1)
   122  }
   123  
   124  func lanczosSumExpgScaled(x float64) float64 {
   125  	return ratevl(x,
   126  		lanczosSumExpgScaledNum[:],
   127  		len(lanczosSumExpgScaledNum)-1,
   128  		lanczosSumExpgScaledDenom[:],
   129  		len(lanczosSumExpgScaledDenom)-1)
   130  }
   131  
   132  func lanczosSumNear1(dx float64) float64 {
   133  	var result float64
   134  
   135  	for i, val := range lanczosSumNear1D {
   136  		k := float64(i + 1)
   137  		result += (-val * dx) / (k*dx + k*k)
   138  	}
   139  
   140  	return result
   141  }
   142  
   143  func lanczosSumNear2(dx float64) float64 {
   144  	var result float64
   145  	x := dx + 2
   146  
   147  	for i, val := range lanczosSumNear2D {
   148  		k := float64(i + 1)
   149  		result += (-val * dx) / (x + k*x + k*k - 1)
   150  	}
   151  
   152  	return result
   153  }