I use formula ((n+b)!/(b!*n!)) * ((n+a)!/(a!*n!))
This is my code :

{$N+}
Var
   D,D1 : Array [1..100] Of Byte;
   Res,Res1 : Extended;
   I,J,N,A,B,Max,T,T1 : Longint;
Begin
Read(N,A,B);
For I:=1 To N+B Do Begin Inc(T);D[T]:=I;End;
For I:=1 To N+A Do Begin Inc(T);D[T]:=I;End;

For I:=1 To B Do Begin Inc(T1);D1[T1]:=I;End;
For I:=1 To A Do Begin Inc(T1);D1[T1]:=I;End;
For I:=1 To N Do Begin Inc(T1);D1[T1]:=I;End;
For I:=1 To N Do Begin Inc(T1);D1[T1]:=I;End;

For I:=1 To 2*N+A+B Do
For J:=1 To 2*N+A+B Do
If D[I]=D1[J] Then Begin D[I]:=0;D1[J]:=0;Break;End;

For I:=1 To 2*N+A+B Do
For J:=1 To 2*N+A+B Do
If (D1[J]<>0) And (D[I]<>0) Then
If D[I] Mod D1[J] = 0 Then Begin D[I]:=D[I] Div D1[J];D1[J]:=0;End;

Res:=1;
For I:=1 To 2*N+A+B Do If D[I]<>0 Then Res:=Res*D[I];

Res1:=1;
For I:=1 To 2*N+A+B Do If D1[I]<>0 Then Res1:=Res1*D1[I];

Writeln(Res/Res1:0:0);
End.

var
 n,a,b:integer;
 p:array[0..20,0..15,0..15]of extended;
procedure init;
var
 i,j:integer;
begin
 read(n);
 read(a);
 read(b);

 fillchar(p,sizeof(p),0);
 for i:=0 to a do
  for j:=0 to b do
   p[0,i,j]:=1;
 for i:=0 to n do p[i,0,0]:=1;

end;
procedure solve(n,a,b:integer);
var
 i,j:integer;
begin
 if n>=1 then
  for i:=0 to a do
   for j:=0 to b do
    begin
     if p[n-1,i,j]=0 then
      solve(n-1,i,j);
     p[n,a,b]:=p[n,a,b]+p[n-1,i,j];
    end;
end;
begin
 init;
 solve(n,a,b);
 writeln(p[n,a,b]:0:0);
end.

#include <stdio.h>

void main()
{
    int n,a,b,i,j;
    double c=1;
    double p[36];

    scanf("%d %d %d",&n,&a,&b);

      if(a>b)
    {
        i=a; a=b; b=i;
    }
    for(i=0;i<36;i++) p[i]=1;
      for(i=1;i<=n+b;i++)
    {
            for(j=i;j>1;j--)
                p[j]=p[j]+p[j-1];
            if(i==n+a) c=p[a+1];
      }
    printf("%.0lf",c*p[b+1]);
}

Please write wha am I getting CE:
// Boxes
#include <iostream.h>
#include <math.h>

int n, a, b;
long num=0;

int min(int a, int b)
{
 if (a<b) return a;
 else return b;
}

int c(int il)
{
 if (min(a, b)>il)
  return il+1;
 else return min(a, b)+1;
}

int count()
{
 int i;
 for (i=a+b; i>=0; i--)
 {
  num += c(min(i, a+b-i)) * long(pow(double(n), double(i)));
 }
 cout << num;
}

int main()
{
 cin >> n >> a >> b;
 count();
 return 0;
}

I can't find my program's mistake, so I check sb's Ac program with
mine.But to my surprise, the answers are the same.So I submit these
Ac programs like destiny's ,Revenger and NSC's,idinbox's

but all get WA ,why????????????????????

I get the correct results at all the tests when the number of balls
is less than 15. Any test over that is plain smaller. I'm using the
following sequence to calculate the result:

x1 = 1 + nb
xn = xn-1 * ( nb + n ) * ( 1 / n )

where nb is the number of boxes. I use this function for both
colours and then multiply the two results.

Is there something wrong with it?

type
 long = record
              a:array[1..250]of byte;
              l:byte;
        end;
procedure add(var a,b:long);
var l,c,k,i:byte;
begin
  k:=0;
  if a.l>b.l then l:=a.l else l:=b.l;
  i:=1;
  for i:=1 to l do
  begin
    c:=a.a[i];
    a.a[i]:=(c+b.a[i]+k)mod 10;
    k:=(c+b.a[i]+k)div 10;
  end;
  a.l:=i;
  if k<>0 then
  begin
    inc(a.l);
    a.a[a.l]:=K;
  end;
end;

procedure mul_i(var a:long;b:longint);
var k:longint;c,i:byte;
begin
  k:=0;
  i:=1;
  for i:=1 to a.l do
  begin
    c:=a.a[i];
    a.a[i]:=(c*b+k)mod 10;
    k:=(c*b+k)div 10;
  end;
  a.l:=i;
  while k<>0 do
  begin
    inc(a.l);
    a.a[a.l]:=K mod 10;
    k:=k div 10;
  end;
end;

procedure mul(var a,b:long);
var tmp,_a:long;
    i,j:byte;
begin
  fillchar(tmp,sizeof(tmp),0);
  for i:=1 to b.l do
  begin
    _a:=a;
    mul_i(_a,b.a[i]);
    for j:=_a.l downto 1 do _a.a[j+i-1]:=_a.a[j];
    for j:=1 to i-1 do _a.a[j]:=0;
    _a.l:=_a.l+i-1;
    add(tmp,_a);
  end;
{  i:=tmp.l;
  for i:=tmp.l downto 1 do if tmp.a[i]<> 0 then break;
  tmp.l:=i;}
  a:=tmp;
end;

var c:array[0..15]of extended;
    N,A,B,MMM:longint;
    i,j:integer;
    tmp:extended;
function cnk(N,K:integer):longint;
var s:longint;
    i:integer;
begin
  s:=1;
  for i:=K+1 to N do
  s:=s*i div (i-K);
  cnk:=s;
end;

var _a,_b:long;
begin
{  cnk(20,10);}
{  assign(input,'1114.dat');reset(input);}
  readln(N,A,B);
  if A>B then mmm:=A else mmm:=B;
  c[0]:=1;
  for i:=1 to mmm do
  begin
    {for j:=i-1 downto 0 do} c[i]:=c[i]+c[i-1];
    for j:=1 to i do
    begin
      c[i]:=c[i]+cnk(N,j)*cnk(i-1,j-1);
    end;
  end;
  tmp:=c[a];
  for i:=1 to 25 do
  begin
    _a.a[i]:=trunc(tmp-int(tmp/10)*10);
    tmp:=int(tmp/10);
  end;
  i:=25;
  for i:=25 downto 1 do if _a.a[i]<>0 then break;
  _a.l:=i;
  tmp:=c[b];
  for i:=1 to 25 do begin _b.a[i]:=trunc(tmp-int(tmp/10)*10);tmp:=int
(tmp/10);end;
  i:=25;
  for i:=25 downto 1 do if _b.a[i]<>0 then break;
  _b.l:=i;
  mul(_a,_b);
  for i:=_a.l downto 1 do write(_a.a[i]);
  writeln;
{  writeln(c[A]*c[B]:0:0);}
end.

All right!I think.
PROGRAM Ural1114;
Var
   S:Extended;
   I,J,K,P,Q,N,A,B:Byte;
   F:Array[1..20,0..15,0..15] Of Extended;
Begin
     ReadLN(N,A,B);
     FillChar(F,SizeOf(F),0);
     For I:=0 To A Do
         For J:=0 To B Do F[1,I,J]:=1;
     For I:=2 To N Do
         For J:=0 To A Do
             For K:=0 To B Do
             Begin
                  For P:=0 To J Do
                      For Q:=0 To K Do F[I,J,K]:=F[I,J,K]+F[I-1,P,Q];
             End;
     S:=0;
     For I:=0 To A Do
         For J:=0 To B Do
             S:=S+F[N,I,J];
     WriteLn(S:0:0);
End.

----------------------MY RROGRAM-----------------
type arr = record
             a:array[1..100]of byte;
             l:integer;
           end;
var N,a,b,i:longint;
    r,_a,_b:arr;
procedure mul_i(var a:arr;x:integer);
var c,k:longint;
    i:integer;
begin
  k:=0;
  for i:=1 to a.l do
  begin
    c:=a.a[i];
    a.a[i]:=(a.a[i]*x+k) mod 10;
    k:=(c*x+k) div 10;
  end;
  while k<>0 do
  begin
    inc(i);
    a.a[i]:=k mod 10;
    k:=k div 10;
  end;
  a.l:=i;
end;

procedure mul_i2(var a:arr;x:integer;var d:arr;p:integer);
var c,k:longint;
    i:integer;
begin
  k:=0;
  for i:=1 to a.l do
  begin
    c:=a.a[i];
    d.a[i+p]:=(a.a[i]*x+k) mod 10;
    k:=(c*x+k) div 10;
  end;
  while k<>0 do
  begin
    inc(i);
    d.a[i+p]:=k mod 10;
    k:=k div 10;
  end;
  d.l:=i+p;
end;

procedure add(var a,b:arr);
var m,i:integer;
    k,c:byte;
begin
  if a.l>b.l then m:=a.l else m:=b.l;
  k:=0;
  for i:=1 to m do
  begin
    c:=a.a[i];
    a.a[i]:=(a.a[i]+b.a[i]+k) mod 10;
    k:=(c+b.a[i]+k) div 10;
  end;
  if k = 1 then
  begin
    inc(i);
    a.a[i]:=1;
  end;
  a.l:=i;
end;

procedure mul(var a,b:arr);
var c,d:arr;
    i:integer;
begin
  d.l:=0;fillchar(d,sizeof(d),0);
  for i:=1 to a.l do
  begin
    c.l:=0;fillchar(c,sizeof(c),0);
    mul_i2(b,a.a[i],c,i-1);
    add(d,c);
  end;
  a:=d;
end;



procedure calc(a:longint;var s:arr);
var i:longint;
    x:arr;
begin
  fillchar(s,sizeof(s),0);
  fillchar(x,sizeof(x),0);
  s.l:=0;
  x.l:=1;
  x.a[1]:=1;
  for i:=0 to a do
  begin
    {s:=s+x;}add(s,x);
    {X:=X*N;}mul_i(x,N);
  end;
end;

begin
  readln(N,A,b);
  calc(a,_a);
  calc(B,_b);
  mul(_a,_b);
  for i:=_a.l downto 1 do write(_a.a[i]);
  writeln;
end.

Here is my solution:

#include <fstream.h>
#include <stdio.h>

int n,a,b;
long double x[21][16][16];
int main()
{
 cin>>n>>a>>b;

 {
  for(int ai=0; ai<=a; ai++)
   for(int bi=0; bi<=b; bi++)
    x[0][ai][bi] = 1;
 }

 for(int in=1; in<=n; in++)
  for(int ia=0; ia<=a; ia++)
   for(int ib=0; ib<=b; ib++)
 {
   long double xx = 0;
   for(int ai=0; ai<=ia; ai++)
    for(int bi=0; bi<=ib; bi++)
     xx += x[in-1][ia-ai][ib-bi];
   x[in][ia][ib] = xx;
 }

/* cout.setf(ios::fixed);
 cout.precision(0);
 cout<<x[n][a][b]<<endl;*/
 printf("%.0Lf\n",x[n][a][b]);
 return 0;
}

There is my solution:

Program t1114;

Const MaxL=15;

Var Mas         :array[1..20,0..MaxL,0..MaxL]of comp;
    Ans         :comp;
    N,i,j,A,B,k :longint;
    u,v         :longint;

begin
Read(N,A,B);
for i:=1 to 20 do
 for j:=1 to MaxL do
  for k:=1 to MaxL do
   mas[i,j,k]:=0;
if A>0 then mas[1,A-1,B]:=1;
if B>0 then mas[1,A,B-1]:=1;
if (A>0)and(B>0) then mas[1,A-1,B-1]:=1;
mas[1,A,B]:=1;
for i:=2 to N do begin
 for j:=0 to MaxL do
  for k:=0 to MaxL do
   if mas[i-1,j,k]<>0 then begin
    for u:=0 to j do
     for v:=0 to k do
      mas[i,j-u,k-v]:=mas[i,j-u,k-v]+mas[i-1,j,k];
   end;
end;
Ans:=0;
for i:=0 to MaxL do
 for j:=0 to MaxL do
  Ans:=Ans+mas[N,i,j];
writeln(Ans:0:0);
end.