Representing small values in C

[hamster_nz]

The range of numbers here is the problem - (2^22)^5 goes from 0 to 1.33e+33. And as one term is the 4th power, and the other is the fifth power, for most of values of 'x' the x^5 term dominates and the other term is insignificantly tiny once x > 70,000 or so (when X is a range of 0 to 4,000,000 or so).

(premature optimization is the root of much evil.)

The full quote is...

"We should forget about small efficiencies, say about 97% of the time: premature optimization is the root of all evil. Yet we should not pass up our opportunities in that critical 3%" - Donald Knuth

On an small PIC this calculation could well be the critical 3%.

I don't think that the original post was optimizing too soon - they were trying to explore their options early enough that they would be able to weight up options.

[Mechatrommer]

and x*x*x*x is probably faster than pow(x,4) anyway.

Code: [Select]

y=x*x
z=y*y
is probably faster

[Kalvin]

If you perform numerical analysis for the typical numbers x and calculate the values, it might be possible that even a simple, piece-wise linear look-up table will give sufficient accuracy. Taking the logarithm of the number x and computing the lookup table for the logarithms would possibly reduce the lookup table size considerably without introducing too much of an error.

Edit: Better yet, you can normalize the x values using simple shifts for some convenient numerical range, compute the lookup table for the normalized values, use the normalized values for indexing the lookup table, and finally denormalize the the lookup values using simple shifts again. Should get pretty accurate results without too much of computational efforts.

[hamster_nz]

So I was a little bored, and wanted to see if custom floating point code would actually work. I tested it on laptop CPU, and not indicative of performance on a PIC, but it at least proves that the math is somewhat sound....

Performance
Calculating approx 1,000,000,000 times (250 times through the 2^22 range)

Hardware Floating point = 5 seconds (200M calcs per second)
My floating point = 60 seconds (16M calcs per second)

Approximately 12 times slower than H/W floating point. Not as bad as I would have guessed, without even optimizing...

If I used 64-bit multiplications, (rather than 4x32-bit ones) performance is then only 7x slower than H/W floating point.

Error vs floating point
Difference between h/w floating point and my floating point:

Max error at 2097333 is 0.000000004570%

Accuracy could be relaxed to improve performance.

Code: [Select]

/**********************************************
* Testing a range limited custom floating point
* for calculating c1*pow(x,5)+c2*pow(x,4) over
* a range of x=0 to x=2^22-1.
* Assumes 32-bit 'signed' and 'unsigned' types.
********************************************/
#include <stdio.h>
#include <time.h>

double c1 = 1.40715e-09;
double c2 = 9.18189e-08;

/* 1.4071499998200009962090462067863e-9 in custom floating point*/
unsigned m_c1 = 3244666990;
signed e_c1 = -29;

/*  9.18189e-08 in custom floating point */
unsigned m_c2 = 3308124510;
signed e_c2 = -23;

/****************************************************************/
/* Conversion for printing                                      */
/****************************************************************/
static double fixed_to_float(unsigned m, int e) {

  double t;
  t = m;
  t /= 65536;
  t /= 65536;
  while(e > 0) {
    t *= 2; e--;
  }
  while(e < 0) {
    t /= 2; e++;
  }
  return t;
}

/****************************************************************/
/* Adjust so MSB of mantissa is set                             */
/****************************************************************/
static void fixed_normalize(unsigned *m, signed *e) {
  /* Adjust to have a one in the MSB of mantissa */
  if(*m == 0) {
    *e = 0;
     return;
  }
  while((*m & 0x80000000) == 0) {
    *m <<= 1;
    *e = *e -1;
  }
}


/****************************************************************/
/* Multiply two numbers together, using only 32bit math         */
/****************************************************************/
static void fixed_mult(unsigned *m1, signed *e1, unsigned *m2, signed *e2) {
  unsigned h1,l1,h2,l2;
  unsigned h1h2, h1l2, l1h2, l1l2;
  unsigned carry,acc;

  if(*m1 == 0 || *m2 == 0) {
    *m1 = 0;
    *e1 = 0;
  }

  h1 = *m1 >> 16;
  l1 = *m1 & 0xFFFF;
  h2 = *m2 >> 16;
  l2 = *m2 & 0xFFFF;

  h1h2 = h1 * h2;
  h1l2 = h1 * l2;
  l1h2 = l1 * h2;
  l1l2 = l1 * l2;

  carry = 0;
  acc = l1l2 >> 16;

  acc += h1l2;
  if(acc < h1l2) carry += 0x10000;

  acc += l1h2;
  if(acc < l1h2) carry += 0x10000;

  acc >>= 16;
  acc += carry + h1h2;
  *m1 = acc;
  *e1 += *e2;
  fixed_normalize(m1,e1);
}

/****************************************************************/
/* Add two numbers together                                     */
/****************************************************************/
void fixed_add(unsigned *m1, signed *e1, unsigned *m2, signed *e2) {
  if(*m2 == 0) {
    return; /* Nothing to add */
  }

  if(*m1 == 0) {
    *m1 = *m2;
    *e1 = *e2;
    return; /* Nothing in m1 at the moment */
  }

  /* Do we need to move the decimal on m1*2^e1 ? */
  while(*e2 > *e1) {
    *m1 >>= 1;
    *e1 = *e1 + 1;
  }

  /* Do we need to move the decimal on m2*2^e2 ? */
  while(*e1 > *e2) {
    *m2 >>= 1;
    *e2 = *e2 + 1;
  }

  /* As e2 = e1 we can now just add */
  *m1 += *m2;

  /* Test for overflow */
  if(*m1 < *m2) {
    *m1 = (*m1>>1) | 0x80000000;
    *e1 = *e1 + 1;
  }
}

/****************************************************************/
/* The main calculation in custom fixed point                   */
/****************************************************************/
double calcFixed(unsigned i, unsigned *m, signed *e) {
  unsigned m_x = i;
  signed e_x = 32;
  unsigned m_t1;
  signed e_t1;
  unsigned m_t2;
  signed e_t2;
  fixed_normalize(&m_x, &e_x);
  m_t2 = m_t1 = m_x;
  e_t2 = e_t1 = e_x;

  fixed_mult(&m_t1, &e_t1, &m_x, &e_x);  /* t1 = x^2 */
  fixed_mult(&m_t1, &e_t1, &m_x, &e_x);  /* t1 = x^3 */
  fixed_mult(&m_t1, &e_t1, &m_x, &e_x);  /* t1 = x^4 */
  fixed_mult(&m_t2, &e_t2, &m_t1, &e_t1);  /* t2 = x^5 */
  fixed_mult(&m_t1, &e_t1, &m_c2, &e_c2); /* t1 = c2 * x ^ 4 */
  fixed_mult(&m_t2, &e_t2, &m_c1, &e_c1); /* t2 = c1 * x ^ 5 */
  fixed_add(&m_t1, &e_t1, &m_t2, &e_t2);
  *m = m_t1;
  *e = e_t1;
  return 0;
}

/****************************************************************/
/* The same calculation in floating point                       */
/****************************************************************/
double calcFloat(unsigned i) {
  double xd = i;
  return c1 * xd*xd*xd*xd*xd + c2* xd*xd*xd*xd;
}

/****************************************************************/
/* Test harness to verify the code                              */
/****************************************************************/
int main(int argc, char *argv[]) {
  int start_time;
  int end_time;
  int j;
  unsigned i;
  double t = 0;
  double max_error=0.0, min_error=0.0;
  int max_i = 0, min_i = 0;

  /* Time the H/W floating point routine */
  start_time = time(NULL);
  for(j = 0; j < 250; j++) {
    for(i = 0; i < 0x400000; i++) {
      t+= calcFloat(i);
    }
  }
  end_time = time(NULL);
  printf("Float took %i seconds\n", end_time-start_time);
  printf("t is %f (to prevent optimization)\n",t);

  /* Time the custom point routine */
  start_time = time(NULL);
  for(j = 0; j < 250; j++) {
    for(i = 0; i < 0x400000; i++) {
      unsigned m;
      signed e;
      calcFixed(i, &m, &e);
    }
  }
  end_time = time(NULL);
  printf("Fixed Took %i seconds\n", end_time-start_time);

  /* Check the accuracy */
  for(i = 1;i < 4*1024*1024; i++) {
    unsigned m;
    signed e;
    double error;
    calcFixed(i, &m, &e);
    error = 1.0 - fixed_to_float(m,e)/calcFloat(i);

    if(max_error < error) {
      max_error = error;
      max_i = i;
    }

    if(min_error > error) {
      min_error = error;
      min_i = i;
    }
  }
  printf("Max error at %i is %3.6f%%\n",max_i, max_error);
  printf("Min error at %i is %3.6f%%\n",min_i, min_error);
  return 0;
}


[dannyf]

"If you perform numerical analysis"

You would quickly come to the realization that for moderately large x, the x5 term dominates, ,,,,

[Kalvin]

"If you perform numerical analysis"

You would quickly come to the realization that for moderately large x, the x5 term dominates, ,,,,

Only if abs(x) > 1.0. If abs(x) < 1.0, the abs(x^5) term will be less than abs(x^4).

[dannyf]

The number to beat is 3000 instructions per calculation.

[dannyf]

And about 600 bytes in code size.

[Howardlong]

I may have missed it, and I apologise if I have, but what performance is the OP looking for?

While we can certainly do many tricks once we understand the inputs and outputs, if, say, 1 or 2 ms is good enough with the compiler's software floating point then it becomes little more than an academic exercise. While I love solving such exercises in performance analysis, without a requirement goal it's a bit meaningless.

I wrote a complete autonomous antenna satellite tracking system about 12 years ago based on a PIC18, it did all the calculations for real time presiction and point antennas for multiple satellites for a given set of orbital parameters (Kelerian elements). That was pretty scary, but it worked pretty well and was remarkably fast, tens of updates/sec, and would predict upcoming AOS and LOS (when each satellite is above the horizon) on the fly too, an even more computationally intensive thing to achieve.

[dannyf]

The number to beat is 3000 instructions per calculation.

If you allow approximation, I can reduce that number to 2000 instructions for 99.98% of the x's values. Flash size no material change.