###################################################################### ##NCF.txt: Save this file as NCF.txt # ## To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read NCF.txt # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Rutgers University , # #zeilberg at math dot rutgers dot edu # ###################################################################### #Created: 2014; Renewed: 2020 print(`Created: 2014`): print(`Renewed: 2020`): print(` This is NCF.txt `): print(`It is one of the packages that accompany the article `): print(` On N-Continued Fractions (temporary title)`): print(`by Robert Dougherty-Bliss and Nathan Fox and Doron Zeilberger`): print(`and also available from Zeilberger's website`): print(``): print(`Please report bugs to zeilberg at math dot rutgers dot edu`): print(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg/ .`): print(`---------------------------------------`): print(`For a list of the Supporting procedures type ezra1();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): print(`---------------------------------------`): print(`For a list of the MAIN procedures type ezra();, for help with`): print(`a specific procedure, type ezra(procedure_name); .`): print(``): print(`---------------------------------------`): with(combinat): ezra1:=proc() if args=NULL then print(` The supporting procedures are: DetectPeriod, IsPer1, OneStep, OneStepPQ, tzug `): print(``): else ezra(args): fi: end: ezra:=proc() if args=NULL then print(`The main procedures are: EvalcfN, EvalPcf, EvalUPcf, Ncf, NcfQ, NcfQg, NcfQplus, PerSeq, PerSeqV, PerSeqP,`): print(` PQseq, PQseqD, PerSeqPv , PQseq, Targem, `): print(` tzug `): elif nops([args])=1 and op(1,[args])=DetectPeriod then print(`DetectPeriod(L,r): inputs a list L, and outputs a pair [a,p] such that starting at the a-th entry`): print(` L is periodic of period p repeated at least r times`): print(`For example, try:`): print(` DetectPeriod([1,2,seq(op([5,6]),i=1..20)],3); `): elif nops([args])=1 and op(1,[args])=EvalcfN then print(`EvalcfN(L,N): evaluates the N-continued fraction L. Try: `): print(` EvalcfN([4,2,5],1); `): elif nops([args])=1 and op(1,[args])=EvalPcf then print(`EvalPcf(L,N): evaluates the purely periodic N-continued fraction L. Try:`): print(`EvalPcf([2,3,5],2);`): elif nops([args])=1 and op(1,[args])=EvalUPcf then print(`EvalUPcf(L,N): evaluates the ultimately-periodic N-continued fraction L=[L1,L2]. Try:`): print(`EvalUPcf([[1],[2,3,5]],2);`): elif nops([args])=1 and op(1,[args])=IsPer1 then print(`IsPer1(L,p,r): Is p a possible period, with at least r repeats`): elif nops([args])=1 and op(1,[args])=Ncf then print(`Ncf(x,N,K): the first K terms of the best N-continued fraction of the real number x, try:`): print(` Ncf(sqrt(2.),2,20); `): elif nops([args])=1 and op(1,[args])=NcfQ then print(`NcfQ(D1,N,K): The periodic N-continued fraction of sqrt(D1) or FAIL, by trying out the first K terms, Try:`): print(` NcfQ(2,2,200); `): elif nops([args])=1 and op(1,[args])=NcfQg then print(`NcfQg(D1,N,K): The periodic N-continued fraction of sqrt(D1), followed by the left-overs`): print(` or FAIL, by trying out the first K terms, Try:`): print(` NcfQg(2,2,200); `): elif nops([args])=1 and op(1,[args])=NcfQplus then print(`NcfQplus(D1,N,K): The periodic N-continued fraction of sqrt(D1) or FAIL, followed by the first K terms of the N-continued fraction`): print(`by trying out the first K terms, Try:`): print(` NcfQplus(2,2,200); `): elif nops([args])=1 and op(1,[args])=OneStep then print(`OneStep(zug,D1,N): one step in the Best N-Continued Fraction expansion of sqrt(D1). Inputs`): print(`a pair of rational numbers zug=[a,b], a positive integer D1 (not a perfect square), and a positive `): print(`integer N, and outputs the pair [x,[a1,b1]] where x is the integer part of a+b*sqrt(D1) and`): print(` a+b*sqrt(D1)-x=a1+b1*sqrt(D1). Try: `): print(` OneStep([0,1],5,1); `): elif nops([args])=1 and op(1,[args])=OneStepPQ then print(`OneStepPQ(zug): one step in generating the simple 1-continued fraction takes [P[i],Q[i]] to [P[i+1],Q[i+1]]`): print(` Try: `): print(` OneStepPQ([0,1],17); `): elif nops([args])=1 and op(1,[args])=PerSeq then print(`PerSeq(N,K,M): the ultimate periodic Best N-continued fraction expansions`): print(`for square-roots of non-perfect squares up to K^2 by looking up to M terms`): print(`It alo returns, the list of the periods, followed by the list of the failures, followed by the pair`): print(`[champion, record], giving the integer whose square-root yields the highest period, followed by that period`): print(`among those that have periods by investigating M terms. `): print(`Try: `): print(`PerSeq(1,20,300); `): elif nops([args])=1 and op(1,[args])=PerSeqP then print(`PerSeqP(N,K,M): the ultimate periodic Best N-continued fraction expansions`): print(`for square-roots of the first K primes up to the by looking up to M terms`): print(`It alo returns, the list of the periods, followed by the list of the failures, followed by the pair`): print(`[champion, record], giving the integer whose square-root yields the highest period, followed by that period`): print(`amont those that have periods by investigating M terms. `): print(`Try: `): print(`PerSeqP(1,20,300); `): elif nops([args])=1 and op(1,[args])=PerSeqPv then print(`PerSeqPv(N,K,M): Verbose version of PerSeqP(N,K,M) (q.v.) `): print(`Try: `): print(`PerSeqPv(1,20,300); `): elif nops([args])=1 and op(1,[args])=PerSeqV then print(`PerSeqV(N,K,M): verbose form of PerSeq(N,K,M) (q.v.)`): print(`Try: `): print(`PerSeqV(1,10,300); `): elif nops([args])=1 and op(1,[args])=PQseq then print(` PQseq(D1,N,K): the sequence of [P_i,Q_i] in the N-continued fraction of sqrt(D1), using K terms. `): print(`If no period is detected, or the propery that is true for N=1 does not hold, then it returns FAIL1, and FAIL2 resp.`): print(` Try: `) : print(` PQseq(109,1,200); `): elif nops([args])=1 and op(1,[args])=PQseqD then print(` PQseqD(D1,K): the sequence of [P_i,Q_i] in the 1-continued fraction of sqrt(D1), using K terms. `): print(`If no period is detected, or the propery that is true for N=1 does not hold, then it returns FAIL1, and FAIL2 resp.`): print(`It is done directly, and only for N=1, so it is not part of the input. `): print(` Try: `) : print(` PQseqD(109,200); `): elif nops([args])=1 and op(1,[args])=Targem then print(`Targem(L,r): inputs a list L, and outputs two lists [L1,L2], such that L is (conjecturally) L1,L2$(infinity),`): print(` with L2 repeated at least r times. Try:`): print(`Targem([1,2,seq(op([5,6]),i=1..20)],2); `): elif nops([args])=1 and op(1,[args])=tzug then print(`tzug(zug): inputs a pair of rational numbers zug=[a,b], outputs the pair [P,Q] such that`): print(` a=P/Q and b=1/Q if possible, or returns FAIL. Try: `): print(` tzug([7/6,1/18]); `): print(``): else print(`There is no ezra for`,args): fi: end: Digits:=1000: #NcfOld(x,N,K): the first K terms of the best N-continued fraction of the real number x, try: #NcfOld(sqrt(2.),2,20); NcfOld:=proc(x,N,K) local gu,x1,lu: if type(x,integer) then RETURN([x]): fi: x1:=trunc(x): if K=1 then RETURN([x1]): fi: lu:=N/(x-x1): gu:=NcfOld(lu,N,K-1): [x1,op(gu)]: end: #Ncf(x,N,K): the first K terms of the best N-continued fraction of the real number x, try: #Ncf(sqrt(2.),2,20); Ncf:=proc(x,N,K) local gu,x1,i: x1:=x: gu:=[]: for i from 1 to K do if type(x1, integer) then RETURN([op(gu),x1]): fi: gu:=[op(gu),trunc(x1)]: x1:=N/(x1-trunc(x1)): od: gu: end: #Targem(L,r): inputs a list L, and outputs two lists [L1,L2], such that L is (conjecturally) L1,L2$(infinity). Try: #L2 must have period<=nops(L)/3 #Targem([1,2,seq(op([5,6]),i=1..20)],r); Targem:=proc(L,r) local L1,L2,P: P:=DetectPeriod(L,r): if P=FAIL then RETURN(FAIL): fi: L1:=[op(1..P[1],L)]: L2:=[op(P[1]+1..P[1]+P[2],L)]: [L1,L2]: end: #EvalcfN(L,N): evaluates the N-continued fraction L. Try: #EvalcfN([4,2,5],1); EvalcfN:=proc(L,N) : if nops(L)=1 then L[1]: else normal(L[1]+N/EvalcfN([op(2..nops(L),L)],N)): fi: end: #EvalPcf(L,N): evaluates the purely periodic N-continued fraction L. Try: #EvalPcf([2,3,5],2); EvalPcf:=proc(L,N) local x,gu,mu: gu:=EvalcfN([op(L),x],N) : gu:=numer(normal(x-gu)): mu:=solve(gu,x): if evalf(mu[1])>0 then RETURN(mu[1]): elif evalf(mu[2])>0 then RETURN(mu[1]): else RETURN(FAIL): fi: end: #EvalUPcf(L,N): evaluates the ultimately-periodic N-continued fraction L=[L1,L2]. Try: #EvalUPcf([[1],[2,3,5]],2); EvalUPcf:=proc(L,N) local L1,L2,gu: if not (type(L,list) and nops(L)=2 and type(L[1],list) and type(L[2],list)) then print(`Bad input`): RETURN(FAIL): fi: L1:=L[1]: L2:=L[2]: gu:=EvalPcf(L2,N): gu:=[op(L1),gu]: EvalcfN(gu,N) : end: #IsPer1(L,p,r): Is p a possible period, with at least r repeats IsPer1:=proc(L,p,r) local i,N,i1,j: N:=nops(L): if N<=p*r then print(`Make L bigger`): RETURN(FAIL): fi: for i from 1 to r do if nops({seq(L[N-p*i1],i1=0..i)})<>1 then RETURN(false): fi: od: if [op(N-2*p+1..N-p,L)]<>[op(N-p+1..N,L)] then RETURN(false): fi: if nops({seq([op(N-j*p+1..N-(j-1)*p,L)],j=1..r)})<>1 then RETURN(false); fi: true: end: #DetectPeriod(L,r): inputs a list L, and outputs a pair [a,p] such that starting at the a-th entry # L is periodic of period p repeated at least r times DetectPeriod:=proc(L,r) local p,N,i: N:=nops(L): for p from 1 to trunc(N/r)-4 do if IsPer1(L,p,r) then for i from 1 to N-p*r while [op(i..i+p-1,L)]<>[op(i+p..i+2*p-1,L)] do od: RETURN(i-1,p): fi: od: FAIL: end: #PerSeqOld(N,K,M): the first terms sequence going over all square-roots of non-perfect squares up to K^2 of periods of best N-continued fractions #Try: #PerSeqOld(1,20,300); PerSeqOld:=proc(N,K,M) local gu,lu,ju,mu,ku: local i: gu:=[]: mu:=[]: ku:=[]: for i from 1 to K^2 do if not type(sqrt(i),integer) then lu:=Ncf(evalf(sqrt(i)),N,M): ju:=Targem(lu,4): if ju=FAIL then ku:=[op(ku),i]: gu:=[op(gu),[i,FAIL]]: mu:=[op(mu),FAIL]: else gu:=[op(gu),[i,ju]]: mu:=[op(mu),nops(ju[2])]: fi: fi: od: gu,mu,ku: end: #PerSeq(N,K,M): the first terms sequence going over all square-roots of non-perfect squares up to K^2 of periods of best N-continued fractions #Try: #PerSeq(1,20,300); PerSeq:=proc(N,K,M) local gu,ju,mu,ku,aluf,si: local i: gu:=[]: mu:=[]: ku:=[]: si:=0: aluf:=0: for i from 1 to K^2 do if not type(sqrt(i),integer) then ju:=NcfQ(i,N,M): if ju=FAIL then ku:=[op(ku),i]: gu:=[op(gu),[i,FAIL]]: mu:=[op(mu),FAIL]: else gu:=[op(gu),[i,ju]]: mu:=[op(mu),nops(ju[2])]: if nops(ju[2])>si then si:=nops(ju[2]): aluf:=i: fi: fi: fi: od: gu,mu,ku,[aluf,si]: end: #OneStep(zug,D1,N): one step in the Best N-Continued Fraction expansion of sqrt(D1). Inputs #a pair of rational numbers zug=[a,b], a positive integer D1 (not a perfect square), and a positive #integer N, and outputs the pair [x,[a1,b1]] where x is the integer part of a+b*sqrt(D1) and #a+b*sqrt(D1)-x=a1+b1*sqrt(D1). Try: #OneStep([0,1],5,1); OneStep:=proc(zug,D1,N) local x,a,b: a:=zug[1]: b:=zug[2]: x:=trunc(a+b*sqrt(D1)): [x,[N*(a-x)/((a-x)^2-b^2*D1),N*b/(b^2*D1-(a-x)^2)]]: end: #NcfQ(D1,N,K): The periodic N-continued fraction of sqrt(D1) or FAIL, by trying out the first K terms, Try: #NcfQ(2,2,200); NcfQ:=proc(D1,N,K) local zug,x,lu,gu,zug1,i,j,k, yofee: zug:=[0,1]: zug1:=OneStep(zug,D1,N): lu:=zug1: gu:=[lu]: zug:=lu: for i from 1 to K do lu:=OneStep(zug[2],D1,N): if member(lu,gu) then for j from 1 to nops(gu) while gu[j]<>lu do od: yofee:=[[seq(gu[k][1],k=1..j-1)],[seq(gu[k][1],k=j..nops(gu))] ]: if simplify(EvalUPcf(yofee,N)-sqrt(D1))<>0 then print(`Something terrible happened`): RETURN(FAIL): else RETURN(yofee): fi: else gu:=[op(gu),lu]: zug:=lu: fi: od: FAIL: end: #PerSeqP(N,K,M): the first terms sequence going over all square-roots of the first K primes of periods of best N-continued fractions #Try: #PerSeqP(1,20,300); PerSeqP:=proc(N,K,M) local gu,ju,mu,ku,aluf,si,p: local i: gu:=[]: mu:=[]: ku:=[]: si:=0: aluf:=0: for i from 1 to K do p:=ithprime(i): ju:=NcfQ(p,N,M): if ju=FAIL then ku:=[op(ku),p]: gu:=[op(gu),[p,FAIL]]: mu:=[op(mu),FAIL]: else gu:=[op(gu),[p,ju]]: mu:=[op(mu),nops(ju[2])]: if nops(ju[2])>si then si:=nops(ju[2]): aluf:=p: fi: fi: od: gu,mu,ku,[aluf,si]: end: #PerSeqV(N,K,M): verbose form of PerSeq(N,K,M) (q.v) #Try: #PerSeqV(1,20,300); PerSeqV:=proc(N,K,M) local gu,t0,i: t0:=time(): gu:=PerSeq(N,K,M): print(`All The `, N, `-continued fractions of Square-Roots of Positive Integers Up to`, K^2-1): print(`By Shalosh B. Ekhad `): print(``): if gu[3]<>[] then print(`The square-roots of the following integers do not seem to have an umtimately periodic,`, N, `-continued fraction `): print(`judging by the first`, M, ` terms `): print(``): lprint(gu[3]): else print(`All of them have periodic continued fractions.`): fi: print(``): print(` Among those that do, the largest period is the square-root of`, gu[4][1], `whose period is`, gu[4][2]): print(`Here are all the successful ones `): for i from 1 to nops(gu[1]) do print(gu[1][i][1], `: `, gu[1][i][2] ): od: print(`Finally, the sequence of periods (for the OEIS) is`): lprint(gu[2]): print(``): print(`This ends this article, that took`, time()-t0, `seconds to generate. `): end: #PerSeqPv(N,K,M): verbose form of PerSeqP(N,K,M) (q.v) #Try: #PerSeqPv(1,20,300); PerSeqPv:=proc(N,K,M) local gu,t0,i: t0:=time(): gu:=PerSeqP(N,K,M): print(`All The `, N, `-continued fractions of the first`, K, ` primes `): print(`By Shalosh B. Ekhad `): print(``): if gu[3]<>[] then print(`The square-roots of the following primes do not seem to have an umtimately periodic,`, N, `-continued fraction `): print(`judging by the first`, M, ` terms `): print(``): lprint(gu[3]): else print(`All of them have periodic continued fractions.`): fi: print(``): print(` Among those that do, the largest period is the square-root of`, gu[4][1], `whose period is`, gu[4][2]): print(`Here are all the successful ones `): for i from 1 to nops(gu[1]) do print(gu[1][i][1], `: `, gu[1][i][2] ): od: print(`Finally, the sequence of periods (for the OEIS) is`): lprint(gu[2]): print(``): print(`This ends this article, that took`, time()-t0, `seconds to generate. `): end: #NcfQplus(D1,N,K): The periodic N-continued fraction of sqrt(D1) or FAIL, by trying out the first K terms, Try: #NcfQplus(2,2,200); NcfQplus:=proc(D1,N,K) local zug,x,lu,gu,zug1,i,j,k, yofee: zug:=[0,1]: zug1:=OneStep(zug,D1,N): lu:=zug1: gu:=[lu]: zug:=lu: for i from 1 to K do lu:=OneStep(zug[2],D1,N): if member(lu,gu) then for j from 1 to nops(gu) while gu[j]<>lu do od: yofee:=[[seq(gu[k][1],k=1..j-1)],[seq(gu[k][1],k=j..nops(gu))] ]: if simplify(EvalUPcf(yofee,N)-sqrt(D1))<>0 then print(`Something terrible happened`): RETURN(FAIL): else RETURN(yofee): fi: else gu:=[op(gu),lu]: zug:=lu: fi: od: gu:=[seq(gu[i][1],i=1..nops(gu))]: FAIL,gu: end: #NcfQg(D1,N,K): The periodic N-continued fraction of sqrt(D1) followed by the corresponding left-overs #or FAIL, by trying out the first K terms, Try: #NcfQg(2,2,200); NcfQg:=proc(D1,N,K) local zug,x,lu,gu,zug1,i,j,k, yofee,tofee: zug:=[0,1]: zug1:=OneStep(zug,D1,N): lu:=zug1: gu:=[lu]: zug:=lu: for i from 1 to K do lu:=OneStep(zug[2],D1,N): if member(lu,gu) then for j from 1 to nops(gu) while gu[j]<>lu do od: yofee:=[[seq(gu[k][1],k=1..j-1)],[seq(gu[k][1],k=j..nops(gu))] ]: tofee:=[[seq(gu[k][2],k=1..j-1)],[seq(gu[k][2],k=j..nops(gu))] ]: if simplify(EvalUPcf(yofee,N)-sqrt(D1))<>0 then print(`Something terrible happened`): RETURN(FAIL): else RETURN(yofee,tofee): fi: else gu:=[op(gu),lu]: zug:=lu: fi: od: FAIL, gu: end: #tzug(zug): inputs a pair of rational numbers zug=[a,b], outputs the pair [P,Q] such that #a=P/Q and b=1/Q if possible, or returns FAIL. Try: #tzug([7/6,1/18]); tzug:=proc(zug) local a,b,c: a:=zug[1]: b:=zug[2]: if numer(b)<>1 then RETURN(FAIL): fi: c:=denom(b)/denom(a): if not type(c,integer) then FAIL: else [numer(a)*c,denom(b)]: fi: end: #PQseq(D1,N,K): the sequence of [P_i,Q_i] in the N-continued fraction of sqrt(D1), using K terms. #If no period is detected, or the propery that is true for N=1 does not hold, then it returns FAIL1, and FAIL2 resp. #Try: #PQseq(109,1,200); PQseq:=proc(D1,N,K) local gu,gu1,gu2,ka,FAIL1, FAIL2, lu1,lu2,i: gu:=NcfQg(D1,N,K): if nops([gu])=2 and gu[1]=FAIL then RETURN(FAIL1): fi: gu:=gu[2]: gu1:=gu[1]: gu2:=gu[2]: lu1:=[]: for i from 1 to nops(gu1) do ka:=tzug(gu1[i]): if ka=FAIL then RETURN(FAIL2): else lu1:=[op(lu1),ka]: fi: od: lu2:=[]: for i from 1 to nops(gu2) do ka:=tzug(gu2[i]): if ka=FAIL then RETURN(FAIL2): else lu2:=[op(lu2),ka]: fi: od: [lu1,lu2]: end: #OneStepPQ(zug): one step in generating the simple 1-continued fraction takes [P[i],Q[i]] to [P[i+1],Q[i+1]] #Try: #OneStepPQ([0,1],17); OneStepPQ:=proc(zug,D1) local a,c: a:=trunc((zug[1]+sqrt(D1))/zug[2]): c:=a*zug[2]-zug[1]: [a,[c,(D1-c^2)/zug[2]]]: end: #PQseqD(D1,K): The periodic 1-continued fraction of sqrt(D1) done via the PQ recurrence #Try: #PQseqD(19,200); PQseqD:=proc(D1,K) local zug,x,lu,gu,zug1,i,j,k, yofee,tofee: zug:=[0,1]: zug1:=OneStepPQ(zug,D1): lu:=zug1: gu:=[lu]: zug:=lu: for i from 1 to K do lu:=OneStepPQ(zug[2],D1): if member(lu,gu) then for j from 1 to nops(gu) while gu[j]<>lu do od: yofee:=[[seq(gu[k][1],k=1..j-1)],[seq(gu[k][1],k=j..nops(gu))] ]: tofee:=[[seq(gu[k][2],k=1..j-1)],[seq(gu[k][2],k=j..nops(gu))] ]: RETURN([yofee,tofee]): else gu:=[op(gu),lu]: zug:=lu: fi: od: FAIL, gu: end: