/*========================================================================
                            ** BuDDy **
                      BINARY DECISION DIAGRAMS
                               with
                  FINITE DOMAIN VARIABLE INTERFACE
                      Library Package ver. 1.9
                        By Jorn Lind-Nielsen

            Copyright (C) 1996-2000 by Jorn Lind-Nielsen

    Permission is hereby granted to reproduce and distribute this
    package by any means and for any fee, whether alone or as part
    of a larger distribution, in source or in binary form, provided
    this notice is included with any such distribution and is visible 
    for the end user, and is not removed from any of its header files. 

      *** I AM NOT RESPONSIBLE FOR ANY KIND OF DAMAGE TO YOUR  ***
      *** FILES, DATA, HARDWARE, LOSS OF MONEY, SYSTEM CRASHES *** 
      *** OR ANY OTHER THING YOU MIGHT COME UP WITH.           ***
      *** - USE THIS PROGRAM OF YOUR OWN FREE WILL !!!         ***

      Happy Hacking
                   Jorn Lind-Nielsen

========================================================================*/

/* 
 *  Copyright (c) 1996 Technical University of Denmark
 *
 * by Henrik Hulgaard
 * Department of Information Technology
 * Building 344
 * Technical University of Denmark
 * DK-2800 Lyngby, Denmark
 * 
 * e-mail: henrik@it.dtu.dk
 *
 * Small mods by Jorn Lind feb 2000
 */

#include <stdlib.h>
#include "prime.h"

typedef int Boolean;

#ifndef false
#define false 0
#endif

#ifndef true
#define true (!false)
#endif



#define WORDSIZE	32		/* Change to 64 on architectures where 
					 * an unsigned int fits in 64 bits.
					 */

#define ODD(n)		((n)&1)
#define EVEN(n)		(!ODD(n))
#define TESTBIT(x,n)	((x) & (1<<(n)))

#define HI(x)	( ((x) & 0xffff0000) >> 16 )
#define LO(x)	  ((x) & 0x0000ffff)

#ifndef min
#define min(a,b)	((a) < (b) ? (a) : (b))
#endif

#ifndef max
#define max(a,b)	((a) > (b) ? (a) : (b))
#endif

#define lrandom(l,u) ( (rand() % ((u)-(l)+1)) + (l) )
#define select_bits(x1,x0,i) ((x0) >> (i)) + ((i) ? ((x1) << (32-(i))) : 0)
#define is_factor(x,n) ( (((x)!=(n)) && ((x)%(n) == 0)) )


static unsigned int number_bits( unsigned int x );
static void subtract( unsigned int x1, unsigned int x0, unsigned int yn,
		      unsigned int i, unsigned int *X1, unsigned int *X0 );


static void dp_mul( unsigned a, unsigned b, unsigned *x1, unsigned *x0 )
{
  unsigned int a_hi = HI(a);
  unsigned int a_lo = LO(a);

  unsigned int b_hi = HI(b);
  unsigned int b_lo = LO(b);
  
  unsigned int ab_hh = a_hi * b_hi;
  unsigned int ab_hl = a_hi * b_lo;
  unsigned int ab_lh = a_lo * b_hi;
  unsigned int ab_ll = a_lo * b_lo;

  unsigned int ab_hl_h = (ab_hl >> 16);
  unsigned int ab_hl_l = ((ab_hl & 0x0000ffff) << 16);

  unsigned int ab_lh_h = (ab_lh >> 16);
  unsigned int ab_lh_l = ((ab_lh & 0x0000ffff) << 16);

  unsigned int carry = 0;

  /* a + b produces a carry iff (a+b) < max(a,b)  */
  unsigned int tmp = ab_hl_l + ab_lh_l; 
  carry += (tmp < max( ab_hl_l, ab_lh_l ) ? 1 : 0);

  *x0 = tmp + ab_ll;
  carry += ( *x0 < max( tmp, ab_ll ) ? 1 : 0);
 
  *x1  = ab_hh + ab_hl_h + ab_lh_h + carry;
}


static void random_init(unsigned int seed)
{
   srand(seed);
}

static unsigned int number_bits( unsigned int x )
{
  int k = (sizeof( unsigned int ) * 8) - 1;

  while ( k && !TESTBIT( x, k ))
    k--;

  if ( k == 0 ) 
    return (TESTBIT( x, k ) ? 1 : 0);
  else
    return k+1;
}

static void subtract( unsigned int x1, unsigned int x0, unsigned int yn,
		      unsigned int i, unsigned int *X1, unsigned int *X0 )
{
  /* subtract (n << i) from (x1,x0) */
  unsigned int x = select_bits( x1, x0, i ) - yn;

  *X1  = i ? (x >> (32-i)) : 0;
  *X0  = (x << i) + (x0 & ((1 << i) - 1));		/* 0b0000111..11 */
}

static unsigned int reminder( unsigned int x1, unsigned int x0, unsigned int n )
{
  int i;

  for ( i = number_bits( x1 ); i >= 0; i-- ) {
    unsigned int k = select_bits( x1, x0, i );	
    unsigned int y = k/n;    

    subtract( x1, x0, y*n, i, &x1, &x0 );
  }    

  return x0;
}

#if WORDSIZE==32
static unsigned int mul_mod( unsigned int a, unsigned int b, unsigned int n )
{
  unsigned int x0, x1;
  dp_mul( a, b, &x1, &x0 );

  return reminder( x1, x0, n );
}
#else
static inline unsigned int mul_mod( unsigned int a, unsigned int b, unsigned int n )
{
  return (a*b) % n;
}
#endif

static Boolean witness( unsigned int a, unsigned int n )
{
  unsigned int b = n - 1;
  unsigned int d = 1;
  unsigned int k = number_bits( b ) - 1;

  int i;
  
  for ( i = k; i >= 0; i-- ) {
    unsigned int x = d;

    d = mul_mod( d, d, n );

    if ( (d == 1) && (x != 1) && (x != b))
      return true;

    if ( TESTBIT( b, i ) )
      d = mul_mod( d, a, n );
  }

  return ( d != 1 );
}


static Boolean miller_rabin_check_prime( unsigned int x, int s )
{
  int j;

  for ( j = 0; j < s; j++ ) {
    unsigned int a = lrandom( 1, x-1 );

    if ( witness( a, x ))
      return false;
  }
  return true;
}


/* fast_check_composite --
 * 	checks if 3, 5, 7, 11 or 13 divides n.
 */

static Boolean fast_check_composite( unsigned int  x )
{
  return ( is_factor( x, 3 ) ||
	   is_factor( x, 5 ) ||
	   is_factor( x, 7 ) ||
	   is_factor( x, 11 ) ||
	   is_factor( x, 13 ));
}

/* is_prime -- return whether x is prime. 
 */

static Boolean is_prime( unsigned int x, unsigned int s )
{
  if ( fast_check_composite( x ) )
    return false;

  return miller_rabin_check_prime( x, s );
}


/* prime_ge -- returns the next prime after x.
 */

unsigned int bdd_prime_gte( unsigned int x )
{
  random_init( 1 );
  
  if ( EVEN(x) )	  /* make sure x is odd. */
    x++;

  while ( !is_prime( x, 20 ))
    x += 2;

  return x; 
}


/* prime_ge -- returns the next prime after x.
 */

unsigned int bdd_prime_lte( unsigned int x )
{
  random_init( 1 );
  
  if ( EVEN(x) )	  /* make sure x is odd. */
    x--;

  while ( !is_prime( x, 20 ))
    x -= 2;

  return x; 
}


#if 0
void main( int argc, char *argv[] )
{ 
  if ( argc > 1 ) {
    unsigned int x = atoi( argv[1] );
    printf( "prime >= %u: %u\n", x, bdd_prime_gte( x ));
    printf( "prime <= %u: %u\n", x, bdd_prime_lte( x ));
  }
}
#endif
