# TMath::Binomial distribution

    // tutorial illustrating the use of TMath::Binomial
    //  can be run with:
    // root > .x binomial.C
    // root > .x binomial.C+ with ACLIC

    #include <TMath.h>
    #include <TRandom.h>

    void binomialSimple() {
      //
      // Simple test for the binomial distribution
      //
      printf("nTMath::Binomial simple testn");
      printf("Build the Tartaglia trianglen");
      printf("============================n");
      const Int_t max=13;
      Int_t j;
      for(Int_t i=0;i<max;i++) {
        printf("n=%2d",i);
        for(j=0;j<(max-i);j++) printf("  ");
        for(j=0;j<i+1;j++) printf("%4d",TMath::Nint(TMath::Binomial(i,j)));
        printf("n");
      }
    }

    void binomialFancy() {
      Double_t x;
      Double_t y;
      Double_t res1;
      Double_t res2;
      Double_t err;
      Double_t serr=0;
      const Int_t nmax=10000;
      printf("nTMath::Binomial fancy testn");
      printf("Verify Newton formula for (x+y)^nn");
      printf("x,y in [-2,2] and n from 0 to 9  n");
      printf("=================================n");
      TRandom r;
      for(Int_t i=0; i<nmax; i++) {
        do {
            x=2*(1-2*r.Rndm());
            y=2*(1-2*r.Rndm());
        } while (TMath::Abs(x+y)<0.75); //Avoid large cancellations
        for(Int_t j=0; j<10; j++) {
           res1=TMath::Power(x+y,j);
           res2=0;
           for(Int_t k=0; k<=j; k++)
              res2+=TMath::Power(x,k)*TMath::Power(y,j-k)*TMath::Binomial(j,k);
           if((err=TMath::Abs(res1-res2)/TMath::Abs(res1))>1e-10)
     	      printf("res1=%e res2=%e x=%e y=%e err=%e j=%dn",res1,res2,x,y,err,j);
           serr +=err;
         }
      }
      printf("Average Error = %en",serr/nmax);
    }

    void binomial () {
       binomialSimple();
       binomialFancy();
    }