bugGNU Octave - Bugs: bug #43492, warning: division by zero for some...

 
 

bug #43492: warning: division by zero for some nthroots of zero

Submitter:  Marco Caliari <caliari>
Submitted:  Wed 29 Oct 2014 09:40:19 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  None
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 3.8.1 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Thu 30 Oct 2014 03:21:29 PM UTC, comment #4: 

The division by zero warning is the title and focus of the bug report.  This warning is only emitted for division by a scalar zero.  Division by a vector or matrix that contains zero does not produce a warning, so there is no code change required.

The existing isfinite code catches any 0/0 divisions in a vector/matrix combination because this operation will produce NaN which is not finite.

For why scalar division produces a warning and vector does not you could ask on the Octave maintainers list.  I suspect the answer is that checking a single scalar for being equal to zero is cheap, but checking an NxN matrix would require N^2 operations which would quickly grow cumbersome.  Octave also probably just forwards on to BLAS any matrix/matrix or matrix/vector operations and BLAS does not emit a warning for division by zero since the result is defined as either Inf (value != 0) or NaN (0 / 0).

If you want to solve the different problem of actually preventing a division by zero entirely, you could search the y vector for non-zeros and then perform the Newton correction only on those entries.  Benchmarking would be required to make sure this isn't slower.

As to why (2^12)^(1/6) does not equal exactly 4, that is probably running into precision issues because the numbers are so different in size.  You can use eps() to see the machine precision for each number, and they differ by about 3 orders of magnitude.


eps (2^12)
ans =    9.0949e-13
eps (1/6)
ans =    2.7756e-17


On the other hand, when calculating 2^(log2(2^12)/6), you are dealing with the numbers 12 and 6 for which the machine precision is just about the same.


eps (12)
ans =    1.7764e-15
eps (6)
ans =    8.8818e-16



Rik <rik5>
Group administrator
Thu 30 Oct 2014 12:33:54 PM UTC, comment #3: 

I found what the correction is: it is one step of Newton method (see http://en.wikipedia.org/wiki/Nth_root_algorithm). I also found that without the correction, nthroot(2^12,6) is not exactly 4 (the error is eps/2), while it is with the correction. Anyway, it is not clear to me why a single step of Newton method should give a better result in general.

What it is less obvious to me is the following: both 2^(log2(2^12)/6) and exp(log(2^12)/6) give exactly 4 but (2^12)^(1/6) not. How is it computed?

Marco

Marco Caliari <caliari>
Group Member
Thu 30 Oct 2014 08:26:31 AM UTC, comment #2: 

Thanks Rik,

your patch does not prevent division by zero (try nthroot([0,0],2)), but [0,0]./[0,0] does not give the warning. The proper fix would be to understand the correction and correct it. I cannot understand it, too.

Marco

Marco Caliari <caliari>
Group Member
Wed 29 Oct 2014 03:23:35 PM UTC, comment #1: 

I pushed a fix for this to the development branch here (http://hg.savannah.gnu.org/hgweb/octave/rev/238522618904).  Thanks for noticing and reporting this.

Rik <rik5>
Group administrator
Wed 29 Oct 2014 09:40:19 AM UTC, original submission:  

Dear all,

I get the warning about division by zero if I compute nthroot(0,n) where n is 2 or integer > 3. The problem is in the final correction. The obvious thing it to do nothing for the elements of x which are zero, but I don't know how to get this result in a nice way.

Best regards,

Marco

Marco Caliari <caliari>
Group Member

 

(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 rik5 (Posted a comment)
  • -email is unavailable- added by caliari (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 group members can vote.

     

    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-10-29 rik5 StatusNone Fixed
        Open/ClosedOpen Closed

    Back to the top

    Powered by Savane 3.13-caa5.
    Corresponding source code