如何求數列的通項公式

在做尤拉計劃427題的時候，雖然在OEIS網站上找到了程式碼，但由於它是遞迴的，不太好總結。 所以自己從數字的因式分解中找b(m,1,n)的規律，但這個效率太低了。一晚上才找到幾個，有什麼好的辦法？
例如

b(1,1,9)=150994944
b(2,1,9)=170496000
b(3,1,9)=20731392
b(4,1,9)=1958976
b(5,1,9)=176256
b(6,1,9)=14976
b(7,1,9)=1152
b(8,1,9)=72
b(9,1,9)=0
72 1152 14976 176256

b(1,1,8)=6588344
b(2,1,8)=7068544
b(3,1,8)=919240
b(4,1,8)=94080
b(5,1,8)=9016
b(6,1,8)=784
b(7,1,8)=56
b(8,1,8)=0
56 784 9016 94080

D:\>a
b(1,1,7)=326592
b(2,1,7)=328104
b(3,1,7)=45612
b(4,1,7)=5040
b(5,1,7)=504
b(6,1,7)=42
b(7,1,7)=0
42 504 5040 45864

D:\>a
b(1,1,6)=18750
b(2,1,6)=17250
b(3,1,6)=2550
b(4,1,6)=300
b(5,1,6)=30
b(6,1,6)=0
30 300 2550 19800

已經發現
b(n,1,n)=0
b(n-1,1,n)=n*(n-1)
b(n-2,1,n)=n*(n-1)*2*(n-1)
b(n-3,1,n)=n*(n-1)*(3*(n-1)+2)
b(n-4,1,n)=n*(n-1)*(n-1)*(2*n-1)*(2*n)) //n=7時已經錯誤，n=6也錯
...
b(1,1,n)=n*(n-1)^(n-1)

只知道隨著第一個引數遞增，函式值遞減，遞減幅度約是除以n-1的數量級，如何求數列的通項公式？
展開式

In[4]:= Expand[(2 + 3 (-1 + n)) (-1 + n)  n]

               2      3
Out[4]= n - 4 n  + 3 n

In[5]:= Expand[n*(n-1)*2*(n-1)]

                 2      3
Out[5]= 2 n - 4 n  + 2 n

In[6]:= Expand[n*(n-1)]

              2
Out[6]= -n + n

In[7]:= Expand[n*(n-1)*  (n-1)*(2*n-1)*(2*n)]

            2      3       4      5
Out[7]= -2 n  + 8 n  - 10 n  + 4 n

假定n-5的通項公式是a + a1*n +a2*n^2 +a3*n^3+a4*n^4+a5*n^5+a6*n^6，用n=10~16代入，聯立方程組。

In[3]:= NSolve[
        {
        a + a1*15 +a2*15^2 +a3*15^3+a4*15^4+a5*15^5+a6*15^6 == 47628000,
        a + a1*14 +a2*14^2 +a3*14^3+a4*14^4+a5*14^5+a6*14^6 == 31070312,
        a + a1*13 +a2*13^2 +a3*13^3+a4*13^4+a5*13^5+a6*13^6 == 19614816,
        a + a1*12 +a2*12^2 +a3*12^3+a4*12^4+a5*12^5+a6*12^6 == 11918016,
        a + a1*11 +a2*11^2 +a3*11^3+a4*11^4+a5*11^5+a6*11^6 == 6921200, 
        a + a1*10 +a2*10^2 +a3*10^3+a4*10^4+a5*10^5+a6*10^6 == 3807000, 
        a + a1*16 +a2*16^2 +a3*16^3+a4*16^4+a5*16^5+a6*16^6 == 70963200
        }, 
        {a, a1, a2, a3, a4, a5, a6}]

                                                                 -6
Out[3]= {{a -> 0.0000378713, a1 -> -0.0000177607, a2 -> 3.4431 10  , a3 -> -3., a4 -> 11., a5 -> -13., 

>     a6 -> 5.}}

排除掉約等於0的結果。f(n)=-3*n^3+11*n^4-13*n^5+5*n^6，代回方程中

D:\>a
...
b(10,1,16)=1356595200
b(11,1,16)=70963200
b(12,1,16)=3571200
b(13,1,16)=169200
b(14,1,16)=7200
b(15,1,16)=240
b(16,1,16)=0
240 7200 169200 3571200 70963200
D:\>g++ p427p1.cpp

D:\>a
...
b(9,1,15)=853335000
b(10,1,15)=47628000
b(11,1,15)=2557800
b(12,1,15)=129360
b(13,1,15)=5880
b(14,1,15)=210
b(15,1,15)=0
210 5880 129360 2557800 47628000

b(3,1,9)=20731392
b(4,1,9)=1958976
b(5,1,9)=176256
b(6,1,9)=14976
b(7,1,9)=1152
b(8,1,9)=72
b(9,1,9)=0
72 1152 14976 176256 1959552

驗算結果是正確的，聯絡到每個函式值都有n(n-1)的因子，所以a,a1等於0也合理。不出意料，n=9時，這個結論是錯的。 對展開式因式分解

In[4]:= Factor[-3*n^3+11*n^4-13*n^5+5*n^6]

                2  3
Out[4]= (-1 + n)  n  (-3 + 5 n)

再計算下一個則要用到n=12~18,n=11及以下都無效了

In[7]:= NSolve[
        {
        a*17^7 + a1*17 +a2*17^2 +a3*17^3+a4*17^4+a5*17^5+a6*17^6 == 2095374848, 
        a*16^7 + a1*16 +a2*16^2 +a3*16^3+a4*16^4+a5*16^5+a6*16^6 == 1356595200,
        a*15^7 + a1*15 +a2*15^2 +a3*15^3+a4*15^4+a5*15^5+a6*15^6 == 853335000,
        a*14^7 + a1*14 +a2*14^2 +a3*14^3+a4*14^4+a5*14^5+a6*14^6 == 519384320,
        a*13^7 + a1*13 +a2*13^2 +a3*13^3+a4*13^4+a5*13^5+a6*13^6 == 304346016,
        a*12^7 + a1*12 +a2*12^2 +a3*12^3+a4*12^4+a5*12^5+a6*12^6 == 170615808,
        a*18^7 + a1*18 +a2*18^2 +a3*18^3+a4*18^4+a5*18^5+a6*18^6 == 3155158656
        }, 
        {a, a1, a2, a3, a4, a5, a6}]

Out[7]= {{a -> 6., a1 -> 0.00152974, a2 -> -0.000637064, a3 -> 0.000110092, a4 -> -4.00001, a5 -> 14., 

>     a6 -> -16.}}

In[8]:= Factor[6*n^7-4*n^4+14*n^5-16*n^6]

                  2  4
Out[8]= 2 (-1 + n)  n  (-2 + 3 n)

b函式的個數太多,為什麼不從計算a函式的每一輪迴圈入手呢,在每一輪輸出b函式的和再除以n，就能發現規律了

D:\>a
256 304 56 8 1 T(5)=5345

D:\>a
3125 3900 655 85 10 1 T(6)=79716

D:\>a
46656 60480 9288 1092 120 12 1 T(7)=1403689

   q: 256 304 56 8 1
2 2 2 2  2 2 2 2
2 2 2 2 19 0 0 0
2 2 2 7  0 0 0 0
2 2 2 0  0 0 0 0
0 0 0 0  0 0 0 0

   q: 3125 3900 655 85 10 1
5   5 5 5 5  0
2   2 3 5 5 13
5 131 0 0 0  0
5  17 0 0 0  0
2   5 0 0 0  0
0   0 0 0 0  0

   q: 46656 60480 9288 1092 120 12 1
2 2 2 2  2 2  3 3 3 3 3 3
2 2 2 2  2 2  3 3 3 5 7 0
2 2 2 3  3 3 43 0 0 0 0 0
2 2 3 7 13 0  0 0 0 0 0 0
2 2 2 3  5 0  0 0 0 0 0 0
2 2 3 0  0 0  0 0 0 0 0 0
0 0 0 0  0 0  0 0 0 0 0 0

後面是對輸出用j語言因式分解的結果，顯然第1行是數(n-1)^(n-1)，而倒數第2行是2*(n-1)，其餘各行也有n-1的因子，倒數第3行是(3n-1)*(n-1)

In[12]:=  NSolve[
                 {
                 a + a1* 5 +a2* 5^2 +a3* 5^3 == 655 /5,
                 a + a1* 6 +a2* 6^2 +a3* 6^3 == 1092/6, 
                 a + a1* 7 +a2* 7^2 +a3* 7^3 == 1680/7, 
                 a + a1* 8 +a2* 8^2 +a3* 8^3 == 2448/8
                 }, 
                 {a, a1, a2, a3}]

Out[12]= {{a -> -54., a1 -> 30.3333, a2 -> 0.5, a3 -> 0.166667}}

In[13]:= -54+30.3333*4+0.5*4^2+0.166667*4^3  

Out[13]= 85.9999

In[14]:=  NSolve[
                 {
                 a + a1* 6 +a2* 6^2 +a3* 6^3 == 1092/6, 
                 a + a1* 7 +a2* 7^2 +a3* 7^3 == 1680/7, 
                 a + a1* 8 +a2* 8^2 +a3* 8^3 == 2448/8,
                 a + a1* 9 +a2* 9^2 +a3* 9^3 == 3420/9
                 }, 
                 {a, a1, a2, a3}] 

Out[14]= {{a -> 2., a1 -> 6., a2 -> 4., a3 -> 0.}}

In[15]:= 2+6*5+4*25

Out[15]= 132

同樣遇到了b函式類似的問題，當倒數的行增大時，較小的n就不符合了.
猜測再歸納總結,得出了suma(x)的公式，pow(n,x-2)*(n-1)*((x+1)*n-x+1)，當x小於n/2-1時成立。下面第一部分輸出是遞迴函式，第2部分是總結的公式。%%%前的都符合

D:\>a
2875106592253 3378689621098 5916951335844 2779275901584 180882453206 11648177463 741121459 46444944 2853760 170716 9828 533 26 1 T(1
4)=52724643726460
26 533 9828 170716 2853760 46444944 %%% 741121472 11648192192 180889572863 2781862370560 42437790094335 643013098011648 9686559885836288

D:\>g++ p427p3.cpp

D:\>a
5340158891349 5821757677568 4391828017082 4747106059637 5597140953796 339667984516 20411999384 1211962500 70875000 4063500 226800 12
180 616 28 1 T(15)=89868341011325
28 616 12180 226799 4063500 %%% 70874999 1211962500 20412000000 339668437499 5597353124999 91495195312499 1485451406250000 2397712324218
7499 385087130859375000
D:\>g++ p427p3.cpp

D:\>a
140380859375 5051409531250 2248001868750 2559926075250 97469505375 4508878967600 612032817375 34477178865 1918894080 105185280 56524
80 295680 14880 705 30 1 T(16)=49535458286096
30 705 14880 295680 5652480 105185280 1918894080%%% 34477178880 612032839680 10758893076480 187604171489280 3249056860078080 5594315162
1242880 958422295699783680 -9223372036854775808