math
math copied to clipboard
racket
Complex special functions
Algorithms for calculating gamma, log-gamma and digamma in the complex field
I'm interested in @soegaard's thoughts here -- I think he knows this code/area the best.
I'm doing some automated checking of my implementation's results against Wolfram-alpha and I found some problems. I will try to fix them next week.
@bdeket If you are looking for way to write some randomized tests, here are some ideas.
For a real nummer kappa, the following holds [1]:
This can be used to test gammain a randomized manner.
Generate random, real k (kappa) and check the results are equal:
(define (h k) (sqr (magnitude (gamma (* 0+i k)))))
(define (j k) (/ pi (* z (sinh (* pi k)))))
Similar identities:
from [2].
/Jens Axel
[1] https://en.wikipedia.org/wiki/Particular_values_of_the_gamma_function [2] https://en.wikipedia.org/wiki/Gamma_function#:~:text=Main%20definition,-The%20notation%20is&text=to%20a%20meromorphic%20function%20defined,to%20as%20the%20gamma%20function.
I tried to get a general idea of the error in https://gist.github.com/bdeket/dd84a9a012dc3e5c7eba45afc43be4b4 the main file submits a calculation to wolfram alpha and stores it since the (free) api only let's you make 2000 calls a month.
Results: The (relative)error is <1e-12 as mentioned in the documentation, however for values where the value should go to 0 or to inf, my implementation often goes to nan, or to 0 to soon.
some figures: green is flulp-error < 10 blue < 100 orange < 1000 red < 50000 red circles are > 50000
digamma: (larger errors are near the negative real axis ±0.1i)
gamma: (errors due to nan instead of 0 or 0 instead of ~1e-300)
erf:
In conclusion, I don't think this should be merged. But if you have any tips on what would be the minimum to be acceptable, please tell.
Really nice diagrams. Your comments about using Wolfram Alpha to compute test values made be resurrect some old code that talks to Maxima. Now Maxima is not as precise as Wolfram, but you can get a quicker turnaround when doing randomized testing. Your data set contains several points, where Wolfram Alpha computes a zero, but Maxima computes an overflow.
https://gist.github.com/soegaard/cdba96801778b37713b4cf3add497316
Could you make a similar diagram comparing your gamma with Maxima? And perhaps Maxima vs Wolfram?
If your gamma and Maxima differs more or less the same from Wolfram, then I'd say that's fine enough for inclusion.
If you need more digits from Maxima, I have just found out how to activate bigfloats:
(%i7) fpprec:32;
(%o7) 32
(%i8) bfloat(gamma(1+%i));
(%o8) 4.980156681183560427136911174622b-1
- 1.5494982830181068512495513048389b-1 %i
I don't really have the math background all loaded in, so this is a kind of a shot in the dark, but in computing m on gamma.rkt:288, you're going to see a lot of error. It looks like you're computing a sum over simple poles, but once you're far from a pole, the alternating signs of the coefficients will lead to a buildup of errors. Unfortunately @ztatlock has my copy of Numerical Recipes, but rewriting this fold to use flsum would help.
How can I run your WolframAlpha tests to see if that works?
Maybe this PR should be changed to draft. I'm rewriting the functions and doing more thorough testing. I currently don't have my computer with me, but next week I can share the Wolfram tests and results.
Well, I'm not totally sure what to do about this PR. IIUC, these implementations have pretty high error on some inputs (bad) and @bdeket was thinking about rewriting them but hasn't. I suppose, given the time elapsed, we shouldn't wait on that.
Perhaps the error is not so bad and we should just merge these implementations—though @bdeket disagrees. Fixing up the remaining errors—including some underflows and possibly some cancellations—might fix it, but to be honest I don't really have the mathematical background to be sure it will.
From what i remember when I started looking into improving it, one hurdle I encountered was: https://github.com/racket/racket/issues/3324.
That led me to start a local branch that offers functions like fl+ fltan, etc but working on complex floats, and I never got far enough to finish that...
I'll spend some time on these complex special functions in the next weeks to better see the error margins, and report back.