bugGNU Scientific Library - Bugs: bug #26807, gsl_sum_levin_u_accel does not...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #26807: gsl_sum_levin_u_accel does not handle zero terms in series

Submitted by:  -Deleted Account- <bjg>
Submitted on:  Mon 15 Jun 2009 10:18:36 AM UTC  
 
Category: Runtime errorSeverity: 3 - Normal
Operating System: Status: Fixed
Assigned to: NoneOpen/Closed: Closed
Release: 1.12

Sun 28 Jun 2009 07:51:43 PM UTC, comment #1:

fixed by commit 29db8472d7faf57f37bf51530a95d7576070b802

-Deleted Account- <bjg>
Mon 15 Jun 2009 10:18:36 AM UTC, original submission:

From: 足立 聡 <adachi@aa.ap.titech.ac.jp>
To: -unavailable-
Cc: 足立 聡 <adachi@aa.ap.titech.ac.jp>
Subject: [Bug-gsl] A Bug of gsl_sum_levin_u_accel()
Date: Fri, 12 Jun 2009 15:21:57 +0900

Dear Developers of GSL,

(1) I found a little problem in GSL (version 1.12),
which concerns the acceleration of summation
based on the Levin u-transformation.

I am using GSL-1.12 on my Machintosh iBook G4 (PowerPC),
which is built from its souce code on my machine.
OS version is MacOS X 10.5.7.
The GNU C compiler is used. Its version is
"powerpc-apple-darwin9-gcc-4.0.1 (GCC) 4.0.1 (Apple Inc. build 5490)".

(1-1) The next sequence (t_{n}) is a trivial sequence for human
for the calculation of summation

t_{0} = 1,
t_{n} = 0 for n >= 1.

The summation sum(t_{n},n,0,inf) becomes

sum(t_{n},n,0,inf) = 1

(1-2) But, gsl_sum_levin_u_accel() seems to fail to calculate the
value of
this summation. The following is a program for demonstration:

-------------------------------------------------------------------------------
/*
* test_sum_accel_trivial.c:
*
* For demonstration of failure of gsl_sum_levin_u_accel() in
GSL-1.12.
* This program is based on gsl-1.12/sum/test.c written by G.
Jungman.
*
* S.Adachi 2009/06/07
*/

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sum.h>

#include <gsl/gsl_ieee_utils.h>

#define N 50

void check_trunc (double * t, double expected, const char * desc);
void check_full (double * t, double expected, const char * desc);

int main(int argc, char** argv)
{
gsl_ieee_env_setup();
gsl_test_verbose( 1 );

{
double t[N];
int n;

const double one = 1.0;

for (n = 0; n < N; n++)
t[n] = ((n == 0) ? 1.0 : 0.0);

check_trunc(t, one, "sequence (1,0,0,...)");
check_full(t, one, "sequence (1,0,0,...)");
}

return gsl_test_summary();
}

void
check_trunc (double * t, double expected, const char * desc)
{
double sum_accel, prec;

gsl_sum_levin_utrunc_workspace * w = gsl_sum_levin_utrunc_alloc (N);

gsl_sum_levin_utrunc_accel (t, N, w, &sum_accel, &prec);
gsl_test_rel (sum_accel, expected, 1e-8, "trunc result, %s", desc);

/* No need to check precision for truncated result since this is not
a meaningful number */

gsl_sum_levin_utrunc_free (w);
}

void
check_full (double * t, double expected, const char * desc)
{
double sum_accel, err_est, sd_actual, sd_est;

gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N);

gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err_est);
gsl_test_rel (sum_accel, expected, 1e-8, "full result, %s", desc);

sd_est = -log10 (err_est/fabs(sum_accel));
sd_actual = -log10 (DBL_EPSILON + fabs ((sum_accel - expected)/
expected));

/* Allow one digit of slop */

gsl_test (sd_est > sd_actual + 1.0, "full significant digits, %s
(%g vs %g)", desc, sd_est, sd_actual);

gsl_sum_levin_u_free (w);
}
-------------------------------------------------------------------------------

(1-3) The result of execution is as follows

-------------------------------------------------------------------------------
PASS: trunc result, sequence (1,0,0,...) (1 observed vs 1 expected)
FAIL: full result, sequence (1,0,0,...) (nan observed vs 1 expected)
PASS: full significant digits, sequence (1,0,0,...) (nan vs nan)
-------------------------------------------------------------------------------

The trivial value 1 can not be obtained by the full calculation.
Its result is Nan.

I think that this is a bug.
Please fix this bug.

(2) I explain that the bug reported above, which may look as to be not
so
serious, causes a serious real problem in programing codes for the
acceleration
of summation of series. I below explain the situation in which I met
the bug.

I need to evaluate a summation of power series in x,
sum(a_{n}x^n, n, 0, inf), where a_{n} does not depend on x.

I write a naive code for the subroutine of this evaluation.

The convergence by the naive code is so slow and I decided to use
gsl_sum_levin_u_accel() to utilize the Levin u-transformation
to accelerate the convergence.

gsl_sum_levin_u_accel() is wonderful.
gsl_sum_levin_u_accel() accelerates the convergence drastically
for general values of x.

The problem occurs when the parameter x = 0 is given to my subroutine.
In this case, the series becomes (a_{0},0,0,...) and the bug reported
in this letter appears.

As a workaround, I add a code gragment to my subroutine at its begining
to treat the case the absolute value of x is small enough separately.
In this case, the subroutine just returns a_{0} as the value.

However, what is the precise condition to judge
the absolute value of x is small enough?

I do not know the answer and I have written the added code fragment
in a ad hoc way.

I think that the situation that I meet the bug of
gsl_sum_levin_u_accel() is a natural one
as the situation in which the technique of the series acceleration is
used.
Many users of gsl_sum_levin_u_accel() may face the same problem.

So, please fix this bug.

Finally, thank you very much for your efforts to develop and maintain
GNU Scientific Library (GSL). GSL is really helpfull for me.

Sincerely yours,

Satoshi Adachi

_____________________________________________
Bug-gsl mailing list
-unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl

-Deleted Account- <bjg>

 

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Please enter the title of George Orwell's famous dystopian book (it's a date):

 

 

Follow 2 latest changes.

Date Changed By Updated Field Previous Value => Replaced By
Sun 28 Jun 2009 07:51:43 PM UTCbjgStatusConfirmed=>Fixed
  Open/ClosedOpen=>Closed

Back to the top


Powered by Savane 3.1-cleanup1