Discussion:
Bessel function scaling + limited range?
(too old to reply)
Kris Kuhlman
2008-02-04 17:32:00 UTC
Permalink
I am trying to understand the scaling that Octave performs with besseli and
besselk for real arguments. I have worked with Amos' Bessel function
routines in fortran, I understand that they return exp(-x)*besseli(n,x) and
exp(x)*besselk(n,x) when you request a scaled Bessel function. I have seen
a couple previous bug reports in the Octave list regarding this topic, but
they did not address my problem.

Why does Octave only return ~ 1.0E+33 as a maximum range? When the argument
is near 1 the exponential scaling has very little effect, so I can't compute
very high orders of Bessel functions. An example of my usage is below.
Basically, for an argument I need a vector of Bessel functions of order 0 to
~35, Octave can only muster 0:23.

for e=1:ne
v(1,e) = sqrt(q)*exp(-eta(e));
v(2,e) = sqrt(q)*exp( eta(e));
besi(1:R+1,e) = besseli(0:R, v(1,e)).';
besk(1:R+1,e) = besselk(0:R, v(2,e)).';

for n=1:N
Ke(n,e) = sum(AEv(1:R,n).* besi(1:R,e).*besk(1:R,e));
Ko(n,e) = sum(AOd(1:R,n).* ...
(besi(1:R,e).*besk(2:R+1,e) + besi(2:R+1,e).*besk(1:R,e)));
end
end

In the case I have been looking at just now, v(2,e) is ~0.6, so it goes like
this: The scaling option (the third argument) doesn't help me here.

GNU Octave, version 3.0.0
-- snip --
Octave was configured for "x86_64-unknown-linux-gnu".
-- snip --
octave:1> besselk(23,0.6)
ans = 5.9453e+32
octave:2> besselk(24,0.6)
ans = Inf + Infi
octave:3> besselk(23,0.6,1)
ans = 1.0833e+33
octave:4>besselk(24,0.6,1)
ans = Inf + Infi

I would appreciate any pointers or advice as to how I can move my matlab
code over to octave (this works fine in Matlab), I guess I will not be able
to have neutral code, I will have to write something Octave-specific into my
code?

Thanks,

Kris
--
Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/
John W. Eaton
2008-02-04 21:08:22 UTC
Permalink
On 4-Feb-2008, Kris Kuhlman wrote:

| I am trying to understand the scaling that Octave performs with besseli and
| besselk for real arguments. I have worked with Amos' Bessel function
| routines in fortran, I understand that they return exp(-x)*besseli(n,x) and
| exp(x)*besselk(n,x) when you request a scaled Bessel function. I have seen
| a couple previous bug reports in the Octave list regarding this topic, but
| they did not address my problem.
|
| Why does Octave only return ~ 1.0E+33 as a maximum range? When the argument
| is near 1 the exponential scaling has very little effect, so I can't compute
| very high orders of Bessel functions. An example of my usage is below.
| Basically, for an argument I need a vector of Bessel functions of order 0 to
| ~35, Octave can only muster 0:23.
|
| for e=1:ne
| v(1,e) = sqrt(q)*exp(-eta(e));
| v(2,e) = sqrt(q)*exp( eta(e));
| besi(1:R+1,e) = besseli(0:R, v(1,e)).';
| besk(1:R+1,e) = besselk(0:R, v(2,e)).';
|
| for n=1:N
| Ke(n,e) = sum(AEv(1:R,n).* besi(1:R,e).*besk(1:R,e));
| Ko(n,e) = sum(AOd(1:R,n).* ...
| (besi(1:R,e).*besk(2:R+1,e) + besi(2:R+1,e).*besk(1:R,e)));
| end
| end

This code snippet doesn't help much since you ahven't sent any of the
data to go along with it (waht are the values for ne, e, eta, q, R,
AEv, AOd, ...)?

| In the case I have been looking at just now, v(2,e) is ~0.6, so it goes like
| this: The scaling option (the third argument) doesn't help me here.
|
| GNU Octave, version 3.0.0
| -- snip --
| Octave was configured for "x86_64-unknown-linux-gnu".
| -- snip --
| octave:1> besselk(23,0.6)
| ans = 5.9453e+32
| octave:2> besselk(24,0.6)
| ans = Inf + Infi
| octave:3> besselk(23,0.6,1)
| ans = 1.0833e+33
| octave:4>besselk(24,0.6,1)
| ans = Inf + Infi
|
| I guess I will not be able
| to have neutral code, I will have to write something Octave-specific into my
| code?

Or, you could help us fix the problem in Octave. You say you have
experience using the Amos functions, so perhaps you can tell us how to
use them properly to solve your problem?

I'm attaching a test program in Fortran that calls ZBESK from Amos
directly. The test program includes ZBESK and all of its
dependencies. It calls ZBESK with KODE = 1 and then 2:

C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
C KODE= 1 RETURNS
C CY(I)=K(FNU+I-1,Z), I=1,...,N
C = 2 RETURNS
C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
C

In both cases, for Z = 0.6 and FNU = 24, ZBESK returns IERR = 2, which
indicates that an overflow has occurred. I don't know how to avoid
that problem.

Is there a newer version of the Amos library around somewhere that
fixes this problem?

jwe
Kris Kuhlman
2008-02-04 21:33:06 UTC
Permalink
I have used a modified fortran90 version of these libraries, as distributed
by Allan Miller; TOMS 644. The fortran code was downloaded from the
address below.

http://users.bigpond.net.au/amiller/toms.html

Kris
Post by John W. Eaton
| I am trying to understand the scaling that Octave performs with besseli and
| besselk for real arguments. I have worked with Amos' Bessel function
| routines in fortran, I understand that they return exp(-x)*besseli(n,x) and
| exp(x)*besselk(n,x) when you request a scaled Bessel function. I have seen
| a couple previous bug reports in the Octave list regarding this topic, but
| they did not address my problem.
|
| Why does Octave only return ~ 1.0E+33 as a maximum range? When the argument
| is near 1 the exponential scaling has very little effect, so I can't compute
| very high orders of Bessel functions. An example of my usage is below.
| Basically, for an argument I need a vector of Bessel functions of order 0 to
| ~35, Octave can only muster 0:23.
|
| for e=1:ne
| v(1,e) = sqrt(q)*exp(-eta(e));
| v(2,e) = sqrt(q)*exp( eta(e));
| besi(1:R+1,e) = besseli(0:R, v(1,e)).';
| besk(1:R+1,e) = besselk(0:R, v(2,e)).';
|
| for n=1:N
| Ke(n,e) = sum(AEv(1:R,n).* besi(1:R,e).*besk(1:R,e));
| Ko(n,e) = sum(AOd(1:R,n).* ...
| (besi(1:R,e).*besk(2:R+1,e) + besi(2:R+1,e).*besk(1:R,e)));
| end
| end
This code snippet doesn't help much since you ahven't sent any of the
data to go along with it (waht are the values for ne, e, eta, q, R,
AEv, AOd, ...)?
| In the case I have been looking at just now, v(2,e) is ~0.6, so it goes like
| this: The scaling option (the third argument) doesn't help me here.
|
| GNU Octave, version 3.0.0
| -- snip --
| Octave was configured for "x86_64-unknown-linux-gnu".
| -- snip --
| octave:1> besselk(23,0.6)
| ans = 5.9453e+32
| octave:2> besselk(24,0.6)
| ans = Inf + Infi
| octave:3> besselk(23,0.6,1)
| ans = 1.0833e+33
| octave:4>besselk(24,0.6,1)
| ans = Inf + Infi
|
| I guess I will not be able
| to have neutral code, I will have to write something Octave-specific into my
| code?
Or, you could help us fix the problem in Octave. You say you have
experience using the Amos functions, so perhaps you can tell us how to
use them properly to solve your problem?
I'm attaching a test program in Fortran that calls ZBESK from Amos
directly. The test program includes ZBESK and all of its
C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
C KODE= 1 RETURNS
C CY(I)=K(FNU+I-1,Z), I=1,...,N
C = 2 RETURNS
C CY(I)=K(FNU+I-1,Z)*EXP(Z), I=1,...,N
C
In both cases, for Z = 0.6 and FNU = 24, ZBESK returns IERR = 2, which
indicates that an overflow has occurred. I don't know how to avoid
that problem.
Is there a newer version of the Amos library around somewhere that
fixes this problem?
jwe
--
Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/ <http://www.u.arizona.edu/%7Ekuhlman/>
--
Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/
John W. Eaton
2008-02-04 22:07:55 UTC
Permalink
On 4-Feb-2008, Kris Kuhlman wrote:

| I have used a modified fortran90 version of these libraries, as distributed
| by Allan Miller; TOMS 644. The fortran code was downloaded from the
| address below.
|
| http://users.bigpond.net.au/amiller/toms.html

OK, looking at this page, it is prefaced with

N.B. I have been asked to provide a link to the copyright policy of
the ACM. Loosely paraphrased, this allows use, and modification, of
the TOMS algorithms for most non-commercial purposes. It also
emphasizes that the ACM accepts no responsibility for the accuracy
of the code.

so if these functions are subject to the ACM copyright, then I don't
think we can use them in Octave since the "non-commercial use only"
clause makes the license incompatible with the GPL.

OK, I see that the toms644.zip file contains a translation to Fortran
90 of an update to the original TOMS 644 algorithm.

The code in Octave comes from netlib.org, which I think is just the
original algorithm from 1986. I assumed (perhaps incorrectly) that it
was OK to use these files, but if anyone knows differently, then we
should probably stop using them. In any case, to avoid using code
with stupid restrictions like this, it would be good to find an
alternative. Maybe we should just use GSL for this, as it doesn't
seem to suffer from this problem? However, switching to GSL may not
be trivial, as it seems to rely on cblas (at least on my system) and I
think it would be best to use ATLAS and not also cblas. If someone
would like to follow up on this, I'd suggest moving the discussion to
hte maintainers list.

Thanks,

jwe
Dmitri A. Sergatskov
2008-02-04 22:38:40 UTC
Permalink
FWIW:

gsl's version of besselk function seams to work:

#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_sf_bessel.h>

#define NPT (20)

int
main(void)
{
int nmin;
int nmax;
int idx;
double x;
double result[NPT+1];

nmin = 15;
nmax = nmin+NPT;

x = 0.6;

gsl_sf_bessel_Kn_scaled_array (nmin, nmax, x, result);

for(idx=0; idx<NPT; idx++) {
printf("result %li = %g \n", nmin+idx, result[idx]);
}
return 0;

[]$ gcc testgsl.c -lgsl -lgslcblas -lm
[]$ ./a.out
result 15 = 5.49978e+18
result 16 = 2.75107e+20
result 17 = 1.46779e+22
result 18 = 8.3202e+23
result 19 = 4.99359e+25
result 20 = 3.16344e+27
result 21 = 2.10946e+29
result 22 = 1.47694e+31
result 23 = 1.0833e+33
result 24 = 8.30676e+34
result 25 = 6.64649e+36
result 26 = 5.53957e+38
result 27 = 4.80163e+40
result 28 = 4.32202e+42
result 29 = 4.03437e+44
result 30 = 3.90032e+46
result 31 = 3.90072e+48
result 32 = 4.03114e+50
result 33 = 4.30027e+52
result 34 = 4.7307e+54

Sincerely,

Dmitri.
--
Kris Kuhlman
2008-02-07 02:02:25 UTC
Permalink
Post by John W. Eaton
so if these functions are subject to the ACM copyright, then I don't
think we can use them in Octave since the "non-commercial use only"
clause makes the license incompatible with the GPL.
So maybe this is an uninformed question or comment, but I will state what
seems to be to be the obvious. Since Matlab uses the same Amos Fortran
library that Octave does (according to the "help besseli" command), can
there really be a non-commercial-use license attached with it?

Secondly, if Octave uses the same Amos Fortran libraries as Matlab, why does
it's precision only go to 1.0E+30, while Matlab's goes much higher? (this
question is sort of the point of my initial question in this thread, but I
suppose moving Bessel function support over to GSL is a good idea anyway.)

Since it sounds like the move to GSL will not happen soon, could it be that
maybe the Amos library is configured or compiled wrong?

I have attached a small f90 program and its output (named screen) that
drives the f90 amos code I posted the earlier link to. It gives me correct
results well beyond where Octave stops. I have tried to do the same with
the amos f77 code in the libcruft directory, but I get junk. I am not sure
why.

Curious,

Kris
--
Kristopher L. Kuhlman
PhD candidate Department of Hydrology & Water Resources
University of Arizona
(520) 621-1380
http://www.u.arizona.edu/~kuhlman/
John W. Eaton
2008-02-07 04:06:03 UTC
Permalink
On 6-Feb-2008, Kris Kuhlman wrote:

| On Feb 4, 2008 3:07 PM, John W. Eaton <***@bevo.che.wisc.edu> wrote:
|
| > so if these functions are subject to the ACM copyright, then I don't
| > think we can use them in Octave since the "non-commercial use only"
| > clause makes the license incompatible with the GPL.
| >
|
| So maybe this is an uninformed question or comment, but I will state what
| seems to be to be the obvious. Since Matlab uses the same Amos Fortran
| library that Octave does (according to the "help besseli" command), can
| there really be a non-commercial-use license attached with it?

Maybe they negotiated some different terms? If you are interested in
having updated versions of these routines in Octave, then please ask
the ACM if they would allow distribution under terms that are
compatible with the GPL.

| Secondly, if Octave uses the same Amos Fortran libraries as Matlab, why does
| it's precision only go to 1.0E+30, while Matlab's goes much higher?

Perhaps they fixed the problem, or used an updated version? The
README file from the version of the Amos library that you linked to
includes the following references:

2. Amos, D. E., Algorithm 644, A Portable Package For Bessel Functions of a
Complex Argument and Nonnegative Order, ACM Transactions on Mathematical
Software, Vol. 12, No. 3, September 1986, Pages 265-273

3. Amos, D. E., Remark on Algorithm 644, ACM Transactions on Mathematical
Software, Vol. 16, No. 4, December 1990, Page 404

4. Amos, D. E., Remark on Algorithm 644 (Improvements in Algorithm 644),
ACM Transactions on Mathematical Software, (Estimated date 1995)

I haven't looked these up, but I guess the later publications include
some bug fixes. As far as I know, the version of the code that we are
using (from netlib.org) is just what was published in 1986.

| Since it sounds like the move to GSL will not happen soon, could it be that
| maybe the Amos library is configured or compiled wrong?

There is nothing to configure.

| I have attached a small f90 program and its output (named screen) that
| drives the f90 amos code I posted the earlier link to. It gives me correct
| results well beyond where Octave stops. I have tried to do the same with
| the amos f77 code in the libcruft directory, but I get junk. I am not sure
| why.

Presumably the f90 version you have includes the updates published in
1990 and 1995. But those are apparently subject to the newer ACM
copyright, so we can't use them.

jwe
Kris Kuhlman
2008-02-07 05:02:05 UTC
Permalink
Post by John W. Eaton
| So maybe this is an uninformed question or comment, but I will state what
| seems to be to be the obvious. Since Matlab uses the same Amos Fortran
| library that Octave does (according to the "help besseli" command), can
| there really be a non-commercial-use license attached with it?
Maybe they negotiated some different terms? If you are interested in
having updated versions of these routines in Octave, then please ask
the ACM if they would allow distribution under terms that are
compatible with the GPL.
Looking at the ACM copyright notice page:
http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html

I guess it is part (b) of #4 that is the reason it can't be in octave?

User may modify the Software and distribute that modified work to third
Post by John W. Eaton
parties provided that: (a) if posted separately, it clearly acknowledges
that it contains material copyrighted by ACM (b) no charge is associated
with such copies, (c) User agrees to notify ACM and the Author(s) of the
distribution, and (d) User clearly notifies secondary users that such
modified work is not the original Software.
Thanks for the info, otherwise. It is all much clearer now.

Kris
John W. Eaton
2008-02-07 06:15:35 UTC
Permalink
On 6-Feb-2008, Kris Kuhlman wrote:

| Looking at the ACM copyright notice page:
| http://www.acm.org/pubs/copyright_policy/softwareCRnotice.html
|
| I guess it is part (b) of #4 that is the reason it can't be in octave?
|
| User may modify the Software and distribute that modified work to third
| > parties provided that: (a) if posted separately, it clearly acknowledges
| > that it contains material copyrighted by ACM (b) no charge is associated
| > with such copies, (c) User agrees to notify ACM and the Author(s) of the
| > distribution, and (d) User clearly notifies secondary users that such
| > modified work is not the original Software.

I think 4b and 4c both conflict with the the terms of the GPL.

jwe
Brian Gough
2008-02-07 10:15:41 UTC
Permalink
At Wed, 06 Feb 2008 23:06:03 -0500,
Post by John W. Eaton
Maybe they negotiated some different terms? If you are interested in
having updated versions of these routines in Octave, then please ask
the ACM if they would allow distribution under terms that are
compatible with the GPL.
I think there are versions of the Amos Bessel functions in SLATEC,
which is public domain.
--
Brian Gough

GNU Scientific Library -
http://www.gnu.org/software/gsl/
John W. Eaton
2008-02-15 22:16:59 UTC
Permalink
On 7-Feb-2008, Brian Gough wrote:

| I think there are versions of the Amos Bessel functions in SLATEC,
| which is public domain.

Thanks. I see that Netlib's slatec directory does appear to have an
updated version of the Amos routines, but unfortunately they don't
solve the overflow problem that was originally reported.

jwe

Dmitri A. Sergatskov
2008-02-07 04:12:32 UTC
Permalink
Unless the Matlab compatibility is a must, you can use
a wrapper to GSL's library from octave-forge:


octave:1> bessel_Kn(24, 0.6)
ans = 4.5588e+34
octave:2> bessel_Kn_scaled(24, 0.6)
ans = 8.3068e+34

octave:5> help bessel_Kn_scaled
-- Loadable Function: Y = bessel_Kn_scaled (N, X)
-- Loadable Function: [Y, ERR] = bessel_Kn_scaled (...)
ERR contains an estimate of the absolute error in the value Y.

This function is from the GNU Scientific Library, see
`http://www.gnu.org/software/gsl/' for documentation.

/usr/libexec/octave/packages/gsl-1.0.3/x86_64-redhat-linux-gnu-api-v32/gsl_sf.oct


Sincerely,

Dmitri.
--
Loading...