atanhf(x) Slight Accuracy Improvements

atanhf(x) Slight Accuracy Improvements

[Damian McGuckin]

Hi all,

Currently this routine does not achieve accuracy < ULP across all of its 
domain. For half of its domain, it worst error exceeds that limit by 
about 10%.

By tweaking some algebra, this can be made more accurate. A comparison of 
the more accurate version and the original is noted below.

 	ATANF ... FROM ... TO     WORST ERROR  % > EPS/2   MEAN ERROR
 	--------------------------------------------------------------
 	Accurate:0.00000..0.17000 1.00023*EPS   0.76907%   0.01967*EPS
 	Original:0.00000..0.17000 1.09363*EPS   0.71616%   0.01941*EPS
 	Accurate:0.17000..0.55000 0.99151*EPS   5.52206%   0.21620*EPS
 	Original:0.17000..0.55000 1.12779*EPS   8.87471%   0.23922*EPS
 	Accurate:0.55000..1.00000 0.68151*EPS   1.34789%   0.19159*EPS
 	Original:0.55000..1.00000 0.68372*EPS   1.40897%   0.19207*EPS

I see a reduction in the worst error across the entire spectrum and reduce 
slightly the percentage exceeding 0.5*ULP in most cases. The mean error is 
much the same. I have yet to rework the double version.

However, across a subset of its argument range, namely

 	[0 ..(sqrt(2)-1)/(sqrt(2)+1)]

I cannot crack the 1.0*ULP barrier if the computation of the argument
reduction

 	f = 2*(y + (y*y)/(1-y))  <---- LITTLE PROBLEM

is done in single precision. The error in the 23rd bit causes me grief.

Doing that sole calculation in double precision and then storing it as a 
float brings the worst error to 0.99*ULP. I want to avoid any extended
precision.

Note that I used the same accuracy tweak for log1p as done in log2 to 
avoid the cancellation error seen in

 	f - (f*f)/2

Any suggestions are welcome on how to get one extra bit of accuracy in my 
calculation of 'LITTLE PROBLEM' above. Using the approach as seen in 'sq() 
in 'hypotf' makes things worse unless I am doing something really wrong.

Thanks - Damian

Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer

Re: atanhf(x) Slight Accuracy Improvements

[Damian McGuckin]

On Wed, 13 Mar 2019, Damian McGuckin wrote:

> By tweaking some algebra, this can be made more accurate. A comparison 
> of the more accurate version and the original is noted below.
>
> 	 ATANF ... FROM ... TO WORST ERROR % > EPS/2 MEAN ERROR
> 	 --------------------------------------------------------------
> 	 Accurate:0.00000..0.17000 1.00023*EPS 0.76907% 0.01967*EPS
> 	 Accurate:0.17000..0.55000 0.99151*EPS 5.52206% 0.21620*EPS

Sorted for the crucial segment.

More algebraic rearrangement and that 1.00023*EPS drops to 0.6667*EPS.

I will work on the second segment.

Regards - Damian

Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer

Re: atanhf(x) Slight Accuracy Improvements

[Damian McGuckin]

Here is the net nesult for the 32-bit version.

To date, the research has only involved FLT_EVAL_METHOD=0.

ATANH(x)  From  .. To      WORST ERROR  % ERROR>EPS/2  MEAN ERROR

Improved:0.00000..0.00012  0.04027*EPS   0.00000%      0.00301*EPS
Original:0.00000..0.00012  0.99741*EPS   0.39818%      0.00581*EPS

Improved:0.00012..0.00100  0.66667*EPS   8.86075%      0.21559*EPS
Original:0.00012..0.00100  0.95814*EPS   7.42942%      0.21211*EPS

Improved:0.00100..0.17157  0.55526*EPS   0.00605%      0.17453*EPS
Original:0.00100..0.17157  1.09363*EPS   8.51956%      0.23306*EPS

Improved:0.17157..0.18248  0.63539*EPS   0.41833%      0.18648*EPS
Original:0.17157..0.18248  1.07990*EPS   5.23646%      0.21706*EPS

Improved:0.18248..0.26667  0.83800*EPS   2.10627%      0.18376*EPS
Original:0.18248..0.26667  1.12385*EPS   8.71681%      0.23764*EPS

Improved:0.26667..0.55000  0.90005*EPS   3.97403%      0.20866*EPS
Original:0.26667..0.55000  1.12779*EPS   9.27245%      0.24201*EPS

Improved:0.55000..1.00000  0.68372*EPS   1.39651%      0.19197*EPS
Original:0.55000..1.00000  0.68372*EPS   1.40897%      0.19207*EPS

Does this sort of gain make it worthwhile pursing a double version?

Regards - Damian

Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer

Re: Re: atanhf(x) Slight Accuracy Improvements

[Szabolcs Nagy]

* Damian McGuckin <damianm@esi.com.au> [2019-03-17 22:41:53 +1100]:
> ATANH(x)  From  .. To      WORST ERROR  % ERROR>EPS/2  MEAN ERROR
> 
> Improved:0.00000..0.00012  0.04027*EPS   0.00000%      0.00301*EPS
> Original:0.00000..0.00012  0.99741*EPS   0.39818%      0.00581*EPS
> 
> Improved:0.00012..0.00100  0.66667*EPS   8.86075%      0.21559*EPS
> Original:0.00012..0.00100  0.95814*EPS   7.42942%      0.21211*EPS
> 
> Improved:0.00100..0.17157  0.55526*EPS   0.00605%      0.17453*EPS
> Original:0.00100..0.17157  1.09363*EPS   8.51956%      0.23306*EPS
> 
> Improved:0.17157..0.18248  0.63539*EPS   0.41833%      0.18648*EPS
> Original:0.17157..0.18248  1.07990*EPS   5.23646%      0.21706*EPS
> 
> Improved:0.18248..0.26667  0.83800*EPS   2.10627%      0.18376*EPS
> Original:0.18248..0.26667  1.12385*EPS   8.71681%      0.23764*EPS
> 
> Improved:0.26667..0.55000  0.90005*EPS   3.97403%      0.20866*EPS
> Original:0.26667..0.55000  1.12779*EPS   9.27245%      0.24201*EPS
> 
> Improved:0.55000..1.00000  0.68372*EPS   1.39651%      0.19197*EPS
> Original:0.55000..1.00000  0.68372*EPS   1.40897%      0.19207*EPS
> 
> Does this sort of gain make it worthwhile pursing a double version?

if you improve accuracy without performance loss that sounds
useful, but i cannot comment on it without seeing the patch.

Re: Re: atanhf(x) Slight Accuracy Improvements

[Damian McGuckin]

On Mon, 18 Mar 2019, Szabolcs Nagy wrote:

> if you improve accuracy without performance loss that sounds useful, but 
> i cannot comment on it without seeing the patch.

The main extra load is rather when computing

 	f = 2 * y + (2 * y * y) / (1 - y)

and
 	hfsq = (0.5 * f) * f

it is computed using a slightly different formula and then dissected like 
(f-hfsq) is currently handled in log2.c because of the cancellation error 
seen around sqrt(2), which in this case is when the absolute value 
of the argument 'x', i.e. 'y' from MUSL's atanh.c,, has a value of

 	sqrt(2) - 1
 	-----------
 	sqrt(2) + 1

I need to work on the code because currently it does not conform to the 
normal way MUSL is coded. That is somewhat because the C code was ported 
from a research project which is using the Chapel language and has some 
design consideration which are irrelevant to MUSL, and others which are 
are more limiting than what MUSL has. It might take me a few weeks because 
I have some other work priorities for the next 3 weeks and also need to 
document the research results for the Chapel project first.

Regards - Damian

Pacific Engineering Systems International, 277-279 Broadway, Glebe NSW 2037
Ph:+61-2-8571-0847 .. Fx:+61-2-9692-9623 | unsolicited email not wanted here
Views & opinions here are mine and not those of any past or present employer