Anonymous user
Addition-chain exponentiation: Difference between revisions
→{{header|C}}
(fixed error in 2↑15; Also: changed M, initial choice was a mistake, particularly as det M!=1, new M has det of 1) |
|||
Line 83:
[[wp:Kudoc|Kudos]] (κῦδος) for providing a routine that generate sequence A003313 in the output.
=={{header|C}}==
Using complex instead of matrix. Requires [[Addition-chain exponentiation/Achain.c|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 main()
{
cplx r1 = {1.0000254989, 0.0000577896},
r2 = {1.0000220632, 0.0000500026};
init();
puts("Precompute chain lengths");
seq_len(31415);
chain_expo(r1, 27182);
chain_expo(r2, 31415);
return 0;
}</lang>output<lang>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</lang>
=={{header|J}}==
| |||