From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 1599 invoked by alias); 22 Nov 2007 16:13:41 -0000 Received: (qmail 1524 invoked by uid 22791); 22 Nov 2007 16:13:41 -0000 X-Spam-Check-By: sourceware.org Received: from qmta09.westchester.pa.mail.comcast.net (HELO QMTA09.westchester.pa.mail.comcast.net) (76.96.62.96) by sourceware.org (qpsmtpd/0.31) with ESMTP; Thu, 22 Nov 2007 16:13:33 +0000 Received: from OMTA02.westchester.pa.mail.comcast.net ([76.96.62.19]) by QMTA09.westchester.pa.mail.comcast.net with comcast id FrVc1Y00C0QuhwU0503A00; Thu, 22 Nov 2007 16:13:31 +0000 Received: from hippogriff.homeunix.org ([75.70.82.180]) by OMTA02.westchester.pa.mail.comcast.net with comcast id FsDW1Y00E3tRyZn0300000; Thu, 22 Nov 2007 16:13:31 +0000 X-Authority-Analysis: v=1.0 c=1 a=j-a4cumVLDQQMNjCRUcA:9 a=LX4ld3Pm7uj8SRR9uJ7qyV9-VL4A:4 a=XF7b4UCPwd8A:10 Received: by hippogriff.homeunix.org (Postfix, from userid 1000) id B2589317C9; Thu, 22 Nov 2007 09:13:31 -0700 (MST) Date: Thu, 22 Nov 2007 16:13:00 -0000 From: Patrick Alken To: Serhiy Lisovenko Cc: gsl-discuss@sourceware.org Subject: Re: Not clear to understand :) Message-ID: <20071122161331.GA17520@hippogriff.physics.drexel.edu> References: <1195731261.5477.14.camel@localhost.localdomain> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <1195731261.5477.14.camel@localhost.localdomain> User-Agent: Mutt/1.4.2.2i Mailing-List: contact gsl-discuss-help@sourceware.org; run by ezmlm Precedence: bulk List-Id: List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sourceware.org X-SW-Source: 2007-q4/txt/msg00037.txt.bz2 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). > > > >