#John Kim #Use whatever you like Help:=proc() : print( `IsConst(L),Di(L),IsPol(L) , GuessPol(L,x) `): end: #Di(L): the sequence of differences of consectutive entries of L #For example,Di([1,4,9]); should be [3,5]. Di:=proc(L) local i: [seq(L[i]-L[i-1],i=2..nops(L))]: end: #IsConst(L): is the sequence of numbers L a constant sequence: IsConst:=proc(L) evalb( nops(convert(L,set))<=1 ): end: #IsPol(L): inputs a list of length L and returns #a conjectured degree if it (conj.) a polynomial #and FAIL otherwise IsPol:=proc(L) local L1, k: L1:=L: for k from 1 to nops(L)-4 do L1:=Di(L1): if L1=[0\$nops(L1)] then RETURN(k-1): fi: od: FAIL: end: #GuessPol(L,x): guesses a polynomial expression P(x) #such that P(i)=L[i] for all i between 1 and nops(L) GuessPol:=proc(L,x) local d,P,var, eq,a,i,c,Q,start,finish: start:=time(): d:=IsPol(L): if d=FAIL then finish:=time(): print(finish-start): RETURN(d): fi: P:=add( a[i]*x^i, i=0..d): var:={ seq( a[i], i=0..d)}: eq:={ seq( subs(x=c,P)=L[c], c=1..d+1)}: var:=solve(eq,var): if var=NULL then finish:=time(): print(finish-start): RETURN(FAIL): fi: Q:=subs( var, P): if convert([seq(subs(x=c,Q)-L[c], c=d+2..nops(L))],set)={0} then finish:=time(): print(finish-start): Q: else finish:=time(): print(finish-start): RETURN(FAIL): fi: end: #GuessRat1(L,x,d1,d2): guesses a rational poly expression P(x)/Q(x) #such that deg(P)=d1, deg(Q)=d2, and P(i)/Q(i)=L[i] for all i between #1 and nops(L). L should be longer than d1+d2+4. GuessRat1:=proc(L,x,d1,d2) local a,b,P,Q,var,eq,i,c: if d1+d2+4>nops(L) then RETURN(FAIL): fi: P:=add(a[i]*x^i, i=0..d1-1)+x^d1: Q:=add(b[i]*x^i, i=0..d2): var:={seq(a[i],i=0..d1),seq(b[i],i=0..d2)}: eq:={seq(subs(x=c,P/Q)=L[c],c=1..nops(L))}: var:=solve(eq,var): if var=NULL then RETURN(FAIL): fi: subs(var,P/Q): end: #GuessRat(L,x): guesses a rational poly expression P(x)/Q(x) #that minimizes deg(P)+deg(Q). Returns FAIL if none is found. GuessRat:=proc(L,x) local d,d1,d2: for d from 0 to nops(L)-4 do for d1 from 0 to d do d2:=d-d1: if GuessRat1(L,x,d1,d2) <> FAIL then RETURN(GuessRat1(L,x,d1,d2)): fi: od: od: RETURN(FAIL): end: #FF(x,k): returns falling factorial of degree k in x. FF:=proc(x,k) local i: if k=0 then 1: else mul(x-i,i=0..k-1): fi: end: #BetterGuessPol(L,x): guesses a polynomial expression P(x) #such that P(i)=L[i] for all i between 1 and nops(L) #and FAIL otherwise BetterGuessPol:=proc(L,x) local L1,k,a,start,finish: start:=time(): L1:=L: a[0]:=L1[1]: for k from 1 to nops(L)-4 do L1:=Di(L1): a[k]:=L1[1]/k!: if L1=[0\$nops(L1)] then finish:=time(): print(finish-start): RETURN(expand(add(a[i]*FF(x-1,i),i=0..k))): fi: od: finish:=time(): print(finish-start): FAIL: end: #For the polynomial sequence x^80+x^23+3*i (100 terms), BetterGuessPol runs in 0.032 seconds. #GuessPol runs in 0.375 seconds. In this case, over 10 times better!