[Bug-gsl] [bug #26807] gsl_sum_levin_u_accel does not handle zero terms in series

[Brian Gough]

URL:
  <http://savannah.gnu.org/bugs/?26807>

                 Summary: gsl_sum_levin_u_accel does not handle zero terms in
series
                 Project: GNU Scientific Library
            Submitted by: bjg
            Submitted on: Mon 15 Jun 2009 11:18:36 AM BST
                Category: Runtime error
                Severity: 3 - Normal
        Operating System: 
                  Status: Confirmed
             Assigned to: None
             Open/Closed: Open
                 Release: 1.12
         Discussion Lock: Any

    _______________________________________________________

Details:

From: 足立 聡 <address@hidden>
To: address@hidden
Cc: 足立 聡 <address@hidden>
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
address@hidden
http://lists.gnu.org/mailman/listinfo/bug-gsl






    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?26807>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/