Complex arg-trig functions

Bruce Evans brde at optusnet.com.au
Wed Aug 15 17:13:32 UTC 2012


On Tue, 14 Aug 2012, Stephen Montgomery-Smith wrote:

> I was looking through the code e_acosh.c, and it made me realize I could get 
> a small fraction more ULP in catrig.c by making the replacements:
>
> 216c216
> < 			*rx = log1p(Am1 + sqrt(Am1*(A+1)));
> ---
>> 			*rx = log1p(Am1 + sqrt(2*Am1 + Am1*Am1));
> 282c282
> < 			*sqrt_A2my2 = sqrt(Amy*(A+y));
> ---
>> 			*sqrt_A2my2 = sqrt(2*y*Amy + Amy*Amy);
>
> I'm not quite sure if the second replacement makes much difference, but the 
> first replacement seemed quite effective.

This seems to be slightly worse.  In my tests, it makes little difference
to the peak error, but unimproves the number of correctly rounded cases
quite often.

1,5c1,4
< amd64 float prec, on 2**12 x 2**12 args:
< rcacos:max_er = 0x58460841 2.7585, avg_er = 0.317, #>=1:0.5 = 29084:255712
< rcacosh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
< rcasin:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38388:382508
< rcasinh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
---
> rcacos:max_er = 0x57352248 2.7252, avg_er = 0.317, #>=1:0.5 = 28694:256172
> rcacosh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.265, #>=1:0.5 = 107904:3459300
> rcasin:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38332:382056
> rcasinh:max_er = 0x5e1e45e6 2.9412, avg_er = 0.265, #>=1:0.5 = 107904:3459300

['<' is the old version, '>' the bew version]

17,20c16,19
< icacos:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
< icacosh:max_er = 0x58460841 2.7585, avg_er = 0.317, #>=1:0.5 = 29084:255712
< icasin:max_er = 0x5e1e45e6 2.9412, avg_er = 0.262, #>=1:0.5 = 85868:3413684
< icasinh:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38388:382508
---
> icacos:max_er = 0x5e1e45e6 2.9412, avg_er = 0.265, #>=1:0.5 = 107904:3459300
> icacosh:max_er = 0x57352248 2.7252, avg_er = 0.317, #>=1:0.5 = 28694:256172
> icasin:max_er = 0x5e1e45e6 2.9412, avg_er = 0.265, #>=1:0.5 = 107904:3459300
> icasinh:max_er = 0x631b8183 3.0971, avg_er = 0.209, #>=1:0.5 = 38332:382056
32,37c31,34
< 
< amd64 double prec, on 2**12 x 2**12 args:
< rcacos:max_er =     0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2394:125988
< rcacosh:max_er =      0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741860
< rcasin:max_er =     0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33296:99152
< rcasinh:max_er =      0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741796
---
> rcacos:max_er =     0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2374:125954
> rcacosh:max_er =      0xf8a 1.9424, avg_er = 0.258, #>=1:0.5 = 8396:2741812
> rcasin:max_er =     0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33312:99184
> rcasinh:max_er =      0xf8a 1.9424, avg_er = 0.258, #>=1:0.5 = 8396:2741748
42,45c39,42
< icacos:max_er =      0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741860
< icacosh:max_er =     0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2394:125988
< icasin:max_er =      0xf7d 1.9360, avg_er = 0.257, #>=1:0.5 = 612:2741796
< icasinh:max_er =     0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33296:99152
---
> icacos:max_er =      0xf8a 1.9424, avg_er = 0.258, #>=1:0.5 = 8396:2741812
> icacosh:max_er =     0x1b5a 3.4189, avg_er = 0.228, #>=1:0.5 = 2374:125954
> icasin:max_er =      0xf8a 1.9424, avg_er = 0.258, #>=1:0.5 = 8396:2741748
> icasinh:max_er =     0x15c5 2.7212, avg_er = 0.113, #>=1:0.5 = 33312:99184

The unimprovement on i386 is similar.  This is surprising for the float case,
since the expressions are evaluated in double precision.

Bruce


More information about the freebsd-numerics mailing list