Not clear to understand :)

Not clear to understand :)

[Serhiy Lisovenko]

in file fit/linear.c

  for (i = 0; i < n; i++)
    {
      m_x += (x[i * xstride] - m_x) / (i + 1.0);
      m_y += (y[i * ystride] - m_y) / (i + 1.0);
    }

and other same...

code equivalent to

  for (i = 0; i < n; i++)
    {
      m_x +=x[i * xstride];
      m_y +=y[i * ystride];
    }
m_x/=(double)n;
m_y/=(double)n;

But the second code is more clear and some faster (no extra divisions).

RE: Not clear to understand :)

[keith.briggs]

These are *not* equivalent in floating point, of course.  The first is better.   
Even better:

  for (i = 1; i <= n; i++)
    {
      m_x += (x[i * xstride] - m_x) / i;
      m_y += (y[i * ystride] - m_y) / i;
    }

Keith

-----Original Message-----
From: gsl-discuss-owner@sourceware.org
[mailto:gsl-discuss-owner@sourceware.org]On Behalf Of Serhiy Lisovenko
Sent: 22 November 2007 11:34
To: gsl-discuss@sourceware.org
Subject: Not clear to understand :)


in file fit/linear.c

  for (i = 0; i < n; i++)
    {
      m_x += (x[i * xstride] - m_x) / (i + 1.0);
      m_y += (y[i * ystride] - m_y) / (i + 1.0);
    }

and other same...

code equivalent to

  for (i = 0; i < n; i++)
    {
      m_x +=x[i * xstride];
      m_y +=y[i * ystride];
    }
m_x/=(double)n;
m_y/=(double)n;

But the second code is more clear and some faster (no extra divisions).

Re: Not clear to understand :)

[Patrick Alken]

The way you propose is not numerically stable. For example suppose
the x vector consists of 2 elements which are both huge numbers
close to the maximum floating point limit. If you just add them
like you propose you will get an overflow, even though their mean
is a representable number.

Computing sums and sums of squares etc can be a challenging numerical
problem. I believe Knuth's book discusses ways to do these things
stably.

On Thu, Nov 22, 2007 at 01:34:21PM +0200, Serhiy Lisovenko wrote:
> in file fit/linear.c
> 
>   for (i = 0; i < n; i++)
>     {
>       m_x += (x[i * xstride] - m_x) / (i + 1.0);
>       m_y += (y[i * ystride] - m_y) / (i + 1.0);
>     }
> 
> and other same...
> 
> code equivalent to
> 
>   for (i = 0; i < n; i++)
>     {
>       m_x +=x[i * xstride];
>       m_y +=y[i * ystride];
>     }
> m_x/=(double)n;
> m_y/=(double)n;
> 
> But the second code is more clear and some faster (no extra divisions).
> 
> 
> 
>