Serving the Quantitative Finance Community

  • 1
  • 3
  • 4
  • 5
  • 6
  • 7
  • 11
 
User avatar
Cuchulainn
Posts: 20253
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Optimized math functions

October 11th, 2010, 6:01 pm

QuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnQuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnDo Microsoft /Intel use Pade, or is it something else? This is how Intel does it: The Computation of Transcendental Functions on the IA-64 ArchitectureI had a look at this. Fair enough.Now play 'advocaat van de duivel'.They wrote x = Nlog(2)/2^K + rand then exp(x) = 2^(N/2)^K * exp(r)Now:1. They they use Remez to approxmate exp(r) by a polynomial of degree 5 (==> Solve a matrix system). 2. They get typically r < log(2)/64 = 0.01 which can be done with Pade, so why use Remez?You could try Intel solution versus the bespoke R(1,1) solution. I miss the rationale for Remez (btw Remez uses rational functions, they mention polynomials which can be notoriously bad). They build table lookups as well.I don't think they solve a matrix system... I expect that they find the coeficient *once*, and cast that in logic. After that they just evaluate the polynomial with the constant coeffiients -at least that's how I would do it-.Would be cool to compare various alternatives, especially for fast evaluation / low precision. Maybe the downside of Pade is *not* the precision / number of terms, but the fact that is has a division? The double precision Intel version described in the papers costs 70 cycles, a single division already 30.According to the Remez algo they do find the polynomial coefficients by solving a linear system. At least from wiki.I just get the feeling that there's a lot going on, so a comparative analysis would indeed be very useful. Maybe renorm would solve it in a day.
 
User avatar
Cuchulainn
Posts: 20253
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Optimized math functions

October 11th, 2010, 6:31 pm

RemezSee Step 1.A good one is to compare Pade approx with Boost::Remez on exp(z).
 
User avatar
Cuchulainn
Posts: 20253
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Optimized math functions

October 11th, 2010, 6:44 pm

QuoteOriginally posted by: outrunQuoteOriginally posted by: CuchulainnRemezSee Step 1.A good one is to compare Pade approx with Boost::Remez on exp(z).Yes I was reading that base on you previous remark about wiki I think that what is does is find a Chebyshev approximation of a function, and that makes sense. Chebyshev approximation generally converge very quickly Compare it to fitting a polynomial to exp(z), finding c0,c1,c2.. in exp(z) = c0 +c1 z + c2 z^2 + ...However, once you have fitted the polynomial approximation and calculates the c0,c1,c2,.., you can evaluate exp(x) as c0 + c1*x + c2*x^2Yes, it is a one-off which is not too bad.
 
User avatar
AVt
Posts: 90
Joined: December 29th, 2001, 8:23 pm

Optimized math functions

October 11th, 2010, 6:52 pm

Here is a brute way, without an attempt for polishing (any much range reductionor bit hacking) for single precision:E0(x) = (8193/8192+8389461/16777216*x)*x+1 ~ exp(x) on +-1/32 (thus enough fordiscounting at 3% for 1 year), relative error are below ~ 1E-6.For |x| <= 1/2 one can use exp(t/2)^2 = exp(t) iteratively with the pseudo code below,resulting in relative errors below ~ 2E-5, covering 'better than 1% in most cases'. Note, that division by pure powers of 2 is 'free'.On +-1/64 a polynomial of degree=6 (~ Intel's mentioned method) is enough tohave (relative) errors below 7E-19 there, thus Tang's table method works (oneof the authors) and one can pass to sound solutions, I would say.left:= -1/32;E := proc(x)local e; 1 if abs(x) <= evalf(abs(left)) then 2 return E0(x) end if; 3 if abs(x) <= 2*evalf(abs(left)) then 4 e := E0(1/2*x); 5 return e*e end if; 6 if abs(x) <= 4*evalf(abs(left)) then 7 e := E0(1/4*x); 8 e := e*e; 9 return e*e end if; 10 if abs(x) <= 8*evalf(abs(left)) then 11 e := E0(1/8*x); 12 e := e*e; 13 e := e*e; 14 return e*e end if; 15 if abs(x) <= 16*evalf(abs(left)) then 16 e := E0(1/16*x); 17 e := e*e; 18 e := e*e; 19 e := e*e; 20 return e*e end if; 21 return exp(x)end proc
 
User avatar
AVt
Posts: 90
Joined: December 29th, 2001, 8:23 pm

Optimized math functions

October 11th, 2010, 7:03 pm

 
User avatar
Cuchulainn
Posts: 20253
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Optimized math functions

October 12th, 2010, 6:05 am

QuoteOriginally posted by: outrunQuoteOriginally posted by: AVtwhat's that?Two emoticons.
 
User avatar
Cuchulainn
Posts: 20253
Joined: July 16th, 2004, 7:38 am
Location: 20, 000

Optimized math functions

October 12th, 2010, 8:28 am

QuoteOriginally posted by: outrunjust found a paper that contains the Ramez approx for exp(z) on [0..1/2] asexp(z) = (32793 +31836 z + 21146 z^2)/32768it should have a 10 bit precision, but I get a 3% error... mmmm?!?!This might be a good experiment; determine the accuracy desired TOL, then choose q and p s.t.f(x) - R(p,q) = O(t^p+q+1) == TOL for t < 1BTW Knuth has 1 page on this..overview
Last edited by Cuchulainn on October 11th, 2010, 10:00 pm, edited 1 time in total.