Addition-chain exponentiation

This page uses content from Wikipedia. The original article was at Addition-chain exponentiation. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

In cases of special objects (such as with matrices) the operation of multiplication can be excessively expensive. In these cases the operation of multiplication should be avoided or reduced to a minimum.

In mathematics and computer science, optimal addition-chain exponentiation is a method of exponentiation by positive integer powers that requires a minimal number of multiplications. It works by creating a shortest addition chain that generates the desired exponent. Each exponentiation in the chain can be evaluated by multiplying two of the earlier exponentiation results. More generally, addition-chain exponentiation may also refer to exponentiation by non-minimal addition chains constructed by a variety of algorithms (since a shortest addition chain is very difficult to find).

The shortest addition-chain algorithm requires no more multiplications than binary exponentiation and usually less. The first example of where it does better is for $a^{15}$ , where the binary method needs six multiplies but a shortest addition chain requires only five:

a^{15}=a\times (a\times [a\times a^{2}]^{2})^{2}\!

(binary, 6 multiplications)

a^{15}=a^{3}\times ([a^{3}]^{2})^{2}\!

(shortest addition chain, 5 multiplications)

On the other hand, the addition-chain method is much more complicated, since the determination of a shortest addition chain seems quite difficult: no efficient optimal methods are currently known for arbitrary exponents, and the related problem of finding a shortest addition chain for a given set of exponents has been proven NP-complete.

Table demonstrating how to do *Exponentiation* using *Addition Chains*
Number of Multiplications	Actual Exponentiation	Specific implementation of Addition Chains to do Exponentiation
0	a¹	a
1	a²	a × a
2	a³	a × a × a
2	a⁴	(a × a→b) × b
3	a⁵	(a × a→b) × b × a
3	a⁶	(a × a→b) × b × b
4	a⁷	(a × a→b) × b × b × a
3	a⁸	((a × a→b) × b→d) × d
4	a⁹	(a × a × a→c) × c × c
4	a¹⁰	((a × a→b) × b→d) × d × b
5	a¹¹	((a × a→b) × b→d) × d × b × a
4	a¹²	((a × a→b) × b→d) × d × d
5	a¹³	((a × a→b) × b→d) × d × d × a
5	a¹⁴	((a × a→b) × b→d) × d × d × b
5	a¹⁵	((a × a→b) × b × a→e) × e × e
4	a¹⁶	(((a × a→b) × b→d) × d→h) × h

The number of multiplications required follows this sequence: 0, 1, 2, 2, 3, 3, 4, 3, 4, 4, 5, 4, 5, 5, 5, 4, 5, 5, 6, 5, 6, 6, 6, 5, 6, 6, 6, 6, 7, 6, 7, 5, 6, 6, 7, 6, 7, 7, 7, 6, 7, 7, 7, 7, 7, 7, 8, 6, 7, 7, 7, 7, 8, 7, 8, 7, 8, 8, 8, 7, 8, 8, 8, 6, 7, 7, 8, 7, 8, 8, 9, 7, 8, 8, 8, 8, 8, 8, 9, 7, 8, 8, 8, 8, 8, 8, 9, 8, 9, 8, 9, 8, 9, 9, 9, 7, 8, 8, 8, 8...

This sequence can be found at: http://oeis.org/A003313

Task requirements:

Using the following values: $A={\begin{bmatrix}{\sqrt {\frac {1}{2}}}&0&{\sqrt {\frac {1}{2}}}&0&0&0&\\0&{\sqrt {\frac {1}{2}}}&0&{\sqrt {\frac {1}{2}}}&0&0&\\0&{\sqrt {\frac {1}{2}}}&0&-{\sqrt {\frac {1}{2}}}&0&0&\\-{\sqrt {\frac {1}{2}}}&0&{\sqrt {\frac {1}{2}}}&0&0&0&\\0&0&0&0&0&1&\\0&0&0&0&1&0&\\\end{bmatrix}}$ $m=31415$ and $n=27182$

Repeat task Matrix-exponentiation operator, except use addition-chain exponentiation to better calculate:

A^{m}

,

A^{n}

and

A^{n\times m}

.

As an easier alternative to doing the matrix manipulation above, generate the addition-chains for 31415 and 27182 and use addition-chain exponentiation to calculate these two equations:

1.00002206445416³¹⁴¹⁵
1.00002550055251²⁷¹⁸²

Also: Display a count of how many multiplications were done in each case.

Note: There are two ways to approach this task:

Brute force - try every permutation possible and pick one with the least number of multiplications. If the brute force is a simpler algorithm, then present it as a subtask under the subtitle "Brute force", eg ===Brute Force===.
Some clever algorithm - the wikipedia page has some hints, subtitle the code with the name of algorithm.

Note: Binary exponentiation does not usually produce the best solution. Provide only optimal solutions.

Kudos (κῦδος) for providing a routine that generate sequence A003313 in the output.

C

Using complex instead of matrix. Requires Achain.c. It takes a long while to compute the shortest addition chains, such that if you don't have the chain lengths precomputed and stored somewhere, you are probably better off with a binary chain (normally not shortest but very simple to calculate) whatever you intend to use the chains for. <lang c>#include <stdio.h>

include "achain.c" /* not common practice */

/* don't have a C99 compiler atm */ typedef struct {double u, v;} cplx;

inline cplx c_mul(cplx a, cplx b) { cplx c; c.u = a.u * b.u - a.v * b.v; c.v = a.u * b.v + a.v * b.u; return c; }

cplx chain_expo(cplx x, int n) { int i, j, k, l, e[32]; cplx v[32];

l = seq(n, 0, e);

puts("Exponents:"); for (i = 0; i <= l; i++) printf("%d%c", e[i], i == l ? '\n' : ' ');

v[0] = x; v[1] = c_mul(x, x); for (i = 2; i <= l; i++) { for (j = i - 1; j; j--) { for (k = j; k >= 0; k--) { if (e[k] + e[j] < e[i]) break; if (e[k] + e[j] > e[i]) continue; v[i] = c_mul(v[j], v[k]); j = 1; break; } } } printf("(%f + i%f)^%d = %f + i%f\n", x.u, x.v, n, v[l].u, v[l].v);

return x; }

int bin_len(int n) { int r, o; for (r = o = -1; n; n >>= 1, r++) if (n & 1) o++; return r + o; }

int main() { cplx r1 = {1.0000254989, 0.0000577896}, r2 = {1.0000220632, 0.0000500026}; int n1 = 27182, n2 = 31415, i;

init(); puts("Precompute chain lengths"); seq_len(n2);

chain_expo(r1, n1); chain_expo(r2, n2); puts("\nchain lengths: shortest binary"); printf("%14d %7d %7d\n", n1, seq_len(n1), bin_len(n1)); printf("%14d %7d %7d\n", n2, seq_len(n2), bin_len(n2)); for (i = 1; i < 100; i++) printf("%14d %7d %7d\n", i, seq_len(i), bin_len(i)); return 0; }</lang> output

...
Exponents:
1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568 27136 27182
(1.000025 + i0.000058)^27182 = -0.000001 + i2.000001
Exponents:
1 2 4 8 16 17 33 49 98 196 392 784 1568 3136 6272 6289 12561 25122 31411 31415
(1.000022 + i0.000050)^31415 = -0.000001 + i2.000000

chain lengths: shortest binary
         27182      18      21
         31415      19      24
             1       0       0
             2       1       1
             3       2       2
             4       2       2
      ...
            89       9       9
            90       8       9
            91       9      10
            92       8       9
            93       9      10
      ...

Go

A non-optimal solution. <lang go>/* Continued fraction addition chains, as described in "Efficient computation of addition chains" by F. Bergeron, J. Berstel, and S. Brlek, published in Journal de théorie des nombres de Bordeaux, 6 no. 1 (1994), p. 21-38, accessed at http://www.numdam.org/item?id=JTNB_1994__6_1_21_0.

/

package main

import (

   "fmt"
   "math"

)

// Representation of addition chains. // Notes: // 1. While an []int might represent addition chains in general, the // techniques here work only with "star" chains, as described in the paper. // Knowledge that the chains are star chains allows certain optimizations. // 2. The paper descibes a linked list representation which encodes both // addends of numbers in the chain. This allows additional optimizations, but // for the purposes of the RC task, this simpler representation is adequate. type starChain []int

// ⊗= operator. modifies receiver. func (s1 *starChain) cMul(s2 starChain) {

   p := *s1
   i := len(p)
   n := p[i-1]
   p = append(p, s2[1:]...)
   for ; i < len(p); i++ {
       p[i] *= n
   }
   *s1 = p

}

// ⊕= operator. modifies receiver. func (p *starChain) cAdd(j int) {

   c := *p
   *p = append(c, c[len(c)-1]+j)

}

// The γ function described in the paper returns a set of numbers in general, // but certain γ functions return only singletons. The dichotomic strategy // is one of these and gives good results so it is the one used for the // RC task. Defining the package variable γ to be a singleton allows some // simplifications in the code. var γ singleton

type singleton func(int) int

func dichotomic(n int) int {

   return n / (1 << uint(λ(n)/2))

}

// integer log base 2 func λ(n int) (a int) {

   for n != 1 {
       a++
       n >>= 1
   }
   return

}

// minChain as described in the paper. func minChain(n int) starChain {

   switch a := λ(n); {
   case n == 1<<uint(a):
       r := make(starChain, a+1)
       for i := range r {
           r[i] = 1 << uint(i)
       }
       return r
   case n == 3:
       return starChain{1, 2, 3}
   }
   return chain(n, γ(n))

}

// chain as described in the paper. func chain(n1, n2 int) starChain {

   q, r := n1/n2, n1%n2
   if r == 0 {
       c := minChain(n2)
       c.cMul(minChain(q))
       return c
   }
   c := chain(n2, r)
   c.cMul(minChain(q))
   c.cAdd(r)
   return c

}

func main() {

   m := 31415
   n := 27182
   show(m)
   show(n)
   show(m * n)
   showEasier(m, 1.00002206445416)
   showEasier(n, 1.00002550055251)

}

func show(e int) {

   fmt.Println("exponent:", e)
   s := math.Sqrt(.5)
   a := matrixFromRows([][]float64{
       {s, 0, s, 0, 0, 0},
       {0, s, 0, s, 0, 0},
       {0, s, 0, -s, 0, 0},
       {-s, 0, s, 0, 0, 0},
       {0, 0, 0, 0, 0, 1},
       {0, 0, 0, 0, 1, 0},
   })
   γ = dichotomic
   sc := minChain(e)
   fmt.Println("addition chain:", sc)
   a.scExp(sc).print("a^e")
   fmt.Println("count of multiplies:", mCount)
   fmt.Println()

}

var mCount int

func showEasier(e int, a float64) {

   fmt.Println("exponent:", e)
   γ = dichotomic
   sc := minChain(e)
   fmt.Printf("%.14f^%d: %.14f\n", a, sc[len(sc)-1], scExp64(a, sc))
   fmt.Println("count of multiplies:", mCount)
   fmt.Println()

}

func scExp64(a float64, sc starChain) float64 {

   mCount = 0
   p := make([]float64, len(sc))
   p[0] = a 
   for i := 1; i < len(p); i++ {
       d := sc[i] - sc[i-1]
       j := i - 1
       for sc[j] != d {
           j--
       }
       p[i] = p[i-1] * p[j]
       mCount++
   }
   return p[len(p)-1]

}

func (m *matrix) scExp(sc starChain) *matrix {

   mCount = 0
   p := make([]*matrix, len(sc))
   p[0] = m.copy()
   for i := 1; i < len(p); i++ {
       d := sc[i] - sc[i-1]
       j := i - 1
       for sc[j] != d {
           j--
       }
       p[i] = p[i-1].multiply(p[j])
       mCount++
   }
   return p[len(p)-1]

}

func (m *matrix) copy() *matrix {

   return &matrix{append([]float64{}, m.ele...), m.stride}

}

// code below copied from matrix multiplication task type matrix struct {

   ele    []float64
   stride int

}

func matrixFromRows(rows [][]float64) *matrix {

   if len(rows) == 0 {
       return &matrix{nil, 0}
   }
   m := &matrix{make([]float64, len(rows)*len(rows[0])), len(rows[0])}
   for rx, row := range rows {
       copy(m.ele[rx*m.stride:(rx+1)*m.stride], row)
   }
   return m

}

func (m *matrix) print(heading string) {

   if heading > "" {
       fmt.Print(heading, "\n")
   }
   for e := 0; e < len(m.ele); e += m.stride {
       fmt.Printf("%6.3f ", m.ele[e:e+m.stride])
       fmt.Println()
   }

}

func (m1 *matrix) multiply(m2 *matrix) (m3 *matrix) {

   m3 = &matrix{make([]float64, (len(m1.ele)/m1.stride)*m2.stride), m2.stride}
   for m1c0, m3x := 0, 0; m1c0 < len(m1.ele); m1c0 += m1.stride {
       for m2r0 := 0; m2r0 < m2.stride; m2r0++ {
           for m1x, m2x := m1c0, m2r0; m2x < len(m2.ele); m2x += m2.stride {
               m3.ele[m3x] += m1.ele[m1x] * m2.ele[m2x]
               m1x++
           }
           m3x++
       }
   }
   return m3

}</lang> Output (manually wrapped at 80 columns.)

exponent: 31415
addition chain: [1 2 4 5 10 20 25 50 55 110 220 245 490 980 1960 3920 7840
15680 31360 31415]
a^e
[ 0.707  0.000  0.000 -0.707  0.000  0.000] 
[ 0.000  0.707  0.707  0.000  0.000  0.000] 
[ 0.707  0.000  0.000  0.707  0.000  0.000] 
[ 0.000  0.707 -0.707  0.000  0.000  0.000] 
[ 0.000  0.000  0.000  0.000  0.000  1.000] 
[ 0.000  0.000  0.000  0.000  1.000  0.000] 
count of multiplies: 19

exponent: 27182
addition chain: [1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568
27136 27182]
a^e
[-0.500 -0.500 -0.500  0.500  0.000  0.000] 
[ 0.500 -0.500 -0.500 -0.500  0.000  0.000] 
[-0.500 -0.500  0.500 -0.500  0.000  0.000] 
[ 0.500 -0.500  0.500  0.500  0.000  0.000] 
[ 0.000  0.000  0.000  0.000  1.000  0.000] 
[ 0.000  0.000  0.000  0.000  0.000  1.000] 
count of multiplies: 18

exponent: 853922530
addition chain: [1 2 3 6 7 10 17 34 68 78 95 173 346 441 614 1055 2110 3165
3779 4834 9668 19336 24170 48340 52119 104238 208476 416952 833904 1667808
3335616 6671232 13342464 26684928 53369856 106739712 213479424 426958848
853917696 853922530]
a^e
[-0.500  0.500 -0.500  0.500  0.000  0.000] 
[-0.500 -0.500 -0.500 -0.500  0.000  0.000] 
[-0.500 -0.500  0.500  0.500  0.000  0.000] 
[ 0.500 -0.500 -0.500  0.500  0.000  0.000] 
[ 0.000  0.000  0.000  0.000  1.000  0.000] 
[ 0.000  0.000  0.000  0.000  0.000  1.000] 
count of multiplies: 39

exponent: 31415
1.00002206445416^31415: 1.99999999989447
count of multiplies: 19

exponent: 27182
1.00002550055251^27182: 1.99999999997876
count of multiplies: 18

MATLAB / Octave

I assume that Matlab and Octave have about such optimization already included. On a single core of an "Pentium(R) Dual-Core CPU E5200 @ 2.50GHz", the computation of 100000 repetitions of the matrix exponentiation A^27182 and A^31415 take about 2 and 2.2 seconds, resp. <lang Matlab>x = sqrt(.5); A = [x,0,x,0,0,0;0,x,0,x,0,0; 0,x,0,-x,0,0;-x,0,x,0,0,0;0,0,0,0,0,1; 0,0,0,0,1,0];A t = cputime(); for k=1:100000,x1=A^27182;end; cputime()-t,x1, t = cputime(); for k=1:100000,x2=A^31415;end; cputime()-t,x2, </lang> Output:

Octave3.4.3>  t = cputime(); for k=1:100000,x1=A^27182;end; cputime()-t,x1,
ans =  1.9900
x1 =

  -0.50000  -0.50000  -0.50000   0.50000   0.00000   0.00000
   0.50000  -0.50000  -0.50000  -0.50000   0.00000   0.00000
  -0.50000  -0.50000   0.50000  -0.50000   0.00000   0.00000
   0.50000  -0.50000   0.50000   0.50000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   1.00000

Octave3.4.3> t = cputime(); for k=1:100000,x2=A^31415;end; cputime()-t,x2,
ans =  2.2000
x2 =

   0.70711   0.00000   0.00000  -0.70711   0.00000   0.00000
   0.00000   0.70711   0.70711   0.00000   0.00000   0.00000
   0.70711   0.00000   0.00000   0.70711   0.00000   0.00000
   0.00000   0.70711  -0.70711   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   1.00000
   0.00000   0.00000   0.00000   0.00000   1.00000   0.00000