An Efficient Simulated Annealing Algorithm for Short Addition Sequence

1. Introduction

Given a set of numbers $N=\left\{n_1,n_2,\dots ,n_k\right\}$ such that ${1<n}_1<n_2<\dots <n_k.$ An addition sequence [1,2] for the set N, denoted by ASeq(N), is an increasing sequence of numbers ${as}_0,{as}_1,{as}_2,\dots ,{as}_l$ such that (1) ${as}_0=1,$ (2) ${as}_l=n_k,$ (3) ${as}_i={as}_j+{as}_h,0\le j,h\le i-1,$ and (4) $N\subseteq \left\{{as}_0,{as}_1,{as}_2,\dots ,{as}_l\right\},$ i.e., each number $n_i$ should appear in the sequence ${as}_0,{as}_1,{as}_2,\dots ,{as}_l.$ The number l is called the length of ASeq(N). The minimal length of an ASeq(N) is denoted by $\mathcal{l}\left(N\right).$ In the case of $k=1,$ the sequence is called addition chain [1,2].

The problem of generating a shortest ASeq(N) is equivalent to the simultaneous evaluation of k power monomials $g^{n_1}, g^{n_2}, \dots ,g^{n_k}$ with a minimum number of multiplications.

For example, let$N=\left\{53, 163, 203, 363\right\}.$ The following are two ASeqs with lengths 15 and 13, respectively.

The elements of the first ASeq are: 1, 2=1+1, 3=2+1, 6=3+3, 12=6+6, 13=12+1, 26=13+13, 39=26+13, 40=39+1, 53=40+13, 106=53+53, 159=106+53, 160=159+1, 163=160+3, 203=163+40, 363=203+160.

The elements of the second ASeq are: 1, 2=1+1, 3=2+1, 5=3+2, 10=5+5, 13=10+3, 20=10+10, 40=20+20, 53=40+13, 80=40+40, 160=80+80, 163=160+3, 203=163+40, 363=203+160

The evaluation of $g^{53},g^{163},g^{203},g^{363},$ using the first sequence is

$$g,g^2,g^3,g^6,g^{12},g^{13},g^{26},g^{39},g^{40},g^{53},g^{106},g^{159},g^{160},g^{163},g^{203},g^{363}$$

while the evaluation of the same powers $g^{53},g^{163},g^{203},g^{363}$ using the second sequence is

$$g,g^2,g^3,g^5,g^{10},g^{13},g^{20},g^{40},g^{53},g^{80},g^{160},g^{163},g^{203},g^{363}$$

A step i is called star if ${as}_i={as}_{i-1}+{as}_h, 0\le h\le i-1;$ and non-star if ${as}_i={as}_j+{as}_h, 0\le j,h\le i-2$. In case of $j=h=i-1, {as}_i=2{as}_{i-1},$ the step is called doubling. If all steps in the sequence are stars, then the sequence is called star. If ${\mathcal{l}}^{*}\left(N\right)$ denotes to the minimal length of star ASeq(N), then we have

$$\mathcal{l}\left(N\right)\le {\mathcal{l}}^{\mathrm{*}}\left(N\right)$$

(1)

Yao [3] showed that:

$$\mathcal{l}\left(N\right)\le \log n_k+(c k)\log n_k/\log \log n_k$$

(2)

for some constant $c=2+4/\sqrt{\log n_k}$.

Bleichenbacher [4] computed the lower bound

$$\mathcal{l}\left(N\cup \left\{n_{k+1}\right\}\right)\ge \mathcal{l}\left(N\right)+\alpha +1,$$

(3)

where $n_{k+1}>{2^{\alpha }n}_k,\alpha \ge 0$

ASeqs have received a lot of consideration among mathematicians and computer scientists for the following reasons:

The first reason is that one of the fundamental operations that play an important role in the efficiency of many public key cryptosystems and protocols is group exponentiation (sometimes it is called multi-modular exponentiation [2,5]), i.e., computing $g^{n_1}, g^{n_2}, \dots ,g^{n_k}$ simultaneously with a minimal number of operations, where g is an element in a group. Designing a fast algorithm for generating a shortest (or short) ASeq increases the efficiency of such public key cryptosystems and protocols since evaluating $g^{n_1}, g^{n_2}, \dots ,g^{n_k}$ with a minimal number of multiplications is equivalent to finding a shortest ASeq(N).

The second reason is that ASeqs (including addition chains) are generalized to the following:

(i)

B-chains [7], where every element in the B-chain has the form $a_i=a_j o a_h, 0\le j,h\le i-1,$ and the binary operation o belongs to a finite set of binary operations over the set of natural numbers B, i.e. $o\in B=\left\{+,-,*,÷\right\}$. Guzmán-Trampe et. al. [8] proposed a method for generating addition-subtraction, $B=\left\{+,-\right\},$ sequence for the Kachisa–Schaefer–Scott family of pairing-friendly elliptic curves.

(ii)

Vectorial addition chain [9,10]: it is a sequence of k-dimensional vectors of nonnegative integers $v_i,-k+1\le i\le l$, such that (1) $v_{-k+1}=[\mathrm{1,0},0\dots ,\mathrm{0,0}]$, $v_{-k+2}=[\mathrm{0,1},0\dots ,\mathrm{0,0}]$, …, $v_0=[\mathrm{0,0},0\dots ,\mathrm{0,1}]$, (2) $v_i=v_j+v_h,1\le i\le l,-k+1\le j,h\le i-1$, and (3) $v_l=\left[n_1,n_2,\dots ,n_k\right].$ Finding a shortest vectorial addition chain is equivalent to evaluating, multi-exponentiation, i.e., the product $\prod _{i=1}^k{g_i}^{n_i}$ with the minimal number of multiplications.

The third reason is that in Internet of Things, IoT, devices with limited resources have a problem when they perform some public-key primitives, such as decryption and signature, because most public-key primitives are (i) time-consuming compared with symmetric-key cryptosystems; and (ii) using private information. One of the common solutions to this problem is to use what is called “server aided secret computation protocols”, denoted by SASCP [11], or sometimes it is called outsourcing protocols [12].

In SASCP, devices with limited power and resources, such as smart cards, can execute public-key primitives efficiently with the aid of an untrusted powerful server without revealing the private information. Examples of such protocols that used ASeqs are [6,11,13]. Other protocols and their security analysis are [12,14]. Another and similar solution to the problem is to define a delegation protocol. It is a protocol that satisfies two security objectives: (i) privacy: the private information should not be recovered by a passive attacker; and (ii) verifiability: the untrusted server should not be able to make the devices accept an invalid value as the result of the delegated computation. Examples of such protocols and their security analysis are [12,15,16,17].

The main challenge of finding a shortest ASeq, i.e., the minimal number of additions needed to compute all elements of N, that it is NP-complete [4]. Additionally, when the size of N is large and the size of exponents is large, the running time for finding ASeq is very large [18]. Therefore, designing a fast algorithm for generating a short (not necessarily shortest) ASeq is interesting using metaheuristics techniques such as simulated annulling, ant colony and evolutionary algorithms.

In this paper, a new metaheuristic algorithm based on simulated annealing strategy is proposed to find a short ASeq for the set N. The proposed simulated annulling algorithm for ASeq has three advantages over the previous ASeq algorithms. The first advantage is the designed algorithm has running time less than the exact algorithm. The second advantage is the length of ASeq generated by the proposed algorithm is shorter than the ASeq generated by previous heuristic algorithms. The third advantage is that there is no comparative study between suboptimal algorithms and the exact algorithm in terms of the length of ASeq.

The remainder of this paper is organized as follows. Section 2 includes the related works of ASeq. In Section 3, the details of the proposed algorithm is given. Section 4 includes the dataset used in the experiments, and the results and analysis of the experimental studies for the proposed simulated annulling algorithm and other algorithms. Finally, Section 5 includes the conclusion of this paper and the future works.

2. Related Works

Algorithms for generating ASeq, i.e. $k\ge 1,$ can be classified into two categories. The first category is to find a shortest ASeq. In fact, there are a few papers that discussed a generation of shortest ASeqs. Bleichenbacher [4] suggested an algorithm to find a shortest ASeq(N) with length r provided that we previously computed $\mathcal{l}\left(y\right)$ for all numbers $y<maximum\left(N\right),$ and $\mathcal{l}\left(y\right)<r.$ He used the suggested algorithm to generate a shortest addition chain up to a certain number.

The authors in [18] suggested a branch and bound depth-first search algorithm to generate a shortest ASeq for any set N. The algorithm starts by computing a lower bound, Eq(3), and looking for an addition chain for the first element $n_1$ in the set N. Then, it extends the chain to addition sequence for ${\{n}_1, n_2\}$, and so on until it generates ASeq(N). The algorithm uses different strategies to speed up the generation as follows. (i) Using bounding sequences to prune some branches in the search tree, which cannot lead to a shortest ASeq. (ii) Determining an upper bound of $\mathcal{l}\left(n_1, n_2, \dots , n_i\right), 1\le i\le k.$ (iii) Using some sufficient conditions for star steps to skip the generation of non-star steps. (vi) If no ASeq(N) of length l is found, then the algorithm increases l by one and repeats the process until either l is equal to the length of the generated short ASeq produced by continued fraction (CF) method [19] or the algorithms finds a shortest ASeq. Recently, the authors in [20] used multicore systems to improve the generation of a shortest ASeq.

The second category is to find a short ASeq. Yao [3] presented an algorithm to compute $g^{n_1}, g^{n_2}, \dots ,g^{n_k}$ in $O(\lg n_k+c\sum _{i=1}^k(\log n_i/loglog(n_i+2))$ multiplications for some constant c. Bos and Coster [21] proposed four methods to generate a short ASeq to use it in the window method [21]. The upper bound of the length of generated an ASeq, by the four methods, could be estimated, experimentally, by

$$\mathcal{l}\left(n_1,n_2,\dots ,n_k\right)\le \frac{3}{2}\log n_k+k+1,$$

(4)

for $n_k\le 1000.$

Bergeron et al. [19] proposed an efficient method based on CF. The suggested method can be considered as an extension and unifying approach of some previously known methods (such as binary and k-ary methods [1]) for generating a short addition chain, i.e., ASeq with k=1.

Recall that the CF expansion of $n/d,$ denoted by $[c_1,$ …, $c_t],$ is

$$\frac{n}{d}=c_t+\frac{1}{c_{t-1}+\frac{1}{\begin{array}{c} \\ \ddots \\ +\frac{1}{\begin{array}{c}{c}_2+\frac{1}{c_1} \end{array}}\end{array}}}$$

(5)

where d is an integer in the range $\left[2..n-1\right].$

Bergeron et al. [19] suggested different strategies for choosing the value of d. One of the efficient strategies that produces a good suboptimal ASeq is dichotomic strategy, where

$$d=⌊\frac{n}{2^{⌈⌊{log}_{2}n⌋/2⌉}}⌋.$$

(6)

Let $N^{\prime }=\left\{n_{\mathrm{k}},n_{\mathrm{k}-1},\dots ,n_1\right\}, \mathrm{a}\mathrm{n}\mathrm{d}\mathcal{L}\left(N\prime \right)$ denotes the length of the ASeq($N^{\prime }$) generated by CF using dichotomic strategy. Then

$$\mathcal{L}\left(N^{\prime }\right)=\left\{\begin{array}{c}L\left(N^{\prime }\setminus \left\{n_k\right\}\right)+l\left(q\right), if rem=0;\\ L\left(N^{\prime }\setminus \left\{n_k\right\}\right)+l\left(q\right)+1, if rem=\mathrm{1,2}\\ L\left(N^{\prime }\cup \left\{rem\right\}\setminus \left\{n_k\right\}\right)+l\left(q\right)+1, otherwise\end{array}\right.;$$

(7)

where $n_k=q{n}_{k-1}+rem,$ and

$$\mathcal{l}\left(n\right)=\left\{\begin{array}{c}\alpha , if n=2^{\alpha };\\ 3, if n=3;\\ L\left(\{n,d\}\right), otherwise\end{array}\right.$$

(8)

where d is defined by Eq. (6).

Enge et. al. [22] proposed a special method to construct a short ASeq to find the first k nonzero terms in the sparse q-series belonging to the Dedekind eta function or the Jacobi theta constants. Nadia and Mourelle [23] used Anti Colony strategy to find a short ASeq. They tested the strategy for a small set of numbers. Abbas and Gustafsson [24] proposed a method based on integer linear programming to generate a short ASeq for a small set of numbers.

In all previous studies, there was not enough experimental study for generating a short ASeq with different sizes of the set N, or with different range values of each number in the set N. Also, there is no comparative study between suboptimal algorithms and the exact algorithm in terms of the length of ASeq.

3. The Proposed Method

In this section, we first present a brief description of the proposed algorithm that is based on simulated annealing strategy to find short ASeq, and then we present its details. The algorithm is named SAAS for simulated annealing addition sequence.

Initially, the algorithm starts by generating the initial state, ${AS}_0$, using the CF method [19], and its energy is equal to the length of ${AS}_0$, $l_{{AS}_0}$. Then, the algorithm assigns these two values to the best state and the best energy, respectively. After that, the algorithm repeats the following steps based on the number of Metropolis cycles, m, for a fixed temperature. In each iteration of this loop, the algorithm performs the following steps:

The first step is generating a new state, ${AS}_{new}$, and its energy, $l_{new}$. The second step is determining whether the algorithm accepts this new state or not. The algorithm accepts the new state and its energy, and then assigns these values to the best state and best energy if either of the following conditions is true. (1) if the energy of the new state is lower than the energy of the best state. (2) If the Boltzmann distribution is greater than a random real number in the range [0,1].

After completing the number of Metropolis cycles for a fixed temperature, the algorithm updates the temperature using the Kirkpatrick quenching method and repeats this process until it reaches the maximum number of annealing iterations.

The details of the algorithm steps are as follows.

Step 1: Generate the initial state, $AS$, using CF method for the set of exponents $N=\left\{n_1,n_2,\dots ,n_k\right\}$, where $AS=\left\{{as}_0,{as}_1,{as}_2,\dots ,{as}_l\right\}$ such that (1) ${as}_0=1$ and ${as}_1=2$, (2) $\exists i,l_i$ s.t. $n_i={as}_{l_i}$ and $1\le i\le k$, (3) $L=\left\{l_1,l_2,\dots ,l_k\right\}$ such that $l_i<l_{i+1}$ and ${l=l}_k$.

Step 2: Repeat the following m times:

Step 2.1: Generate a random integer number, r, from the range $\left[0,k-1\right]$. This number will be used as a start point of mutation based on the elements of N.

Step 2.2: Generate a new state, ${AS}_{new}$, by mutate $AS$, from the location $l_r$. If r=0, then the algorithm mutates $AS$ from the element ${as}_1=2$. Otherwise, the algorithm mutates the state $AS$ from ${as}_{l_r}=n_r$ to $n_k$. The process of generating the new elements from $n_i$ to $n_{i+1}$ is based on the following rules.

• Rule # 1: Doubling the current element, ${as}_{j+1}=2 {as}_j$.
• Rule # 2: Summing the last two elements, ${as}_{j+1}={as}_j+{as}_{j-1}$.
• Rule # 3: Summing the last element with any other random element in the sequence, ${as}_{j+1}={as}_j$+${as}_h$, $0\le h<j$.

This step can be done as follows (Steps 2.2.1-2.2.3).

Step 2.2.1 (Generate one element in the sequence): If the current goal is $n_{i+1}$ and the current ASeq is $\left\{{as}_0, {as}_1, \dots , {{as}_{l_i}=n}_i,{as}_{l_i+1},\dots ,{as}_{l_i+j}\right\}, j\ge 0$, then the steps of generating a new element in the chain are as follows.

•  $d=n_{i+1}-{as}_{l_i+j}$
•  If $d={as}_{l_i+j}$ then apply rule # 1
•  Else if $d={as}_{l_i+j-1}$ then apply rule # 2
•    Else if $d>{as}_{l_i+j}$ then
•       Generate a random real number $\alpha \in \left[\mathrm{0,1}\right]$
•       If $\alpha \ge 0.5$ then apply the rule # 1
•       Else
•        Generate a random real number $\alpha \in \left[\mathrm{0,1}\right]$
•        If $\alpha \ge 0.5$ then apply the rule # 2
•        Else
•         Generate a random integer number $r\in [0,l_i+j-2]$
•         Apply the rule # 3, where h=r.
•      Else // $d<{as}_{l_i+j}$
•     Generate a random integer number $r\in [0,l_i+j-2]$
•     Apply the rule # 3, where h=r.
•     If the new element is less than or equal to $n_{i+1}$ then the element is
•     accepted. Otherwise, decrease the value of r and apply rule #3 until
•     we found a certain value of h such that the new element is less than or
•     equal to $n_{i+1}$.

Step 2.2.2 (Generate all elements between $n_i$ and $n_{i+1}$): Repeat Step 2.2.1 starting from $j=0$, and ${{as}_{l_i}=n}_i$, until the algorithm finds ${{as}_{l_i+j_i}=n}_{i+1}$. In this case, the algorithm updates the value of $l_{i+1}=l_i+j_i, 1\le j_i.$

Step 2.3.3 (Generate the ASeq from $n_r$ to $n_k$): Repeat Steps 2.2.1 and 2.2.2 until generate $n_k$. Therefore,${AS}_{new}=\left\{{as}_0,{as}_1,\dots ,n_r={as}_{l_r},{as}_{l_r+1}^{\mathrm{’}},\dots ,n_{i+1}={as}_{l_i+j_i}^{\mathrm{’}},\dots ,n_k={as}_{l_{k-1}+j_{k-1}}^{\mathrm{’}}={as}_{j_{k-1}}^{\mathrm{’}}\right\},j_i\ge 1$; and $L_{new}=\left\{l_1,l_2,\dots ,l_r,l_{r+1}^{\mathrm{’}},l_{r+2}^{\mathrm{’}},\dots ,l_k^{\mathrm{’}}\right\}$.

Step 3: Test the acceptance of the new state by the following substeps.

• If $l_k^{\mathrm{’}}<l$ then $AS={AS}_{new}$ and $l=l_k^{\mathrm{’}}$
• Else generate a random real number $r$.
•  $d_e=l_k^{\mathrm{’}}-l$
•  If $e^{-d_e/T}>r$ then $AS={AS}_{new}$ and $l=l_k^{\mathrm{’}}$

Step 4: Decrease the temperature using Kirkpatrick quenching method: $T=\gamma T$, where $\gamma =0.99$.

Step 5: Repeat Steps 2, 3, and 4 until reach the maximum number of annealing iterations.

The complete pseudocodes for the new proposed SAAS is given in Algorithms 1 and 2.

Algorithm 1: SAAS

Algorithm 2: MutateAS

3.1. Data Generation

The data used in the experimental study is based on two factors. The first factor is the number of elements k in the set of exponents $N$. The experimental values of k are 2, 4, 6, 8, and 10. The second factor is the domain of each exponent in the set N. According to the window method and its variations, the range of exponents is the integer interval $\left[1, 2^e-1\right]$, where e is the window length (of size e-bits). Also, according to the performance of the window method, the value of each exponent should be odd. The experimental values of e are equal to 7, 8, 9, and 10. The reason for starting the values of e with 7, the running times for all compared algorithms are fast when $e<7.$

The methodology of generating the ASeq is based on fixing the size of the window, i.e., e-bits, say e=7, and then generating different sets $N_{k,e}$ with lengths k= 2, 4, 6, 8, and 10. For each value of k, the algorithm generates 25 sets of exponents in the range$\left[1, 2^e-1\right]$. The process of generating different sets of exponents is as follows.

• Set e to the maximum number of bits in the exponents, i.e., the window size.
• Set the set $N_{0,e}=\mathrm{\varnothing }$ and i=2.
• While $i\le k=10$ do the following
• Construct a new set $N_{i,e},$ by adding two randomly odd numbers, in the range $\left[1..2^e-1\right],$ to the set $N_{i-2,e},i.e.,N_{i,e}=N_{i-2,e}\cup${ the two generated randomly odd numbers}
• Set i=i+2.
• Make sure that $N_{i,e}$ is sorted.
• Repeat Steps 2-4, 25 times to generate 25 sets of exponents with at most e-bits.
• Repeat Steps 1-5 for different size of exponents e=7, 8, 9, and 10.

The following example illustrates the generation of five sets with different values of k and fixed size of exponents e=8.

• $N_{\mathrm{2,8}}=\left\{177,241\right\}.$
• $N_{\mathrm{4,8}}=\left\{65,125,177,241\right\}.$
• $N_{\mathrm{6,8}}=\left\{65,89,125,177,189,241\right\}.$
• $N_{\mathrm{8,8}}=\left\{43,65,89,125,177,189,221,241\right\}.$
• $N_{\mathrm{10,8}}=\left\{43,65,89,103,125,177,189,203,221,241\right\}$

3.2. Results

The results of implementing the three algorithms on the generated data in terms of the length of the output are shown in Table 1. The first two columns represent the two factors e and k, while the three last columns represent the percentage of differences in the lengths of the output for the following cases: (1) ExAS and SAAS algorithms, (2) ExAS and CFAS algorithms, and (3) SAAS and CFAS algorithms. Since the exact algorithm always produces the shortest ASeq, the methodology for analyzing the results is computing the number of cases in which the lengths of ASeqs generated by the SAAS and CFAS algorithms are greater than the shortest ASeq generated by the ExAs algorithm. The percentages of these cases represent the third and fourth columns. Also, Table 1 presents the difference between the lengths of the ASeqs generated by the SAAS algorithm and those generated by the CFAS algorithm, see the, the last column in Table 1.

The analysis of data results shows the following observations.

First, as in Table 1, the percentage of the difference between the lengths of the ASeq generated by the exact algorithm, ExAS, and the heuristic algorithms, SAAS and CFAS, increases with the increase in the number of elements in the set N. For example, for fixed e=7 and k=2, 4, 6, 8, and 10, the percentages of cases that the exact algorithm generates ASeq with a length less than that generated by the SAAS algorithm are 12%, 40%, 56%, 76%, and 88%. Similarly, for the CFAS algorithm, the differences are 28%, 64%, 76%, 84%, and 92%.

Second, as in Table 1, the comparison between the lengths of ASeqs generated by the SAAS and CFAS algorithms, independent of the ExAS algorithm, , is presented in the last column. The data shows that the SAAS algorithm outperforms the CFAS algorithm in terms of the short ASeq for all studied cases.

Third, the length of the output generated by the SAAS algorithm is near to the shortest length compared to that generated by the CFAS algorithm. Figure 1 shows the distribution of difference between the length of the output for the SAAS algorithm (similarly the CFAS algorithm) and the output of the exact algorithm. It is clear that the SAAS algorithm generates short ASeq with lengths that are near to the shortest ASeq than that generated by the CFAS algorithm. For example, when e=8 and k=2, there are 20% of the instances where the length of ASeq generated by the SAAS algorithm is greater by one than the length of ASeq generated by the ExAS algorithm. On the other side, using the CFAS algorithm, there are 44% and 20% of instances have lengths greater than the shortest by one and two, respectively.

The comparison between the three algorithms, ExAS, SAAS, and CFAS, in terms of execution time is shown in Table 2. The analysis of data in the table demonstrates the following observations. (1) The fastest running time for all compared algorithms is CFAS algorithm. (2) The CFAS algorithm is not affected by the values of e and k in general. On the other side, the SAAS algorithm is slightly affected by increasing e and k, whereas the ExAS algorithm is significantly affected by increasing e and k. (3) The running time for the SAAS algorithm is affected by the two parameters, succNo and metropolis. The increase in the values of two parameters leads to a slight increase in the running time. (4) The running time for SAAS algorithm is faster than the exact algorithm, and the difference between the two algorithms in running time increases with increase in e and k. (5) The last column of Table 2 shows the percentage improvement for the SAAS algorithm compared to the exact algorithm.

5. Conclusion and Future Works

In this paper, finding a short addition sequence for a set of positive integers was studied. A new metaheuristic algorithm was proposed to find an addition sequence with short length. The proposed algorithm starts with generating sequence using continued fraction and then apply the simulated annealing strategy to improve the length of the sequence. The proposed algorithm is fast compared to the exact algorithm and able to generate addition sequence with length less than the previous heuristic algorithm.

The efficiency of the proposed algorithm was conducted with considering different parameters such as the number of elements in the set and the size of positive integer.

There are many research directions related to addition sequence such as (1) extend the concept of B-chains and verctorial chain to ASeq, (2) use high-performance system to accelerate the computation of ASeq, and (3) accelerating the multi-modular exponentiation used ASeq.