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

Submitter:  -Deleted Account- <bjg>
Submitted:  Mon 15 Jun 2009 10:18:36 AM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Status:  Fixed
Assigned to:  None Open/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: -email is 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
-email is unavailable-
http://lists.gnu.org/mailman/listinfo/bug-gsl


-Deleted Account- <bjg>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

CC list is empty

 

Follow 2 latest changes.

Date Changed by Updated Field Previous Value => Replaced by
2009-06-28 bjg StatusConfirmed Fixed
    Open/ClosedOpen Closed

Back to the top

Powered by Savane 3.13-758e.
Corresponding source code