| This page uses content fromWikipedia. The original article was atAddition-chain exponentiation. The list of authors can be seen in thepage history.As with Rosetta Code, the text of Wikipediais available under theGNU FDL. (See links for details on variance) |
In cases of special objects (such as withmatrices) the operation of multiplication can be excessively expensive. In these cases the operation of multiplication should be avoided or reduced to a minimum.
Inmathematics andcomputer science, optimaladdition-chain exponentiation is a method ofexponentiation by positiveinteger powers that requires a minimal number of multiplications. It works by creating a shortestaddition 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-chainalgorithm requires no more multiplications thanbinary exponentiation and usually less. The first example of where it does better is for, where the binary method needs six multiplies but a shortest addition chain requires only five:
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 provenNP-complete.
| Number of Multiplications | Actual Exponentiation | Specific implementation of Addition Chains to do Exponentiation |
|---|---|---|
| 0 | a1 | a |
| 1 | a2 | a × a |
| 2 | a3 | a × a × a |
| 2 | a4 | (a × a→b) × b |
| 3 | a5 | (a × a→b) × b × a |
| 3 | a6 | (a × a→b) × b × b |
| 4 | a7 | (a × a→b) × b × b × a |
| 3 | a8 | ((a × a→b) × b→d) × d |
| 4 | a9 | (a × a × a→c) × c × c |
| 4 | a10 | ((a × a→b) × b→d) × d × b |
| 5 | a11 | ((a × a→b) × b→d) × d × b × a |
| 4 | a12 | ((a × a→b) × b→d) × d × d |
| 5 | a13 | ((a × a→b) × b→d) × d × d × a |
| 5 | a14 | ((a × a→b) × b→d) × d × d × b |
| 5 | a15 | ((a × a→b) × b × a→e) × e × e |
| 4 | a16 | (((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: and
Repeat taskMatrix-exponentiation operator, except useaddition-chain exponentiation to better calculate:
As an easier alternative to doing the matrix manipulation above, generate theaddition-chains for 12509, 31415 and 27182 and useaddition-chain exponentiation to calculate these two equations:
Also: Display a count of how many multiplications were done in each case.
Note: There are two ways to approach this task:
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.
Also, see the Rosetta Code task: [http://rosettacode.org/wiki/Knuth%27s_power_treeKnuth's power tree].
Using complex instead of matrix. RequiresAchain.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.
#include<stdio.h>#include"achain.c" /* not common practice *//* don't have a C99 compiler atm */typedefstruct{doubleu,v;}cplx;inlinecplxc_mul(cplxa,cplxb){cplxc;c.u=a.u*b.u-a.v*b.v;c.v=a.u*b.v+a.v*b.u;returnc;}cplxchain_expo(cplxx,intn){inti,j,k,l,e[32];cplxv[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);returnx;}intbin_len(intn){intr,o;for(r=o=-1;n;n>>=1,r++)if(n&1)o++;returnr+o;}intmain(){cplxr1={1.0000254989,0.0000577896},r2={1.0000220632,0.0000500026};intn1=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));return0;}
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.000001Exponents: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.000000chain 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 ...
ConstNAsInteger=32ConstNMAXAsInteger=40000TypeSaveTypepAsIntegerPtrvAsIntegerEndType' ArraysDimSharedAsIntegeru(N-1)={1,2}' upper boundsDimSharedAsIntegerl(N-1)={1,2}' lower boundsDimSharedAsIntegerout_(N-1),sum(N-1),tail(N-1),cache(NMAX)cache(2)=1DimSharedAsIntegerstack=0,known=2DimSharedAsSaveTypeundo(N*N-1)DeclareFunctionseq(nAsInteger,leAsInteger,bufAsIntegerPtr)AsIntegerDeclareFunctionseqRecur(leAsInteger)AsBooleanSubreplace(x()AsInteger,iAsInteger,tAsInteger)undo(stack).p=@x(i)undo(stack).v=x(i)x(i)=tstack+=1EndSubSubrestore_(tAsInteger)Whilestack>tstack-=1*undo(stack).p=undo(stack).vWendEndSub' lower and upper boundsFunctionlower(tAsInteger,upAsIntegerPtr)AsIntegerIft<=2Orelse(t<=NMAXAndalsocache(t)<>0)ThenIfup<>0Then*up=cache(t)Returncache(t)EndIfDimAsIntegeri=-1,o=0,tmp=tWhiletmp<>0If(tmpAnd1)<>0Theno+=1tmpShr=1i+=1WendIfup<>0Theni-=1*up=o+iEndIfDoi+=1oShr=1LoopUntilo=0ReturniEndFunctionFunctioninsert(xAsInteger,posicAsInteger)AsBooleanDimAsIntegersave=stackDimAsIntegeri,tIfl(posic)>xOrelseu(posic)<xThenReturnFalseIfl(posic)=xThenGotoreplUreplace(l(),posic,x)i=posic-1Whileu(i)*2<u(i+1)t=l(i+1)+1Ift*2>u(i)ThenGotobailreplace(l(),i,t)i-=1Wendi=posic+1Whilel(i)<=l(i-1)t=l(i-1)+1Ift>u(i)ThenGotobailreplace(l(),i,t)i+=1WendreplU:Ifu(posic)=xThenReturnTruereplace(u(),posic,x)i=posic-1Whileu(i)>=u(i+1)t=u(i+1)-1Ift<l(i)ThenGotobailreplace(u(),i,t)i-=1Wendi=posic+1Whileu(i)>u(i-1)*2t=u(i-1)*2Ift<l(i)ThenGotobailreplace(u(),i,t)i+=1WendReturnTruebail:restore_(save)ReturnFalseEndFunctionFunctiontry(pAsInteger,qAsInteger,leAsInteger)AsBooleanDimAsIntegerpl=cache(p)Ifpl>=leThenReturnFalseDimAsIntegerql=cache(q)Ifql>=leThenReturnFalseDimAsIntegerpu,quWhilepl<leAndalsou(pl)<ppl+=1Wendpu=pl-1Whilepu<le-1Andalsou(pu+1)>=ppu+=1WendWhileql<leAndalsou(ql)<qql+=1Wendqu=ql-1Whilequ<le-1Andalsou(qu+1)>=qqu+=1WendIfp<>qAndalsopl<=qlThenpl+=1Ifpl>puOrelseql>quOrelseql>puThenReturnFalseIfout_(le)=0Thenpu=le-1pl=puEndIfDimAsIntegerps=stackWhilepu>=plIfinsert(p,pu)Thenout_(pu)+=1sum(pu)+=leIfp<>qThenDimAsIntegerqs=stackDimAsIntegerj=quIfj>=puThenj=pu-1Whilej>=qlIfinsert(q,j)Thenout_(j)+=1sum(j)+=letail(le)=qIfseqRecur(le-1)ThenReturnTruerestore_(qs)out_(j)-=1sum(j)-=leEndIfj-=1WendElseout_(pu)+=1sum(pu)+=letail(le)=pIfseqRecur(le-1)ThenReturnTrueout_(pu)-=1sum(pu)-=leEndIfout_(pu)-=1sum(pu)-=lerestore_(ps)EndIfpu-=1WendReturnFalseEndFunctionFunctionseqRecur(leAsInteger)AsBooleanDimAsIntegert=l(le)Ifle<2ThenReturnTrueDimAsIntegerlimit=t-1Ifout_(le)=1Thenlimit=t-tail(sum(le))Iflimit>u(le-1)Thenlimit=u(le-1)DimAsIntegerq,p=limitWhilep>=(t+1)\2q=t-pIftry(p,q,le)ThenReturnTruep-=1WendReturnFalseEndFunctionFunctionseq(tAsInteger,leAsInteger,bufAsIntegerPtr)AsIntegerDimAsIntegeriIfle=0Thenle=seqLen(t)stack=0l(le)=t:u(le)=tFori=0Toleout_(i)=0:sum(i)=0NextFori=2Tole-1l(i)=l(i-1)+1u(i)=u(i-1)*2NextFori=le-1To3Step-1Ifl(i)*2<l(i+1)Thenl(i)=(1+l(i+1))\2Ifu(i)>=u(i+1)Thenu(i)=u(i+1)-1NextIfNotseqRecur(le)ThenReturn0Ifbuf<>0ThenFori=0Tolebuf[i]=u(i)NextEndIfReturnleEndFunctionFunctionseqLen(tAsInteger)AsIntegerIft<=knownThenReturncache(t)' Need all lower t to compute sequenceWhileknown+1<tseqLen(known+1)WendDimAsIntegerub,lb=lower(t,@ub)Whilelb<ubAndalsoseq(t,lb,0)=0lb+=1Wendknown=tIf(tAnd1023)=0ThenPrint"Cached";knowncache(t)=lbReturnlbEndFunctionFunctionbinLen(tAsInteger)AsIntegerDimAsIntegerr=-1,o=-1,tmp=tWhiletmp<>0If(tmpAnd1)<>0Theno+=1tmpShr=1r+=1WendReturnr+oEndFunctionSubprintMatrix(mAsDoublePtr)DimAsIntegeri,jFori=0To5Print"[";Forj=0To5PrintUsing"##.###### ";m[i*6+j];NextPrintChr(8);"]"NextPrintEndSubSubmatrixMul(m1AsDoublePtr,m2AsDoublePtr,resultAsDoublePtr)DimAsIntegeri,j,kFori=0To5Forj=0To5result[i*6+j]=0Fork=0To5result[i*6+j]+=m1[i*6+k]*m2[k*6+j]NextNextNextEndSubSubmatrixPow(mAsDoublePtr,tAsInteger,resultAsDoublePtr,printoutAsBoolean)DimAsIntegeri,j,kDimAsIntegere(N-1)DimAsDoublev(N-1,5,5)DimAsIntegerle=seq(t,0,@e(0))IfprintoutThenPrint"Addition chain:"Fori=0TolePrintStr(e(i));Iif(i=le,""," ");NextEndIf' Copy initial matrix to v(0)Fori=0To5Forj=0To5v(0,i,j)=m[i*6+j]NextNext' Calculate v(1) = m * mmatrixMul(@v(0,0,0),m,@v(1,0,0))Fori=2ToleForj=i-1To1Step-1Fork=jTo0Step-1Ife(k)+e(j)<e(i)ThenExitForIfe(k)+e(j)>e(i)ThenContinueFormatrixMul(@v(j,0,0),@v(k,0,0),@v(i,0,0))j=1ExitForNextNextNext' Copy resultFori=0To5Forj=0To5result[i*6+j]=v(le,i,j)NextNextEndSubSubmain()DimAsIntegeri,q=27182,r=31415DimAsDoublerh=Sqr(0.5)DimAsDoublemx(5,5)mx(0,0)=rh:mx(0,2)=rhmx(1,1)=rh:mx(1,3)=rhmx(2,1)=rh:mx(2,3)=-rhmx(3,0)=-rh:mx(3,2)=rhmx(4,5)=1mx(5,4)=1Print"Precompute chain lengths:"seqLen(r)Print!"\nThe first 100 terms of A003313 are:"Fori=1To100PrintStr(seqLen(i));" ";IfiMod10=0ThenPrintNextDimAsIntegerexs(1)={q,r}DimAsDoublemxs(1,5,5)' Store resultsFori=0To1Print!"\nExponent:"&exs(i)matrixPow(@mx(0,0),exs(i),@mxs(i,0,0),True)Print!"\nA ^ "&exs(i)&!":-\n"printMatrix(@mxs(i,0,0))Print"Number of A/C multiplies:";seqLen(exs(i))Print" c.f. Binary multiplies:";binLen(exs(i))NextPrint!"\nExponent: "&q&" x "&r&" = "&q*rPrint"A ^ "&q*r&" = (A ^ "&q&") ^ "&r&":-"&!"\n"DimAsDoublemx2(5,5)matrixPow(@mxs(0,0,0),r,@mx2(0,0),False)printMatrix(@mx2(0,0))EndSubmain()Sleep
Precompute Chain lengths:Cached 1024Cached 2048Cached 3072Cached 4096Cached 5120Cached 6144Cached 7168Cached 8192Cached 9216Cached 10240Cached 11264Cached 12288Cached 13312Cached 14336Cached 15360Cached 16384Cached 17408Cached 18432Cached 19456Cached 20480Cached 21504Cached 22528Cached 23552Cached 24576Cached 25600Cached 26624Cached 27648Cached 28672Cached 29696Cached 30720The first 100 terms of A003313 are:0 1 2 2 3 3 4 3 4 45 4 5 5 5 4 5 5 6 56 6 6 5 6 6 6 6 7 67 5 6 6 7 6 7 7 7 67 7 7 7 7 7 8 6 7 77 7 8 7 8 7 8 8 8 78 8 8 6 7 7 8 7 8 89 7 8 8 8 8 8 8 9 78 8 8 8 8 8 9 8 9 89 8 9 9 9 7 8 8 8 8Exponent:27182Addition Chain:1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568 27136 27182A ^27182:-[-0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 -0.500000 0.000000 0.000000][ 0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]Number of A/C multiplies: 18c.f. Binary multiplies: 21Exponent:31415Addition Chain:1 2 4 8 16 17 33 49 98 196 392 784 1568 3136 6272 6289 12561 25122 31411 31415A ^31415:-[ 0.707107 0.000000 0.000000 -0.707107 0.000000 0.000000][ 0.000000 0.707107 0.707107 0.000000 0.000000 0.000000][ 0.707107 0.000000 0.000000 0.707107 0.000000 0.000000][ 0.000000 0.707107 -0.707107 0.000000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000]Number of A/C multiplies: 19c.f. Binary multiplies: 24Exponent: 27182 x 31415 = 853922530A ^ 853922530 = (A ^ 27182) ^ 31415:-[-0.500000 0.500000 -0.500000 0.500000 0.000000 0.000000][-0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]
Though adjusted to deal with matrices rather than complex numbers.
Calculating A ^ (m * n) from scratch using this method would take 'for ever' so I've calculated it instead as (A ^ m) ^ n.
packagemainimport("fmt""math")const(N=32NMAX=40000)var(u=[N]int{0:1,1:2}// upper boundsl=[N]int{0:1,1:2}// lower boundsout=[N]int{}sum=[N]int{}tail=[N]int{}cache=[NMAX+1]int{2:1}known=2stack=0undo=[N*N]save{})typesavestruct{p*intvint}funcreplace(x*[N]int,i,nint){undo[stack].p=&x[i]undo[stack].v=x[i]x[i]=nstack++}funcrestore(nint){forstack>n{stack--*undo[stack].p=undo[stack].v}}/* lower and upper bounds */funclower(nint,up*int)int{ifn<=2||(n<=NMAX&&cache[n]!=0){ifup!=nil{*up=cache[n]}returncache[n]}i,o:=-1,0for;n!=0;n,i=n>>1,i+1{ifn&1!=0{o++}}ifup!=nil{i--*up=o+i}for{i++o>>=1ifo==0{break}}ifup==nil{returni}foro=2;o*o<n;o++{ifn%o!=0{continue}q:=cache[o]+cache[n/o]ifq<*up{*up=qifq==i{break}}}ifn>2{if*up>cache[n-2]+1{*up=cache[n-1]+1}if*up>cache[n-2]+1{*up=cache[n-2]+1}}returni}funcinsert(x,posint)bool{save:=stackifl[pos]>x||u[pos]<x{returnfalse}ifl[pos]==x{gotoreplU}replace(&l,pos,x)fori:=pos-1;u[i]*2<u[i+1];i--{t:=l[i+1]+1ift*2>u[i]{gotobail}replace(&l,i,t)}fori:=pos+1;l[i]<=l[i-1];i++{t:=l[i-1]+1ift>u[i]{gotobail}replace(&l,i,t)}replU:ifu[pos]==x{returntrue}replace(&u,pos,x)fori:=pos-1;u[i]>=u[i+1];i--{t:=u[i+1]-1ift<l[i]{gotobail}replace(&u,i,t)}fori:=pos+1;u[i]>u[i-1]*2;i++{t:=u[i-1]*2ift<l[i]{gotobail}replace(&u,i,t)}returntruebail:restore(save)returnfalse}functry(p,q,leint)bool{pl:=cache[p]ifpl>=le{returnfalse}ql:=cache[q]ifql>=le{returnfalse}varpu,quintforpl<le&&u[pl]<p{pl++}forpu=pl-1;pu<le-1&&u[pu+1]>=p;pu++{}forql<le&&u[ql]<q{ql++}forqu=ql-1;qu<le-1&&u[qu+1]>=q;qu++{}ifp!=q&&pl<=ql{pl=ql+1}ifpl>pu||ql>qu||ql>pu{returnfalse}ifout[le]==0{pu=le-1pl=pu}ps:=stackfor;pu>=pl;pu--{if!insert(p,pu){continue}out[pu]++sum[pu]+=leifp!=q{qs:=stackj:=quifj>=pu{j=pu-1}for;j>=ql;j--{if!insert(q,j){continue}out[j]++sum[j]+=letail[le]=qifseqRecur(le-1){returntrue}restore(qs)out[j]--sum[j]-=le}}else{out[pu]++sum[pu]+=letail[le]=pifseqRecur(le-1){returntrue}out[pu]--sum[pu]-=le}out[pu]--sum[pu]-=lerestore(ps)}returnfalse}funcseqRecur(leint)bool{n:=l[le]ifle<2{returntrue}limit:=n-1ifout[le]==1{limit=n-tail[sum[le]]}iflimit>u[le-1]{limit=u[le-1]}// Try to break n into p + q, and see if we can insert p, q into// list while satisfying bounds.p:=limitforq:=n-p;q<=p;q,p=q+1,p-1{iftry(p,q,le){returntrue}}returnfalse}funcseq(n,leint,buf[]int)int{ifle==0{le=seqLen(n)}stack=0l[le],u[le]=n,nfori:=0;i<=le;i++{out[i],sum[i]=0,0}fori:=2;i<le;i++{l[i]=l[i-1]+1u[i]=u[i-1]*2}fori:=le-1;i>2;i--{ifl[i]*2<l[i+1]{l[i]=(1+l[i+1])/2}ifu[i]>=u[i+1]{u[i]=u[i+1]-1}}if!seqRecur(le){return0}ifbuf!=nil{fori:=0;i<=le;i++{buf[i]=u[i]}}returnle}funcseqLen(nint)int{ifn<=known{returncache[n]}// Need all lower n to compute sequence.forknown+1<n{seqLen(known+1)}varubintlb:=lower(n,&ub)forlb<ub&&seq(n,lb,nil)==0{lb++}known=nifn&1023==0{fmt.Printf("Cached %d\n",known)}cache[n]=lbreturnlb}funcbinLen(nint)int{r,o:=-1,-1for;n!=0;n,r=n>>1,r+1{ifn&1!=0{o++}}returnr+o}type(vector=[]float64matrix[]vector)func(m1matrix)mul(m2matrix)matrix{rows1,cols1:=len(m1),len(m1[0])rows2,cols2:=len(m2),len(m2[0])ifcols1!=rows2{panic("Matrices cannot be multiplied.")}result:=make(matrix,rows1)fori:=0;i<rows1;i++{result[i]=make(vector,cols2)forj:=0;j<cols2;j++{fork:=0;k<rows2;k++{result[i][j]+=m1[i][k]*m2[k][j]}}}returnresult}func(mmatrix)pow(nint,printoutbool)matrix{e:=make([]int,N)varv[N]matrixle:=seq(n,0,e)ifprintout{fmt.Println("Addition chain:")fori:=0;i<=le;i++{c:=' 'ifi==le{c='\n'}fmt.Printf("%d%c",e[i],c)}}v[0]=mv[1]=m.mul(m)fori:=2;i<=le;i++{forj:=i-1;j!=0;j--{fork:=j;k>=0;k--{ife[k]+e[j]<e[i]{break}ife[k]+e[j]>e[i]{continue}v[i]=v[j].mul(v[k])j=1break}}}returnv[le]}func(mmatrix)print(){for_,v:=rangem{fmt.Printf("% f\n",v)}fmt.Println()}funcmain(){m:=27182n:=31415fmt.Println("Precompute chain lengths:")seqLen(n)rh:=math.Sqrt(0.5)mx:=matrix{{rh,0,rh,0,0,0},{0,rh,0,rh,0,0},{0,rh,0,-rh,0,0},{-rh,0,rh,0,0,0},{0,0,0,0,0,1},{0,0,0,0,1,0},}fmt.Println("\nThe first 100 terms of A003313 are:")fori:=1;i<=100;i++{fmt.Printf("%d ",seqLen(i))ifi%10==0{fmt.Println()}}exs:=[2]int{m,n}mxs:=[2]matrix{}fori,ex:=rangeexs{fmt.Println("\nExponent:",ex)mxs[i]=mx.pow(ex,true)fmt.Printf("A ^ %d:-\n\n",ex)mxs[i].print()fmt.Println("Number of A/C multiplies:",seqLen(ex))fmt.Println(" c.f. Binary multiplies:",binLen(ex))}fmt.Printf("\nExponent: %d x %d = %d\n",m,n,m*n)fmt.Printf("A ^ %d = (A ^ %d) ^ %d:-\n\n",m*n,m,n)mx2:=mxs[0].pow(n,false)mx2.print()}
Precompute chain lengths:Cached 1024Cached 2048Cached 3072....Cached 28672Cached 29696Cached 30720The first 100 terms of A003313 are: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 Exponent: 27182Addition chain:1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568 27136 27182A ^ 27182:-[-0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 -0.500000 0.000000 0.000000][ 0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]Number of A/C multiplies: 18 c.f. Binary multiplies: 21Exponent: 31415Addition chain:1 2 4 8 16 17 33 49 98 196 392 784 1568 3136 6272 6289 12561 25122 31411 31415A ^ 31415:-[ 0.707107 0.000000 0.000000 -0.707107 0.000000 0.000000][ 0.000000 0.707107 0.707107 0.000000 0.000000 0.000000][ 0.707107 0.000000 0.000000 0.707107 0.000000 0.000000][ 0.000000 0.707107 -0.707107 0.000000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000]Number of A/C multiplies: 19 c.f. Binary multiplies: 24Exponent: 27182 x 31415 = 853922530A ^ 853922530 = (A ^ 27182) ^ 31415:-[-0.500000 0.500000 -0.500000 0.500000 0.000000 0.000000][-0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]
Below is the original solution which was ruled inadmissible because it uses 'star chains' and is therefore non-optimal.
I think it should nevertheless be retained as it is an interesting approach and there are other solutions to this task which are based on it.
/*Continued fraction addition chains, as described in "Efficient computationof addition chains" by F. Bergeron, J. Berstel, and S. Brlek, published inJournal 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.*/packagemainimport("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.typestarChain[]int// ⊗= operator. modifies receiver.func(s1*starChain)cMul(s2starChain){p:=*s1i:=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(jint){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γsingletontypesingletonfunc(int)intfuncdichotomic(nint)int{returnn/(1<<uint((λ(n)+1)/2))}// integer log base 2funcλ(nint)(aint){forn!=1{a++n>>=1}return}// minChain as described in the paper.funcminChain(nint)starChain{switcha:=λ(n);{casen==1<<uint(a):r:=make(starChain,a+1)fori:=ranger{r[i]=1<<uint(i)}returnrcasen==3:returnstarChain{1,2,3}}returnchain(n,γ(n))}// chain as described in the paper.funcchain(n1,n2int)starChain{q,r:=n1/n2,n1%n2ifr==0{c:=minChain(n2)c.cMul(minChain(q))returnc}c:=chain(n2,r)c.cMul(minChain(q))c.cAdd(r)returnc}funcmain(){m:=31415n:=27182show(m)show(n)show(m*n)showEasier(m,1.00002206445416)showEasier(n,1.00002550055251)}funcshow(eint){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},})γ=dichotomicsc:=minChain(e)fmt.Println("addition chain:",sc)a.scExp(sc).print("a^e")fmt.Println("count of multiplies:",mCount)fmt.Println()}varmCountintfuncshowEasier(eint,afloat64){fmt.Println("exponent:",e)γ=dichotomicsc:=minChain(e)fmt.Printf("%.14f^%d: %.14f\n",a,sc[len(sc)-1],scExp64(a,sc))fmt.Println("count of multiplies:",mCount)fmt.Println()}funcscExp64(afloat64,scstarChain)float64{mCount=0p:=make([]float64,len(sc))p[0]=afori:=1;i<len(p);i++{d:=sc[i]-sc[i-1]j:=i-1forsc[j]!=d{j--}p[i]=p[i-1]*p[j]mCount++}returnp[len(p)-1]}func(m*matrix)scExp(scstarChain)*matrix{mCount=0p:=make([]*matrix,len(sc))p[0]=m.copy()fori:=1;i<len(p);i++{d:=sc[i]-sc[i-1]j:=i-1forsc[j]!=d{j--}p[i]=p[i-1].multiply(p[j])mCount++}returnp[len(p)-1]}func(m*matrix)copy()*matrix{return&matrix{append([]float64{},m.ele...),m.stride}}// code below copied from matrix multiplication tasktypematrixstruct{ele[]float64strideint}funcmatrixFromRows(rows[][]float64)*matrix{iflen(rows)==0{return&matrix{nil,0}}m:=&matrix{make([]float64,len(rows)*len(rows[0])),len(rows[0])}forrx,row:=rangerows{copy(m.ele[rx*m.stride:(rx+1)*m.stride],row)}returnm}func(m*matrix)print(headingstring){ifheading>""{fmt.Print(heading,"\n")}fore:=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}form1c0,m3x:=0,0;m1c0<len(m1.ele);m1c0+=m1.stride{form2r0:=0;m2r0<m2.stride;m2r0++{form1x,m2x:=m1c0,m2r0;m2x<len(m2.ele);m2x+=m2.stride{m3.ele[m3x]+=m1.ele[m1x]*m2.ele[m2x]m1x++}m3x++}}returnm3}
Output (manually wrapped at 80 columns.)
exponent: 31415addition chain: [1 2 4 5 10 20 25 50 55 110 220 245 490 980 1960 3920 784015680 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: 19exponent: 27182addition chain: [1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 1356827136 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: 18exponent: 853922530addition chain: [1 2 4 5 7 12 24 48 96 103 206 309 412 721 1133 1854 3708 48419682 19364 21218 26059 52118 104236 208472 416944 833888 1667776 3335552 667110413342208 26684416 53368832 106737664 213475328 426950656 853901312 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: 37exponent: 314151.00002206445416^31415: 1.99999999989447count of multiplies: 19exponent: 271821.00002550055251^27182: 1.99999999997876count of multiplies: 18
Generators of a nearly-optimal addition chains for a given number.
dichotomicChain::Integrala=>a->[a]dichotomicChainn|n==3=[3,2,1]|n==2^log2n=takeWhile(>0)$iterate(`div`2)n|otherwise=letk=n`div`(2^((log2n+1)`div`2))inchainnkwherechainn1n2|n2<=1=dichotomicChainn1|otherwise=casen1`divMod`n2of(q,0)->dichotomicChainq`mul`dichotomicChainn2(q,r)->[r]`add`(dichotomicChainq`mul`chainn2r)c1`mul`c2=map(headc2*)c1++tailc2c1`add`c2=map(headc2+)c1++c2log2=floor.logBase2.fromIntegralbinaryChain::Integrala=>a->[a]binaryChain1=[1]binaryChainn|evenn=n:binaryChain(n`div`2)|oddn=n:binaryChain(n-1)
λ> dichotomicChain 31415 [31415,31360,15680,7840,3920,1960,980,490,245,220,110,55,50,25,20,10,5,4,2,1]λ> length $ dichotomicChain 3141520λ> length $ binaryChain 3141525λ> length $ dichotomicChain (31415*27182)38λ> length $ binaryChain (31415*27182)45
Universal monoid multiplication that uses additive chain
times::(Monoidp,Integrala)=>a->p->p0`times`_=memptyn`times`x=reswhere(res:_,_,_)=foldlf([x],1,0)$tailchch=reverse$dichotomicChainnf(p:ps,c1,i)c2=letJustj=elemIndex(c2-c1)chin((p<>((p:ps)!!(i-j))):p:ps,c2,i+1)
λ> 31 `times` "a""aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa"λ> getSum $ 31 `times` 262λ> getProduct $ 31 `times` 22147483648
Implementation of matrix multiplication as monoidal operation.
dataMa=M[[a]]|IderivingShowinstanceNuma=>Semigroup(Ma)whereI<>m=mm<>I=mMm1<>Mm2=M$map(\r->map(\c->r`dot`c)(transposem2))m1wheredotab=sum$zipWith(*)abinstanceNuma=>Monoid(Ma)wheremempty=I
-- matrix multiplicationλ> M [[2,4],[5,7]] <> M [[4,1],[6,2]]M [[32,10],[62,19]]-- matrix exponentiationλ> 2 `times` M [[2,4],[5,7]]M [[24,36],[45,69]]-- calculation of large Fibonacci numbersλ> 2 `times` M [[1,1],[1,0]]M [[2,1],[1,1]]λ> 3 `times` M [[1,1],[1,0]]M [[3,2],[2,1]]λ> 4 `times` M [[1,1],[1,0]]M [[5,3],[3,2]]λ> 20 `times` M [[1,1],[1,0]]M [[10946,6765],[6765,4181]]λ> 200 `times` M [[1,1],[1,0]]M [[453973694165307953197296969697410619233826,280571172992510140037611932413038677189525],[280571172992510140037611932413038677189525,173402521172797813159685037284371942044301]]
The task
λ> :{Main:> | let a = a = let s = sqrt (1/2) in M [[s,0,s,0,0,0]Main:> | ,[0,s,0,s,0,0]Main:> | ,[0,s,0,-s,0,0]Main:> | ,[-s,0,s,0,0,0]Main:> | ,[0,0,0,0,0,1]Main:> | ,[0,0,0,0,1,0]] Main:> | :}λ> 31415 `times` aM [[0.7071067811883202,0.0,0.0,-0.7071067811883202,0.0,0.0],[0.0,0.7071067811883202,0.7071067811883202,0.0,0.0,0.0],[0.7071067811883202,0.0,0.0,0.7071067811883202,0.0,0.0],[0.0,0.7071067811883202,-0.7071067811883202,0.0,0.0,0.0],[0.0,0.0,0.0,0.0,0.0,1.0],[0.0,0.0,0.0,0.0,1.0,0.0]]λ> 27182 `times` aM [[-0.5000000000010233,-0.5000000000010233,-0.5000000000010233,0.5000000000010233,0.0,0.0],[0.5000000000010233,-0.5000000000010233,-0.5000000000010233,-0.5000000000010233,0.0,0.0],[-0.5000000000010233,-0.5000000000010233,0.5000000000010233,-0.5000000000010233,0.0,0.0],[0.5000000000010233,-0.5000000000010233,0.5000000000010233,0.5000000000010233,0.0,0.0],[0.0,0.0,0.0,0.0,1.0,0.0],[0.0,0.0,0.0,0.0,0.0,1.0]]λ> (31415*27182) `times` aM [[-0.500000030038797,0.5000000300387971,-0.5000000300387971,0.5000000300387973,0.0,0.0],[-0.5000000300387971,-0.5000000300387972,-0.5000000300387969,-0.5000000300387969,0.0,0.0],[-0.5000000300387972,-0.5000000300387971,0.5000000300387969,0.5000000300387969,0.0,0.0],[0.5000000300387971,-0.500000030038797,-0.5000000300387973,0.5000000300387971,0.0,0.0],[0.0,0.0,0.0,0.0,1.0,0.0],[0.0,0.0,0.0,0.0,0.0,1.0]]importjava.util.Arrays;publicclassMain{// ==== toggle this if you want exact minimal chains (much slower) ====staticfinalbooleanFAST=true;staticfinalintN=32;staticfinalintNMAX=40000;staticint[]u=newint[N],l=newint[N],out=newint[N],sum=newint[N],tail=newint[N];staticint[]cache=newint[NMAX+1];staticintknown=2,stack=0;staticSave[]undo=newSave[N*N];static{u[0]=1;u[1]=2;l[0]=1;l[1]=2;cache[2]=1;for(inti=0;i<undo.length;i++)undo[i]=newSave();}staticfinalclassIntBox{intvalue;IntBox(){}IntBox(intv){value=v;}}staticfinalclassSave{int[]arr;intidx;intval;}staticvoidreplace(int[]arr,inti,intn){undo[stack].arr=arr;undo[stack].idx=i;undo[stack].val=arr[i];arr[i]=n;stack++;}staticvoidrestore(intn){while(stack>n){Saves=undo[--stack];s.arr[s.idx]=s.val;}}staticintlower(intn,IntBoxup){if(n<=2||(n<=NMAX&&cache[n]!=0)){if(up!=null)up.value=cache[n];returncache[n];}inti=-1,o=0,tmp=n;for(;tmp!=0;tmp>>=1,i++)if((tmp&1)!=0)o++;if(up!=null){i--;up.value=o+i;}while(true){i++;o>>=1;if(o==0)break;}if(up==null)returni;for(intd=2;d*d<n;d++){if(n%d!=0)continue;intq=cache[d]+cache[n/d];if(q<up.value){up.value=q;if(q==i)break;}}if(n>2){if(up.value>cache[n-1]+1)up.value=cache[n-1]+1;if(up.value>cache[n-2]+1)up.value=cache[n-2]+1;}returni;}staticbooleaninsert(intx,intpos){intsave=stack;if(l[pos]>x||u[pos]<x)returnfalse;if(l[pos]!=x){replace(l,pos,x);for(inti=pos-1;i>=0&&u[i]*2<u[i+1];i--){intt=l[i+1]+1;if(t*2>u[i]){restore(save);returnfalse;}replace(l,i,t);}for(inti=pos+1;i<N&&l[i]<=l[i-1];i++){intt=l[i-1]+1;if(t>u[i]){restore(save);returnfalse;}replace(l,i,t);}}if(u[pos]==x)returntrue;replace(u,pos,x);for(inti=pos-1;i>=0&&u[i]>=u[i+1];i--){intt=u[i+1]-1;if(t<l[i]){restore(save);returnfalse;}replace(u,i,t);}for(inti=pos+1;i<N&&u[i]>u[i-1]*2;i++){intt=u[i-1]*2;if(t<l[i]){restore(save);returnfalse;}replace(u,i,t);}returntrue;}staticbooleantryPQ(intp,intq,intle){intpl=cache[p];if(pl>=le)returnfalse;intql=cache[q];if(ql>=le)returnfalse;intpu,qu;while(pl<le&&u[pl]<p)pl++;for(pu=pl-1;pu<le-1&&u[pu+1]>=p;pu++){}while(ql<le&&u[ql]<q)ql++;for(qu=ql-1;qu<le-1&&u[qu+1]>=q;qu++){}if(p!=q&&pl<=ql)pl=ql+1;if(pl>pu||ql>qu||ql>pu)returnfalse;if(out[le]==0){pu=le-1;pl=pu;}intps=stack;for(;pu>=pl;pu--){if(!insert(p,pu))continue;out[pu]++;sum[pu]+=le;if(p!=q){intqs=stack;intj=qu;if(j>=pu)j=pu-1;for(;j>=ql;j--){if(!insert(q,j))continue;out[j]++;sum[j]+=le;tail[le]=q;if(seqRecur(le-1))returntrue;restore(qs);out[j]--;sum[j]-=le;}}else{out[pu]++;sum[pu]+=le;tail[le]=p;if(seqRecur(le-1))returntrue;out[pu]--;sum[pu]-=le;}out[pu]--;sum[pu]-=le;restore(ps);}returnfalse;}staticbooleanseqRecur(intle){intn=l[le];if(le<2)returntrue;intlimit=n-1;if(out[le]==1)limit=n-tail[sum[le]];if(limit>u[le-1])limit=u[le-1];intp=limit;for(intq=n-p;q<=p;q++,p--)if(tryPQ(p,q,le))returntrue;returnfalse;}staticintseq(intn,intle,int[]buf){if(le==0)le=seqLen(n);stack=0;l[le]=n;u[le]=n;for(inti=0;i<=le;i++){out[i]=0;sum[i]=0;}for(inti=2;i<le;i++){l[i]=l[i-1]+1;u[i]=u[i-1]*2;}for(inti=le-1;i>2;i--){if(l[i]*2<l[i+1])l[i]=(1+l[i+1])/2;if(u[i]>=u[i+1])u[i]=u[i+1]-1;}if(!seqRecur(le))return0;if(buf!=null)for(inti=0;i<=le;i++)buf[i]=u[i];returnle;}staticintseqLen(intn){if(FAST){// fast fallback: binary upper bound onlyreturnbinLen(n);}if(n<=known)returncache[n];while(known+1<n)seqLen(known+1);IntBoxub=newIntBox();intlb=lower(n,ub);while(lb<ub.value&&seq(n,lb,null)==0)lb++;known=n;if((n&1023)==0)System.out.printf("Cached %d%n",known);cache[n]=lb;returnlb;}staticintbinLen(intn){intr=-1,o=-1;for(;n!=0;n>>=1,r++)if((n&1)!=0)o++;returnr+o;}staticfinalclassMatrix{finaldouble[][]a;Matrix(double[][]d){a=d;}Matrix(intr,intc){a=newdouble[r][c];}Matrixmul(Matrixm2){intr=a.length,c=a[0].length,r2=m2.a.length,c2=m2.a[0].length;if(c!=r2)thrownewIllegalArgumentException("Matrices cannot be multiplied.");Matrixres=newMatrix(r,c2);for(inti=0;i<r;i++)for(intj=0;j<c2;j++){doubles=0;for(intk=0;k<r2;k++)s+=a[i][k]*m2.a[k][j];res.a[i][j]=s;}returnres;}// exact chain powering (slow) – used only when FAST=falseMatrixpowByChain(intn,booleanprintout){int[]e=newint[N];Matrix[]v=newMatrix[N];intle=seq(n,0,e);if(printout){System.out.println("Addition chain:");for(inti=0;i<=le;i++){charc=(i==le)?'\n':' ';System.out.printf("%d%c",e[i],c);}}v[0]=this;v[1]=this.mul(this);for(inti=2;i<=le;i++){for(intj=i-1;j!=0;j--){for(intk=j;k>=0;k--){if(e[k]+e[j]<e[i])break;if(e[k]+e[j]>e[i])continue;v[i]=v[j].mul(v[k]);j=1;break;}}}returnv[le];}// fast square-and-multiplyMatrixpowBinary(intexp,booleanprintout){if(printout){System.out.println("Binary chain (exponents where we multiply):");booleanfirst=true;for(inti=31-Integer.numberOfLeadingZeros(exp);i>=0;i--){if(((exp>>i)&1)==1){if(!first)System.out.print(' ');System.out.print(1<<i);first=false;}}System.out.println();}Matrixbase=this;Matrixresult=identity(a.length);inte=exp;while(e>0){if((e&1)==1)result=result.mul(base);e>>=1;if(e!=0)base=base.mul(base);}returnresult;}Matrixpow(intn,booleanprintout){returnFAST?powBinary(n,printout):powByChain(n,printout);}staticMatrixidentity(intn){MatrixI=newMatrix(n,n);for(inti=0;i<n;i++)I.a[i][i]=1.0;returnI;}voidprint(){for(double[]row:a){StringBuildersb=newStringBuilder();for(inti=0;i<row.length;i++){if(i>0)sb.append(' ');sb.append(String.format("% f",row[i]));}System.out.println(sb);}System.out.println();}}publicstaticvoidmain(String[]args){intm=27182,n=31415;System.out.println("Precompute chain lengths:");// In FAST mode, seqLen uses binLen and returns immediately.seqLen(n);doublerh=Math.sqrt(0.5);Matrixmx=newMatrix(newdouble[][]{{rh,0,rh,0,0,0},{0,rh,0,rh,0,0},{0,rh,0,-rh,0,0},{-rh,0,rh,0,0,0},{0,0,0,0,0,1},{0,0,0,0,1,0},});System.out.println("\nThe first 100 terms of A003313 are:");for(inti=1;i<=100;i++){System.out.printf("%d ",seqLen(i));if(i%10==0)System.out.println();}int[]exs={m,n};Matrix[]mxs=newMatrix[2];for(inti=0;i<exs.length;i++){intex=exs[i];System.out.printf("%nExponent: %d%n",ex);mxs[i]=mx.pow(ex,true);System.out.printf("A ^ %d:-%n%n",ex);mxs[i].print();System.out.println("Number of A/C multiplies: "+(FAST?binLen(ex):seqLen(ex)));System.out.println(" c.f. Binary multiplies: "+binLen(ex));}System.out.printf("%nExponent: %d x %d = %d%n",m,n,m*n);System.out.printf("A ^ %d = (A ^ %d) ^ %d:-%n%n",m*n,m,n);Matrixmx2=(FAST?mxs[0].powBinary(n,false):mxs[0].powByChain(n,false));mx2.print();}}
Precompute chain lengths:The first 100 terms of A003313 are:0 1 2 2 3 3 4 3 4 4 5 4 5 5 6 4 5 5 6 5 6 6 7 5 6 6 7 6 7 7 8 5 6 6 7 6 7 7 8 6 7 7 8 7 8 8 9 6 7 7 8 7 8 8 9 7 8 8 9 8 9 9 10 6 7 7 8 7 8 8 9 7 8 8 9 8 9 9 10 7 8 8 9 8 9 9 10 8 9 9 10 9 10 10 11 7 8 8 9 8 Exponent: 27182Binary chain (exponents where we multiply):16384 8192 2048 512 32 8 4 2A ^ 27182:--0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000 0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000-0.500000 -0.500000 0.500000 -0.500000 0.000000 0.000000 0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000Number of A/C multiplies: 21 c.f. Binary multiplies: 21Exponent: 31415Binary chain (exponents where we multiply):16384 8192 4096 2048 512 128 32 16 4 2 1A ^ 31415:- 0.707107 0.000000 0.000000 -0.707107 0.000000 0.000000 0.000000 0.707107 0.707107 0.000000 0.000000 0.000000 0.707107 0.000000 0.000000 0.707107 0.000000 0.000000 0.000000 0.707107 -0.707107 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000Number of A/C multiplies: 24 c.f. Binary multiplies: 24Exponent: 27182 x 31415 = 853922530A ^ 853922530 = (A ^ 27182) ^ 31415:--0.500000 0.500000 -0.500000 0.500000 0.000000 0.000000-0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000-0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000 0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000
Uses an iterator to generate A003313 (the first 100 values). Uses the Knuth path method for the larger integers. This gives a chain length of 18 for both 27182 and 31415.
importBase.iteratemutablestructAdditionChains{T}chains::Vector{Vector{T}}work_chain::Intwork_element::IntAdditionChains{T}()whereT=new{T}([[one(T)]],1,1)endfunctionBase.iterate(acs::AdditionChains,state=1)i,j=acs.work_chain,acs.work_elementnewchain=[acs.chains[i];acs.chains[i][end]+acs.chains[i][j]]push!(acs.chains,newchain)ifj==length(acs.chains[i])acs.work_chain+=1acs.work_element=1elseacs.work_element+=1endreturnnewchain,state+1endfunctionfindchain!(acs::AdditionChains,n)@assertn>0n==1&&return[one(eltype(first(acs.chains)))]idx=findfirst(a->a[end]==n,acs.chains)ifidx==nothingfor(i,chain)inenumerate(acs)chain[end]==n&&returnchainendendreturnacs.chains[idx]end""" memoization for knuth_path """constknuth_path_p,knuth_path_lvl=Dict(1=>0),[[1]]""" knuth path method for addition chains """functionknuth_path(n)::Vector{Int}iszero(n)&&returnInt[]while!haskey(knuth_path_p,n)q=Int[]forxinfirst(knuth_path_lvl),yinknuth_path(x)if!haskey(knuth_path_p,x+y)knuth_path_p[x+y]=xpush!(q,x+y)endendknuth_path_lvl[begin]=qendreturnpush!(knuth_path(knuth_path_p[n]),n)endfunctionpow(x,chain)p,products=0,Dict{Int,typeof(x)}(0=>one(x),1=>x)foriinchainproducts[i]=products[p]*products[i-p]p=iendreturnproducts[chain[end]]endfunctiontest_addition_chains()additionchain=AdditionChains{Int}()println("First one hundred addition chain lengths:")foriin1:100print(rpad(length(findchain!(additionchain,i))-1,3),i%10==0?"\n":"")endprintln("\nKnuth chains for addition chains of 31415 and 27182:")expchains=Dict(i=>knuth_path(i)foriin[31415,27182])for(n,chn)inexpchainsprintln("Exponent: ",rpad(n,10),"\n Addition Chain:$(chn[begin:end-1]))")endprintln("\n1.00002206445416^31415 = ",pow(1.00002206445416,expchains[31415]))println("1.00002550055251^27182 = ",pow(1.00002550055251,expchains[27182]))println("1.00002550055251^(27182 * 31415) = ",pow(BigFloat(pow(1.00002550055251,expchains[27182])),expchains[31415]))println("(1.000025 + 0.000058i)^27182 = ",pow(Complex(1.000025,0.000058),expchains[27182]))println("(1.000022 + 0.000050i)^31415 = ",pow(Complex(1.000022,0.000050),expchains[31415]))x=sqrt(1/2)matrixA=[x0x000;0x0x00;0x0-x00;-x0x000;000001;000010]println("matrix A ^ 27182 = ")display(pow(matrixA,expchains[27182]))println("matrix A ^ 31415 = ")display(round.(pow(matrixA,expchains[31415]),digits=6))println("(matrix A ^ 27182) ^ 31415 = ")display(pow(pow(matrixA,expchains[27182]),expchains[31415]))endtest_addition_chains()
First one hundred addition chain lengths:0 1 2 2 3 3 4 3 4 45 4 5 5 5 4 5 5 6 56 6 6 5 6 6 6 6 7 67 5 6 6 7 6 7 7 7 67 7 7 7 7 7 8 6 7 77 7 8 7 8 7 8 8 8 78 8 8 6 7 7 8 7 8 89 7 8 8 8 8 8 8 9 78 8 8 8 8 8 9 8 9 89 8 9 9 9 7 8 8 8 8Knuth chains for addition chains of 31415 and 27182:Exponent: 27182 Addition Chain: [1, 2, 3, 5, 7, 14, 21, 35, 70, 140, 143, 283, 566, 849, 1698, 3396, 6792, 6799, 13591])Exponent: 31415 Addition Chain: [1, 2, 3, 5, 7, 14, 28, 56, 61, 122, 244, 488, 976, 1952, 3904, 7808, 15616, 15677, 31293])1.00002206445416^31415 = 1.99999999989246381.00002550055251^27182 = 1.99999999997926881.00002550055251^(27182 * 31415) = 7.199687435551025768365237017391630520475805934721292725697031724530209692195819e+9456(1.000025 + 0.000058i)^27182 = -0.01128636963542673 + 1.9730308496660347im(1.000022 + 0.000050i)^31415 = 0.00016144681325535107 + 1.9960329014194498immatrix A ^ 27182 =6×6 Matrix{Float64}: -0.5 -0.5 -0.5 0.5 0.0 0.0 0.5 -0.5 -0.5 -0.5 0.0 0.0 -0.5 -0.5 0.5 -0.5 0.0 0.0 0.5 -0.5 0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0matrix A ^ 31415 =6×6 Matrix{Float64}: 0.707107 0.0 -0.0 -0.707107 0.0 0.0 -0.0 0.707107 0.707107 -0.0 0.0 0.0 0.707107 -0.0 -0.0 0.707107 0.0 0.0 0.0 0.707107 -0.707107 -0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0(matrix A ^ 27182) ^ 31415 =6×6 Matrix{Float64}: -0.5 0.5 -0.5 0.5 0.0 0.0 -0.5 -0.5 -0.5 -0.5 0.0 0.0 -0.5 -0.5 0.5 0.5 0.0 0.0 0.5 -0.5 -0.5 0.5 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0importmath,sequtils,strutilsconstN=32NMax=40_000typeSave=tuple[p:ptrint;v:int]varu:array[N,int]# Upper bounds.l:array[N,int]# Lower bounds.u[0]=1;u[1]=2l[0]=1;l[1]=2varoutp,sum,tail:array[N,int]varcache:array[NMax+1,int]cache[2]=1varknown=2stack=0undo:array[N*N,Save]procreplace(x:varopenArray[int];i,n:int)=undo[stack]=(x[i].addr,x[i])x[i]=nincstackprocrestore(n:int)=whilestack>n:decstackundo[stack].p[]=undo[stack].vprocbounds(n:int):tuple[lower,upper:int]=# Return lower and upper bounds.ifn<=2orn<=NMaxandcache[n]!=0:return(cache[n],cache[n])vari=-1o=0n=nwhilen!=0:if(nand1)!=0:incon=nshr1incideciresult.upper=o+iwhiletrue:incio=oshr1ifo==0:breako=2whileo*o<n:ifnmodo==0:letq=cache[o]+cache[ndivo]ifq<result.upper:result.upper=qifq==i:breakincoifn>2:ifresult.upper>cache[n-2]+1:result.upper=cache[n-1]+1ifresult.upper>cache[n-2]+1:result.upper=cache[n-2]+1result.lower=iprocinsert(x,pos:int):bool=letsave=stackifl[pos]>xoru[pos]<x:returnfalseifl[pos]!=x:l.replace(pos,x)vari=pos-1whileu[i]*2<u[i+1]:lett=l[i+1]+1ift*2>u[i]:restore(save)returnfalsel.replace(i,t)decii=pos+1whilel[i]<=l[i-1]:lett=l[i-1]+1ift>u[i]:restore(save)returnfalsel.replace(i,t)inciifu[pos]==x:returntrueu.replace(pos,x)vari=pos-1whileu[i]>=u[i+1]:lett=u[i+1]-1ift<l[i]:restore(save)returnfalseu.replace(i,t)decii=pos+1whileu[i]>u[i-1]*2:lett=u[i-1]*2ift<l[i]:restore(save)returnfalseu.replace(i,t)inciresult=true# Forward reference.procseqRecur(le:int):boolproc`try`(p,q,le:int):bool=varpl=cache[p]ifpl>=le:returnfalsevarql=cache[q]ifql>=le:returnfalsewhilepl<leandu[pl]<p:incplvarpu=pl-1whilepu<le-1andu[pu+1]>=p:incpuwhileql<leandu[ql]<q:incqlvarqu=ql-1whilequ<le-1andu[qu+1]>=q:incquifp!=qandpl<=ql:pl=ql+1ifpl>puorql>quorql>pu:returnfalseifoutp[le]==0:pu=le-1pl=puletps=stackwhilepu>=pl:ifinsert(p,pu):incoutp[pu]incsum[pu],leifp!=q:letqs=stackvarj=quifj>=pu:j=pu-1whilej>=ql:ifinsert(q,j):incoutp[j]incsum[j],letail[le]=qifseqRecur(le-1):returntruerestore(qs)decoutp[j]decsum[j],ledecjelse:incoutp[pu]incsum[pu],letail[le]=pifseqRecur(le-1):returntruedecoutp[pu]decsum[pu],ledecoutp[pu]decsum[pu],lerestore(ps)decpuprocseqRecur(le:int):bool=letn=l[le]ifle<2:returntruevarlimit=n-1ifoutp[le]==1:limit=n-tail[sum[le]]iflimit>u[le-1]:limit=u[le-1]# Try to break n into p + q, and see if we can insert p, q into# list while satisfying bounds.varp=limitvarq=n-pwhileq<=p:if`try`(p,q,le):returntruedecpincq# Forward reference.procsequence(n,le:int):intprocseqLen(n:int):int=ifn<=known:returncache[n]# Need all lower n to compute sequence.whileknown+1<n:discardseqLen(known+1)var(lb,ub)=bounds(n)whilelb<ubandsequence(n,lb)==0:inclbknown=nif(nand1023)==0:echo"Cached: ",knowncache[n]=lbresult=lbprocsequence(n,le:int):int=letle=ifle!=0:leelse:seqLen(n)stack=0l[le]=nu[le]=nforiin0..le:outp[i]=0sum[i]=0foriin2..<le:l[i]=l[i-1]+1u[i]=u[i-1]*2foriincountdown(le-1,3):ifl[i]*2<l[i+1]:l[i]=(1+l[i+1])div2ifu[i]>=u[i+1]:u[i]=u[i+1]-1ifnotseqRecur(le):return0result=leprocsequence(n,le:int;buf:varopenArray[int]):int=letle=sequence(n,le)foriin0..le:buf[i]=u[i]result=leprocbinLen(n:int):int=varr,o=-1varn=nwhilen!=0:if(nand1)!=0:incon=nshr1incrresult=r+otypeVector=seq[float]Matrix=seq[Vector]func`*`(m1,m2:Matrix):Matrix=letrows1=m1.lencols1=m1[0].lenrows2=m2.lencols2=m2[0].lenifcols1!=rows2:raisenewException(ValueError,"Matrices cannot be multiplied.")result=newSeqWith(rows1,newSeq[float](cols2))foriin0..<rows1:forjin0..<cols2:forkin0..<rows2:result[i][j]+=m1[i][k]*m2[k][j]procpow(m:Matrix;n:int;printout:bool):Matrix=vare=newSeq[int](N)varv:array[N,Matrix]letle=sequence(n,0,e)ifprintout:echo"Addition chain:"echoe[0..le].join(" ")v[0]=mv[1]=m*mforiin2..le:blockloop2:forjincountdown(i-1,0):forkincountdown(j,0):ife[k]+e[j]<e[i]:breakife[k]+e[j]>e[i]:continuev[i]=v[j]*v[k]breakloop2result=v[le]func`$`(m:Matrix):string=forvinm:result.add'['letstart=result.lenforxinv:result.addSep(" ",start)result.addx.formatFloat(ffDecimal,precision=6).align(9)result.add"]\n"varm=27182varn=31415echo"Precompute chain lengths:"discardseqlen(n)letrh=sqrt(0.5)letmx:Matrix=@[@[rh,0.0,rh,0.0,0.0,0.0],@[0.0,rh,0.0,rh,0.0,0.0],@[0.0,rh,0.0,-rh,0.0,0.0],@[-rh,0.0,rh,0.0,0.0,0.0],@[0.0,0.0,0.0,0.0,0.0,1.0],@[0.0,0.0,0.0,0.0,1.0,0.0]]echo"\nThe first 100 terms of A003313 are:"foriin1..100:stdout.writeseqLen(i),' 'ifimod10==0:echo()letexs=[m,n]varmxs:array[2,Matrix]fori,exinexs:echo()echo"Exponent: ",exmxs[i]=mx.pow(ex,true)echo"A ^$#:\n".format(ex)echomxs[i]echo"Number of A/C multiplies: ",seqLen(ex)echo" c.f. Binary multiplies: ",binLen(ex)echo()echo"Exponent:$1 x$2 =$3".format(m,n,m*n)echo"A ^$1 = (A ^$2) ^$3:\n".format(m*n,m,n)echomxs[0].pow(n,false)
Precompute chain lengths:Cached: 1024Cached: 2048Cached: 3072Cached: 4096Cached: 5120Cached: 6144Cached: 7168Cached: 8192Cached: 9216Cached: 10240Cached: 11264Cached: 12288Cached: 13312Cached: 14336Cached: 15360Cached: 16384Cached: 17408Cached: 18432Cached: 19456Cached: 20480Cached: 21504Cached: 22528Cached: 23552Cached: 24576Cached: 25600Cached: 26624Cached: 27648Cached: 28672Cached: 29696Cached: 30720The first 100 terms of A003313 are: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 Exponent: 27182Addition chain:1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568 27136 27182A ^ 27182:[-0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 -0.500000 0.000000 0.000000][ 0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]Number of A/C multiplies: 18 c.f. Binary multiplies: 21Exponent: 31415Addition chain:1 2 4 8 16 17 33 49 98 196 392 784 1568 3136 6272 6289 12561 25122 31411 31415A ^ 31415:[ 0.707107 0.000000 0.000000 -0.707107 0.000000 0.000000][ 0.000000 0.707107 0.707107 0.000000 0.000000 0.000000][ 0.707107 0.000000 0.000000 0.707107 0.000000 0.000000][ 0.000000 0.707107 -0.707107 0.000000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000]Number of A/C multiplies: 19 c.f. Binary multiplies: 24Exponent: 27182 x 31415 = 853922530A ^ 853922530 = (A ^ 27182) ^ 31415:[-0.500000 0.500000 -0.500000 0.500000 0.000000 0.000000][-0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]
Naieve brute force search, no attempt to optimise, manages about 4 million checks/s.
Replacing the recursion with an internal stack and chosen with a fixed length array might help, but otherwise I got no good ideas at all for trimming the search space.
Giving it the same length, I think, yields the same result asKnuth's_power_tree#Phix, and at leastin the cases that I tried, somewhat faster than the method on that page.
If you know the A003313 number, you can throw that at it and wait (for several billion years) or get the power tree length and loop trying to find a path one shorter (and wait several trillion years). The path() routine is a direct copy from the link above.
Note that "tries" overflows (crashes) at 1073741824, which I kept in as a deliberate limiter.
withjavascript_semanticsconstantp=new_dict({{1,0}})sequencelvl={1}functionpath(integern)ifn=0thenreturn{}endifwhilegetd_index(n,p)=NULLdosequenceq={}fori=1tolength(lvl)dointegerx=lvl[i]sequencepx=path(x)forj=1tolength(px)dointegery=x+px[j]ifgetd_index(y,p)!=NULLthenexitendifsetd(y,x,p)q&=yendforendforlvl=qendwhilereturnpath(getd(n,p))&nendfunctionatomt1=time()+1integertries=0functionaddition_chain(integertarget,len,sequencechosen={1})-- target and len must be >=2tries+=1integerl=length(chosen),last=chosen[$]iflast=targetthenreturnchosenendififl=lenthenifplatform()!=JSandtime()>t1then?{"addition_chain",chosen,tries}t1=time()+1endifelsefori=lto1by-1dointegernext=last+chosen[i]ifnext<=targetthensequenceres=addition_chain(target,len,deep_copy(chosen)&next)iflength(res)thenreturnresendifendifendforendifreturn{}endfunction-- first, some proof of correctness at the lower/trivial end:sequenceres=repeat(0,120)res[2]=1forn=3tolength(res)dointegerl=length(path(n))whiletruedosequenceac=addition_chain(n,l)iflength(ac)=0thenexitendifl=length(ac)-1endwhileres[n]=lendforputs(1,"The first 120 members of A003313:\n")pp(res)printf(1,"addition_chain(31,8):%s\n",{sprint(addition_chain(31,8))})
The first 120 members of A003313:{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,9,8,9,8,9,9,9,8,9,9,9, 8,9,9,9,9,9,9,9,8}addition_chain(31,8):{1,2,4,8,10,20,30,31}On the task numbers, however, as to be expected, it struggles but probably would eventually get there if the overflow on tries were removed:
--?path(12509) -- {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,12288,12416,12480,12496,12504,12508,12509}--?addition_chain(12509,21) -- 0s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,12288,12416,12480,12496,12504,12508,12509}--?addition_chain(12509,20) -- 12.3s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,4160,8320,12480,12496,12504,12508,12509}--?addition_chain(12509,19) -- 1.1s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,4160,4168,8336,12504,12508,12509}--?addition_chain(12509,18) -- bust--?addition_chain(12509,17) -- bust--?path(31415) -- {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,28672,30720,31232,31360,31392,31408,31412,31414,31415}--?addition_chain(31415,25) -- 0s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,28672,30720,31232,31360,31392,31408,31412,31414,31415}--?addition_chain(31415,24) -- bust--?addition_chain(31415,23) -- bust--?addition_chain(31415,22) -- 137s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,10240,10368,10384,10386,20772,31158,31414,31415}--?addition_chain(31415,21) -- 116s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,6144,6272,6288,6290,12562,25124,31414,31415}--?addition_chain(31415,20) -- bust --?addition_chain(31415,19) -- bust--?path(27182) -- {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,26624,27136,27168,27176,27180,27182}--?addition_chain(27182,22) -- 0s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,26624,27136,27168,27176,27180,27182}--?addition_chain(27182,21) -- 19.4s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,8704,8712,8714,17428,26142,27166,27182}--?addition_chain(27182,20) -- 14.2s {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,5120,5128,5640,5642,10770,21540,27182}--?addition_chain(27182,19) -- bust--?addition_chain(27182,18) -- bustOnce you have an addition chain, of course, apply it as perKnuth's_power_tree#Phix, and of course length(pn) is the number of multiplies that will perform.
--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,28672,30720,31232,31360,31392,31408,31412,31414,31415}--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,10240,10368,10384,10386,20772,31158,31414,31415}--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,6144,6272,6288,6290,12562,25124,31414,31415}--sequence pn = {1,2,4,8,16,17,33,49,98,196,392,784,1568,3136,6272,6289,12561,25122,31411,31415} -- (from C)--printf(1,"%3g ^ %d (%d) = %s\n", {1.00002206445416,31415,length(pn),treepow(1.00002206445416,31415,pn)})--1.00002 ^ 31415 (25) = 1.9999999998949602994638558--1.00002 ^ 31415 (22) = 1.9999999998949602994638556--1.00002 ^ 31415 (21) = 1.9999999998949602994638552--1.00002 ^ 31415 (20) = 1.9999999998949602994638291--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,24576,26624,27136,27168,27176,27180,27182}--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,8192,8704,8712,8714,17428,26142,27166,27182}--sequence pn = {1,2,4,8,16,32,64,128,256,512,1024,2048,4096,5120,5128,5640,5642,10770,21540,27182}--sequence pn = {1,2,4,8,10,18,28,46,92,184,212,424,848,1696,3392,6784,13568,27136,27182} -- (from C)--printf(1,"%3g ^ %d (%d) = %s\n", {1.00002550055251,27182,length(pn),treepow(1.00002550055251,27182,pn)})--1.00003 ^ 27182 (22) = 1.9999999999727968238298861--1.00003 ^ 27182 (21) = 1.9999999999727968238298859--1.00003 ^ 27182 (20) = 1.9999999999727968238298855--1.00003 ^ 27182 (19) = 1.9999999999727968238298527''' Rosetta Code task Addition-chain_exponentiation '''frommathimportsqrtfromnumpyimportarrayfrommpmathimportmpfclassAdditionChains:''' two methods of calculating addition chains '''def__init__(self):''' memoization for knuth_path '''self.chains,self.idx,self.pos=[[1]],0,0self.pat,self.lvl={1:0},[[1]]defadd_chain(self):''' method 1: add chains depth then breadth first until done '''newchain=self.chains[self.idx].copy()newchain.append(self.chains[self.idx][-1]+self.chains[self.idx][self.pos])self.chains.append(newchain)ifself.pos==len(self.chains[self.idx])-1:self.idx+=1self.pos=0else:self.pos+=1returnnewchaindeffind_chain(self,nexp):''' method 1 interface: search for chain ending with n, adding more as needed '''assertnexp>0ifnexp==1:return[1]chn=next((aforainself.chainsifa[-1]==nexp),None)ifchnisNone:whileTrue:chn=self.add_chain()ifchn[-1]==nexp:breakreturnchndefknuth_path(self,ngoal):''' method 2: knuth method, uses memoization to search for a shorter chain '''ifngoal<1:return[]whilenotngoalinself.pat:new_lvl=[]foriinself.lvl[0]:forjinself.knuth_path(i):ifnoti+jinself.pat:self.pat[i+j]=inew_lvl.append(i+j)self.lvl[0]=new_lvlreturnpath=self.knuth_path(self.pat[ngoal])returnpath.append(ngoal)returnreturnpathdefcpow(xbase,chain):''' raise xbase by an addition exponentiation chain for what becomes x**chain[-1] '''pows,products=0,{0:1,1:xbase}foriinchain:products[i]=products[pows]*products[i-pows]pows=ireturnproducts[chain[-1]]if__name__=='__main__':# test both addition chain methodsacs=AdditionChains()print('First one hundred addition chain lengths:')forkinrange(1,101):print(f'{len(acs.find_chain(k))-1:3}',end='\n'ifk%10==0else'')print('\nKnuth chains for addition chains of 31415 and 27182:')chns={m:acs.knuth_path(m)formin[31415,27182]}for(num,cha)inchns.items():print(f'Exponent:{num:10}\n Addition Chain:{cha[:-1]}')print('\n1.00002206445416^31415 =',cpow(1.00002206445416,chns[31415]))print('1.00002550055251^27182 =',cpow(1.00002550055251,chns[27182]))print('1.000025 + 0.000058i)^27182 =',cpow(complex(1.000025,0.000058),chns[27182]))print('1.000022 + 0.000050i)^31415 =',cpow(complex(1.000022,0.000050),chns[31415]))sq05=mpf(sqrt(0.5))mat=array([[sq05,0,sq05,0,0,0],[0,sq05,0,sq05,0,0],[0,sq05,0,-sq05,0,0],[-sq05,0,sq05,0,0,0],[0,0,0,0,0,1],[0,0,0,0,1,0]])print('matrix A ^ 27182 =')print(cpow(mat,chns[27182]))print('matrix A ^ 31415 =')print(cpow(mat,chns[31415]))print('(matrix A ** 27182) ** 31415 =')print(cpow(cpow(mat,chns[27182]),chns[31415]))
First one hundred addition chain lengths: 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 8Knuth chains for addition chains of 31415 and 27182:Exponent: 31415 Addition Chain: [1, 2, 3, 5, 7, 14, 28, 56, 61, 122, 244, 488, 976, 1952, 3904, 7808, 15616, 15677, 31293]Exponent: 27182 Addition Chain: [1, 2, 3, 5, 7, 14, 21, 35, 70, 140, 143, 283, 566, 849, 1698, 3396, 6792, 6799, 13591]1.00002206445416^31415 = 1.99999999989246381.00002550055251^27182 = 1.99999999997926881.000025 + 0.000058i)^27182 = (-0.01128636963542673+1.9730308496660347j)1.000022 + 0.000050i)^31415 = (0.00016144681325535107+1.9960329014194498j)matrix A ^ 27182 =[[mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0 0] [0 mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0] [0 mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0] [mpf('5.0272320351359433e-4092') 0 mpf('5.0272320351359433e-4092') 0 0 0] [0 0 0 0 0 1] [0 0 0 0 1 0]]matrix A ^ 31415 =[[mpf('3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0 0] [0 mpf('3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0] [0 mpf('3.7268602513250562e-4729') 0 mpf('-3.7268602513250562e-4729') 0 0] [mpf('-3.7268602513250562e-4729') 0 mpf('3.7268602513250562e-4729') 0 0 0] [0 0 0 0 0 1] [0 0 0 0 1 0]](matrix A ** 27182) ** 31415 =[[mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148') 0 0 0] [0 mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148') 0 0] [0 mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148') 0 0] [mpf('1.77158544475749e-128528148') 0 mpf('1.77158544475749e-128528148') 0 0 0] [0 0 0 0 0 1] [0 0 0 0 1 0]]This example isincorrect. Please fix the code and remove this message.
|
The addition chains correspond to binary exponentiation.
#langracket(define(chainn); computes a simple addition chain for n(cond[(=n1)'()][(even?n)(definen/2(/n2))(cons(listnn/2n/2)(chainn/2))][(odd?n)(definen-1(-n1))(cons(listnn-11)(chain(-n1)))]))(definemult(let([n0])(λxs(cond[(equal?xs(list'count))n][(equal?xs(list'reset))(set!n0)][else(set!n(+n1))(apply*xs)]))))(define(expt/chainxnchain); computes x^n using the addition chain(defineht(make-hash))(hash-set!ht1x)(define(expt1n)(or(hash-refhtn#f)(let()(definex^n(match(assocnchain)[(list_st)(mult(expt1s)(expt1t))]))(hash-set!htnx^n)x^n)))(expt1n))(define(testxn)(displayln(~a"Chain for "n"\n"(chainn)))(mult'reset)(displayln(~ax" ^ "n" = "(expt/chainxn(chainn))))(displayln(~a"Multiplications: "(mult'count)))(newline))(test1.0000220644541631415)(test1.0000255005525127182)
Output:
Chainfor31415((31415314141)(314141570715707)(15707157061)(1570678537853)(785378521)(785239263926)(392619631963)(196319621)(1962981981)(9819801)(980490490)(490245245)(2452441)(244122122)(1226161)(61601)(603030)(301515)(15141)(1477)(761)(633)(321)(211))1.00002206445416^31415=1.9999999998913485Multiplications:24Chainfor27182((271821359113591)(13591135901)(1359067956795)(679567941)(679433973397)(339733961)(339616981698)(1698849849)(8498481)(848424424)(424212212)(212106106)(1065353)(53521)(522626)(261313)(13121)(1266)(633)(321)(211))1.00002550055251^27182=1.9999999999774538Multiplications:21
# 20230327 Raku programming solutionuseMath::Matrix;classAdditionChains {has ($.idx,$.pos,@.chains,@.lvl,%.pat )isrw;methodadd_chain {# method 1: add chains depth then breadth first until donereturngathergivenself {takemy$newchain = .chains[.idx].clone.append: .chains[.idx][*-1] + .chains[.idx][.pos]; .chains.append:$newchain; .pos == +.chains[.idx]-1 ?? ( .idx +=1 && .pos =0 ) !! .pos +=1; } }methodfind_chain(\nexpwherenexp >0) {# method 1 interface: search for chain ending with n, adding more as neededreturn ([1],)ifnexp ==1;my@chn =self.chains.grep: *.[*-1] ==nexp // [];unless@chn.Bool {repeat {@chn =self.add_chain }until@chn[*-1][*-1] ==nexp }return@chn }methodknuth_path(\ngoal) {# method 2: knuth method, uses memoization to search for a shorter chainreturn []ifngoal <1;untilself.pat{ngoal}:exists {self.lvl[0] = [gatherforself.lvl[0].Array -> \i {forself.knuth_path(i).Array -> \j {unlessself.pat{i +j}:exists {self.pat{i +j} =iandtakei +j } } } ] }returnself.knuth_path(self.pat{ngoal}).append:ngoal }multimethodcpow(\xbase, \chain) {# raise xbase by an addition exponentiation chain for what becomes x**chain[-1]my ($pows,%products) =0,1 =>xbase;%products{0} =xbase ~~Math::Matrix ??Math::Matrix.new([ [1xxxbase.size[1] ]xxxbase.size[0] ]) !!1;forchain.Array -> \i {%products{i} =%products{$pows} *%products{i -$pows};$pows =i }return%products{chain[*-1] } }}my$chn =AdditionChains.new(idx =>0,pos =>0,chains => ([1],),lvl => ([1],),pat => {1=>0} );say'First one hundred addition chain lengths:';.Str.sayfor ( (1..100).map: { +$chn.find_chain($_)[0] -1 } ).rotor:10;my%chns = (31415,27182).map: {$_ =>$chn.knuth_path:$_ };say"\nKnuth chains for addition chains of 31415 and 27182:";say"Exponent: $_\n Addition Chain: %chns{$_}[0..*-2]"for%chns.keys;say'1.00002206445416^31415 = ',$chn.cpow(1.00002206445416,%chns{31415});say'1.00002550055251^27182 = ',$chn.cpow(1.00002550055251,%chns{27182});say'(1.000025 + 0.000058i)^27182 = ',$chn.cpow:1.000025+0.000058i,%chns{27182};say'(1.000022 + 0.000050i)^31415 = ',$chn.cpow:1.000022+0.000050i,%chns{31415};my \sq05 =0.5.sqrt;my \mat =Math::Matrix.new( [[sq05,0,sq05,0,0,0], [0,sq05,0,sq05,0,0], [0,sq05,0, -sq05,0,0], [-sq05,0,sq05,0,0,0], [0,0,0,0,0,1], [0,0,0,0,1,0]] );say'matrix A ^ 27182 =';saymy$res27182 =$chn.cpow(mat,%chns{27182});say'matrix A ^ 31415 =';say$chn.cpow(mat,%chns{31415});say'(matrix A ** 27182) ** 31415 =';say$chn.cpow($res27182,%chns{31415});
First one hundred addition chain lengths:0 1 2 2 3 3 4 3 4 45 4 5 5 5 4 5 5 6 56 6 6 5 6 6 6 6 7 67 5 6 6 7 6 7 7 7 67 7 7 7 7 7 8 6 7 77 7 8 7 8 7 8 8 8 78 8 8 6 7 7 8 7 8 89 7 8 8 8 8 8 8 9 78 8 8 8 8 8 9 8 9 89 8 9 9 9 7 8 8 8 8Knuth chains for addition chains of 31415 and 27182:Exponent: 27182 Addition Chain: 1 2 3 5 7 14 21 35 70 140 143 283 566 849 1698 3396 6792 6799 13591Exponent: 31415 Addition Chain: 1 2 3 5 7 14 28 56 61 122 244 488 976 1952 3904 7808 15616 15677 312931.00002206445416^31415 = 1.99999999989246381.00002550055251^27182 = 1.999999999974053(1.000025 + 0.000058i)^27182 = -0.01128636963542673+1.9730308496660347i(1.000022 + 0.000050i)^31415 = 0.00016144681325535107+1.9960329014194498imatrix A ^ 27182 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0matrix A ^ 31415 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0 0 0 -0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0(matrix A ** 27182) ** 31415 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0
This example isincorrect. Please fix the code and remove this message.
|
Using code atMatrix multiplication#Tcl andMatrix Transpose#Tcl (not shown here).
# 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.## Uses the dichotomic strategy, which produces good results with simpler# coding than for a pluggable non-deterministic strategy.packagerequireTcl8.5namespacepath{::tcl::mathop::tcl::mathfunc}procminchain{n}{if{!($n&($n-1))}{for{seti1}{$i<=$n}{incri$i}{lappendc$i}return$c}elseif{$n==3}{return{123}}return[chain$n[expr{$n>>int(ceil(floor(log($n)/log(2))/2))}]]}procchain{n1n2}{setq[expr{$n1/$n2}]setr[expr{$n1%$n2}]if{$r==0}{return[chain.*[minchain$n2][minchain$q]]}else{return[chain.+[chain.*[chain$n2$r][minchain$q]]$r]}}procchain.+{nsk}{return[lappendns[expr{[lindex$nsend]+$k}]]}procchain.*{nsms}{setn_k[lindex$nsend]foreachm_i$ms{if{$m_i==1}continuelappendns[expr{$n_k*$m_i}]}return$ns}# Generate a lambda term to do exponentiation with a given multiplier command.# Works by extracting information from the addition chain; the lambda term# generated is minimalprocmakeExponentiationLambda{nmulfunc}{setchain[minchain$n]setcmd{seta0}setidxes0foreachc0[lrange$chain0end-1]c1[lrange$chain1end]{lappendidxes[lsearch$chain[expr{$c1-$c0}]]}for{seti1}{$i<[llength$chain]}{incri}{setcmd"$mulfunc \[$cmd\] \$a[lindex $idxes $i]"if{$iin$idxes}{setcmd"set a$i \[$cmd\]"}}lista0$cmd}# Demonstrating application of problem to matrix exponentiationproccount_mult{ab}{incr::countMult;matrix_multiply$a$b}setm31415setn27182setmn[expr{$m*$n}]setpow_m[makeExponentiationLambda$mcount_mult]setpow_n[makeExponentiationLambda$ncount_mult]setpow_mn[makeExponentiationLambda$mncount_mult]setrh[expr{sqrt(0.5)}]setmrh[expr{-$rh}]setA[subst{{$rh0$rh000}{0$rh0$rh00}{0$rh0$mrh00}{$mrh0$rh000}{000001}{000010}}]puts"A**$m";setcountMult0print_matrix[apply$pow_m$A]%6.3fputs"$countMult matrix multiplies"puts"A**$n";setcountMult0print_matrix[apply$pow_n$A]%6.3fputs"$countMult matrix multiplies"puts"A**$mn";setcountMult0print_matrix[apply$pow_mn$A]%6.3fputs"$countMult matrix multiplies"
A**31415 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 19 matrix multipliesA**27182-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 18 matrix multipliesA**853922530-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 37 matrix multiplies
This is very slow (12 mins 50 secs) compared to Go (25 seconds) but, given the amount of calculation involved, not too bad for the Wren interpreter.
import"./dynamic"forStructimport"./fmt"forFmtvarSave=Struct.create("Save",["p","i","v"])varN=32varNMAX=40000varu=List.filled(N,0)// upper boundsvarl=List.filled(N,0)// lower boundsvarout=List.filled(N,0)varsum=List.filled(N,0)vartail=List.filled(N,0)varcache=List.filled(NMAX+1,0)varknown=2varstack=0varundo=List.filled(N*N,null)for(iin0..1)l[i]=u[i]=i+1for(iin0...N*N)undo[i]=Save.new(null,0,0)cache[2]=1varreplace=Fn.new{|x,i,n|undo[stack].p=xundo[stack].i=iundo[stack].v=x[i]x[i]=nstack=stack+1}varrestore=Fn.new{|n|while(stack>n){stack=stack-1undo[stack].p[undo[stack].i]=undo[stack].v}}/* lower and upper bounds */varlower=Fn.new{|n,up|if(n<=2||(n<=NMAX&&cache[n]!=0)){if(up.count>0)up[0]=cache[n]returncache[n]}vari=-1varo=0while(n!=0){if((n&1)!=0)o=o+1n=n>>1i=i+1}if(up.count>0){i=i-1up[0]=o+i}while(true){i=i+1o=o>>1if(o==0)break}o=2while(o*o<n){if(n%o!=0){o=o+1continue}varq=cache[o]+cache[(n/o).floor]if(q<up[0]){up[0]=qif(q==i)break}o=o+1}if(n>2){if(up[0]>cache[n-2]+1)up[0]=cache[n-1]+1if(up[0]>cache[n-2]+1)up[0]=cache[n-2]+1}returni}varinsert=Fn.new{|x,pos|varsave=stackif(l[pos]>x||u[pos]<x)returnfalseif(l[pos]!=x){replace.call(l,pos,x)vari=pos-1while(u[i]*2<u[i+1]){vart=l[i+1]+1if(t*2>u[i]){restore.call(save)returnfalse}replace.call(l,i,t)i=i-1}i=pos+1while(l[i]<=l[i-1]){vart=l[i-1]+1if(t>u[i]){restore.call(save)returnfalse}replace.call(l,i,t)i=i+1}}if(u[pos]==x)returntruereplace.call(u,pos,x)vari=pos-1while(u[i]>=u[i+1]){vart=u[i+1]-1if(t<l[i]){restore.call(save)returnfalse}replace.call(u,i,t)i=i-1}i=pos+1while(u[i]>u[i-1]*2){vart=u[i-1]*2if(t<l[i]){restore.call(save)returnfalse}replace.call(u,i,t)i=i+1}returntrue}vartry// forward declarationvarseqRecur=Fn.new{|le|varn=l[le]if(le<2)returntruevarlimit=n-1if(out[le]==1)limit=n-tail[sum[le]]if(limit>u[le-1])limit=u[le-1]// Try to break n into p + q, and see if we can insert p, q into// list while satisfying bounds.varp=limitvarq=n-pwhile(q<=p){if(try.call(p,q,le))returntrueq=q+1p=p-1}returnfalse}try=Fn.new{|p,q,le|varpl=cache[p]if(pl>=le)returnfalsevarql=cache[q]if(ql>=le)returnfalsewhile(pl<le&&u[pl]<p)pl=pl+1varpu=pl-1while(pu<le-1&&u[pu+1]>=p)pu=pu+1while(ql<le&&u[ql]<q)ql=ql+1varqu=ql-1while(qu<le-1&&u[qu+1]>=q)qu=qu+1if(p!=q&&pl<=ql)pl=ql+1if(pl>pu||ql>qu||ql>pu)returnfalseif(out[le]==0){pu=le-1pl=pu}varps=stackwhile(pu>=pl){if(!insert.call(p,pu)){pu=pu-1continue}out[pu]=out[pu]+1sum[pu]=sum[pu]+leif(p!=q){varqs=stackvarj=quif(j>=pu)j=pu-1while(j>=ql){if(!insert.call(q,j)){j=j-1continue}out[j]=out[j]+1sum[j]=sum[j]+letail[le]=qif(seqRecur.call(le-1))returntruerestore.call(qs)out[j]=out[j]-1sum[j]=sum[j]-lej=j-1}}else{out[pu]=out[pu]+1sum[pu]=sum[pu]+letail[le]=pif(seqRecur.call(le-1))returntrueout[pu]=out[pu]-1sum[pu]=sum[pu]-le}out[pu]=out[pu]-1sum[pu]=sum[pu]-lerestore.call(ps)pu=pu-1}returnfalse}varseq// forward declarationvarseqLen// recursive functionseqLen=Fn.new{|n|if(n<=known)returncache[n]// Need all lower n to compute sequence.while(known+1<n)seqLen.call(known+1)varub=0varpub=[ub]varlb=lower.call(n,pub)ub=pub[0]while(lb<ub&&seq.call(n,lb,[])==0){lb=lb+1}known=nif(n&1023==0)System.print("Cached%(known)")cache[n]=lbreturnlb}seq=Fn.new{|n,le,buf|if(le==0)le=seqLen.call(n)stack=0l[le]=u[le]=nfor(iin0..le)out[i]=sum[i]=0vari=2while(i<le){l[i]=l[i-1]+1u[i]=u[i-1]*2i=i+1}i=le-1while(i>2){if(l[i]*2<l[i+1]){l[i]=((1+l[i+1])/2).floor}if(u[i]>=u[i+1]){u[i]=u[i+1]-1}i=i-1}if(!seqRecur.call(le))return0if(buf.count>0){for(iin0..le)buf[i]=u[i]}returnle}varbinLen=Fn.new{|n|varr=-1varo=-1while(n!=0){if(n&1!=0)o=o+1n=n>>1r=r+1}returnr+o}varmul=Fn.new{|m1,m2|varrows1=m1.countvarrows2=m2.countvarcols1=m1[0].countvarcols2=m2[0].countif(cols1!=rows2)Fiber.abort("Matrices cannot be multiplied.")varresult=List.filled(rows1,null)for(iin0...rows1){result[i]=List.filled(cols2,0)for(jin0...cols2){for(kin0...rows2){result[i][j]=result[i][j]+m1[i][k]*m2[k][j]}}}returnresult}varpow=Fn.new{|m,n,printout|vare=List.filled(N,0)varv=List.filled(N,null)for(iin0...N)v[i]=List.filled(N,0)varle=seq.call(n,0,e)if(printout){System.print("Addition chain:")for(iin0..le){varc=(i==le)?"\n":" "System.write("%(e[i])%(c)")}}v[0]=mv[1]=mul.call(m,m)vari=2while(i<=le){varj=i-1while(j!=0){for(kinj..0){if(e[k]+e[j]<e[i])breakif(e[k]+e[j]>e[i])continuev[i]=mul.call(v[j],v[k])j=1break}j=j-1}i=i+1}returnv[le]}varprint=Fn.new{|m|for(vinm)Fmt.print("$9.6f",v)System.print()}varm=27182varn=31415System.print("Precompute chain lengths:")seqLen.call(n)varrh=(0.5).sqrtvarmx=[[rh,0,rh,0,0,0],[0,rh,0,rh,0,0],[0,rh,0,-rh,0,0],[-rh,0,rh,0,0,0],[0,0,0,0,0,1],[0,0,0,0,1,0]]System.print("\nThe first 100 terms of A003313 are:")for(iin1..100){Fmt.write("$d ",seqLen.call(i))if(i%10==0)System.print()}varexs=[m,n]varmxs=List.filled(2,null)vari=0for(exinexs){System.print("\nExponent:%(ex)")mxs[i]=pow.call(mx,ex,true)Fmt.print("A ^ $d:-\n",ex)print.call(mxs[i])System.print("Number of A/C multiplies:%(seqLen.call(ex))")System.print(" c.f. Binary multiplies:%(binLen.call(ex))")i=i+1}Fmt.print("\nExponent: $d x $d = $d",m,n,m*n)Fmt.print("A ^ $d = (A ^ $d) ^ $d:-\n",m*n,m,n)varmx2=pow.call(mxs[0],n,false)print.call(mx2)
Precompute chain lengths:Cached 1024Cached 2048Cached 3072....Cached 28672Cached 29696Cached 30720The first 100 terms of A003313 are: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 Exponent: 27182Addition chain:1 2 4 8 10 18 28 46 92 184 212 424 848 1696 3392 6784 13568 27136 27182A ^ 27182:-[-0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 -0.500000 0.000000 0.000000][ 0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]Number of A/C multiplies: 18 c.f. Binary multiplies: 21Exponent: 31415Addition chain:1 2 4 8 16 17 33 49 98 196 392 784 1568 3136 6272 6289 12561 25122 31411 31415A ^ 31415:-[ 0.707107 0.000000 0.000000 -0.707107 0.000000 0.000000][ 0.000000 0.707107 0.707107 0.000000 0.000000 0.000000][ 0.707107 0.000000 0.000000 0.707107 0.000000 0.000000][ 0.000000 0.707107 -0.707107 0.000000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000]Number of A/C multiplies: 19 c.f. Binary multiplies: 24Exponent: 27182 x 31415 = 853922530A ^ 853922530 = (A ^ 27182) ^ 31415:-[-0.500000 0.500000 -0.500000 0.500000 0.000000 0.000000][-0.500000 -0.500000 -0.500000 -0.500000 0.000000 0.000000][-0.500000 -0.500000 0.500000 0.500000 0.000000 0.000000][ 0.500000 -0.500000 -0.500000 0.500000 0.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 1.000000 0.000000][ 0.000000 0.000000 0.000000 0.000000 0.000000 1.000000]