Re: left-right (up-down) parsing may lead to over/underflow of prod

[David Bateman]

Rolf Fabian wrote:
> I consider the following 'feature' as a
> major drawback of Octave because
> it infringes the expected commutativity
> and associativity while building products.
>
> Parsing products from left to right (for vectors)
> or up-to-down (columnwise product of matrices)
> may lead to avoidable over /underflow if
> some involved factors are extraordinarily
> large and/or small
>
> :> computer, version
> i686-pc-msdosmsvc
> ans = 3.0.0
>
> EXAMPLE 1 (prod of a vector)
> :> x1 = [1e150, 1e160, 1e150, 1e160, 1e-310, 1e-310]
> It is obvious that the correct product of x is equal to 1,
> but Octave either overflows if the product is
> evaluated from left-to-right:
> :> prod (x)
> ans = Inf
>
> or underflows (without any warning)
> if the product is evaluated from right-to-left
> :> prod  (fliplr (x))
> ans = 0
>
> Thus simply rearranging the order of factors changes
> Octave's result of the product from 'Inf' to '0', which
> is unexpected because  generally a product should
> not depend on the ordering of its factors.
>
> Expected correct answer is returned by
> :> prod_via_sum_of_logs = exp (sum (log (x)))
> ans =  1
>
> Similar artefacts occur if input to prod is a matrix like:
> :> y = [1e300, 1e-200; 1e300, 1e-200; 1e-200, 1e300; 1e-200, 1e300]
>
> :> prod(y)
> ans =
>    Inf     0
>
> :> prod (flipud (y))
> ans =
>      0   Inf
>
> The correct result can again be obtained via sums of logarithms:
> :> prod_via_sum_of_logs  = exp (sum (log (y)))
> ans =
>    1.0000e+200   1.0000e+200
>
> Of course, some spurious complex noise might appear if some factors
> are negative while using sums of logarithms, but this can be adjusted
> easily.
>
> In order to avoid similar artefacts as outlined above I propose to implement
> the 'prod' function via (column) sums of logarithms of factors.
>
> Rolf
>
>
> -----
> Rolf Fabian
> <r dot fabian at jacobs-university dot de>
>
>   
This is consistent with MatlabR2007b. See

>> x = [1e150, 1e160, 1e150, 1e160, 1e-310, 1e-310]

x =

  1.0e+160 *

    0.0000    1.0000    0.0000    1.0000    0.0000    0.0000

>> prod(x)

ans =

   Inf

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary