# Finding polynomial relations between power series and sequences

In my research I often need to find polynomial relations between power series and number sequences. For example, in my essay "A Multisection of q-Series", I trisect a specific power series $$\, A = A_0 + A_1 + A_2 \,$$ and find a homogeneous cubic polynomial relation $$\, 0 = A_2A_0^2 + A_0A_1^2 + A_1A_2^2\,$$ between the three sections. Another example is in my note "A Remarkable eta-product Identity" which is one of the many "Dedekind eta-product identities" that I collected.

Many years ago I wrote code in PARI/GP which can find such relations automatically, essentially by finding the kernel of a constructed matrix. I could just translate my gp code into Wolfram code but I wonder if the functionality already exists and I just don't know about it. I do essentially the same thing in PARI/GP with sequences instead of power series. Again, I wonder if the functionality already exists since the two functions I wrote are so useful for my own research purposes and it is likely that others may find the functionality useful also.

As a specific, concrete example for sequences, consider the sequence OEIS A006720 given by the initial values $$a(0)=a(1)=a(2)=a(3)=1$$ and the recursion $$a(n) = (a(n-1)a(n-3) + a(n-2)^2) / a(n-4).$$ Define the shifted sequences $$a_k(n) = a(n+k).$$ In Wolfram code define some functions and lists:

ClearAll[a, c1,c2, a0,a1,a2,a3,a4,a5, x1,x2,x3,x4,x5,x6];
(* Define the sequence A006720 *)
a[0] = a[1] = a[2] = a[3] = 1;
a[n_] := (a[n] = (a[n-1]*a[n-3] + a[n-2]^2)/a[n-4]);
(* Initialize the finite shifted sequences *)
a0=Array[a, 26, 0]; a1=Drop[a0, 1]; a2=Drop[a1, 1];
a3=Drop[a2, 1]; a4=Drop[a3, 1]; a5=Drop[a4, 1];
(* Solve for the coefficients of some quadratic homogeneous relations *)
Solve[Table[a0[[n]]*a4[[n]] == c1*a1[[n]]*a3[[n]]+c2*a2[[n]]^2, {n,21}],{c1,c2}]
(* {{c1->1,c2->1}} *)
Solve[Table[a1[[n]]*a5[[n]] == c1*a2[[n]]*a4[[n]]+c2*a3[[n]]^2, {n,21}],{c1,c2}]
(* {{c1->1,c2->1}} *)
Solve[Table[a0[[n]]*a5[[n]] == c1*a1[[n]]*a4[[n]]+c2*a2[[n]]*a3[[n]], {n,21}],{c1,c2}]
(* {{c1->-1,c2->5}} *)


A proposed Wolfram equivalent of my gp function might behave like this:

findSequenceRelation[2, 21, {a0, a1, a2, a3, a4, a5}, {x1, x2, x3, x4, x5, x6}]
(* {x3^2 + x4*x2 - x5*x1, 5*x4*x3 - x5*x2 - x6*x1, x4^2 + x5*x3 - x6*x2} *)


which returns three linearly independent degree $$2$$ polynomial relations (when equated to zero) between the six finite sequences {a0, a1, a2, a3, a4, a5} using the first $$21$$ terms of the sequences and where, in the output, the six sequences are referred to as {x1, x2, x3, x4, x5, x6} similar to the way that the Wolfram function SymmetricReduction uses its optional last argument.

A previous question Efficient way of finding poynomial relation is another concrete example.

Notice the following details:

• The first expression in the returned list, namely x3^2 + x4*x2 - x5*x1 (when equated to zero), corresponds to the first polynomial relation $$a(n)a(n+4) = a(n+1)a(n+3) + a(n+2)^2$$.
• The total degree of the polynomial relations sought is given as the first argument.
• The number of terms of each sequence to use is given as the second argument.
• The list of the sequences themselves is given as the third argument.
• The optional fourth argument is a list of symbols to use to refer to the sequences in the output. The default is to use {x1, x2, x3, ...}.
• If one of the sequences is a constant sequence, then the polynomial relation need not be homogeneous.
• A proposed very similar function findSeriesRelation would find polynomial relations between given power series instead of sequences.
• The Wolfram function FindLinearRecurrence almost solves the degree 1 case for sequences (but see next item) and perhaps a higher degree extension may exist already.
• In general, the sequences and power series need not have any obvious connection (that is, the sequences need not be shifted copies of the same sequence), but if there exists polynomial relations of the given degree, then the function will find a basis for them.
• Really difficult (for me, impossible) to answer without a concrete example. Mar 29, 2022 at 0:22
• @DanielLichtblau I hope you like my concrete example. Mar 29, 2022 at 1:32
• Thank you, it's a good example. I have follow-up questions. Your example shows a quadratic relation while the note mentions cubic. Is it the case that the degree will be known in advance? Also will the relations always be homogeneous? Mar 29, 2022 at 14:52
• @DanielLichtblau Thanks for your comment. I have added further information. Mar 29, 2022 at 16:20

Here is a fairly blunt approach. Given the number of terms seqlen and the number of shifts n to use, and also the degree of desired relations we proceed as follows.

(1) Compute seqlen terms of the sequence.

(2) Compute the n shifted sequences and truncate to get a square array of values.

(3) Form all products of that degree using appropriate subsequences of the shifted sequences. (In the example below I do not use a parameter but just brute-force the quadratic case.)

(4) Set up an appropriate matrix using these values, and find the null space generators.

(5) Use this to recover the polynomial relations.

I'll illustrate with the quadratic example provided. Here the degree is 2 and the number of shifts from the sequence is 5. I use more terms than necessary to make sure the null space is minimal; for this purpose 30 is fine.

a[0] = a[1] = a[2] = a[3] = 1;
a[n_] := (a[n] = (a[n - 1]*a[n - 3] + a[n - 2]^2)/a[n - 4]);
seqlen = 30;
a0 = Array[a, seqlen, 0];


Get the subsequences.

n = 5;
aList = NestList[RotateLeft, a0, n][[All, 1 ;; Length[a0] - n]];
Dimensions[aList]

(* Out[668]= {6, 25} *)


Form the matrix of pairwise list products (including multiplying each by itself, in 0order to allow for square terms in the result). Here I also create a list of "tag" variables, that indicate which list pairs get multiplied. This is useful for going from the nulls to the desired polynomial relations.

products =
Flatten[Table[
aList[[i]]*aList[[j]], {i, Length[aList]}, {j, i, Length[aList]}],
1];
Dimensions[products]
vars = Flatten[
Table[c[i, j], {i, Length[aList]}, {j, i, Length[aList]}], 1];


Find the null space of the transpose of this matrix (that's what will give the needed relations *)

nulls = NullSpace[Transpose@products]

(* Out[687]= {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0}, {0, 0, 0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 5, 0, 0, 0, 0, 0, 0,
0, 0}, {0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0}} *)


Recast this as polynomials in some new variable.

nulls . vars /. c[i_, j_] :> w[i]*w[j]

(* Out[700]= {w[4]^2 + w[3] w[5] - w[2] w[6],
5 w[3] w[4] - w[2] w[5] - w[1] w[6], w[3]^2 + w[2] w[4] - w[1] w[5]} *)


I'll give some thought as to whether I can find a different method. Not optimistic though...

• Thanks for your response. What you have done is essentially what I do in my PARI/GP function but for general degree and without some of the symbolic Wolfram functionality. I was hoping that I would not need to break it into the individual steps as you have done. Mar 30, 2022 at 1:52

Here is a different method, using, in effect, a module Groebner basis. I will skip the machinery I sometimes use for that purpose and just show how to apply to this example. A benefit is we do not need to know what is the degree of the relations we seek. Drawbacks include (1) Getting more relations than we want, some of which are seemingly incidental artifacts that depend on how long are the subsequences we use. (2) We still need to guess at how long should be the subsequences in order to get all the desired relations.

Again I'll set up the subsequences.

a[0] = a[1] = a[2] = a[3] = 1;
a[n_] := (a[n] = (a[n - 1]*a[n - 3] + a[n - 2]^2)/a[n - 4]);
seqlen = 25;
a0 = Array[a, seqlen, 0];
n = 5;
aList = NestList[RotateLeft, a0, n][[All, 1 ;; Length[a0] - n]];


We need "tag" variables that encode the separate sequences, "column" variables across the sequences for purposes of emulating a module GB.

relVars = Array[x, n + 1];
colVars = Array[c, Length[aList[[1]]]];
polys = relVars - Map[Apply[Plus, colVars*#] &, aList];


Just to see what these are I'll show the first one.

In[194]:= First[polys]

(* Out[194]= -c[1] - c[2] - c[3] - c[4] - 2 c[5] - 3 c[6] - 7 c[7] -
23 c[8] - 59 c[9] - 314 c[10] - 1529 c[11] - 8209 c[12] -
83313 c[13] - 620297 c[14] - 7869898 c[15] - 126742987 c[16] -
1687054711 c[17] - 47301104551 c[18] - 1123424582771 c[19] -
32606721084786 c[20] + x[1] *)


To emulate a module GB we use relations that state, in effect, that multiplying elements from different columns gives zero, and multiplying a column by itself does not change the power of the column variable.

rels1 = Flatten[
Table[c[i]*c[j], {i, Length[colVars] - 1},
{j, i + 1, Length[colVars]}]];
rels2 = Map[#^2 - # &, colVars];
rels = Join[rels1, rels2];


Now compute a basis that eliminates the column variables. The "tag" variables record the relations between the columns, that is, the subsequences. I omit most of the superfluous elements. The first should give an idea of what they look like.

Timing[
gb25 = GroebnerBasis[Join[polys, rels], relVars, colVars,
MonomialOrder -> EliminationOrder];]
Length[gb25]
gb25[[1 ;; 4]]

(* Out[247]= {0.901119, Null}

Out[248]= 29

Out[249]= {x[4]^2 + x[3] x[5] - x[2] x[6],
5 x[3] x[4] - x[2] x[5] - x[1] x[6],
x[3]^2 + x[2] x[4] -
x[1] x[5], \
-236798963718008626811297046384258003877373596813360675521791948854292\
0593761146250696187977588738769004474502561260601746120346834022631850\
8657798784033783215934631494973910662634543674003457093143207571646681\
1628550590528055616212109846768476600 x[1] -
27044491370810331188290762241379707799311052794931895815749778699411\
3773466231363304784977294975491642199373990827222182048498982471410553\
9391412559652594614397215830526475519679415048139613260841207050512760\
42269375218380992664891247461283157800 x[2] -
96483469443915695262002387991838761459280695037160264688362434202614\
7005003040465822127678385886609409076273302050297742285317588754658221\
5546889164823078336881862189394489897209816892719788859984735297262476\
53177548982220960502769875411047449400 x[3] +
13967714878116049536262013559513486932993109936369632349923083743651\
4847979224395906260523811184252441991937887170760036003478127779342411\
7038451437053351136275315847563809155754424886698755435321902134468284\
977720968965097348449450452356710881370 x[1] x[3] +
12785009281822849568160605497688142919680039748367652862531380504452\
4569156060799963047387409932708926472127185698993775885025224269959771\
3594754250416793144352570981339591789975874185074537638922127201208366\
65570732766244509885671789873700637395 x[2] x[3] +
35707880370204474429180377777973593921701036566292122548093083954221\
1774054100276567906459144171741953009528946622076437824248381858745855\
1622982738172828492247934204817741508686738386763044122767869381234794\
00774926183759673405232034462138442200 x[4] -
28767587096600675884934225211543018655528512347596301538156254405462\
2106401288834066055276617256499703158061024499659510327811092917199957\
5213098740528808455370109203591096492952747769416664988514430292504712\
55177909578341666896055867553076590715 x[1] x[4] +
11904801793275392074448688666406695404801675377419749175247777691890\
2255297291908203875850692929087965875836824524976462211727372564776476\
5808376560513432686330512524547516798014542045329886802069364662463936\
96007161559398304537920157017753339620 x[2] x[4] +
84707970736924496449297712926091034468365367728227722210476661419308\
4629801878973483616247352763830277204405102742570559965972190952752152\
6695090216791074394040167659011710081681438317899462239480391233751680\
188535917204366029978758151484846600 x[5] -
90047957796272621987209918007244639984486352544932812429620057035991\
3017947713811419707297071366038198903385530976659937691807486478301330\
1772696357876836914514752809083800641826984754032588292987591957517855\
7275797104362370803281696057436372825 x[1] x[5] +
10872881947079824857595225981910246921224910661160068600433183596397\
8861779006517629834232734435022215908022098562411776962814762822639358\
0778475854054619794434660139097196865929142706991302274199893949926110\
365584804636623732354173029658348188 x[2] x[5] +
81049807089672747090493167146861810036031305694682623649307055275900\
0471761530060680251943036052287625235977462847565926180373004783867934\
0771871022974583313967483035164533904011162427261728253012147144589782\
60890009151471170063553316606016920 x[3] x[5] +
57577102137556673011513987121691325004522558731166445296680942796429\
1820167550229927130898590547907017316398315142297645305964295513593400\
4499939530916067123729251115987332915421789632985150262951843166423432\
3661003319153495759595200048826560 x[4] x[5] -
40377621574896386012665234808165860959176132619520488206583741927966\
8639736505110142192177913721045241273601493237616434769860731028262482\
0665854729606836742073684093042510360926341717794044265668222076579087\
61144244065604790027324888950505 x[5]^2 -
30599136804669661040390410080933537937920145496894894152303103478743\
4677252017781535215677319691148737352026162759317756742638828284105905\
6577027142963369841806906222677139201596536750775511189925836570572614\
7801961450600558802835374525485059200 x[6] +
45113044878639624716186466005020331240302287790637181948074456099610\
8213396540915264538997837945010566355512442154803915835718221443264706\
9868476807314385269391491520255971913958584218053787834997569680888955\
840100250463699468154781898973725913 x[1] x[6] -
11598207815834465187770516285943264493967300169949906289280258292253\
0248596009660326308865200517518721257544146861513392309325147384637817\
5385489032233259770330480069487204338352594432329477441811371299069743\
584994427854002495654730945050707460 x[2] x[6] -
17084299214340410360372258186766596107378189342172679478940655860696\
0919185995376937977950649962080102385923383372144118024399831763067653\
8756660102917296731281245739354135077447337334516265141997088143352868\
746926411693696410547185660264210 x[3] x[6] -
21659207927713713170084856896067943385422862631520245274193794859222\
5136327783035493666013859858527682878889200210137557017370694491663244\
1068506307973623595361005046071586847573790381431320955519586429739030\
06625474259411656557590811821800 x[4] x[6] -
26313082045547384747002304826778720131988800798728942137395153274322\
7952980549764009548620219432382573887630425995431657470017140098952260\
5233748863492040569181810978050823155092435590949536922096498614566881\
591063972581411662808306448285 x[5] x[6] +
58911237691744677415173766734047425424318250406032330336390497044461\
8790706374603710922107974257209571702220265534548976027299478059102226\
2644504473168462072822216891133988583460091381877975913956935960511192\
2007099430527030052738435005 x[6]^2} *)

• From your answer it seems to me that this approach has a tradeoff. That is, in exchange for "we do not need to know what is the degree of the relations we seek" we get spurious relations with large coefficients which will be eliminated if more sequence terms are used only to be replaced with other spurious relations with even larger coefficients. Apr 1, 2022 at 23:27
• Because of this tradeoff, the only practical way to use this seems to be by doing it with two different seqlen values and taking the intersection of the Groebner basis results. For example, use seqlen=23 and then seqlen=24 and the intersection will be the three degree 2 relations sought. Apr 1, 2022 at 23:53
• Yes, except there is a question of when the stabilization occurs.Perhaps with say seqlen=18 only two of the three would appear. (I didn't try that specifically, I'm just pointing out a possible further issue.) Apr 2, 2022 at 3:24
• I think that the issue about stabilization is already covered in your "need to guess at how long should be the subsequences". Apr 2, 2022 at 3:29
• Small improvement: We do not need the relations of the form c^2-c for the column marker variables (that is, rels2 in the code). Without them we get a stabilization at seqlen=23` wherein all but the first three found relations in the Groebner basis have degree larger than 2. Apr 2, 2022 at 15:25