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

Discussion in 'Ruby' started by Mike Hall, Jul 25, 2003.

  1. Mike Hall

    Mike Hall Guest

    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 ----
    Mike Hall, Jul 25, 2003
    #1
    1. Advertising

  2. Hi,

    In message "Patches to 1.8.0p4 to add Bessel functions for those that have 'em"
    on 03/07/25, Mike Hall <> writes:

    |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.
    Yukihiro Matsumoto, Jul 25, 2003
    #2
    1. Advertising

  3. 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
    --
    N'attribuez jamais à la malice ce que l'incompétence explique !
    -- Napoléon
    Josef 'Jupp' Schugt, Jul 25, 2003
    #3
  4. Mike Hall

    Martin Weber Guest

    On Sat, Jul 26, 2003 at 01:41:35AM +0900, Yukihiro Matsumoto wrote:
    > Hi,
    >
    > In message "Patches to 1.8.0p4 to add Bessel functions for those that have 'em"
    > on 03/07/25, Mike Hall <> writes:
    >
    > |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
    Martin Weber, Jul 25, 2003
    #4
  5. 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
    --
    N'attribuez jamais à la malice ce que l'incompétence explique !
    -- Napoléon
    Josef 'Jupp' Schugt, Jul 26, 2003
    #5
  6. Mike Hall

    Martin Weber Guest

    On Sun, Jul 27, 2003 at 07:08:09AM +0900, Josef 'Jupp' Schugt wrote:
    > 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?


    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
    Martin Weber, Jul 27, 2003
    #6
  7. 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
    --
    N'attribuez jamais à la malice ce que l'incompétence explique !
    -- Napoléon
    Josef 'Jupp' Schugt, Jul 29, 2003
    #7
  8. Saluton!

    * Guillaume Marcais; 2003-07-30, 13:57 UTC:
    > >In my opinion we
    > >should rather strive for including GSL (GNU Scientific Library)
    > >support in core Ruby 2.0.
    > >

    > 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
    --
    N'attribuez jamais à la malice ce que l'incompétence explique !
    -- Napoléon
    Josef 'Jupp' Schugt, Jul 30, 2003
    #8
  9. Mike Hall

    Mike Hall Guest

    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! :)
    Mike Hall, Aug 1, 2003
    #9
  10. ----- Original Message -----
    From: "Mike Hall" <>
    Newsgroups: comp.lang.ruby
    To: "ruby-talk ML" <>
    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

    --
    Hal Fulton
    Hal E. Fulton, Aug 1, 2003
    #10
  11. 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
    --
    N'attribuez jamais à la malice ce que l'incompétence explique !
    -- Napoléon
    Josef 'Jupp' Schugt, Aug 2, 2003
    #11
  12. Mike Hall

    daz Guest

    Float#to_f (was Re: Patches to 1.8.0p4 to add Bessel functions ...)

    "Josef 'Jupp' Schugt" <> wrote:
    >
    > [...] 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
    daz, Aug 2, 2003
    #12
  13. furlan primus, Aug 3, 2003
    #13
  14. <> wrote in message news:<bggl6e$d94$>...
    >
    > - - 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.
    Jean-Claude Arbaut, Aug 4, 2003
    #14
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Karl Seguin
    Replies:
    0
    Views:
    403
    Karl Seguin
    Dec 8, 2005
  2. BinnuChowdary
    Replies:
    1
    Views:
    522
    Swanand Mokashi
    May 1, 2006
  3. BinnuChowdary
    Replies:
    0
    Views:
    404
    BinnuChowdary
    May 2, 2006
  4. BinnuChowdary
    Replies:
    1
    Views:
    539
    =?UTF-8?B?R8O2cmFuIEFuZGVyc3Nvbg==?=
    May 2, 2006
  5. Replies:
    0
    Views:
    257
Loading...

Share This Page