title: Fastest pow() replacement via modified exp. by squaring when lower powers are already calculated
body: <p><em><strong>EDIT:</em></strong> <br/></p>
  
  <p><em>Goal :</em> <br/>
  Generate a ubiquitous method for deriving a custom power function that outperforms the built-in <code>pow(double, uint)</code> by reusing precalculated/cached powers from power calculations on common variables.</p>
  
  <p><em>What's already been done:</em><br/>
  I've already derived such a function that's roughly 40% faster than the built-in, however this is a brute-force hand-derived function -- I want a method for autogenerating such a power function block for an arbitrary <code>uint</code> power.</p>
  
  <hr>
  
  <p><strong>KNOWNS</strong></p>
  
  <p>To derive an optimal custom <code>pow(double, uint)</code> you need some knowns.  For this question the knowns (to clarify) are:</p>
  
  <ol>
  <li>The power will be an integer.</li>
  <li>The maximum the power can be is known (<code>N_MAX</code>).</li>
  <li>The precalculated powers that can be (re)used are known
  at compile time (e.g. in my example <code>r2</code>, <code>r4</code>, and <code>r6</code>).</li>
  <li>The square <code>r2</code> can be assumed to always have been calculated regardless
  of the other precalculated powers.</li>
  </ol>
  
  <hr>
  
  <p><strong>SOLUTION REQUIREMENTS</strong></p>
  
  <p>An optimal solution requiring a separate program to write a <code>case</code> lookup table or preprocessor logic to generate such a table is acceptable, however, non-optimal solutions using hand-generated (i.e. brute force derived) lookup tables using the powers on hand will not be accepted (as I have that already, and show that in my example... the idea is to get away from this).</p>
  
  <hr>
  
  <p><strong>POSSIBLE SOLUTION ROUTE</strong></p>
  
  <p>As a suggestion, you know <code>N_MAX</code> and a set of powers that are precalculated <code>B</code> (<code>B={2,4,6}</code> for my example).  You can produce either in a separate program or in the preprocessor a table of all squares of <code>Sq(Bi, x</code>) &lt;= N_MAX<code>.  You can use this to form a basis set</code>A<code>, which you then search somehow to determine the least number of terms that can be summed to produce an arbitrary exponent of</code>n>>1<code>, where</code>n&lt;=N_MAX` (the shift is due to that we take care of the odd case by checking the LSB and multiplying by the sqrt(r2)).</p>
  
  <hr>
  
  <p><strong>THEORETICAL BACKGROUND</strong></p>
  
  <p>I believe formally the below method is a modified version of exponentations by squaring:</p>
  
  <p><a href="http://en.wikipedia.org/wiki/Exponentiation_by_squaring" rel="nofollow">http://en.wikipedia.org/wiki/Exponentiation_by_squaring</a></p>
  
  <p>....which takes advantage of the fact that certain lower order powers are already by necessity precalculated, hence it shifts the optimal set of multiplications from a vanilla exponentation by squaring (which I assume <code>pow(double, int)</code> uses).</p>
  
  <p>However there are significant savings by using the stored small power intermediates instead of simple exp. by squares on the <code>r2</code>.</p>
  
  <hr>
  
  <p><strong>THEORETICAL PERFORMANCE</strong></p>
  
  <p>For example, for one set of objects <code>n=14</code>.... in this scenario exp. by powers gives</p>
  
  <pre><code>double r4 = Sq(r2), r14=Sq(r4)*r4*r2; //4 op.
  </code></pre>
  
  <p>... which takes <strong>4 FP multiplications</strong>..... but using the <code>r2</code> and <code>r6</code> we have</p>
  
  <pre><code>double r14=Sq(r6)*r2; //2 op.
  </code></pre>
  
  <p>.... <strong>2 FP multiplications</strong>.... in other words, by going from "dumb" exponentation by squares to my modified exp. by squares using the common exponent precaching, I've cut my cost of calculations for 50% in terms of multiplications ... at least until memory costs are considered.</p>
  
  <p><strong>REAL PERFORMANCE</strong></p>
  
  <p>With my current method (compiled with <code>gcc -O3</code>) I get <strong>35.1 sec.</strong> to run 1 million cycles of my program, versus (w/ no other modifications) <strong>56.6 s</strong> using the built int <code>pow(double, int)</code>.... so almost the theoretical speedup.</p>
  
  <p>At this point you may be scratching your head at how a 50% cut in multiplications on a single instruction line can deliver a ~40% speedup.  But basically this line of code is called 1,000+ times per cycle and is by far the most evaluated/most expensive line of code in the entire program.  Hence the program appears highly sensitive to a small optimization/improvement in this chunk.</p>
  
  <hr>
  
  <p><strong>ORIGINAL POST and EXAMPLE CODE</strong></p>
  
  <p>I need to replace the <code>pow(double, int)</code> function as I already have calculated a 6th power term and have 2nd, 4th power intermediates saved, all of which can be used to reduce multiplications in the second <code>pow</code> call, which uses the same <code>double</code> base.</p>
  
  <p>More specifically, in my c++ code I have a performance critical calculation snippet of code where I raise the reciprocal of the distance between 3D points to the 6th power and nth power.  e.g.:</p>
  
  <pre><code>double distSq = CalcDist(p1,p2), r2 = a/distSq, r6 = r2 * r2 * r2;
  results += m*(pow(sqrt(r2), n) - r6);
  </code></pre>
  
  <p>Where <code>m</code> and <code>a</code> are constants related to the fitted equation and <code>n</code> is the arbitrary power.</p>
  
  <p>A slightly more efficient form is:</p>
  
  <pre><code>double distSq = CalcDist(p1,p2), r2 = a/distSq, r6 = r2 * r2 * r2;
  results += m*(pow(r2, n)*(n&amp;0x1?sqrt(r2):1.0) - r6);
  </code></pre>
  
  <p>However, this is also not optimal.  What I've found to be significantly faster is to have a custom <code>pow</code> function that uses the multiples r2, r4, and r6, which I have to calculate already anyways for the second term.</p>
  
  <p>e.g.:</p>
  
  <pre><code>double distSq = CalcDist(p1,p2), r2 = a/distSq, r4 = r2 * r2, r6 = r4 * r2;
  results += m*(POW(r2, r4, r6 n) - r6);
  </code></pre>
  
  <p>Inside the function:</p>
  
  <pre><code>double POW(double r2, double r4, double r6, uint n)
  {
     double results = (n&amp;0x1 : sqrt(r2) : 1.0);
     n &gt;&gt;= 1;
     switch (n)
     {
       case 1:
       ....
       case 12:
          Sq(Sq(r6));
  
     }
     return result;
  }
  </code></pre>
  
  <p>The good thing is that my function appears fast in preliminary testing.  The bad news is that it's not very ubiquitous and is very long as I need <code>case</code> statements for <code>int</code> powers from <code>8</code> to <code>50</code> or so (potentially even higher in the future).  Further each case I had to examine and try different combinations to find by brute force derivation which combination of <code>r2</code>, <code>r4</code>, and <code>r6</code> yielded the least multiplications</p>
  
  <p>Does anyone have a more ubiquitous solution for a <code>pow(double, int)</code> replacement that uses precalculated powers of the base to cut the number of necessary multiplications, and/or have a a ubiquitous theory of how you can determine the ideal combination to produce the least multiplications for an arbitrary <code>n</code> and some set of precalculated multiples??</p>