Patches to 1.8.0p4 to add Bessel functions for those that have 'em

[Mike Hall]

Here's some simple patches to configure.in, configure and math.c
to add the Bessel functions on those systems that have them
in their math library. They were made against 1.8.0 preview 4.

It's probably better to have one #define for all the Bessels,
and avoid the pollution of all the little #define's.

$ diff -p configure.in configure.in.orig
*** configure.in Fri Jul 25 01:03:32 2003
--- configure.in.orig Thu Jul 24 02:37:56 2003
*************** AC_CHECK_FUNCS(fmod killpg wait4 waitpid
*** 373,380 ****
setrgid setegid setregid setresgid issetugid pause lchown lchmod\
getpgrp setpgrp getpgid setpgid initgroups getgroups setgroups\
getpriority getrlimit dlopen sigprocmask sigaction _setjmp\
! setsid telldir seekdir fchmod mktime timegm cosh sinh tanh \
! j0 j1 y0 y1 jn yn)
AC_ARG_ENABLE(setreuid,
[ --enable-setreuid use setreuid()/setregid() according to need ev
en if obsolete.],
[use_setreuid=$enableval])
--- 373,379 ----
setrgid setegid setregid setresgid issetugid pause lchown lchmod\
getpgrp setpgrp getpgid setpgid initgroups getgroups setgroups\
getpriority getrlimit dlopen sigprocmask sigaction _setjmp\
! setsid telldir seekdir fchmod mktime timegm cosh sinh tanh)
AC_ARG_ENABLE(setreuid,
[ --enable-setreuid use setreuid()/setregid() according to need ev
en if obsolete.],
[use_setreuid=$enableval])


$ diff -p configure configure.orig
*** configure Fri Jul 25 01:02:15 2003
--- configure.orig Thu Jul 24 02:53:33 2003
*************** for ac_func in fmod killpg wait4 waitpid
*** 10311,10318 ****
setrgid setegid setregid setresgid issetugid pause lchown lchmod\
getpgrp setpgrp getpgid setpgid initgroups getgroups setgroups\
getpriority getrlimit dlopen sigprocmask sigaction _setjmp\
! setsid telldir seekdir fchmod mktime timegm cosh sinh tanh\
! j0 j1 jn y0 y1 yn
do
as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
echo "$as_me:$LINENO: checking for $ac_func" >&5
--- 10311,10317 ----
setrgid setegid setregid setresgid issetugid pause lchown lchmod\
getpgrp setpgrp getpgid setpgid initgroups getgroups setgroups\
getpriority getrlimit dlopen sigprocmask sigaction _setjmp\
! setsid telldir seekdir fchmod mktime timegm cosh sinh tanh
do
as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
echo "$as_me:$LINENO: checking for $ac_func" >&5


$ diff -p math.c math.c.orig
*** math.c Fri Jul 25 01:01:47 2003
--- math.c.orig Thu Jul 10 02:36:46 2003
*************** math_erfc(obj, x)
*** 291,347 ****
return rb_float_new(erfc(RFLOAT(x)->value));
}

- #ifdef HAVE_J0
- static VALUE math_j0(VALUE obj, VALUE x)
- {
- Need_Float(x);
- return rb_float_new(j0(RFLOAT(x)->value));
- }
- #endif
-
- #ifdef HAVE_J1
- static VALUE math_j1(VALUE obj, VALUE x)
- {
- Need_Float(x);
- return rb_float_new(j1(RFLOAT(x)->value));
- }
- #endif
-
- #ifdef HAVE_JN
- static VALUE math_jn(VALUE obj, VALUE n, VALUE x)
- {
- Need_Float(x);
- n = rb_Integer(n);
- return rb_float_new(jn(FIX2LONG(n), RFLOAT(x)->value));
- }
- #endif
-
- #ifdef HAVE_Y0
- static VALUE math_y0(VALUE obj, VALUE x)
- {
- Need_Float(x);
- return rb_float_new(y0(RFLOAT(x)->value));
- }
- #endif
-
- #ifdef HAVE_Y1
- static VALUE math_y1(VALUE obj, VALUE x)
- {
- Need_Float(x);
- return rb_float_new(y1(RFLOAT(x)->value));
- }
- #endif
-
- #ifdef HAVE_YN
- static VALUE math_yn(VALUE obj, VALUE n, VALUE x)
- {
- Need_Float(x);
- n = rb_Integer(n);
- return rb_float_new(yn(FIX2LONG(n), RFLOAT(x)->value));
- }
- #endif
-
void
Init_Math()
{
--- 291,296 ----
*************** Init_Math()
*** 388,410 ****

rb_define_module_function(rb_mMath, "erf", math_erf, 1);
rb_define_module_function(rb_mMath, "erfc", math_erfc, 1);
-
- #ifdef HAVE_J0
- rb_define_module_function(rb_mMath, "j0", math_j0, 1);
- #endif
- #ifdef HAVE_J1
- rb_define_module_function(rb_mMath, "j1", math_j1, 1);
- #endif
- #ifdef HAVE_JN
- rb_define_module_function(rb_mMath, "jn", math_jn, 2);
- #endif
- #ifdef HAVE_Y0
- rb_define_module_function(rb_mMath, "y0", math_y0, 1);
- #endif
- #ifdef HAVE_Y1
- rb_define_module_function(rb_mMath, "y1", math_y1, 1);
- #endif
- #ifdef HAVE_YN
- rb_define_module_function(rb_mMath, "yn", math_yn, 2);
- #endif
}
--- 337,340 ----

 

[Yukihiro Matsumoto]

Hi,

In message "Patches to 1.8.0p4 to add Bessel functions for those that have 'em"

|Here's some simple patches to configure.in, configure and math.c
|to add the Bessel functions on those systems that have them
|in their math library. They were made against 1.8.0 preview 4.
|
|It's probably better to have one #define for all the Bessels,
|and avoid the pollution of all the little #define's.

I want Ruby to be as portable as possible among platforms, so that I
decided when we add Bessel functions to Math module, we will prepare
substitution functions for those who don't have them.

If you (or others) come up with your version of Bessel functions (that
need not to be elegant) in a few days, I can merge them before 1.8.0.

matz.

 

[Josef 'Jupp' Schugt]

Saluton!

* Yukihiro Matsumoto; 2003-07-25, 18:41 UTC:

If you (or others) come up with your version of Bessel functions
(that need not to be elegant) in a few days, I can merge them
before 1.8.0.

Numerical Recipes in... (e.g. ...in C) has implementations of
different kinds of Bessel functions. I see no technical problem
porting them to Ruby but up to now I did not do that because I am
usure about wether such code can be released under Ruby's license or
not. Besides that it's nothing but a quick hack.

Gis,

Josef 'Jupp' Schugt

 

[Martin Weber]

Hi,

In message "Patches to 1.8.0p4 to add Bessel functions for those that have 'em"

|Here's some simple patches to configure.in, configure and math.c
|to add the Bessel functions on those systems that have them
|in their math library. They were made against 1.8.0 preview 4.
|
|It's probably better to have one #define for all the Bessels,
|and avoid the pollution of all the little #define's.

I want Ruby to be as portable as possible among platforms, so that I
decided when we add Bessel functions to Math module, we will prepare
substitution functions for those who don't have them.

If you (or others) come up with your version of Bessel functions (that
need not to be elegant) in a few days, I can merge them before 1.8.0.

Why not simply use those of BSD ? Okay, the originate from Sun (due
to the comment at the start of the files) but you can basically do
what you want as long as the comment stays in the source file, so
could include those in the ruby distribution ? *shrug* I mean why
reinvent the wheel when usable code is out there.

-martin

 

[Josef 'Jupp' Schugt]

Saluton!

* Martin Weber; 2003-07-25, 22:53 UTC:

Why not simply use those of BSD ? Okay, the originate from Sun (due
to the comment at the start of the files) but you can basically do
what you want as long as the comment stays in the source file, so
could include those in the ruby distribution ? *shrug* I mean why
reinvent the wheel when usable code is out there.

Out where exactly?

Gis,

Josef 'Jupp' Schugt

 

[Martin Weber]

Saluton!

* Martin Weber; 2003-07-25, 22:53 UTC:

Out where exactly?

E.g. in the netbsd source tree, e.g. j0 in src/lib/libm/src/ej0.c.
The other bessel funcs are lurking there (in this dir), too.

-Martin

 

[Josef 'Jupp' Schugt]

Saluton!

* Yukihiro Matsumoto; 2003-07-25, 18:41 UTC:

If you (or others) come up with your version of Bessel functions
(that need not to be elegant) in a few days, I can merge them
before 1.8.0.

When implementating special functions like Bessel ones, care must be
taken about limited precision of floating point numbers - a problem
alien to most programmers. Moreover the algorithms used to evaluate
special functions often depend on quite a number of constants each of
wich may for some reason or another be imprecise. For these reasons I
suggest *not* adding Bessel functions to Ruby 1.8. In my opinion we
should rather strive for including GSL (GNU Scientific Library)
support in core Ruby 2.0.

Gis,

Josef 'Jupp' Schugt

 

[Josef 'Jupp' Schugt]

Saluton!

* Guillaume Marcais; 2003-07-30, 13:57 UTC:

Isn't GSL GPL licensed? And that could cause problem to include it
into the core of Ruby.

I did not write 'include GSL' but include 'GSL support'. The GPL does
only apply to software released under it and software based on it
where 'based on it' means that you took the GPLd software and did
modify it.

What are you doing when you create Ruby bindings for GSL? You write a
wrapper for some library that implements functionality and API of
GSL. By *no* means do you require that GSL is the library that is
actually being used.

Therefore Ruby bindings to GSL are completely independent of GSL.
Once we have a Ruby way of using GSL we can still think about
creating a RLSL (Ruby's License Scientific Library).

Taking into account the typical use of special functions I don't see
the need for reinventing the wheel.

Is there already a Ruby binding to GSL?

Yes.

Gis,

Josef 'Jupp' Schugt

 

[Mike Hall]

I reworked the patches for the Bessel et.al. math functions into an extension.
Direct link: http://users.rcn.com/m3ha11/ruby/mathext-0.1.tar.gz
Tested on 1.8.0preview6 on OS X last night.

I don't know why this bugs me so much. I guess it's just a matter of:
if ruby has access to _some_ functions in 'libm', then it should have
access to all of them. That doesn't really make sense, I know.

I agree with Josef that if we were to rewrite our own Bessels (etc.),
we'd want someone who really knows what they are doing.
I re-read 'What Every Computer Scientist Should Know About Floating Point"
last night -- very scarey stuff!

 

[Hal E. Fulton]

----- Original Message -----
From: "Mike Hall" <[email protected]>
Newsgroups: comp.lang.ruby
To: "ruby-talk ML" <[email protected]>
Sent: Friday, August 01, 2003 12:13 PM
Subject: Re: Patches to 1.8.0p4 to add Bessel functions for those that have
'em

I agree with Josef that if we were to rewrite our own Bessels (etc.),
we'd want someone who really knows what they are doing.
I re-read 'What Every Computer Scientist Should Know About Floating Point"
last night -- very scarey stuff!

Hmm, I don't know that. Is it on the web?

Hal

 

[Josef 'Jupp' Schugt]

Saluton!

* Mike Hall; 2003-08-01, 20:19 UTC:

I agree with Josef that if we were to rewrite our own Bessels
(etc.), we'd want someone who really knows what they are doing. I
re-read 'What Every Computer Scientist Should Know About Floating
Point" last night -- very scarey stuff!

There is reason to scare. Even the most trivial computations are
nontrivial when one has a closer look at them:

1.upto(16) { |i|
puts 1.0 + (0.1 ** i) - 1.0 - (0.1 ** i)
}

on my machine results in:

8.32667268468867e-17
6.93889390390723e-18
-1.10371781159024e-16
-1.10317571050400e-17
6.55095281406476e-17
-8.22670162108899e-17
5.83866825033984e-17
-6.07747148815308e-17
8.27403705232185e-17
8.27403704456703e-18
8.27403704133586e-19
8.89005823404259e-17
-7.99277837359775e-17
-7.99277837359901e-18
1.10223024625156e-16
-1 e-16

The mathematical results of course all are 0. I did stop at 16
because we have reached what is called 'machine's epsilon'. This
value is defined as the largest number that can be added to 1 without
changing the value. This means that at this point

1.0 + (0.1 ** i) - 1.0 - (0.1 ** i)

becomes equivalent to

- (0.1 ** i)

Gis,

Josef 'Jupp' Schugt

 

[daz]

[...] Even the most trivial computations are
nontrivial when one has a closer look at them:

1.upto(16) { |i|
puts 1.0 + (0.1 ** i) - 1.0 - (0.1 ** i)
}

on my machine results in:

8.32667268468867e-17
[snipped some]
-1 e-16

The mathematical results of course all are 0.

Josef 'Jupp' Schugt

Just trollin', but is it right that Float#to_f
should be a no-op when it could have saved
doing 'to_s.to_f' (or other) below ?:


expect = [ nil,
8.32667268468867e-17,
6.93889390390723e-18,
-1.10371781159024e-16,
-1.10317571050400e-17,
6.55095281406476e-17,
]

1.upto(5) do |i|
ans = 1.0 + (0.1 ** i) - 1.0 - (0.1 ** i)
unless ans == expect
p [ans, expect]
puts 'Oops' unless Float === ans and Float === expect
if ans.to_s.to_f == expect # <<----#####
# ans.to_f short-circuits (just returns self)
puts "Not equal but equal"; puts
end
end
end


#-> [8.32667268468867e-17, 8.32667268468867e-17]
#-> Not equal but equal
#->
#-> [6.93889390390723e-18, 6.93889390390723e-18]
#-> Not equal but equal
#->
#-> [-1.10371781159024e-16, -1.10371781159024e-16]
#-> Not equal but equal
#->
#-> [-1.103175710504e-17, -1.103175710504e-17]
#-> Not equal but equal
#->
#-> [6.55095281406476e-17, 6.55095281406476e-17]
#-> Not equal but equal


daz

 

[furlan primus]

Hmm, I don't know that. Is it on the web?

yes, here is at least one source:

http://docs.sun.com/source/806-3568/ncg_goldberg.html

ciao,
furlan

--

http://thispaceavailable.uxb.net/blog/index.html

Q: How many surrealists does it take to screw in a lightbulb?
A: Fish

 

[Jean-Claude Arbaut]

- - The computational results are all "0" as well, they way to test
floating point computations for equality is to see if the
absolute value of their difference is less than EPS. Frankly,
I don't understand how this comes as a surprise to people.
You should not be allowed to use the / operator unless you
understand this.

- - Booker C. Bense

I don't understand the reason why you show such vanity, or why you
need to make fun of previous posts. It's a fact that numerical analysis
is a difficult domain, and unless you are able to explain each line
of code in Netlib (for example , you shouldn't "be allowed to"
post such things.
In fact it is even funny to read your answer. I don't think it's
a good idea to test equality that way. It will work in some cases,
whene your numbers are near 1. I'd rather use a "relative" error,
and sometimes it may be more complicated. It depends on what you
need to compare. Here it's stupid to test absolute difference, since
the goal was to show why you must not use "a!=b" with floating point
numbers.


To answer another post, "What every computer scientist..." is a chapter
of the "Sun Computation Guide", which can be found at docs.sun.com.
Those who are interested can visit William Kahan's home page, or
search "ieee754" with google. Other sources of information are Netlib,
at www.netlib.org (or www.netlib.no, much faster in Europe), or
GAMS at NIST (gams.nist.gov).
You could also have a look at Numerical Recipes (www.nr.com).
You will find good sources for special functions and many other at Netlib.