# Finding recurrence formulas from procedural code and output lists associated with integer sequences

This code outputs 37 sequences and takes about 15-minutes to run. I would like to be able to see how many and which variables  { i, j, k, m, n }  each of the 37 sequences corresponds to. The 37 sequences are the output of recurrence formulas, but I am not sure how to retrieve the formulas from the code that generated these sequences in a good way. For example for the output list { 1232, 192192, 52276224, 17878468608, 74764505088, 1238954655744, 1152227829841920, 1534767469349437440, 2517018649733077401600 } I have found more than one formula that returns it and would like to see what mathematical representation or recurrence formula is working here from the code. The code took about 15-minutes in average to run on my machine.

In[1]:=
(* Why these sequences in particular? *)
A000004={0,0,0,0,0,0,0,0,0,0,0};
A000012={1,1,1,1,1,1,1,1,1,1,1};
A000040={5,7,11,13,17,19,23,29,31,37,41};
A038110={1,4,8,16,192,3072,55296,110592,442368,13271040,477757440};
(* What exactly offset1 means? Are you adding to the sequence +1 to each term or are you skipping one term or are you dropping the first element or the last element? Why? *)
A038110offset1={4,8,16,192,3072,55296,110592,442368,13271040,477757440,19110297600};
A038111={15,105,385,1001,17017,323323,7436429,19605131,86822723,3212440751,131710070791};
A060753={3,15,35,77,1001,17017,323323,676039,2800733,86822723,3212440751};
A060753offset1={15,35,77,1001,17017,323323,676039,2800733,86822723,3212440751,131710070791};
A161527={11,27,61,809,13945,268027,565447,2358365,73551683,2734683311,112599773191};
A161527offset1={2,11,27,61,809,13945,268027,565447,2358365,73551683,2734683311};
A072044={25,1225,29645,715715,206841635,14933966047,718188003533,86285158710179,82920037520482019,5974606913975783369,10043314222393291843289};
A072044offset1={1225,29645,715715,206841635,14933966047,718188003533,86285158710179,82920037520482019,5974606913975783369,10043314222393291843289,1688189817927745147112851};
A072045={16,768,18432,442368,127401984,9172942848,440301256704,52836150804480,50722704772300800,3652034743605657600,6135418369257504768000};
A072045offset1={768,18432,442368,127401984,9172942848,440301256704,52836150804480,50722704772300800,3652034743605657600,6135418369257504768000,1030750286035260801024000};
A236435={2,12,96,1152,2304,41472,165888,3981312,119439360,3822059520,7644119040};
A236435offset1={12,96,1152,2304,41472,165888,3981312,119439360,3822059520,7644119040,321052999680};
A236436={1,5,35,385,715,12155,46189,1062347,30808063,955049953,1859834119};
A236436offset1={5,35,385,715,12155,46189,1062347,30808063,955049953,1859834119,76253198879};
A002110={6,30,210,2310,30030,510510,9699690,223092870,6469693230,200560490130,7420738134810};
A002110offset1={30,210,2310,30030,510510,9699690,223092870,6469693230,200560490130,7420738134810,304250263527210};
A002110offset2={210,2310,30030,510510,9699690,223092870,6469693230,200560490130,7420738134810,304250263527210,13082761331670030};
A005867={2,8,48,480,5760,92160,1658880,36495360,1021870080,30656102400,1103619686400};
A005867offset1={8,48,480,5760,92160,1658880,36495360,1021870080,30656102400,1103619686400,44144787456000};
A005867offset2={48,480,5760,92160,1658880,36495360,1021870080,30656102400,1103619686400,44144787456000,1854081073152000};
(*Dropped first two values in each list*)

listofSequences={A000012,A000040,A038110,A038110offset1,A038111,A060753,A060753offset1,A161527,A161527offset1,A161527offset1,A072044,A072044offset1,A072045,A072045offset1,A236435,A236435offset1,A236436,A236436offset1,A002110,A002110offset1,A002110offset2,A005867,A005867offset1,A005867offset2};
(* Why do combine them here to refer to them by indices later instead of using something like A000004[[ nth_term_of_sequence ]] ? *)

x=Length[listofSequences];

(* {var1_, var2_, ...} = Table[{}, number_of_empty_lists]; *)
formulas={};
term1={};
term2={};
term3={};
term4={};
term5={};
term6={};
term7={};
term8={};
term9={};
term10={};
term11={};
term12={};
uniqueterms1={};(* Why are these terms unique? *)
For[i=1,i<=x,i++,Print[StringForm["i=",i]]; (* Why you print to the console the intermediate value of i? *)
(* Break your computation into smaller manageable steps and use the name of each series. Check out this link: https://en.wikipedia.org/wiki/Naming_convention_(programming) . Explain what this code is doing, what input it takes and what output it gives. *)
For[j=1,j<=x,j++,For[k=1,k<=x,k++,For[m=1,m<=x,m++,For[n=1,n<=x,n++,If[(listofSequences[[i]]*listofSequences[[j]])+((listofSequences[[k]]*listofSequences[[m]]*listofSequences[[n]])/(listofSequences[[m]]-listofSequences[[n]]))==A000004,temp1=listofSequences[[i]]*listofSequences[[j]];
temp2=listofSequences[[k]]*listofSequences[[m]]*listofSequences[[n]];
temp3=listofSequences[[m]]-listofSequences[[n]];
temp4=(listofSequences[[k]]*listofSequences[[m]]*listofSequences[[n]])/(listofSequences[[m]]-listofSequences[[n]]);
AppendTo[formulas,{i,j,k,m,n}];
AppendTo[term1,temp1];
AppendTo[term2,temp2];
AppendTo[term3,temp3];
AppendTo[term4,temp4];
AppendTo[term5,Numerator[temp1]];
AppendTo[term6,Numerator[temp2]];
AppendTo[term7,Numerator[temp3]];
AppendTo[term8,Numerator[temp4]];
AppendTo[term9,Denominator[temp1]];
AppendTo[term10,Denominator[temp2]];
AppendTo[term11,Denominator[temp3]];
AppendTo[term12,Denominator[temp4]];]]]]]]
(* These are not formulas these are procedures over lists of integers you could choose a more descriptive name. On integers you can use Mod[m,n] to represent integers so they never become too large. *)
formulas

allterms={term1,term2,term3,term4,term5,term6,term7,term8,term9,term10,term11,term12};
Length[allterms]
DeleteDuplicates[allterms];
Length[allterms]
Length[DeleteDuplicates[allterms]]
xyz1=DeleteDuplicates@Flatten[allterms,1]
xyz=Length[DeleteDuplicates@Flatten[allterms,1]]

For[i=1,i<Length[xyz1],i++,If[Max[xyz1[[i]]]<0,(*Abs if all values are<0*)xyz1[[i]]=Abs[xyz1[[i]]];]]
firsttwoelementsdropped={};
For[i=1,i<=Length[xyz1],i++,AppendTo[firsttwoelementsdropped,Drop[xyz1[[i]],2]]]
firsttwoelementsdropped
Length[firsttwoelementsdropped]
(*gives 51 lists*)
DeleteDuplicates[firsttwoelementsdropped]
Length[DeleteDuplicates[firsttwoelementsdropped]]
(*gives 37 lists*)


Below is the compressed output of the code shown above. Uncomment and evaluate this expression to see it if it takes too long to evaluate on your machine.

(* output = Uncompress[ "1:eJztW33InWUZP/vSqdurm5vOOXXmNmfwyv39QasGRXlE1EqoEE1GiE8MlTSpEbIEzcBgaER/JKyVFSTkKGr9ES0aZCE6gtBU6M2yVFBcWurmdvr9rvu88zzPObp31nrPu/f8c+19nvM89319/K7P5975m276+PUntVqtW+aDXFbdcuv1tnZ1Su9VdTL+qBLJaXNBrlsBsvX2Oa2qtWsliPrI3Kp13fIPtvmW1SHGHI1Lqr0Y11pbZ6KPKSWl2mO44zx+dCpqZ5xSja3mcZcNJCvmgOzFD9Xul+ZilwfPmScb6JhiciEFleQ6uhicV16lVDY0NmXvgvfRufZpwoI3xsRkcnI6G9U+nTe9xYvRhWxddhbPNlk5lVzcQ3IDpd72AbLypAYre+z7iqzRJmO88dhoEa9Vxp0QbLZBRLWU02N5k1Rw7SWU3gRrnMlROe21cqq9DHeDyVglBOOsNWDJqz52TiAnXyTZdgKNsI1G2P9JsLN+yzWyvbfBeqd8gl5le+NCNAYMxcyNqAtc6ZzIcoQKXdGFCUEHF/AblJGsN7DUGWTfO2+1SR5qi7jp+EaTr3PJ0iWwW6UuAZm46kTwtePxBaIhyJmhWuxVOIJxbNAQ0mbs0qaKo/YuB6ttyEYHJUpKKRrcNinDarB2URLMl1TWBsz6lLIifNoruCrY9HxSGaMgnM55ELamCVaiSuN1VBr6A0KsitEpyNpkcBis2z6LeyXYSFuVIL11KsGIIYPfgZh8gGQnAfDppSAfu3pRq9L71oyJGsGezyZb2LO4CBwmOGMcMCB6DUF57VVIsLcXIYJNPvtoM96D32QRAsbxSiWtMgCkVPBBN3iZLjbaZwoCk1WIayonZYFom6HL1GTx3MMAvHnsHV2lBE6jYWEfES+Kp2jAW1vEE+yCaHoMPeVS8rgOIKluXAcnuf0nn5/TPll8VWcFiPnCovE5R0R2rYPtsogtXdTZKeAtFxbhJtgSionOhsR3lwknEQoPGmaxiJiQA49TlUYZKjdFazyQy6d81qq9kmZzxidEzRAAaHiCDx6SgKumALTxMEoxKCKJw6/eiXS68Z9PtYRDg629jUCYoDVli/WdVTTjqZJ4AEJcI6J4PrSU7yCwAh1wcI14ZJhtlwsjiEMBsQyQzclYrSaxgNjlACXAFF5glEW4Q7TCj2dzB8AkQsXKGY1QlfBXCFl+b8gw09hv0wvxKx1BJ+Ri+G7OPsFEtKzql3A9rXQlwbR9GcPyGlY+z50J8vLEnKq1eqHpOnAAPq2hu1JsYARpTjMFdAsflXDDJAQHiic5BG6KYgFB+O3j6zdYfdx9DpGyFjvv/vKHT5L1I7CZAVIyLYu5iEvkKQWoBadDt95A9WUDnCYxIaEwkZyAXJtVBMtICiHTDALb5KJDQgs0D95A8PNIzKq9iishTSAEIbaEDMgrT30iArLIafD9U/KdNGuCS+eB5cf+0hIVeSiHARSBtwAB9rFwGWUVC4MCBERqBTUxHcE48HTh2Cd4Z8QlWNUBDqpNyWKQTyESgBM4mlfaRDJMEGjUUh7VlEJSj4jKqLl8yg4qUu3zKCwiYISDAkEme2s13BTQRAE2KDmvEBScTLBfOB9gH1s4X8ygoSTk+ix5lUKRw0xwIrkSb+Qe2jQI1CgGAdRYIEHuXUC6QAjHP5pw0VnwvVICDgIG+EXcYZB3LqXCV0GwY/Hhkcci8jgiS7SoMHFTCY7aF4iHAGk6aAlN0CVyOGRFHeDxbDADfPkGgvzKC4Gy7Ts3Aemfum1lyUWeGIHhUHEUrEGJUdOPgSI4c4mGKQp6wIRHYolF05IewSNg5lGXIDVGYEYnrnSWFFrIuMi+WSE5oX6JyEgip4gJmyfNjIqMB5CgSjYI3VFqn66YqKr4DJxIMUMiHSpW1qxuoIxBbjUTPKp9vtRzFl0CFmXFGKgEpDqAOUDDyvdj9HiDZ3u97MSuAQbXGQhAM4kSi8U0vNui1czWyLZNXeylkSf2sNC/8bYSfvCKRj6XnCMBOSro3asu3iT6oIHTDloWfUjdJGkILR1ALDsmhPQsyW1FKRMhOvo1Gg1sIYGV8H4ON0SBjuDG/gZSIitCtYB6AaVYmLULfk5wJc0CEMY2ASpl2WnFPuuKfRxaBxR1CE9Bo8cO2DVzq0CHioMUsH4B4+9OauH695TEjA7TAQtpsoOA62DL1M2qS6WgQYGJCJrRxqcsXfyygiJ0QuhZmUmlM0gFIWh1IAyKMngEq+nA6NsNv6iekXnRDdPtkYqVSZFepEr4hW96ujfyu/RdDhiBv2GZKOysEbDCNQA8vJWYJ7Ge9QQp7qHvsQOcW0q/GS57Q6TtlGbj/fPZEk4sLg0MPBoeIfVKqSRYNGFVhAeXJv06oYVAhgXLrG2j4IlBCRtHeJQHiqBFODyaAf62qpS0qKsSnBUxAUy5wLgVygOrZWv4AIIt6iuEJbTyjCModHwXhu21XAXJ24tPo6CDUuFA6JGxY2aIRrksT15EVWSCGRUAagE0z1CSZ9FGlrB/SKzy1QB4z3p9tC8m5xZ8COfo5iTIo1hGx4+cgziDijoxeHc5aarwzU6nU+0n+TfJqyQvk7xE8gLJ30ieIXma5AmQxiKHeP9fJM+S/IbkcaxcffV3IAueP9Cp8mde61Tfnfd0p/rlQ/cfkFA8bqjwxKKrsZ5w8X2SHVi6+iyXuvnXILd/AkvdeNt9nWrDe/98qFr54sH9ZSkoHKImuJgUK+MwPAM0NThQYhHsWpInXwcx3wH5w/1vdKpvXn5np9ryq4lXy7roJyJbjm6fOK5DzozzMHKJIeOAAgI4sIGWpbnVf6OS7n66tKZIu00lHSMRxHPGWf8qjpBQb5s+uTrcVfT4+mG8bCRpvUqynuTKQyAbwU5rqwdp3fXXQXY4mvfFycehCNZLsb/yeY4rXEbyyGsgu+4DWf2j1zvVnodvAYonPvpKUQgibMDbrMCKQpCUUV0hHndrgnGmbA2HQaZOkwpBgqZzebYlfVv/n3aVFAVsI9h4VoHsWuDkDWbEkR8WxR4E2f4PkK1foToveRZk7y8OVq3tL+wu6uyd1Y8VxPUO61ktjr/jtP4gN9pD8hx9VR2gLZccxEZXPfNm12Q9c1W50TtYLZvWJqtLCh+N0eoyuXvEkf0+snIrybcOToaNjWsfBT8bHtnTlblnZi+2qA3tReTG1F703j+2X14Me8S5vcTXu0lufoMWuYkWWfBD8PTEnTsKC7XRrrDQmO2KUvqGu8uKs/ZPd88UKaYyvJfIvouhae/PJ7G79Zr1+4uuekeSp3ZN1Tu+F5g25/eirv6xpKhr4FzyLFl5aoPJ6cRaUeuUBvlDY25p1canOs4XqF5NcgUh8YMXQX78vVc61WNjT+0rGu0dpBf36Z2kl4xSG6UvKUmmMUsXWY48059WXqS9G5/aZP+Zw7jctu+d/agbZ3tn++JGzeH+MXWjn5HTJ58HufdP8KA7Lvv6oZKoekbjYwXtvbPxwmhjOF7SUt90fHlhZ8B4XNT6dlP+s8WOUx3zS604tMIMil0SFCYuR0LevXhNp1v3vTUuFxjX5uWnlaRVH5ifLq/1T8zPKNwMGpkLPN52Zr5KdpnyzH+GSiGt9/hRzf7n4I93QXiaYhZ8J2n1Sih/XP6FzZsbZ0lq50yOxQmS+gbTcG6kzsA0nRapMzHdZ0Tq3AzfyZA6f0NzHqTO1tCeAqmzOR1nP+ocDNGJjzpjQ3fOo87e0JyLONLpjjrbQ3amo87cDErxDc5nBtNTqUvqch37Uxv1/Yb8y3IzN8/gExoNzzt+Pnw3EXx8ncZolKpD6ylTPIPRiKLHBQCP7uRFXQOz5rxFXewZedJg6qcsGllj5p6tqAsym08Q1FUxW7Xwbs5R1DEk/z9Kzh1RYxXxXhHEFeUup/7pZnKQuPEqAV29n+QidtFfCiAHzkAd89vtC1vVNW9e0Ko+dMWm+aUsf+vIRH2Razny+va8ybHTHQ/i/Xu/trlV7fnj6jnV35fOXSDv95yTkPzUe0yivuAOMrTuRJBHrwW5eNMJrepzO7e0qrs2nneKrNVzsEDW6j1XIN117WREffWjFrds0XsaotEP/s84lTzdPP9Q32zhYTvvJmktInmC5CF2GLvBQGvr70FaW1Y1Z59Te6vkzp6jDvVV3EL+h6vNk4OFDfEeAO28XYtE2p6TBiJt70GDxVJe9Z4zEGkbhxvqmx2zfUoD3neaYTQsHg2LR8Pi0bB4NCweDYtHw+JhmrvOkmHx0X/mX9v7/qiAGRUwowJmVMCMCphRATMqYIaqFpglBczoa/foa3drZn1sHH3tHgZPGX3tHn3tHn3tHn3tPl6/846+do++do++do++dg/t1+6jGj0zPPwHnmrEcQ==" ]*)


Bellow is a shortened version, so you know what to expect.

In:= output //Short[#, 10]& (* I used the Out[numbers_] of my local session to condense the output into one dissectable chunk *)

Out:={{{10,56,528,6240,97920,1751040,38154240,1058365440,31677972480,1134275788800,45248407142400},{3,60,280,1232,192192,52276224,17878468608,74764505088,1238954655744,1152227829841920,1534767469349437440},{15,140,616,16016,3267264,993248256,37382252544,309738663936,38407594328064,42632429704151040,62925466243326935040},<<45>>,{-60,-1680,-110880,-14414400,-2940537600,-893923430400,-370084300185600,-236113783518412800,-204946764093982310400,-227490908144320364544000,-335776580421016858066944000},{-1680,-110880,-14414400,-2940537600,-893923430400,-370084300185600,-236113783518412800,-204946764093982310400,-227490908144320364544000,-335776580421016858066944000,-577535718324148995875143680000},{1,1,1,1,1,1,1,1,1,1,1}},51,<<5>>,{{528,6240,<<6>>,45248407142400},<<36>>},37}


I will try to replace the Append lines and also here is another formula that can be used in the For loop:

AppendTo[term1,temp1]; (*change to append the i,j,k,m,n Anumber indexes*)
AppendTo[term1,{temp1,{i,j,k,m,n}}];

(* alternate formula requires additional more nested For loop *)
listofSequences[[i]]/listofSequences[[j]])/((listofSequences[[i]]/listofSequences[[k]])-(listofSequences[[i]]/listofSequences[[m]]))==listofSequences[[n]]/listofSequences[[p]],
AppendTo[formulas, {i, j, k, m, n, p}]
(* Consider using sequence name A0000XX[[ nth_term_of_XX ]]. Write in Latex what is the mathematical background of your effort to give context. When writting about math consider using this method to explain your problem and efforts to solve it before asking for help: https://en.wikipedia.org/wiki/Toulmin_method . When solving your problem try the following heuristics: https://en.wikipedia.org/wiki/How_to_Solve_It *)

• Use Drop[#, 2] & /@ { list_1, list_2, list_3, ... } to drop the first two values in each list. Commented Sep 10, 2019 at 6:41
• Use Select[ expr_, # < n & ] or Cases instead of an If statement. Also use DeleteDuplicates[Flatten[Table[expr, {i, o, n}]]] instead of AppendTo or check the documentation for Sow and Reap for colecting intermediate results. Commented Sep 10, 2019 at 6:55
• If you are trying to process the n_th element of each list you can do it simply with an expression like this var = ( {list_1, list_2, list_3, ... } [[ All, n_th ]] ) instead of using C like indices that are hard to track. If this is part of come complex C experiment then you could try the book Numerical Recipies in C numerical.recipes and generate some of the code in Mathematica using CForm you could also leave them in C and execute them from the front-end via a link. Check documentation for details. Commented Sep 10, 2019 at 7:18
• Use MapIndexed to keep track of the index of each item. Commented Sep 10, 2019 at 7:46

Use FindGeneratingFunction to derive the generating function for a sequence. Recall that the power series of a generating function encodes the sequence in its coefficients. Also FindGeneratingFunction, FindSequenceFunction, InterpolatingPolynomial can be used with Series but this would be wasteful if these series are well-kown. Also instead of using a For you can use Table. The ambiguity of my answer is related to the lack of context/use in the original question, so comment your code if you can, also post sample output as well using //Short.Timing,AbsoluteTiming are useful to show the processing time in your sample output. On the other hand, if you just want to debug this code I suggest you try the Wolfram Workbench debug and profiling functionality. You need to analyze the following expression trees to understand what is going on. Below is a diagram of the problem as it was originally. The problem seems to be asking us to explain how each element of the output tree results from navegating the input tree. To account for each element in the output tree you have to traverse the input tree a very large number of times creating a slightly different path.

Finding a formula in this sense of the word would be similar to finding the shortest distance path between the first node of the input tree where an input first appears and the final destination node of the output tree where a correspondent transformed input ends each time. For example:

I will assume that the user, in this case, wants to answer questions about well-known sequences. A database is the right place to start so you do not waste time in any case. If the query returns an output with Head Missing then a call to FindGeneratingFunction is warranted in my opinion.

Here you have all related series in the https://oeis.org/ database for the sample input sequence shown. ToExpression could be use to evaluate the mathematica implementations of the sequence in question obtained from the database directly on your machine.

A very approachable book on generating functions is Generating- functionology by Herbert S. Wilf (A K Peters) https://www.math.upenn.edu/~wilf/gfology2.pdf. You can use an expression like the one below as template or any of the methods shown above to generate a function associated with a particular sequence or to infer a function that will produce the sequence for successive integers.

The image contains a table with the name of each sequence, a brief description of what it is and several Mathematica implementations of the formulas for each sequence.

The iteration variables  { i, j, k, m, n }  in the code that you provided refer to each element and to each list as many times as there are elements in all the lists that you provided in your code combined at different times during execution, multiplying, dividing, adding, comparing, appending and creating even larger integers (or divisions by zero) many, many times over that is why is so slow. If you want to know the resulting series from this computation query the database only with parts of the integer sequence associated with a series. It is better than reverse engineer hellishly nested for loops an inhuman number of times which kinda defeats the purpose of having computers in the first place. An essential issue with functions defined over the set of Integers for computer science, in general, is the amount of memory that it takes to represent large integers. You should encode them with some other mathematical representation that is practical given a set of constraints. Mathematica uses PackedArrays and other methods to tackle some of these issues however the truly clever and powerful bits lie in the mathematics, not in programming tricks. Design a sound mathematical strategy to tackle obstacles ahead of time and share what you learn.

• Thank you for all the edits, I wasn't expecting that. I will try to modify the code to maintain the formulas, and either edit it here or post it as a new thread later. Thanks! Commented Sep 12, 2019 at 1:13
• Hi, traversing a tree or graph for shortest path is an interesting idea, for now I just can use the output to manually create a formula since programming it to be more user friendly is not easy for me, but yesterday I made this formula with this tool: A060753(n+1) = (A161527(n)*A002110(n))/(A002110(n)-A005867(n)). Commented Sep 20, 2019 at 15:19