/*-----------------------------------------------------------*
 *
 * File : fonc.c
 * Author : NTM
 * Created : 30/09/03
 *
 *
 * Copyright (C) Nicolas Thierry-Mieg, 2006.
 *
 *
 * This file is part of InterPool, written by 
 * Nicolas Thierry-Mieg (CNRS, France) Nicolas.Thierry-Mieg@imag.fr
 *
 * InterPool is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * InterPool is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with InterPool; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *------------------------------------------------------------*/

#include <stdio.h> // printf
#include <stdlib.h> // exit
#include <string.h> // memset

#include "fonc.h"

/************************************************************************
 ******************* EXPORTED FUNCTIONS *********************************
 ************************************************************************/

/////////////////////////////////////////////////////////////////////
// firstDivisor: return the smallest divisor of nn, or 0 if nn is prime.
// given by JTM
int firstDivisor(int nn)
{
#define NN 10000              /* max size of prime table */
  static int primeList[NN] ; /* lazy table of known primes */
  static int maxKnown = 0 ;
  int i, p, q ;               /* running integers */
  int testing = 0 ;           /* 1: show prime list when new number computed */

  if (!maxKnown) /* initialise */
    {
      memset (primeList, 0, NN * sizeof (int)) ;
      maxKnown = 1 ;
      primeList[0] = 2 ;
    }
  while (1) /* lazy fill of primeList */
    /* NB if p is our largest known prime, there
       is always another prime q < p*p
    */
    {
      p = primeList[maxKnown - 1] ;
      if (p*p > nn) 
	break ;
      for (q = p + 1 ; firstDivisor (q) ; q++)
	;
      if (testing)
	printf ("%6d is prime\n", q) ;
      primeList [maxKnown++] = q ;
      if (maxKnown == NN)
	{ 
	  fprintf (stderr,"Overflow in firstDivisor nn=%d, maxKnown=%d q = %d\n", nn, maxKnown, q) ;
	  exit (1) ;
	}
    }
  /* testing */
  if (nn<0)
    { // nn is negative, should not happen
      fprintf(stderr, "in firstDivisor, nn is negative (value %d)\n", nn) ;
      exit(1) ;
    }
  for (i = 0 ; i < maxKnown ; i++)
    {
      p = primeList[i] ;
      if (p*p > nn)
	return 0 ; /* nn is prime, return 0 */
      if ( nn == p * (nn/p))
	return p ; /* not a prime */
    }
  return 0 ; /* nn is prime, return 0 */
}


///////////////////////////////////////////////////////////////
// compInts: compare 2 ints (for qsort)
int compInts(const void* a, const void* b)
{
  return (*(int*)a - *(int*)b) ;
}


///////////////////////////////////////////////////////////////
// choose:
// calculate C(t,k) = k!/(t!*(k-t)!)
unsigned long int choose(int t, int k)
{
  if (t>k)
    return 0 ;

  // C(t,k) = C(k-t,k), we want t as small as possible
  if (t > (k-t))
    t = k-t ;

  {
    unsigned long int result = 1 ;
    unsigned long int divideby = 1 ;

    while (t > 0)
      {
	result *= k-- ;
	divideby *= t-- ;
      }

    result /= divideby ;
    return result ;
  }
}


