[Octave-bug-tracker] [bug #62321] Wrong result for Airy funtion

[Arun Giridhar]

Follow-up Comment #3, bug #62321 (project octave):

The bug has been isolated to the scaling code in lo-specfun.cc here:
https://hg.savannah.gnu.org/hgweb/octave/file/30f7f409861a/liboctave/numeric/lo-specfun.cc#l146.
I added print statements as follows:


      if (! scaled)
        {
          std::cout << "\tBefore scaling: ar = " << ar << "  ai = " << ai <<
'\n';
          Complex expz = exp (- 2.0 / 3.0 * z * sqrt (z));

          double rexpz = expz.real ();
          double iexpz = expz.imag ();

          double tmp = ar*rexpz - ai*iexpz;

          ai = ar*iexpz + ai*rexpz;
          ar = tmp;
          std::cout << "\tAfter scaling: ar = " << ar << "  ai = " << ai <<
'\n';
        }


If a scaling option is passed, that if-block is skipped and the two results
agree, but if the scaling option is absent in the original call, then the
function changes things inconsistently:


octave:7> P = airy (0, V);
        Before scaling: ar = 0.458007  ai = -0.0982673
        After scaling: ar = 0.427202  ai = -0.192161

        Before scaling: ar = 0.296766  ai = 0.037377
        After scaling: ar = 0.301429  ai = 0.00181715

        Before scaling: ar = 0.431173  ai = -0.0506924
        After scaling: ar = 0.422305  ai = -0.100691

        Before scaling: ar = 0.222491  ai = 0.00567896
        After scaling: ar = 0.22802  ai = 0.00557015

        Before scaling: ar = 0.409248  ai = -0.0273306
        After scaling: ar = 0.406518  ai = -0.0545397

        Before scaling: ar = 0.113542  ai = -0.019165
        After scaling: ar = 0.122791  ai = 0.0179474

octave:8> Q = zeros(size(V));
octave:9> for t = 1:numel(V)
>   Q(t) = airy (0, V(t));
> end

        Before scaling: ar = 0.458007  ai = -0.0982673
        After scaling: ar = 0.468431  ai = 2.93374e-18

        Before scaling: ar = 0.296766  ai = 0.037377
        After scaling: ar = 0.301429  ai = 0.00181715

        Before scaling: ar = 0.431173  ai = -0.0506924
        After scaling: ar = 0.434143  ai = -5.71596e-18

        Before scaling: ar = 0.222491  ai = 0.00567896
        After scaling: ar = 0.22802  ai = 0.00557015

        Before scaling: ar = 0.409248  ai = -0.0273306
        After scaling: ar = 0.41016  ai = 1.93428e-18

        Before scaling: ar = 0.113542  ai = -0.019165
        After scaling: ar = 0.122791  ai = 0.0179474


The two values ar and ai are equal going into the block, and get changed
differently in there. I do not understand the math well enough to say what the
error is -- please comment on what is wrong there. If you have a known good
scaled value, it would help say whether P is in error or Q is in error.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62321>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/