My b^x approximation

[juha_p]

I'm not sure what you're referring to as "extension of the base function", maybe because I didn't read the entire thread, also not sure to which phenomenon you're referring to. However, since your relative error seems to increase towards the right, I would assume that you did the basic Remez with equal absolute error. You can do a equal relative error Remez, for functions passing through 0, e.g. for 2^x-1 by approximatig f(x)/x with an equal absolute error and then multiplying the result by x (that's IIRC what 2DaT did).

With "extension" I mean as like
- mystran's tanh(x) approximation suggestion (tanh(sqrt(x))/sqrt(x)) few pages earlier here,
- 2Dat's tanh(0.5*ln(2)*x) in his tanh() -thread,
- using sqrt() in my example here and
-OP's suggestion to use 0.25x as input for pade approximation with result raised to ^4 on page 2 ...
(and many more which ones I have bumped to).

and with "base function" the standard function found in math library ... exp(x) as for an example.

Many of those "extensions" may manipulate just the input value but as my example shows, taking sqrt() from exp() before approximation improves the accuracy (inside selected range) a bit. I've tested this with many approximation methods.

What comes to error measures , sollya gives these :

Code: Select all

f=exp(x);
p=remez(f,5,[-0.5;0.5]);
accurateinfnorm(p-f,[-0.5;0.5],20); 6.85118720866739749908447265625e-7
dirtyinfnorm(f-p,[-0.5;0.5]); 6.8511871736519530567224004924030262426642803963291e-7

Code: Select all

f=sqrt(exp(x));
p=remez(f,5,[-0.5;0.5]);
accurateinfnorm(p^2-exp(x),[-0.5;0.5],20); 2.72815441348939202725887298583984375e-8
dirtyinfnorm(exp(x)-p^2,[-0.5;0.5]); 2.72815398477066347610026190423265159037885082468976e-8

[Z1202]

IIRC I tried doing the ^2 trick when building a smooth approximation in another thread (or maybe when doing Remez, not sure anymore). That is, I was building an approximation with the intention of squaring it afterwards. I ended up with having pretty much the same precision for the entire thing for the same number of operations (approximately), so I gave up on the idea. Maybe you have better luck.

[juha_p]

Took further this sqrt() idea:

Code: Select all

f=sqrt(sqrt(exp(x)));
p=remez(f,4,[-0.5;0.5]);
accurateinfnorm(p^2^2-exp(x),[-0.5;0.5],20); 9.25762151382514275610446929931640625e-8
dirtyinfnorm(exp(x)-p^2^2,[-0.5;0.5]); 9.2576145658543045940128221683879222288763569773956e-8

Test functions (iPowxF() given earlier):
Degree 4:

Code: Select all

float rem4_expx4(float x){
	float z = floor(fmaf(x, 0x1.715476p+0, 0x1p-1));
	x -= z * 0x1.62e43p-1;
	int n = (int)z;
	float i = iPowxF(n);
	z =             0x1.55a501da39dbep-13;
	z = fmaf (z, x, 0x1.55aab322542bep-9);
	z = fmaf (z, x, 0x1.ffffed807e51p-6);
	z = fmaf (z, x, 0x1.ffffeaa7df52p-3);
	z = fmaf (z, x, 0x1.000000016c4b8p+0 );
	z *= z;
	return z * z;
}

and for degree 5 :

Code: Select all

f=sqrt(sqrt(exp(x)));
p=remez(f,5,[-0.5;0.5]);
accurateinfnorm(p^2^2-exp(x),[-0.5;0.5],20); 9.642189269243317539803683757781982421875e-10
dirtyinfnorm(exp(x)-p^2^2,[-0.5;0.5]); 9.6421744607584857063507686998172439094192596972843e-10

Test function (iPowxF() given earlier):

Code: Select all

float rem5_expx3(float x){
	float z = floor(fmaf(x, 0x1.715476p+0, 0x1p-1));
	x -= z * 0x1.62e43p-1;
	int n = (int)z;
	float i = iPowxF(n);
	z =             0x1.1145187b2556cp-17;
	z = fmal (z, x, 0x1.55999f42cb656p-13);
	z = fmal (z, x, 0x1.55554d9c05f35p-9);
	z = fmal (z, x, 0x1.fffff331cac03p-6);
	z = fmal (z, x, 0x1.000000009c1cbp-2);
	z = fmal (z, x, 0x1.00000000b621bp+0 );
	z *= z;
	return i * (z * z);
}

For comparison sake, plain exp(x) approximations:

Code: Select all

f=exp(x);
p=remez(f,4,[-0.5;0.5]);
accurateinfnorm(p-f,[-0.5;0.5],20); 1.647486351430416107177734375e-5
dirtyinfnorm(f-p,[-0.5;0.5]); 1.64748398340958573461873447323545505310736597160375e-5

1.00000136866487945431860995861874027267927111112956 + x * (0.999835902322769575719961607133190530306073158231 + x * (0.49992890794832408296427273635089544087696275115881 + x * (0.16928703594022782324130826690835672248935682575477 + x * 4.227791287072496266475229913932607438936510686987e-2)))

p=remez(f,5,[-0.5;0.5]);
accurateinfnorm(p-f,[-0.5;0.5],20); 6.85118720866739749908447265625e-7
dirtyinfnorm(f-p,[-0.5;0.5]); 6.8511871736519530567224004924030262426642803963291e-7

1.0000006833767667815164248059635618527935473425416 + x * (1.0000005857099931532709398461242779713113216546197 + x * (0.49995083513505615854282041532403502889200053619223 + x * (0.16665184176491087696465278790969569234340066213279 + x * (4.2190206834121582060871386383950326075133965985979e-2 + x * 8.433037380381611723971355540499581034978321121093e-3))))

By sollya, degree 8 gives accuracy of around 5e-16 when using this double sqrt() method. Am I in the forest with this?

[2DaT]

You need to consider that multiply is equivalent to the fma operation, therefore you need to compare it to +1 polynomial degree.

Also, e^x core polynomial is range-reduced to [-ln(2)/2;ln(2)/2] not [-0.5,0.5].

[mystran]

Also, e^x core polynomial is range-reduced to [-ln(2)/2;ln(2)/2] not [-0.5,0.5].

...which is equivalent to range reducing 2^x to [-0.5, 0.5] and you really always(!) want to approximate 2^x instead because you really want to handle the integer part by bumping the floating point exponent...

...and then you probably want to have this exp2() function (ie. base-2 exponential) available directly (rather than buried inside a function that only does base-e), because whenever you compute powers of anything other than base-e you waste a multiply converting to base-e only to convert it back to base-2... and in the common case of compute 2^x you waste two multiplies mapping back and forth... so there's very few (approximately sqrt(-1) as whatever reason you think you have is probably imaginary) good reasons to even bother having an exp() approximation if you can have a base-2 exponential instead, because you can only ever just lose cycles.

[Z1202]

I actually wouldn't mind having (2^x-1)/x, as this kind of expression pops up occasionally and it's a waste to do a comparison to zero (or actually to a range around zero) plus a division if you have a routine which computes 2^x-1 or (yuck!) 2^x instead

OTOH, if you're not limited to a small range, then you do need to actually compute (2^x-1)/x.

[2DaT]

Also, e^x core polynomial is range-reduced to [-ln(2)/2;ln(2)/2] not [-0.5,0.5].

...which is equivalent to range reducing 2^x to [-0.5, 0.5] and you really always(!) want to approximate 2^x instead because you really want to handle the integer part by bumping the floating point exponent...

...and then you probably want to have this exp2() function (ie. base-2 exponential) available directly (rather than buried inside a function that only does base-e), because whenever you compute powers of anything other than base-e you waste a multiply converting to base-e only to convert it back to base-2... and in the common case of compute 2^x you waste two multiplies mapping back and forth... so there's very few (approximately sqrt(-1) as whatever reason you think you have is probably imaginary) good reasons to even bother having an exp() approximation if you can have a base-2 exponential instead, because you can only ever just lose cycles.

For high-precision evaluation (<1ulp), e^x core is better, because it has 3 nice coefficients (multiplying by those is error-free) 1+x+x^2/2, and efficient Cody-Waite reduction can be used to compute a reduced argument. However, accurate range reduction rarely matters, because if is matters, computation is ill-conditioned anyway.

I actually wouldn't mind having (2^x-1)/x, as this kind of expression pops up occasionally and it's a waste to do a comparison to zero (or actually to a range around zero) plus a division if you have a routine which computes 2^x-1 or (yuck!) 2^x instead

OTOH, if you're not limited to a small range, then you do need to actually compute (2^x-1)/x.

A rational "core" can be used for 2^x-1. That way the division will be still there, but more powerful rational approximations can be used.

[juha_p]

Hmm... tried this double sqrt() with 2^x as well and here are the coeffs and other data (couldn't test this now because I have it just in range -1,1 implemented):

Code: Select all

p=remez(f,4,dom);
accurateinfnorm(p^2^2-2^x,[-0.5;0.5],20); 1.31971802375119295902550220489501953125e-8
dirtyinfnorm(2^x-p^2^2,[-0.5;0.5]); 1.3197171501518192197598758444488459621739396608267e-8

1.0000000000367348607289336412036915102213520379882 + x * (0.173286769701991415626795268443594646743710730346953 + x * (1.50141547747944792917778046149247090542777402712443e-2 + x * (8.6765868155922171692784981446128844754131880183393e-4 + x * 3.758727249543877436852672005625649211951091899402e-5)))

Code: Select all

p=remez(f,5,dom);
accurateinfnorm(p^2^2-2^x,[-0.5;0.5],20); 9.52817824639851096435450017452239990234375e-11
dirtyinfnorm(2^x-p^2^2,[-0.5;0.5]); 9.5281739758430379946129042834601272507663795286692e-11

double remez5_powx4(double x){

	double r = 1.0000000000183666595120045519598613442637657650793 + x * (0.1732867951427143495324664921016693522491084565535 + 
			x * (1.50141553625752641632557409228264081801811626212897e-2 + x * (8.672516287797738395395898209753499641154743364598e-4 + 
				x * (3.7584921374457713589931932314789806691527044082933e-5 + x * 1.3025708341515378827425544057794225636314496429767e-6))));
	return (r*r)*(r*r);			

}

[mystran]

For high-precision evaluation (<1ulp), e^x core is better, because it has 3 nice coefficients (multiplying by those is error-free) 1+x+x^2/2, and efficient Cody-Waite reduction can be used to compute a reduced argument. However, accurate range reduction rarely matters, because if is matters, computation is ill-conditioned anyway.

https://github.com/evanphx/ulysses-libc ... math/exp.c has a pretty cool approach that claims less than 1 ulp for doubles.

That said, it's hard to imagine a situation in audio DSP where such an accuracy would really matter, because as you note, if it does matter then you probably have bigger problems to deal with.

[juha_p]

Here are plots for two low order Pade approximations done using double sqrt() method :
https://www.desmos.com/calculator/dvbpscgctf

and triple sqrt() method:
https://www.desmos.com/calculator/wfxtkwgkhb

[juha_p]

Is it ("phenomenon") just that taking n times sqrt() from 2^x (or multiplying x with 1/2^2n) before approximation makes it easier to approximate ... as triple sqrt() makes the response of original function ( 2^x in this case) close to straight line? Are there any drawbacks in doing this way?

[juha_p]

Noticed that response of exp() approximation polynomial (taylor) 1+x(1+x/2) can be improved to 8e-5 level by adding x^pi*0.21 and to 2.6e-5 level by adding 6.6*10^-5*sin(x*8.0). Few degrees higher Taylor poly is needed to get as low error.
Here are couple plots:
https://www.desmos.com/calculator/79ydvvefee
https://tinyurl.com/y2deq8vj (WolframAlpha)

Problem is, how those two fixes can be "joined to the base polynomial"/"included in function to be approximated" without loosing the effects of them ... ? Sollya's composepolynomials(poly1, poly2) didn't give good answer. Should one 1st approximate those fixes and try join them after that?

[2DaT]

Noticed that response of exp() approximation polynomial (taylor) 1+x(1+x/2) can be improved to 8e-5 level by adding x^pi*0.21 and to 2.6e-5 level by adding 6.6*10^-5*sin(x*8.0). Few degrees higher Taylor poly is needed to get as low error.
Here are couple plots:
https://www.desmos.com/calculator/79ydvvefee
https://tinyurl.com/y2deq8vj (WolframAlpha)

Problem is, how those two fixes can be "joined to the base polynomial"/"included in function to be approximated" without loosing the effects of them ... ? Sollya's composepolynomials(poly1, poly2) didn't give good answer. Should one 1st approximate those fixes and try join them after that?

This does not make any sense. x^pi and sin(x) are much more expensive than 1 polynomial term. It's not that functions can't be approximated in arbitrary bases, in fact, Remez algorithm can be applied to any so called "Chebyshev system" of functions. It's just not worth it from cost-benefit analysis.

[juha_p]

This does not make any sense. x^pi and sin(x) are much more expensive than 1 polynomial term. It's not that functions can't be approximated in arbitrary bases, in fact, Remez algorithm can be applied to any so called "Chebyshev system" of functions. It's just not worth it from cost-benefit analysis.

... I understand your point (3rd "order" remez has about same error level) but, I quess, remez and chebyshev routines are not an easy task compared to Taylor series polynomial + x^y + sin(x)?
I didn't try approximate those fixes yet but, IIRC, sin(x) is quite easy case since range is not wide (x: 0-0.5) ... . But, as degree 5 Taylor reaches about same error level, it is needed to get one degree less polynomial from degree 2 poly + one or 2 fixes for it ... is this technically possible?

[2DaT]

This does not make any sense. x^pi and sin(x) are much more expensive than 1 polynomial term. It's not that functions can't be approximated in arbitrary bases, in fact, Remez algorithm can be applied to any so called "Chebyshev system" of functions. It's just not worth it from cost-benefit analysis.

... I understand your point (3rd "order" remez has about same error level) but, I quess, remez and chebyshev routines are not an easy task compared to Taylor series polynomial + x^y + sin(x)?
I didn't try approximate those fixes yet but, IIRC, sin(x) is quite easy case since range is not wide (x: 0-0.5) ... . But, as degree 5 Taylor reaches about same error level, it is needed to get one degree less polynomial from degree 2 poly + one or 2 fixes for it ... is this technically possible?

Use Chebyshev polynomials, those are easy to construct.