Approximate concentric square to disc

January 11th, 2017

This post is pretty much automatically generated + copy-n-paste.

In my last post I more or less assumed everyone is at least passingly familar with polymomial approximations of sine and cosine on a narrow range. And approximating concentric square to disc is only interesting if you don’t have fast hardware sine and cosine.

Example of approximate concentric

A straight forward implementation of concentric with approximations might look like this:

void foo(vec2_t* d, vec2_t* s) 
{
  float x = s->x;
  float y = s->y;

  // capt. obvious says (n=num of terms):
  // sin_n -> (n+1) products, n constant loads, (n-1) adds
  //   and
  // cos_n -> (n-1) products, n constant loads, (n-1) adds
  //   or 
  // sqrt  -> seq after sin complete

  if (fabsf(x) > fabsf(y)) {
    float t = y*recip(x);
    float s = sin_n(t);
    d->x = x*cos_n(t); // or x*sqrt(1.f-s*s) & longer dep-chain
    d->y = x*s;
  }
  else {
    float r = x*recip(y+EPS);
    float s = sin_n(t);
    d->x = y*s;
    d->y = y*cos_n(t);  // ditto
  }
}

The area distortion of the approximate equal area method from the previous post:

Which almost all outside the plot range here. Now concentric with sine and cosine approximated using 2 terms, minimizing relative error and clamped to correct at end-points (sin2c and cos2c below):

and now using square root instead of cosine:

Which is a nice improvement (the approx sin squared + cos squared != 1 is corrected). Change the sine approximation to 3-terms (still root for cosine):

Using the scale of the previous plots this would be all white.

Spewed out approximations

Here are some low rank approximations of sine and cosine. The non ‘c’ variant are minimax polynomials and the ‘c’ variants are sloppy end-point contrained. These were generated with this (script). All of the approximations are for minimizing relative error, but the error plots below are abs error (my sollya base scripts need some work). The awkward form of the functions is to be bit exact singles.

\[\begin{align} sin_{2c} \left(\frac{\pi x}{4}\right) & =13171413 \times 2^{-24} x-5232519 \times 2^{-26} x^{3}\\ sin_{2} \left(\frac{\pi x}{4}\right) & =13171413 \times 2^{-24} x-10503801 \times 2^{-27} x^{3}\\ sin_{3c} \left(\frac{\pi x}{4}\right) & =13176775 \times 2^{-24} x-2708675 \times 2^{-25} x^{3}+2614125 \times 2^{-30} x^{5}\\ sin_{3} \left(\frac{\pi x}{4}\right) & =13176775 \times 2^{-24} x-2708675 \times 2^{-25} x^{3}+5230543 \times 2^{-31} x^{5}\\ sin_{4c} \left(\frac{\pi x}{4}\right) & =13176795 \times 2^{-24} x-5418747 \times 2^{-26} x^{3}+668505 \times 2^{-28} x^{5}-4969345 \times 2^{-37} x^{7}\\ sin_{4} \left(\frac{\pi x}{4}\right) & =13176795 \times 2^{-24} x-5418747 \times 2^{-26} x^{3}+668505 \times 2^{-28} x^{5}-9940169 \times 2^{-38} x^{7}\\ \end{align}\]

\[\begin{align} cos_{2c} \left(\frac{\pi x}{4}\right) & =16739113 \times 2^{-24}-2437915 \times 2^{-23} x^{2}\\ cos_{2} \left(\frac{\pi x}{4}\right) & =16739113 \times 2^{-24}-9805545 \times 2^{-25} x^{2}\\ cos_{3c} \left(\frac{\pi x}{4}\right) & =8388509 \times 2^{-23}-5171259 \times 2^{-24} x^{2}+16481549 \times 2^{-30} x^{4}\\ cos_{3} \left(\frac{\pi x}{4}\right) & =8388509 \times 2^{-23}-5171259 \times 2^{-24} x^{2}+4122623 \times 2^{-28} x^{4}\\ cos_{4c} \left(\frac{\pi x}{4}\right) & =16777215 \times 2^{-24}-10348991 \times 2^{-25} x^{2}+4254515 \times 2^{-28} x^{4}-341983 \times 2^{-30} x^{6}\\ cos_{4} \left(\frac{\pi x}{4}\right) & =16777215 \times 2^{-24}-10348991 \times 2^{-25} x^{2}+4254515 \times 2^{-28} x^{4}-10944157 \times 2^{-35} x^{6}\\ \end{align}\]

The inline function defs generated by the script:

static inline float sin2(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -7.8259415924549102783203125e-2f);
  r  = x *(r+0.785077393054962158203125f);
  return r;
}

static inline float sin2c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -7.797060906887054443359375e-2f);
  r  = x *(r+0.785077393054962158203125f);
  return r;
}

static inline float sin3(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( +2.4356613866984844207763671875e-3f);
  r  = x2*(r-8.07248055934906005859375e-2f);
  r  = x *(r+0.785396993160247802734375f);
  return r;
}

static inline float sin3c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( +2.434593625366687774658203125e-3f);
  r  = x2*(r-8.07248055934906005859375e-2f);
  r  = x *(r+0.785396993160247802734375f);
  return r;
}

static inline float sin4(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -3.616212416091002523899078369140625e-5f);
  r  = x2*(r+2.4903751909732818603515625e-3f);
  r  = x2*(r-8.074562251567840576171875e-2f);
  r  = x *(r+0.785398185253143310546875f);
  return r;
}

static inline float sin4c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -3.61567435902543365955352783203125e-5f);
  r  = x2*(r+2.4903751909732818603515625e-3f);
  r  = x2*(r-8.074562251567840576171875e-2f);
  r  = x *(r+0.785398185253143310546875f);
  return r;
}

static inline float cos2(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -0.2922280132770538330078125f);
  r  =    (r+0.997728884220123291015625f);
  return r;
}

static inline float cos2c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -0.29062211513519287109375f);
  r  =    (r+0.997728884220123291015625f);
  return r;
}

static inline float cos3(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( +1.53579674661159515380859375e-2f);
  r  = x2*(r-0.308231055736541748046875f);
  r  =    (r+0.99998819828033447265625f);
  return r;
}

static inline float cos3c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( +1.5349638648331165313720703125e-2f);
  r  = x2*(r-0.308231055736541748046875f);
  r  =    (r+0.99998819828033447265625f);
  return r;
}

static inline float cos4(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -3.1851688981987535953521728515625e-4f);
  r  = x2*(r+1.58493034541606903076171875e-2f);
  r  = x2*(r-0.3084239661693572998046875f);
  r  =    (r+0.999999940395355224609375f);
  return r;
}

static inline float cos4c(float x) {
  float r;
  float x2 = x*x;
  r  = x2*( -3.18496488034725189208984375e-4f);
  r  = x2*(r+1.58493034541606903076171875e-2f);
  r  = x2*(r-0.3084239661693572998046875f);
  r  =    (r+0.999999940395355224609375f);
  return r;
}

Comments

math (33)