Which function in libm are portable enough to assume that they will be on all platforms?
I'm thinking about writing improved complex versions of existing Scheme functions, but some of these use floating-point versions of functions that are not now in Scheme, or in Gambit, e.g., sinh.
So I decided to look through what's declared by <math.h> on Ubuntu 13.10 (which I happen to be sitting in front of) and I include a test file with all the declared "double" functions (there are also "float" versions of these functions and "long double" versions of these functions).
So, my question is, would it be reasonable to assume that a good quality sinh (or other routines) is available on all platforms we want Gambit to run on, so we could define flsinh and use it elsehwere? (I wouldn't even want a general, complex sinh at first.)
This would really make the job of writing complex library code in Scheme easier.
Brad
Afficher les réponses par date
On Oct 30, 2013, at 4:04 PM, Bradley Lucier lucier@math.purdue.edu wrote:
I'm thinking about writing improved complex versions of existing Scheme functions, but some of these use floating-point versions of functions that are not now in Scheme, or in Gambit, e.g., sinh.
So I decided to look through what's declared by <math.h> on Ubuntu 13.10 (which I happen to be sitting in front of) and I include a test file with all the declared "double" functions (there are also "float" versions of these functions and "long double" versions of these functions).
So, my question is, would it be reasonable to assume that a good quality sinh (or other routines) is available on all platforms we want Gambit to run on, so we could define flsinh and use it elsehwere? (I wouldn't even want a general, complex sinh at first.)
This would really make the job of writing complex library code in Scheme easier.
Brad
It is pretty easy to add flsinh, etc to Gambit and I can do that if needed. Do you think the quality of the C sinh is higher than the alternative (i.e. using the existing transcendental functions to do the same computation)? Can you show a few examples where using sinh is a win?
Marc
On Oct 30, 2013, at 5:32 PM, Marc Feeley wrote:
It is pretty easy to add flsinh, etc to Gambit and I can do that if needed. Do you think the quality of the C sinh is higher than the alternative (i.e. using the existing transcendental functions to do the same computation)? Can you show a few examples where using sinh is a win?
Our code fails some fairly simple tests:
;(asin 1234000000.+0.i) got 0.+inf.0i, but expected 1.5707963267949+21.62667394298955i ;(asin -1234000000.-0.i) got -0.-inf.0i, but expected -1.5707963267949+21.62667394298955i ;(acos 1234000000.+0.i) got 1.5707963267948966-inf.0i, but expected 0.-21.62667394298955i ;(acos -1234000000.-0.i) got 1.5707963267948966+inf.0i, but expected 3.14159265358979-21.62667394298955i
This is from a modified R6RS test file on Clinger's web site. Looking at the complex asin and acos code shows that it is naive, does not try to avoid spurious underflows, etc.
Kahan in his paper "Much ado about the sign of nothing" gives algorithms for accurate and stable computations of complex functions (this was the basis for the routines in sbcl) but they use functions listed in math.h but not in standard scheme.
So my question is---is it reasonable to import these floating-point functions into Gambit to make better versions of the standard functions?
Brad
On Oct 30, 2013, at 11:30 PM, Bradley Lucier lucier@math.purdue.edu wrote:
On Oct 30, 2013, at 5:32 PM, Marc Feeley wrote:
It is pretty easy to add flsinh, etc to Gambit and I can do that if needed. Do you think the quality of the C sinh is higher than the alternative (i.e. using the existing transcendental functions to do the same computation)? Can you show a few examples where using sinh is a win?
Our code fails some fairly simple tests:
;(asin 1234000000.+0.i) got 0.+inf.0i, but expected 1.5707963267949+21.62667394298955i ;(asin -1234000000.-0.i) got -0.-inf.0i, but expected -1.5707963267949+21.62667394298955i ;(acos 1234000000.+0.i) got 1.5707963267948966-inf.0i, but expected 0.-21.62667394298955i ;(acos -1234000000.-0.i) got 1.5707963267948966+inf.0i, but expected 3.14159265358979-21.62667394298955i
This is from a modified R6RS test file on Clinger's web site. Looking at the complex asin and acos code shows that it is naive, does not try to avoid spurious underflows, etc.
Kahan in his paper "Much ado about the sign of nothing" gives algorithms for accurate and stable computations of complex functions (this was the basis for the routines in sbcl) but they use functions listed in math.h but not in standard scheme.
So my question is---is it reasonable to import these floating-point functions into Gambit to make better versions of the standard functions?
Brad
If you define
(define (flsinh x) (fl* 0.5 (fl- (flexp x) (flexp (fl- x)))))
(define (flcosh x) (fl* 0.5 (fl+ (flexp x) (flexp (fl- x)))))
(define (fltanh x) (let ((a (flexp (fl* 2.0 x)))) (fl/ (fl- a 1.0) (fl+ a 1.0))))
can you use those functions to improve the precision of complex asin and acos? Can these Scheme definitions be improved? And how good is the precision when compared to the C sinh and cosh?
Marc
On Oct 31, 2013, at 8:16 AM, Marc Feeley wrote:
If you define
(define (flsinh x) (fl* 0.5 (fl- (flexp x) (flexp (fl- x)))))
(define (flcosh x) (fl* 0.5 (fl+ (flexp x) (flexp (fl- x)))))
(define (fltanh x) (let ((a (flexp (fl* 2.0 x)))) (fl/ (fl- a 1.0) (fl+ a 1.0))))
can you use those functions to improve the precision of complex asin and acos?
No
Can these Scheme definitions be improved?
Yes
And how good is the precision when compared to the C sinh and cosh?
Crap.
Marc, I'm too busy right now to prepare a lecture on this stuff.
Brad
On Oct 31, 2013, at 9:40 AM, Bradley Lucier lucier@math.purdue.edu wrote:
On Oct 31, 2013, at 8:16 AM, Marc Feeley wrote:
If you define
(define (flsinh x) (fl* 0.5 (fl- (flexp x) (flexp (fl- x)))))
(define (flcosh x) (fl* 0.5 (fl+ (flexp x) (flexp (fl- x)))))
(define (fltanh x) (let ((a (flexp (fl* 2.0 x)))) (fl/ (fl- a 1.0) (fl+ a 1.0))))
can you use those functions to improve the precision of complex asin and acos?
No
Can these Scheme definitions be improved?
Yes
And how good is the precision when compared to the C sinh and cosh?
Crap.
Marc, I'm too busy right now to prepare a lecture on this stuff.
Brad
I'm asking these questions to respond to your initial question:
On Oct 30, 2013, at 11:30 PM, Bradley Lucier lucier@math.purdue.edu wrote:
So my question is---is it reasonable to import these floating-point functions into Gambit to make better versions of the standard functions?
I think it is "reasonable" if both
1) most C libm implementations have good precision sinh, etc (where "good" means better than the above portable definitions)
2) in case the C libm sinh, etc have poor precision, the use of the above portable definitions in the definition of complex asin and acos gives no worse precision than the current implementation in Gambit.
Marc
On Oct 31, 2013, at 10:40 AM, Marc Feeley wrote:
On Oct 31, 2013, at 9:40 AM, Bradley Lucier lucier@math.purdue.edu wrote:
On Oct 31, 2013, at 8:16 AM, Marc Feeley wrote:
If you define
(define (flsinh x) (fl* 0.5 (fl- (flexp x) (flexp (fl- x)))))
(define (flcosh x) (fl* 0.5 (fl+ (flexp x) (flexp (fl- x)))))
(define (fltanh x) (let ((a (flexp (fl* 2.0 x)))) (fl/ (fl- a 1.0) (fl+ a 1.0))))
can you use those functions to improve the precision of complex asin and acos?
No
Can these Scheme definitions be improved?
Yes
Here's an example for sinh:
If x < 0, return -sinh(-x)
There are numbers 0 < x_0 < x_1 < x_2 < x_3 < +inf.0 such that
If 0<x<x_0 then use a power series about 0 (avoids catastrophic cancelation)
If x_0 < x < x_1 return (e^x-e^{-x})/2
If x_1 < x < x_2 return e^x/2
If x_2< x < x_3 return e^{x-1}*(e/2) (avoids spurious overflow)
If x_3 < x return +inf.0
And that's if you're satisfied with an error of up to about 2 ulps.
And you have to figure out what x_i, i=0,1,2,3 are and what the power series is.
This is somewhat what fdlibm, glibc, eglibc, ..., Mac OS's library, ... do. Only they probably use better formulas than I just put up there.
I don't know that all platforms do this. But I think "most" do. And even the ones who don't, probably do as well as a quick and dirty definition like you gave above for sinh does.
And how good is the precision when compared to the C sinh and cosh?
Crap.
Marc, I'm too busy right now to prepare a lecture on this stuff.
Brad
I'm asking these questions to respond to your initial question:
On Oct 30, 2013, at 11:30 PM, Bradley Lucier lucier@math.purdue.edu wrote:
So my question is---is it reasonable to import these floating-point functions into Gambit to make better versions of the standard functions?
I think it is "reasonable" if both
- most C libm implementations have good precision sinh, etc (where "good" means better than the above portable definitions)
Yes, they do. "Most" includes linux, mac os x. I don't know what else it includes. I'd like to know about iOS and Android.
- in case the C libm sinh, etc have poor precision, the use of the above portable definitions in the definition of complex asin and acos gives no worse precision than the current implementation in Gambit.
I'm sorry, I don't understand this sentence.
On Oct 31, 2013, at 11:01 AM, Bradley Lucier wrote:
On Oct 31, 2013, at 10:40 AM, Marc Feeley wrote:
- most C libm implementations have good precision sinh, etc (where "good" means better than the above portable definitions)
Yes, they do. "Most" includes linux, mac os x. I don't know what else it includes. I'd like to know about iOS and Android.
And the BSDs, of course.
Brad
On Oct 31, 2013, at 11:01 AM, Bradley Lucier lucier@math.purdue.edu wrote:
- in case the C libm sinh, etc have poor precision, the use of the above portable definitions in the definition of complex asin and acos gives no worse precision than the current implementation in Gambit.
I'm sorry, I don't understand this sentence.
In your (future) implementation of complex asin and acos which is based on flsinh, etc, I am asking if the precision of complex asin and acos will be no worse than what it is currently if flsinh, etc are defined using the portable naive implementation.
Marc
On 10/31/2013 11:08 AM, Marc Feeley wrote:
In your (future) implementation of complex asin and acos which is based on flsinh, etc, I am asking if the precision of complex asin and acos will be no worse than what it is currently if flsinh, etc are defined using the portable naive implementation. Marc
Marc:
I still can't parse the sentence, sorry, my broken leg isn't healing properly and I can't always think straight right now.
To be more concrete, let me list the floating-point functions that Kahan uses in his implementations:
The following we have:
sqrt arctan ln tan
The following are trivial (maybe we already have them, I haven't checked).
copysign logb scaleb
The following we don't have:
ln1p ; ln1p(x)=ln(1+x) arcsinh sinh
These floating-point functions need to be computed accurately for the complex functions to be computed accurately.
There is already an algorithm for something like ln1p in the last case of the cond in exact-log in the definition of ##log, and we may need to have Scheme versions of (ln1p, arcsinh, and sinh) to calculate accurate values of the complex functions for exact arguments (Kahan considers only floating-point arguments, so there are limits on how big and how small they can be).
So I need only flln1p, flasinh, and flsinh. Which is not so much. (But which may be more than what Windows provides, according to Joe's email.)
Brad
On Thu, Oct 31, 2013 at 04:21:27PM -0400, Bradley Lucier wrote:
The following we don't have:
ln1p ; ln1p(x)=ln(1+x) arcsinh sinh
So I need only flln1p, flasinh, and flsinh. Which is not so much. (But which may be more than what Windows provides, according to Joe's email.)
Well, according to that page, Windows has flsinh, but lacks ln1p and asinh. flasinh can be computed with a logarithm, and ln1p can be implemented fairly simply as well (see http://www.johndcook.com/cpp_log_one_plus_x.html).
-Joe
On 10/31/2013 04:34 PM, Joe Doyle Ardent wrote:
On Thu, Oct 31, 2013 at 04:21:27PM -0400, Bradley Lucier wrote:
The following we don't have:
ln1p ; ln1p(x)=ln(1+x) arcsinh sinh
So I need only flln1p, flasinh, and flsinh. Which is not so much. (But which may be more than what Windows provides, according to Joe's email.)
Well, according to that page, Windows has flsinh, but lacks ln1p and asinh. flasinh can be computed with a logarithm, and ln1p can be implemented fairly simply as well (see http://www.johndcook.com/cpp_log_one_plus_x.html).
That code has a relative error of about 1e-12, which is not good enough for our purposes. I suppose we could write a few helper functions for Windows that may not have optimal accuracy or range, but that's another pain ...
Brad
On Oct 31, 2013, at 4:21 PM, Bradley Lucier lucier@math.purdue.edu wrote:
On 10/31/2013 11:08 AM, Marc Feeley wrote:
In your (future) implementation of complex asin and acos which is based on flsinh, etc, I am asking if the precision of complex asin and acos will be no worse than what it is currently if flsinh, etc are defined using the portable naive implementation. Marc
Marc:
I still can't parse the sentence, sorry, my broken leg isn't healing properly and I can't always think straight right now.
To be more concrete, let me list the floating-point functions that Kahan uses in his implementations:
The following we have:
sqrt arctan ln tan
The following are trivial (maybe we already have them, I haven't checked).
copysign logb scaleb
The following we don't have:
ln1p ; ln1p(x)=ln(1+x) arcsinh sinh
These floating-point functions need to be computed accurately for the complex functions to be computed accurately.
There is already an algorithm for something like ln1p in the last case of the cond in exact-log in the definition of ##log, and we may need to have Scheme versions of (ln1p, arcsinh, and sinh) to calculate accurate values of the complex functions for exact arguments (Kahan considers only floating-point arguments, so there are limits on how big and how small they can be).
So I need only flln1p, flasinh, and flsinh. Which is not so much. (But which may be more than what Windows provides, according to Joe's email.)
Brad
Given that the hyperbolic functions are interesting on their own, I'm enclined to add them to Gambit regardless of their use for implementing complex asin and acos. Then we can experiment to see if it is worthwhile using them for complex asin and acos.
Can you give me a consistent set of the math functions I should add? If sinh is added then cosh and tanh should be added as well, etc
Marc
On Wed, Oct 30, 2013 at 11:30:47PM -0400, Bradley Lucier wrote:
So my question is---is it reasonable to import these floating-point functions into Gambit to make better versions of the standard functions?
Here's the POSIX standard for math.h, which includes 32 and 64-bit versions of the trigonometric functions, including the hyperbolic ones:
http://pubs.opengroup.org/onlinepubs/9699919799/basedefs/math.h.html
The Android version of that C header file is compliant (http://mobilepearls.com/labs/native-android-api/include/math.h). I can't find the iOS version, but I assume it's also compliant, as it'll likely be the same as the system one for OS X. Windows seems to lack the inverse hyperbolic functions, though:
http://www.johndcook.com/math_h.html
I hope this helps.
-Joe
-
Bradley Lucier -
Joe Doyle Ardent -
Marc Feeley