Small challenge

[robve]

(04-23-2023 11:10 AM)robve Wrote:  Indeed! It is common to find exponentiation by squaring for integer powers, which tends to be more accurate than repeated multiplication and log*exp closed forms.

It's pretty straightforward to do.

We can use the Chains of Recurrences algebra to optimize sequence evaluations such as x^10. We only need 10 additions to produce the next value in the sequence.

To expand \( x^{10} \) to CR form, let \( x \) be the integer sequence CR \( x =\{1,+,1\} \) expand to the CR \( \{1,+,1\}\times\{1,+,1\}\times\{1,+,1\}\times\{1,+,1\}\times\{1,+,1\}\times\{1,+​,1\}\times\{1,+,1\}\times\{1,+,1\}\times\{1,+,1\}\times\{1,+,1\} \). Simplify using the CR rule \( \{\phi_0,+,f_1\} \times \{\psi_0,+,g_1\} \Rightarrow
\{\phi_0\times\psi_0,+,\{\phi_0,+,f_1\}\times g_1+\{\psi_0,+,g_1\}\times f_1+f_1\times g_1\} \) giving the polynomial CR \( \{1,+,1023,+, 57002,+, 874500,+, 5921520,+, 21538440,+, 46070640,+, 59875200,+, 46569600,+, 19958400,+, 3628800 \} \). Note that each value is produced with just 10 additions.

Written in C:

Code:

#include <stdio.h>
#include <stdlib.h>
int main()
{
  uint64_t r[11] = {1, 1023, 57002, 874500, 5921520, 21538440, 46070640, 59875200, 46569600, 19958400, 3628800 };
  int i, j;
  for (i = 0; i < 24; ++i)
  {
    printf("%llu\n", r[0]);
    for (j = 0; j < 10; ++j)
      r[j] += r[j+1];
  }
}

Output:
1
1024
59049
1048576
9765625
60466176
282475249
1073741824
3486784401
10000000000
25937424601
61917364224
137858491849
289254654976
576650390625
1099511627776
2015993900449
3570467226624
6131066257801
10240000000000
16679880978201
26559922791424
41426511213649
63403380965376

With floating point we need to be a bit careful with error accumulation in long CR iterations. The polynomial case, such as \( x^{10} \), is stable. The product case such as \( 10^x \) may not be stable. The above code uses integers, which is always stable (but may overflow eventually).

- Rob