CODEWITHSUNDEEP
cgccfloating-point
camelccc
camelccc

Reputation: 2992

accuracy of pow vs repeated multiplications in c

I have some floating point code that under code profiling gets bogged down in repeated calls to pow. Since the number of values that is being powered is very limited, and the exponents can be factored to be all integers, I want to replace these calls with pre-computed arrays.

Code performance wise, pre-computing the arrays is unlikely to be very significant, however what are the accuracy implications of :

array[0]= 1 ;
for (int n=1 ; n<100 ; n++)
    array[n]=array[n-1]*some_double ;

vs

for (int n=0 ; n<100 ; n++)
    array[n]=pow(some_double , n) ;

assume n does not exceed 100, probably less. Answers for long double also appreciated if different.

Upvotes: 1

Views: 561

Answers (4)

chux
chux

Reputation: 153338

code profiling gets bogged down in repeated calls to pow.

To improve overall speed performance, rather than look how to do xy faster, post the larger code. Best optimizations involve a wider view of code.

* vs. pow() performance.

Repeated multications * incur up to 1/2 ULP in the product each time. With 10100 that is 99 multiplications and an error of √99 * 0.5 ULP can be expected.

The compiler, depending on FLT_EVAL_METHOD may use wider math for * and result in negligible error.

pow() is a tricky library function to implement well. The error in z = xy with simplistic implementations is proportional to ln(|z|) (IIRC) - this is quite bad for large z. A good pow() will use extended math and other crafted code to minimize the error, potentially down to 1.0 ULP. This comes at a modest time performance cost.

For a quality answer, use pow() unless raising to a simplistic small integer power like 2,3,4.

Some alternative code for double raise to an positive integer power. More time efficient than a loop of multiplications once n is maybe 5 or more. Likely more accurate.

double pow_n(double base, unsigned exp) {
  double y = 1.0;
  while (1) {
    if (exp % 2) y *= base;
    exp /= 2;
    if (exp == 0) break;  
    base *= base;
  }
  return y;
}

Upvotes: 1

gnasher729
gnasher729

Reputation: 52530

Your consecutive multiplications lead to 100 rounding errors being added up. You can do that in more clever ways: Store x^0 and x^1, then use x^(2n) = (x^n)^2, and x^(2n+1) = x^(2n) * x. Much less build-up on rounding errors.

Now compare to pow(x, n): A good library implementation will realise that n is an integer and calculate the result using multiplications, and it will use the same multiplications that my code above did, or it will do something more clever. A more clever implementation would use long double for the intermediate results, or would use quad precision arithmetic to get the results with the smallest possible rounding errors.

Since you only calculate 100 powers, if you are really interested in keeping the rounding error as low as possible, I'd be tempted to use pow().

Upvotes: 3

phoenixstudio
phoenixstudio

Reputation: 2598

As far as I know there is no circuit to calculate Exponentiation, which means that the pow() function must be an algorithm (can't be done by hardware at single clock cycle). There are many papers about this and the main concern is precision and not performance (each time you do a multiplication you lose precision) some algorithms compensate after each multiplication and this will require more computation. My assumptions is not about your case, in worst case scenario the algorithm for pow() will be the same as doing it manually, but taking into account all CPUs, operating systems and compilers.

//this should be faster
array[0]= 1 ;
for (int n=1 ; n<100 ; n++)
array[n]=array[n-1]*some_double ;

//this should be more accurate 
for (int n=0 ; n<100 ; n++)
    array[n]=pow(some_double , n) ;

Upvotes: 0

Jabberwocky
Jabberwocky

Reputation: 50778

This should give you a vague idea:

#include <stdio.h>
#include <math.h>

double mypow(double d, int n)
{
  double result = d;
  for (int i = 1; i < n; i++)
    result *= d;

  return result;
}

int main(int argc, char* argv[])
{
  for (double d = 0; d < 10; d += 0.1)
  {
    for (int n = 1; n < 100; n++)
    {
      double r1 = pow(d, n);
      double r2 = mypow(d, n);    
      printf("d = %f, n = %d %f %f %f\n", d,n,r1, r2, r1-r2);
    }
  }
}

Upvotes: 3

Related Questions