[Octave-bug-tracker] [bug #60682] betainc is inaccurate

[Rik]

Update of bug #60682 (project octave):

                  Status:                    None => Fixed                  
             Open/Closed:                    Open => Closed                 
                 Release:                   6.1.0 => dev                    

    _______________________________________________________

Follow-up Comment #6:

I made some modifications for speed and then checked in Michele's code here:
http://hg.savannah.gnu.org/hgweb/octave/rev/ad6a57b215e8.

Michael was in favor of adopting code for these special (although important)
cases if the overhead was not too large.

I benchmarked using a 1M point vector generated with "rand (1e6, 1)".  Hence,
no value of input 'a' or 'b' was actually equal to 1.  Thus, original code and
new code were both using continued fraction expansion for calculation of
output, and differences in timing were solely down to the extra input
processing.

For original code, mean run time (N = 25) was 380 milliseconds.

For Michele's code shown below, mean run time was 396 milliseconds.


  a_b_one = (a == 1) & (b == 1);
  a_one = (a == 1) & (b != 1);
  b_one = (a != 1) & (b == 1);
  non_trivial = (a != 1) & (b != 1);


An increase of 4.2% in run time is probably okay, but I decided to re-code to
trade memory for speed since most machines have copious memory these days.

New code:


  ## Trivial cases (long code here trades memory for speed)
  a_one = (a == 1);
  b_one = (b == 1);
  a_b_one = a_one & b_one;
  a_not_one = ! a_one;
  b_not_one = ! b_one;
  non_trivial = a_not_one & b_not_one;
  a_one &= b_not_one;
  b_one &= a_not_one;


Mean run time is now 382 milliseconds.  Percentage increase is just 0.5% so I
checked in the change.

Marking as Fixed and closing report.

    _______________________________________________________

Reply to this item at:

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

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