Skip to content

Conversation

@ZERICO2005
Copy link
Contributor

@ZERICO2005 ZERICO2005 commented Nov 8, 2024

The accuracy of log1pf(float x) has been fixed. It will return x when x is small for 12.0bits of precision, or it can return x - 0.5f * (x * x) when x is small for 16.0bits of precision. It falls-back to the naive return logf(x + 1.0f) when x is large.

…re provided, one with 12.0bits of precision, and another with 16.0bits of precision. Tested across all inputs on x86_64
@adriweb adriweb linked an issue Nov 8, 2024 that may be closed by this pull request
#include <math.h>

float log1pf(float x)
float log1pf(const float x)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not needed or in the standard.

@parisseb
Copy link
Contributor

parisseb commented Nov 12, 2024

This fix (862e260) should be improved, because it generates a relative error that is smaller than 2^(-13), and that's far from 2^-24. It could be improved like this

	if (fabsf(x) < 2.44140625e-4f /* 0x1.0p-12f */) {
		return x*(1-ld_exp(x,-1));
	}

but for values just above 2^-12, using log(1+x) will still lead to low relative precision (about 2^-12), because computing 1+x will loose the last 12 bits of x.

Here is an improved log1p version, using Pade approximant, 4 multiplications and one division, relative error is bounded by about 2^-21=5e-7

float log1pf(float x){
  if (fabs(x)<=0.125){
    // pade(series(ln(1+x),x=0,6,polynom),x,5,3)
    // (-57*x**2-90*x)/(x**3-21*x**2-102*x-90)
    // relative error less than 1e-7
       return x*(57*x+90)/(((21-x)*x+102)*x+90);
  }
  // relative error about 2^-21 if abs(x) is just above 0.125
  return logf(1 + x);
}

FYI, I have also re-implemented exp/ln/sin/cos/tan/atan in order to speed up a little bit plotting inside KhiCAS (between 10% and 20%), because the toolchain source (from Zilog) seems to be adapted for 14 or 15 digits relative precision, about twice the 24 bits relative precision of floats. The new versions are available in the archive containing all changes I had to make on the toolchain for optimal compilation of KhiCAS:
https://www-fourier.univ-grenoble-alpes.fr/~parisse/ti/celibc.tgz
(files should be copied in the toolchain src/libc directory, beware that it also contains my modified version of malloc to access 2 heaps, one of them being the upper LCD part).

@adriweb
Copy link
Member

adriweb commented Nov 12, 2024

Indeed, very nice!
I'll commit that version, then.

CleanShot 2024-11-12 at 14 18 56@2x

adriweb pushed a commit that referenced this pull request Nov 12, 2024
Improved log1p version, using Pade approximant,
4 multiplications and one division, relative error is bounded
by about 2^-21=5e-7.

From #516 (comment)
@adriweb
Copy link
Member

adriweb commented Nov 12, 2024

I'm closing this PR since it's now obsolete.

@adriweb adriweb closed this Nov 12, 2024
@ZERICO2005
Copy link
Contributor Author

ZERICO2005 commented Nov 12, 2024

because the toolchain source (from Zilog) seems to be adapted for 14 or 15 digits relative precision

Yeah I noticed that too. We could keep those routines around should we decide on implementing Float64 (long double)

@ZERICO2005 ZERICO2005 deleted the fixed_log1pf branch April 16, 2025 04:14
@ZERICO2005 ZERICO2005 restored the fixed_log1pf branch April 16, 2025 04:14
@ZERICO2005 ZERICO2005 deleted the fixed_log1pf branch April 16, 2025 04:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

log1pf is inaccurate for small x

3 participants