Complex special functions by bdeket · Pull Request #30 · racket/math

bdeket · 2019-10-11T16:15:44Z

Algorithms for calculating gamma, log-gamma and digamma in the complex field

bdeket · 2019-10-15T10:58:18Z

and complex erf and Fresnel-integrals

jbclements · 2020-01-19T22:28:50Z

math-doc/math/scribblings/math-special-functions.scrbl

 this the only type documented for @racket[lambert].

-@section[#:tag "real-functions"]{Real Functions}
+@section[#:tag "complex-real-functions"]{Complex and real Functions}


Capitalize Real in title?

jbclements

I'm largely unqualified to review this...

jbclements · 2020-01-19T22:30:02Z

math-doc/math/scribblings/math-special-functions.scrbl

+@defproc[(gamma [x Number]) Number]{
 Computes the @hyperlink["http://en.wikipedia.org/wiki/Gamma_function"]{gamma function},
-a generalization of the factorial function to the entire real line, except nonpositive integers.
+a generalization of the factorial function to the entire complex plane, except nonpositive exact integers.


This suggests that (gamma -13.0) works; is this correct?

this was already the case previously. (gamma -13.0) -> +nan.0

jbclements · 2020-01-19T22:31:41Z

math-doc/math/scribblings/math-special-functions.scrbl

-Error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
-ulps.
+On the real line the error is no more than 10 @tech{ulps} everywhere that has been tested, and is usually no more than 4
+ulps. In the rest of the complex plane the relative error is smaller than 1e-13 (@tech{ulps}=450).


I'm confused by this comment. Does the parenthetical mean "smaller than 450 ulps"?

I don't remember why I wrote that. I think I should just drop the part about ulps.
Also I did some extra testing and in the very close region around negative integers the error will be bigger.
I need to think on how to write this down

It may be better to leave a vaguer "Outside the real line, errors are significantly higher". If we're not certain of the guarantees we shouldn't give them, and should warn users that they are higher.

jbclements · 2020-01-19T22:33:48Z

math-doc/math/scribblings/math-special-functions.scrbl

+@section[#:tag "complex-real-functions"]{Complex and real Functions}

-@defproc[(gamma [x Real]) (U Positive-Integer Flonum)]{
+@defproc[(gamma [x Number]) Number]{


Could the result type be more specific? For instance, something including Float-Complex but not all complex numbers?

Is this a question about the documentation, or making the type of my function more specific?

jbclements · 2020-01-19T22:34:24Z

math-doc/math/scribblings/math-special-functions.scrbl

 }

-@defproc[(psi0 [x Real]) Flonum]{
+@defproc[(psi0 [x Number]) Number]{


cf earlier comment about output types

jbclements · 2020-01-19T22:35:28Z

math-doc/math/scribblings/math-special-functions.scrbl

-For @racket[erf], error is no greater than 2 @tech{ulps} everywhere that has been tested, and
-is almost always no greater than 1. For @racket[erfc], observed error is no greater than 4 ulps,
+For @racket[erf], on the real line the error is no greater than 2 @tech{ulps} everywhere that has been tested, and
+is almost always no greater than 1. In the complex plane the relative error is smaller 1e-12. For @racket[erfc], observed error is no greater than 4 ulps,


ITYM "smaller than 1e-12".

jbclements · 2020-01-19T22:37:27Z

math-lib/math/private/functions/lanczos.rkt

 (require "../../flonum.rkt")

-(provide lanczos-sum lanczos-g)
+(provide (all-defined-out))


I hate to see a specific provide turn into a generic (all-defined-out), unless... does it really make sense here?

The module had only 2 definitions, which were provided. I added 2 more, which will also be provided. Instead of writing all 4, I thought to use the generic. If this is not the right way, I see no problem to spelling everything out.

We should spell them out. (all-defined-out) is a code smell, since it fails to even describe the interface provided.

Algorithms for calculating gamma, log-gamma and digamma in the complex field

samth · 2020-07-02T02:45:08Z

I'm interested in @soegaard's thoughts here -- I think he knows this code/area the best.

bdeket · 2020-07-03T00:14:26Z

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.

soegaard · 2020-07-04T11:38:58Z

@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.

bdeket · 2020-07-05T22:24:02Z

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.

soegaard · 2020-07-06T15:03:49Z

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.

soegaard · 2020-07-06T15:27:20Z

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

pavpanchekha · 2020-08-05T19:40:02Z

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?

bdeket · 2020-08-06T07:12:02Z

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.

pavpanchekha · 2020-08-06T23:32:28Z

Sounds good, looking forward to it!

soegaard · 2023-03-26T16:04:56Z

@bdeket @pavpanchekha

Ping.

pavpanchekha · 2023-03-26T22:32:50Z

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.

bdeket · 2023-03-27T07:34:47Z

From what i remember when I started looking into improving it, one hurdle I encountered was: racket/racket#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.

jbclements reviewed Jan 19, 2020

View reviewed changes

bdeket added 6 commits January 25, 2020 21:45

Complex special functions

d0bde91

Algorithms for calculating gamma, log-gamma and digamma in the complex field

Complex erf, start fresnel

3f3cd78

Complex erf and Fresnel functions

6d959a9

separate bigfloat-fresnel into it's own file

fe463f5

Update 'type' for untyped code

9363d58

Update according to comments

91856a8

bdeket force-pushed the complex-gamma branch from eba68fb to 91856a8 Compare January 25, 2020 21:01

right sign for values with negative real-part

77c855e

complex-gamma: improve accuracy and sign for large (magnitude of) z

db737eb

Uh oh!

Conversation

bdeket commented Oct 11, 2019

Uh oh!

bdeket commented Oct 15, 2019

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jbclements left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

samth commented Jul 2, 2020

Uh oh!

bdeket commented Jul 3, 2020

Uh oh!

soegaard commented Jul 4, 2020

Uh oh!

bdeket commented Jul 5, 2020

Uh oh!

soegaard commented Jul 6, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

soegaard commented Jul 6, 2020

Uh oh!

pavpanchekha commented Aug 5, 2020

Uh oh!

bdeket commented Aug 6, 2020

Uh oh!

pavpanchekha commented Aug 6, 2020

Uh oh!

soegaard commented Mar 26, 2023

Uh oh!

pavpanchekha commented Mar 26, 2023

Uh oh!

bdeket commented Mar 27, 2023

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

soegaard commented Jul 6, 2020 •

edited

Loading