bugGNU Octave - Bugs: bug #64801, besselj(1,z) is not pure...

 
 

bug #64801: besselj(1,z) is not pure imaginary for pure imaginary inputs

Submitter:  Lapresté Jean-Thierry <jtlap>
Submitted:  Sat 21 Oct 2023 12:54:53 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * 8.3.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 22 Oct 2023 03:39:26 AM UTC, comment #9: 

Are other cases than this zero-real case that should produce specific / exact answers? may be worth building up a set of tests to see how big a problem it is. if this is the only specific case, it might not be that onerous to correct it. But a wide range of cases could tedious to address.

I'm still surprised that the real and imaginary parts don't have separate eps values, since my understanding is that the real and imag parts are actually stored as two separate floats that could each have independent precision. maybe the calculation done in polar form which forces equal precision?

Nicholas Jankowski <nrjank>
Group Member
Sun 22 Oct 2023 12:30:31 AM UTC, comment #8: 

I'm still very ambivalent about whether it is a good idea to rewrite a numerical answer with known floating point characteristics to its mathematical counterpart. From a mathematical purity pov, that's OK, but it would hide any signs of e.g. badly conditioned inputs. Tough to assess that tradeoff.

Arun Giridhar <arungiridhar>
Group Member
Sun 22 Oct 2023 12:24:26 AM UTC, comment #7: 

In general with the amos Fortran code, Octave calls the external library, then passes the results back with very little reprocessing, especially if the library can already do what the user wants. But Octave can always rewrite the answers based on intent, such as overwriting the real part with true zero instead of only eps times the imaginary value, with all the positives and negatives that go with that decision. The last time I edited that side of the code was something to do with the Airy functions where we ended up calling the Fortran without extra intervention. But the code allows for it. See besselj.cc for examples of how the upper level functions look (take airy as an example), then lo-specfun.cc for how the C++ calls Fortran after reshaping values etc, then the amos directory for the low level Fortran itself. Any change like what OP asks for is best made in lo-specfun.cc after the Fortran call but before returning.

Arun Giridhar <arungiridhar>
Group Member
Sat 21 Oct 2023 11:14:27 PM UTC, comment #6: 

and maybe a follow up - in general how is the precision of complex numbers handled? are the real and imaginary parts stored in such a way that the two parts share floating point limits? or are they independently handled. Is this besselj issue inherent to how octave handles complex values or is it result of the way the calculation's being done?  (e.g., a vector magnitude calculation being later deconstructed into components.)

Nicholas Jankowski <nrjank>
Group Member
Sat 21 Oct 2023 09:54:07 PM UTC, comment #5: 

Thanks for confirming.  I don’t know much about the library itself, but the fact that the error for one thing that should be identically zero grows as it does is a bit concerning.  There is a compatibility interest to determining if there’s a reasonable way to independently address the accuracy of the two parts.

Nicholas Jankowski <nrjank>
Group Member
Sat 21 Oct 2023 05:35:01 PM UTC, comment #4: 

I don't know either. But you are right a floating point problem.

In fact I am writing the kyosu library based on eve which implements in scalar and simd complex and caley-dickson (quaternions octonions x-nions) and classical functions on them.

I used octave to verify some of my ouputs and I spotted this problem.

What I think is that a complex function is meant to return pure imaginary or real values in quite obvious an general circumstances it is better that it does.

However the problem can be complicated to fix.

I just wanted to signal the issue. Sorry for the noise.

Lapresté Jean-Thierry <jtlap>
Sat 21 Oct 2023 05:14:29 PM UTC, comment #3: 

The amos libraries are indeed used, yes.

The ratio of the real value divided by the imaginary value is indeed close to eps, so this is a limitation of floating point precision. The question is whether to return the raw values to the user (and deem them to be knowledgeable enough to expect the factor of eps) or whether Octave should (artificially) force the real value to zero in such cases. I don't have an opinion myself.

Arun Giridhar <arungiridhar>
Group Member
Sat 21 Oct 2023 02:45:23 PM UTC, comment #2: 

it's surely a bug as the error on the real part seems to tend to infinity as the i coef increase.

 besselj(1,1.e-300+100i)
ans =  6.541875770923201e+25 + 1.068369390338163e+42i

>>  besseli(1,1.e-300i+100)

ans =  1.068369390338163e+42 - 1.489643709947372e-45i

also the relation between besseli(1,z) and besselj(1,iz) is broken even for non pure imaginaries.

perhaps as besseli seems correct a workaround would be for n=1 to define bessej(1,...) in term of besseli ?

In fact do you use amos fortran here ? Is it the extern implementation which has problems ?


Lapresté Jean-Thierry <jtlap>
Sat 21 Oct 2023 02:26:07 PM UTC, comment #1: 

I can confirm the result with v8.3.0, and the result is incompatible with matlab.  Just a wild guess, but based on the magnitude of the error, could it be the real part is showing float precision limits based on the magnitude of eps for the imaginary part?:

Octave 8.3.0:

>> besseli(1,40)
ans = 1.470739616325935e+16

>> besselj(1,40*i)
ans =  9.005385394855129e-01 + 1.470739616325935e+16i

>> eps(besseli(1,40))
ans = 2


Matlab 2023b:

>> besseli(1,40)

ans =

     1.470739616325935e+16

>> besselj(1,40*i)

ans =

      0.000000000000000e+00 + 1.470739616325935e+16i



Nicholas Jankowski <nrjank>
Group Member
Sat 21 Oct 2023 12:54:53 PM UTC, original submission:  

the power serie defining beselj(1,z) in the complex plane is of the form:
 
   z\sum_0^{\infty}(a_i (z^2)^i) (DLMF 10.2.2 for details)

for pure imaginary z = 0+ib (b real) the result is obviously pure imaginary as the sum is real and mutiplied by z

>>  besselj(1,40*i)

ans =  9.005682817563816e-01 + 1.470739616325935e+16i

The imaginary part is ok but not the real part

The imaginary is also correctly equal to besseli(1,40)

40 is just an example. For smaller values of |z| the real part of this implementation is small but not zero as I think it should.

The error also affect scilab. I have no matlab to test with.

Wofram confirms my opinion see:

https://www.wolframalpha.com/input?i=Bessel+function+calculator&assumption=%7B%22F%22%2C+%22BesselFunctionsCalculator%22%2C+%22n%22%7D+-%3E%221%22&assumption=%7B%22FP%22%2C+%22BesselFunctionsCalculator%22%2C+%22besselfunction%22%7D+-%3E+%22BesselJ%22&assumption=%7B%22F%22%2C+%22BesselFunctionsCalculator%22%2C+%22z%22%7D+-%3E%2240i%22

Lapresté Jean-Thierry <jtlap>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by arungiridhar (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by jtlap (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.

     

    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2023-10-21 nrjank Summarybessej(1,z) is not pure imaginary for pure imaginary inputs besselj(1,z) is not pure imaginary for pure imaginary inputs
    2023-10-21 nrjank StatusNone Confirmed
        Release6.4.0 8.3.0
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.12