complex arc-trig etc
Stephen Montgomery-Smith
stephen at missouri.edu
Mon Dec 17 13:30:49 UTC 2012
On 12/17/2012 02:38 AM, Bruce Evans wrote:
> On Sun, 16 Dec 2012, Stephen Montgomery-Smith wrote:
>
>> Hey guys, my complex arc-trig functions at
>> http://people.freebsd.org/~stephen/ have been sitting there a long time.
>> Anyone want to commit them?
>
> I was sort of waiting for the next rounds of changes:
>
> 1. update the template to give similar code for float and long double
> 2. make it work for i386 using ENTERI()
> 3. test (2)
> 4. fix some style bugs. It would be easy for me to run it through
> indent(1) and get large diffs, but I want you to do it.
> 5. decide what to do for atanhl(). I would prefer a version that is a
> direct translation of e_atanh.c including style bugs and the filename
> prefix but not the comment, unless you change the algorithm.
> 6. decide what to do with the template and/or the generation from it
> that strips comments.
>
> I haven't done (2), (3) or (6) for clog*() either. For (6), instead of
> a template and stripped comments, I have all the functions in the same
> file and all comments duplicated, and have to split the file manually
> to compare the functions.
>
> I only have these local patches. They correspond to (1) and (5). Also,
> keep up with the API change in LD80C().
>
> % --- /home/stephen/public_html/catrigf.c 2012-09-22
> 21:13:50.000000000 +0000
> % +++ catrigf.c 2012-09-22 21:35:51.287614000 +0000
> % @@ -353,12 +353,7 @@
> % }
> % % - if (ax == 1 && ay < FLT_EPSILON) {
> % -#if 0
> % - if (ay > 2*FLT_MIN)
> % - rx = - logf(ay/2) / 2;
> % - else
> % -#endif
> % - rx = - (logf(ay) - m_ln2) / 2;
> % - } else
> % + if (ax == 1 && ay < FLT_EPSILON)
> % + rx = - (logf(ay) - m_ln2) / 2;
> % + else
> % rx = log1pf(4*ax / sum_squares(ax-1, ay)) / 4;
> % % --- /home/stephen/public_html/catrigl.c 2012-09-22
> 21:14:24.000000000 +0000
> % +++ catrigl.c 2012-09-28 19:39:47.799585000 +0000
> % @@ -60,6 +60,6 @@
> % #if LDBL_MANT_DIG == 64
> % static const union IEEEl2bits
> % -um_e = LD80C(0xadf85458a2bb4a9b, 1, 0,
> 2.71828182845904523536e0L),
> % -um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 0,
> 6.93147180559945309417e-1L);
> % +um_e = LD80C(0xadf85458a2bb4a9b, 1, 2.71828182845904523536e+0L),
> % +um_ln2 = LD80C(0xb17217f7d1cf79ac, -1, 6.93147180559945309417e-1L);
> % #define m_e um_e.e
> % #define m_ln2 um_ln2.e
> % @@ -348,5 +348,5 @@
> % % if (y == 0 && ax <= 1)
> % - return (cpackl(atanhl(x), y)); /* XXX need atanhl() */
> % + return (cpackl(atanh(x), y)); /* XXX need atanhl() */
> % % if (x == 0)
> % @@ -369,12 +369,7 @@
> % }
> % % - if (ax == 1 && ay < LDBL_EPSILON) {
> % -#if 0
> % - if (ay > 2*LDBL_MIN)
> % - rx = - logl(ay/2) / 2;
> % - else
> % -#endif
> % - rx = - (logl(ay) - m_ln2) / 2;
> % - } else
> % + if (ax == 1 && ay < LDBL_EPSILON)
> % + rx = - (logl(ay) - m_ln2) / 2;
> % + else
> % rx = log1pl(4*ax / sum_squares(ax-1, ay)) / 4;
> %
>
> You also need log*l() be committed before catrigl() can work.
>
> Bruce
>
>
OK, I didn't realize it needed this much work (which I admit isn't that
much but...). So maybe I'll sit on these things for a while. I am not
swamped like Steve is. But I do have other projects that should take
priority for now.
Thanks, Stephen
More information about the freebsd-numerics
mailing list