COME ON CODE ON

Question : What is the number of rooted plane binary trees with n nodes and height h?

Solution :
Let t_n,h denote the number of binary trees with n nodes and height h.

Lets take a binary search tree of n nodes with height equal to h. Let the number written at its root be the m^thnumber in sorted order, 1 ≤ m ≤ n. Therefore, the left subtree is a binary search tree on m-1 nodes, and the right subtree is a binary search tree on n-m nodes. The maximum height of the left and the right subtrees must be equal to h-1. Now there are two cases :

1. The height of left subtree is h-1 and the height of right subtree is less than equal to h-1 i.e from 0 to h-1.
2. The height of right subtree is h-1 and the height of left subtree is less than equal to h-2 i.e from 0 to h-2, since we have considered the case when the left and right subtrees have the same height h-1 in case 1.

Therefore t_n,h is equal to the sum of number of trees in case 1 and case 2. Let’s find the number of trees in case 1 and case 2.

1. The height of the left subtree is equal to h-1. There are t_m-1,h-1 such trees. The right subtree can have any height from 0 to h-1, so there are such trees. Therefore the total number of such tress are .
2. The height of the right subtree is equal to h-1. There are t_n-m,h-1 such trees. The left subtree can have any height from 0 to h-2, so there are such trees. Therefore the total number of such tress are .

Hence we get the recurrent relation

or

where t_0,0=1 and t_0,i=t_i,0=0 for i>0.
Here’s a sample C++ code.

#include<iostream>
using namespace std;

#define MAXN 35

int main()
{
    long long t[MAXN+1][MAXN+1]={0},n,h,m,i;
    
    t[0][0]=1;
    for (n=1;n<=MAXN;n++)
        for (h=1;h<=n;h++)
            for (m=1;m<=n;m++)
            {
                for (i=0;i<h;i++)
                    t[n][h]+=t[m-1][h-1]*t[n-m][i];
                
                for (i=0;i<h-1;i++)
                    t[n][h]+=t[n-m][h-1]*t[m-1][i];
            }
            
    while (scanf("%lld%lld",&n,&h))
          printf("%lld\n",t[n][h]);
}

Note :

1. The total number of binary trees with n nodes is , also known as Catalan numbers or Segner numbers.
2. t_n,n = 2^n-1.
3. The number of trees with height not lower than h is .
4. There are other recurrence relation as well such as
   .

NJOY!
-fR0DDY

Convert a given decmal number to any other base (either positive or negative).

For example, 100 in Base 2 is 1100100 and in Base -2 is 110100100.

Here’s a simple algorithm to convert numbers from Decimal to Negative Bases :

def tonegativeBase(N,B):
   digits = []
   while i != 0:
	i, remainder = divmod (i, B)
	if (remainder < 0):
	    i, remainder = i + 1, remainder + B*-1
	digits.insert (0, str (remainder))
   return ''.join (digits)

We can just tweak the above algorithm a bit to convert a decimal to any Base. Here’s a sample code :

#include<iostream>
using namespace std;

void convert10tob(int N,int b)
{
     if (N==0)
        return;
     
     int x = N%b;
     N/=b;
     if (x<0)
        N+=1;
        
     convert10tob(N,b);
     printf("%d",x<0?x+(b*-1):x);
     return;
}


int main()
{
    int N,b;
    while (scanf("%d%d",&N,&b)==2)
    {
          if (N!=0)
          {
              convert10tob(N,b);
              printf("\n");
          }
          else
              printf("0\n");
    }
}

NJOY!
-fR0DDY

Question : What is the number of trailing zeroes and the number of digits in N factorial in Base B.

For example,
20! can be written as 2432902008176640000 in decimal number system while it is equal to “207033167620255000000″ in octal number system and “21C3677C82B40000″ in hexadecimal number system. That means that 10 factorial has 4 trailing zeroes in Base 10 while it has 6 trailing zeroes in Base 8 and 4 trailing zeroes in Base 16. Also 10 factorial has 19 digits in Base 10 while it has 21 digits in Base 8 and 16 digits in Base 16. Now the question remains how to find it?

Now we can break the Base B as a product of primes :
B = a^p1 * b^p2 * c^p3 * …
Then the number of trailing zeroes in N factorial in Base B is given by the formulae
min{1/p1(n/a + n/(a*a) + ….), 1/p2(n/b + n/(b*b) + ..), ….}.

And the number of digits in N factorial is :
floor (ln(n!)/ln(B) + 1)

Here’s a sample C++ code :

#include<iostream>
using namespace std;
#include<math.h>

int main()
{
    int N,B,i,j,p,c,noz,k;
    while (scanf("%d%d",&N,&B)!=EOF)
    {
          noz=N;
          j=B;
          for (i=2;i<=B;i++)
          {
              if (j%i==0)
              {   
                 p=0;
                 while (j%i==0)
                 {
                       p++;
                       j/=i;
                 }
                 c=0;
                 k=N;
                 while (k/i>0)
                 {
                       c+=k/i;
                       k/=i;
                 }
                 noz=min(noz,c/p);
              }
          }
          double ans=0;
          for (i=1;i<=N;i++)
          {
              ans+=log(i);
          }
          ans/=log(B);ans+=1.0;
          ans=floor(ans);
          printf("%d %.0lf\n",noz,ans);
    }
}

NJOY!
-fR0DDY

In number theory, a multiplicative function is an arithmetic function f(n) of the positive integer n with the property that f(1) = 1 and whenever a and b are coprime, then f(ab) = f(a)f(b).

This type of functions can be programmed in quick time using the multiplicative formmula. Examples of multiplicative functions include Euler’s totient function, Möbius function and divisor function.

A multiplicative function is completely determined by its values at the powers of prime numbers, a consequence of the fundamental theorem of arithmetic. Thus, if n is a product of powers of distinct primes, say n = p^a q^b …, then f(n) = f(p^a) f(q^b) …

Lets take an example. If we need to calculate the divisor of 12 then div(12) = div(2²3¹) or div(2²)div(3¹). If we are calculating it for a range then it becomes even more easier. We only need to calculte the div(2²)part because the div(3¹) part will already be present in the array. So lets look at an sample code for calculating the divisor of all numbers from 1 to 10000000.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int A[10000000]={0};
int main()
{ 

    int isprime[3163],d,n,e,p,k,i;

    for (n=2;n<3163;n++)
        isprime[n]=1;

    //Sieve for Eratosthenes for Prime
    //Storing the smallest prime which divides n.
    //If A[n]=0 it means it is prime number.
    for(d=2;d<3163;d++)
    {
        if(isprime[d])
        {
            for (n=d*d;n<3163;n+=d)
            {
                isprime[n]=0;
                A[n]=d;
            }
            for (;n<=10000000;n+=d)
                A[n]=d;
        }
    }

    //Applying the formula
    //Divisor(N)=Divisors(N/p^f(N,p))*(f(N,p)+1)
    A[1]=1;
    for(n=2;n<=10000000;n++)
    {
        if (A[n]==0)
           A[n]=2;
        else
        {
           p=A[n],k=n/p,e=2;
           while (k%p==0)
                 k/=p,e++;
           A[n]=A[k]*e;
        }
    }
    printf("time=%.3lf sec.\n",
    (double) (clock())/CLOCKS_PER_SEC);
    while (scanf("%d",&i),i)
    {
          printf("%d\n",A[i]);
    }
    return 0;
}

If you run the above program you will be able to see how fast have we made this program. Similar techniques can be used to calculate the divisor function or the sum of squares of divisors which is also called sigma2. The general formula for calculating sigma function is :

or

-fR0DDY

Pell’s equation is any Diophantine equation of the form
                                     x² – Dy² = 1
where D is a nonsquare integer and x and y are integers. Lagrange proved that for any natural number D that is not a perfect square there are x and y > 0 that satisfy Pell’s equation. Moreover, infinitely many such solutions of this equation exist.

There are many ways to get the solutiont to these equations some of which are given on Wikipedia, MathWorld and probably a good explanation here.

Though at first read it may seem difficult to understand the solution. The basic aim here is to get the convergents of Sqrt(D).  Sqrt(D) = [a₀,a₁,a₂,....,a_r] such that a_r+1=2a₀.
Very important point is that if r is even then solution is p_2r+1 where p_n is the numerator in nth convergent.
The following function is a direct application of the logic given on Mathworld and works for all D<=1000.

double PellSolution(double D)
{
       double x1,y1;
       double Pn[100],Qn[100],a[100],p[100],q[100];
       long long n,r;
       if (sqrt(D)!=floor(sqrt(D)))
       {
          //Initialization
          Pn[0]=0;Qn[0]=1;
          a[0]=floor(sqrt(D));
          p[0]=a[0];q[0]=1;

          n=1;
          Pn[1]=a[0];Qn[1]=D-a[0]*a[0];
          a[1]=floor((a[0]+Pn[1])/Qn[1]);
          p[1]=a[0]*a[1]+1.0;
          q[1]=a[1];

          while (a[n]!=2.0*a[0])
          {
                n=n+1;
                Pn[n]=a[n-1]*Qn[n-1]-Pn[n-1];
                Qn[n]=(D-Pn[n]*Pn[n])/Qn[n-1];
                a[n]=floor((a[0]+Pn[n])/Qn[n]);
                p[n]=a[n]*p[n-1]+p[n-2];
                q[n]=a[n]*q[n-1]+q[n-2];
          }
          r=n-1;
          if (r%2==0)
          {
             for (n=r+2;n<=2*r+1;n++)
             {
                 Pn[n]=a[n-1]*Qn[n-1]-Pn[n-1];
                 Qn[n]=(D-Pn[n]*Pn[n])/Qn[n-1];
                 a[n]=floor((a[0]+Pn[n])/Qn[n]);
                 p[n]=a[n]*p[n-1]+p[n-2];
                 q[n]=a[n]*q[n-1]+q[n-2];
             }
             x1=p[2*r+1];
             y1=q[2*r+1];
          }
          else
          {
              x1=p[r];
              y1=q[r];
          }
          return x1;
       }
       else
           return 0;
}

Note that the above program only returns the value of x for fundamental solution. To get fundamental y you can return y1. Also to get the nth solution u can use the following recurrence relation
x_nth=((x1+y1*Sqrt_D)^nth + (x1-y1*Sqrt_D)^nth)/2
y_nth=((x1+y1*Sqrt_D)^nth – (x1-y1*Sqrt_D)^nth)/(2*Sqrt_D)

-fR0D

Suppose that you have to find all prime numbers till say some value N or let’s say sum of all prime numbers till N. How would you proceed. Let’s play around a bit. A grade 5 student would tell you that a prime has only two divisors 1 and itself and would code it like this :

void basic1(int N)
{
     int i,j,k,c;
     long long s=0;
     for (i=1;i<=N;i++)
     {
         c=0;
         for (j=1;j<=i && c<3;j++)
             if (i%j==0)
                c++;
         if (c==2)
            s+=i;
     }
     cout<<s<<endl;
}

A little smarter kid would tell you that a prime number has no divisor between 2 and N/2 and would code it like this :

void basic2(int N)
{
     int i,j,k,F;
     long long s=0;
     for (i=2;i<=N;i++)
     {
         F=1;
         for (j=2;j<=i/2 && F==1;j++)
             if (i%j==0)
                F=0;
         if (F)
            s+=i;
     }
     cout<<s<<endl;
}

An even smarter kid will tell you that a prime number has no divisor between 2 and its root and would code it like this :

void moderate1(int N)
{
     vector<int> p(N/2,0);
     int i,j,k,F,c=0;
     long long s=0;
     for (i=2;i<=N;i++)
     {
         F=1;
         for (j=2;j*j<=i && F==1;j++)
             if (i%j==0)
                F=0;
         if (F)
         {
            p[c++]=i;
            s+=i;
         }
     }
     cout<<s<<endl;
}

A good maths student can tell you that a prime number has no prime divisors between 2 and its root.

void moderate2(int N)
{
     vector<int> p(N/2,0);
     p[0]=2;
     int i,j,F,c=1;
     long long s=2;
     for (i=3;i<=N;i+=2)
     {
         F=1;
         for (j=0; p[j]*p[j]<=i && F==1; j++)
             if (i%p[j] == 0)
                F=0;
         if (F)
         {
            p[c++]=i;
            s+=i;
         }
     }
     cout<<s<<endl;
}

A good programmer will tell you that Sieve of Eratosthenes is the best way to find the list of prime numbers.
Its algorithm is as follows :

1. Create a contiguous list of numbers from two to some highest number n.
2. Strike out from the list all multiples of two (4, 6, 8 etc.).
3. The list’s next number that has not been struck out is a prime number.
4. Strike out from the list all multiples of the number you identified in the previous step.
5. Repeat steps 3 and 4 until you reach a number that is greater than the square root of n (the highest number in the list).
6. All the remaining numbers in the list are prime.

And the C++ implementation would be :

void sieve(int N)
{
     int x=sqrt(N),i,j;
     vector<bool> S(N+1,0);
     for (i=4;i<=N;i+=2)
         S[i]=1;
     for (i=3;i<=x;i+=2)
         if (!S[i])
            for (j=i*i;j<=N;j+=2*i)
                S[j]=1;

     long long s=0;
     for (i=2;i<=N;i++)
         if (!S[i])
            s+=i;

     printf("%lld\n",s);
}

But then we can optimize even the Sieve of Eratosthenes. If you look closer you will realise that apart from 2 all the other even numbers are of no use and just waste our time and memory. So we should get rid of them. We can do that by sieving only odd numbers. Let an element with index i correspond to number 2*i+1. So in the sieve we only need to go upto (N-1)/2. Also if p = 2*i+1 p² = 4*i^2 + 4*i + 1 which is represented by index 2*i*(i+1). Also the step will be of order 2*i+1. So the code would look something like this:

void improved_sieve(int N)
{
     int M=(N-1)/2;
     int x=(floor(sqrt(N))-1)/2,i,j;
     vector<bool> S(M+1,0);
     for (i=1;i<=x;i++)
         if (!S[i])
            for (j=2*i*(i+1);j<=M;j+=(2*i+1))
                S[j]=1;

     long long s=2;
     for (i=1;i<=M;i++)
         if (!S[i])
            s+=(2*i+1);

     printf("%lld\n",s);
}

and since i am learning Python. Here’s the python code as well :

import math
def improved_sieve(N):
    M=(N-1)/2
    x=(int(math.sqrt(N))-1)/2
    S=[]
    for i in range(M+1):
        S.append(False);

    for i in range (1,x+1):
        if (S[i]==False):
            for j in range (2*i*(i+1),M+1,2*i+1):
                S[j]=True

    s=2;
    for i in range (1,M+1):
        if (S[i]==False):
            s+=(2*i+1)

    print s

N=input()
improved_sieve(N)

Lets compare all these codes on runtime :

Limit10^N where N=	basic1	basic2	moderate1	moderate2	sieve	improved_sieve
4	0.079	0.026	0.003	0.007	0.002	0.001
5	7.562	2.611	0.050	0.117	0.020	0.011
6	-	-	1.160	2.218	0.208	0.121
7	-	-	26.350	44.946	2.183	1.269

All time are in seconds. If you notice that though algorithm 4 seems to be better than 3 it takes more time because of the array refrences that are made several times. Choose the best.

Happy Coding!
-fR0D

This is a very interesting problem with a lots of history. Anyways we will not wonder into it. We shall see how fast can we calculate xⁿ, given x and n where n is a positive integer. The brute force method would be to run a loop from 2 to n and calculate in n-1 steps. We will discuss three methods to do it quickly.

Binary Method
This the most common method used in programs today. It is also called the Dynamic Programming method. Here is an C++ implementation :

int binary(int x,int n)
{
    if (n==1)
       return x;
    if (n%2==0)
    {
       int t=binary(x,n/2);
       return t*t;
    }
    else
        return x*binary(x,n-1);
}

But does this method give the minimum number of multiplications. The answer is NO.
The smallest counterexample is n=15. The binary method takes six multiplications but the best method can do it in five multiplications. We can calculate y = x³ in two multiplications and x¹⁵ = y⁵ in three more, needing only five multiplications in total.

Factor Method
Let n=pq where p is smallest prime factor of n. We can calculate n by first calculating p and then raising this quantity to q-th power. If n is prime we calculate x^n-1 and multiply by x. For example to calculate x⁵⁵, we first calculate y = x⁵ = x⁴x = (x²)²x; then we form y¹¹ = y¹⁰y = (y²)⁵y. The whole process takes eight multiplications.
But does this method give the minimum number of multiplications. The answer is NO.
The smallest counterexample is n=33. The factor method takes seven multiplications but the best method(binary method) can do it in six multiplications.

The “Power Tree Method”
Lets first see the power tree.

Figure above shows the first few levels of the “power tree.” The (k + 1)-st level of this tree is defined as follows, assuming that the first k levels have been constructed: Take each node n of the lath level, from left to right in turn, and attach below it the nodes
                              n + 1, n + a₁, n + a₂, . . . , n + a_k-1 = 2n
(in this order), where 1, a₁, a₂, . . . , a₂ is the path from the root of the tree to n;but discard any node that duplicates a number that has already appeared in the tree. The C++ implementation of the above can be :

/** LINKU[j], LINKR[j] for 0 <= j <= 2^r;  
   point upwards and to the right, respectively,  
   if j is a number in the tree  
**/  
int LINKU[2050]={0},k=0,LINKR[2050]={0},q,s,nm,m,i,n;   
LINKR[0]=1;   
LINKR[1]=0;   
//11 being number of level   
while (k < 11)   
{   
    n=LINKR[0];m=0;   
    do  
    {   
        q=0,s=n;   
        do  
        {   
            if (LINKU[n+s]==0)   
            {   
                if (q==0)   
                   nm=n+s;   
                LINKR[n+s]=q;   
                LINKU[n+s]=n;   
                q=n+s;   
            }   
            s=LINKU[s];   
        }while (s!=0);   
        if (q!=0)   
        {   
            LINKR[m]=q;   
            m=nm;   
        }   
        n=LINKR[n];   
    }while (n!=0);   
    LINKR[m]=0;   
    k=k+1;   
}

But does this method give the minimum number of multiplications. The answer is NO.
The smallest counterexample is n=77.Other examples are 154 and 233.

Some Analysis

• The first case for which the power tree is superior to both the binary method and the factor method is 23.
• The first case for which the factor method beats the power tree method is 19879 = 103*193. Such cases are rare, only 6 for n<10⁵.
• Power tree never gives more multiplications for the computation of xⁿ than the binary method.

---

If You are not interested in Maths stop reading here.

Let l(n) denote the minimum number of multiplications for a given number n.
λ(n) = floor(lg n), where floor(x) is the largest integer not greater than x.
v(n) = number of 1s in the binary representation of n.

They follow the following recurrence relations :

λ(1)=0, λ(2n) = λ(2n+1) = λ(n) + 1;
v(1)=1, v(2n)=v(n), v(2n+1) = v(n)+1.

The binary method requires λ(n) + v(n) – 1  steps.

So we have following theorems

• l(n)<= λ(n) + v(n) -1.
• l(n)>= ceiling(lg n), where ceiling(x) is the smallest integer not less than x.
• l(2^A) = A.
• l(2^A + 2^B) = A+1 if A>B.
• l(2^A + 2^B + 2^C) = A+2 if A>B>C.

Conjecture

• l(2n) = l(n) +1; It fails since l(191) = l(382) = 11.
• The smallest four values for which l(2n) = l(n) are n = 191,701,743,1111.

Conclusion

The table of l(n) may be prepared for 2<=n<=1000 by using the formula l(n) = min(l(n-1)+1 , l_n) – δ_n,
where l_n = ∞ if n is prime, otherwise l_n = l(p)+l(n/p) if p is the smallest prime dividing n; and δ_n = 1 for n in Table below and 0 otherwise.

23
163
229
319
371
413
453
553
599
645
707
741
813
849
903


43
165
233
323
373
419
455
557
611
659
709
749
825
863
905


59
179
281
347
377
421
457
561
619
667
711
759
835
869
923


77
203
283
349
381
423
479
569
623
669
713
779
837
887
941


83
211
293
355
382
429
503
571
631
677
715
787
839
893
947


107
213
311
359
395
437
509
573
637
683
717
803
841
899
955


149
227
317
367
403
451
551
581
643
691
739
809
845
901
983

Phew!!! That was long.

-fR0D
(notes from TAOCP written by Donald E Knuth)

---

Links :
http://www.research.att.com/~njas/sequences/a003313.txt
http://wwwhomes.uni-bielefeld.de/achim/addition_chain.html

Farey Sequence of order N is the ascending sequnce of all reduced fractions between 0 and 1 that have denominators <=N. For example, the Farey Series of order 4 is

⁰⁄₁, ¹⁄₄, ¹⁄₃, ¹⁄₂, ²⁄₃, ³⁄₄, ¹⁄₁

This series can be generated using the following recurrence relation

x₀=0     y₀=1
x₁=1     y₁=N
x_k+2= floor((y_k+N)/y_k+1)x_k+1- x_k
y_k+2= floor((y_k+N)/y_k+1)y_k+1- y_k

The C++ Implementation would look something like this :

x1=0,y1=1,x2=1,y2=N;
x=1;y=N;
printf("%.0lf/%.0lf, %.0lf/%.0lf",x1,y1,x2,y2);
while (y!=1.0)
{
    x=floor((y1+N)/(y2))*x2-x1;
    y=floor((y1+N)/(y2))*y2-y1;
    printf(", %.0lf/%.0lf",x,y);
    x1=x2,x2=x,y1=y2,y2=y;
}
printf("\n");

Note that the number of terms in the Farey Series is given by :

where phi(m) is the Euler Toteint Function.

---

Did You Know?
77 x 999 = 777 x 99 and
112 x 112 = 12544 and its back order
211 x 211 = 44521

NJOY!!!
-fR0D

Most of us know that to find the GCD of two numbers, Euclid Algorithm is the best algorithm. But thanks to Josef Stein, there is the Binary GCD Algorithm which is touch faster than the Euclid’s Algorithm.

The Euclids Algorithm implemented is :

long long gcd(long long a, long long b)
{
    if(a==0) return(b);
    return(gcd(b%a, a));
}

The Binary GCD Algorithm as given in the The Art of Computer Programming written by D.E.Knuth is as follows :

Given positive integers u and v,
this algorithm finds their greatest common divisor.
A1. [Find power of 2.] Set k←0, and then repeatedly set k← k + 1, 
u ← u/2,v ← v/2, zero or more times until u and v are not both even.
A2. [Initialize.] (Now the original values of u and v have 
been divided by 2^k,and at least one of their present values is odd.)
If u is odd, set t  ← -v, and go to A4. Otherwise set t ← u.
A3. [Halve t.] (At this point, t is even, and nonzero.) Set t ← t/2.
A4. [Is t even?] If t is even, go back to A3.
A5. [Reset max(u, v).] If t > 0, set u ← t; otherwise set v←  -t.
(The larger of u and v has been replaced by |t|,
except perhaps during the first time this step is performed.)
A6. [Subtract.] Set t ← u-v. If t != 0, go back to A3.
Otherwise the algorithm terminates with u* 2^k as the output.

Implementation in C++ would look something like this :

int binarygcd(int u,int v)
{
    int k=0,t=0,i;
    while (!(u&1) && !(v&1))
    {
          k++;
          u>>=1;
          v>>=1;
    }
    if (u&1)
       t=u;
    else
        t=-v;
    do
    {
        while (!(t&1))
              t>>=1;
        if (t>0)
           u=t;
        else
            v=-t;
        t=u-v;
    }while (t);
    for (i=0;i<k;i++)
        u<<=1;
    return(u);
}

Though this program may seem to be longer and have more iterations but still its faster because it requires no division instruction, it depends only on subtraction,parity checking and halving using bit shifting.

-fR0D

This is going to be long even without the proof.
The sum of the first n m-th powers is given by the formula :

where S2 is Stirling numbers of the Second Kind :

m\k| 1    2    3    4    5    6
——————————-
 1  | 1
 2  | 1    1
 3  | 1    3    1
 4  | 1    7    6    1
 5  | 1   15   25   10    1
 5  | 1   31   90   65   15    1

given by the recurrence relation :

for m >= 1,

S2(m,0) = 0,
S2(m,1) = 1,
S2(m,k) = 0 for all all k > m,
S2(m+1,k) = S2(m,k-1) + k*S2(m,k).

and f(x,k) = x*(x-1)*(x-2)*…*(x-k+1).

Let us take two examples :

T(2,n) = S2(2,1)*f(n+1,2)/2 + S2(2,2)*f(n+1,3)/3,
          = 1*(n+1)*n/2 + 1*(n+1)*n*(n-1)/3,
          = n*(n+1)*(1/2 + [n-1]/3),
          = n*(n+1)*(1/2 + n/3 – 1/3),
          = n*(n+1)*(2*n+1)/6.
and

T(4,n) = S2(4,1)*f(n+1,2)/2 + S2(4,2)*f(n+1,3)/3 +
              S2(4,3)*f(n+1,4)/4 + S2(4,4)*f(n+1,5)/5,
          = 1*(n+1)*n/2 + 7*(n+1)*n*(n-1)/3 +
              6*(n+1)*n*(n-1)*(n-2)/4 +
              1*(n+1)*n*(n-1)*(n-2)*(n-3)/5,
          = n*(n+1)*(1/2 + 7*[n-1]/3 + 3*[n-1]*[n-2]/2 +
              [n-1]*[n-2]*[n-3]/5),
          = n*(n+1)*(1/2 + 7*[n-1]/3 + 3*[n^2-3*n+2]/2 +
              [n^3-6*n^2+11*n-6]/5),
          = n*(n+1)*(1/2 + 7*n/3 – 7/3 + 3*n^2/2 – 9*n/2 + 3 +
              n^3/5 – 6*n^2/5 + 11*n/5 – 6/5),
          = n*(n+1)*(n^3/5 + 3*n^2/10 + n/30 – 1/30),
          = n*(n+1)*(6*n^3+9*n^2+n-1)/30,
          = n*(n+1)*(2*n+1)*(3*n^2+3*n-1)/30.

Try higher powers for yourself.
-fR0D