Complex arg-trig functions

Stephen Montgomery-Smith stephen at missouri.edu
Sun Sep 16 21:00:19 UTC 2012


On 09/16/2012 02:53 PM, Bruce Evans wrote:
> On Mon, 17 Sep 2012, Bruce Evans wrote:
>
>> On Sun, 16 Sep 2012, Stephen Montgomery-Smith wrote:
>>
>>> On 09/16/2012 12:14 AM, Bruce Evans wrote:
>>>> On Sat, 15 Sep 2012, Stephen Montgomery-Smith wrote:
>>>>
>>>>> One more thing I would like an opinion on.
>>>>>
>>>>> In my code I check for |z| being small, and then use the
>>>>> approximations:
>>>>> casinh(z) = z
>>>>> cacos(z) = Pi - z
>>>>
>>>> Actually Pi/2 - z.
>>>>
>>>>> catanh(z) = z

>>> So all things being said and done, I am going to remove the use of
>>> these approximations.

>
> It gives the expected pessimizations, and unexpected accuracy improvements
> and unimprovements.  On amd64:

I got unexpected accuracy improvements as well!  I thought it might just 
be a coincidence, so I ignored it.

> @ @@ -313,5 +323,6 @@
> @              return (cpackf(copysignf(0, x), y+y));
> @          if (isinf(y))
> @ -            return (cpackf(copysignf(0, x), copysignf(m_pi_2, y)));
> @ +            return (cpackf(copysignf(0, x),
> @ +                copysignf(m_pi_2 + tiny, y)));
> @          if (x == 0)
> @              return (cpackf(x, y+y));
> @ @@ -320,9 +331,16 @@
> @ @      if (isinf(x) || isinf(y))
> @ -        return (cpackf(copysignf(0, x), copysignf(m_pi_2, y)));
> @ +        return (cpackf(copysignf(0, x), copysignf(m_pi_2 + tiny, y)));
> @ @      if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
> @          if ((int)(1+tiny)==1)
> @ -            return (cpackf(copysignf(real_part_reciprocal(ax, ay),
> x), copysignf(m_pi_2, y)));
> @ +            return (cpackf(
> @ +                copysignf(real_part_reciprocal(ax, ay), x),
> @ +                copysignf(m_pi_2 + tiny, y)));
> @ +
> @ +    /* XXX the numbers are related to sqrt(6 * FLT_EPSILON). */
> @ +    if (ax < 2048 * FLT_EPSILON && ay < 2048 * FLT_EPSILON)
> @ +        if ((int)ax==0 && (int)ay==0)
> @ +            return (z);
> @ @      if (ax == 1 && ay < FLT_EPSILON) {

I implemented all the m_pi_2 + tiny changes.

Let me still ponder the |z| being small issue.  Or you can put that code 
back in when it is committed.




More information about the freebsd-numerics mailing list