#include <stdio.h>
#include <math.h>

/* 
Demonstration of control structures and loops in C.

Compiles on a UNIX system using 

  gcc -lm -Wall -o C-examples-2 C-examples-2.c

Execute from the shell using "C-examples-2" or "./C-examples-2".
*/



/*  
Let C be the circle of radius R with center located at coordinates
(X, Y).  This function returns -1 if C is completely contained
within the unit circle, returns 1 if C is completely outside of
the unit circle, and returns 0 otherwise.
*/
int circle(const double R,
	   const double X,
	   const double Y)
{
  /* The distance from the center of C to the origin. */
  double L = sqrt(X*X + Y*Y);

  /* The greatest distance from C to the origin. */
  double L1 = L + R;

  /* The least distance from C to the origin. */
  double L2 = L - R;

  /* Determine which of the three situations we are in. */
  if (L1 < 1) return -1;
  else if (L2 > 1) return 1;
  else return 0;
}



/* 
Count the number of quadrants that intersect with the circle C having
radius R and center at coordinates (X, Y).
*/
int count_quadrants(const double R,
		    const double X,
		    const double Y)
{
  /* Start at 1 to account for the quadrant where the center lies. */
  int N = 1;
  
  /* The distance from the center of C to the origin. */
  double L = sqrt(X*X + Y*Y);

  /* Check if C intersects the adjacent quadrant in the X direction. */
  if (R >= fabs(X)) ++N;

  /* Check if C intersects the adjacent quadrant in the Y direction. */
  if (R >= fabs(Y)) ++N;

  /* Check if C intersects the quadrant diagonal to the center. */
  if (R >= L) ++N;

  return N;
}


/* 
Evaluates a step function at X, where the step function is defined to
return A when X is less than or equal to C and B when X is greater
than C.
*/
double stepfun(const double X, 
	       const double A, 
	       const double B, 
	       const double C)
{
  return (X <= C) ? A : B;
}


/* 
Returns 1 if b is strictly beetween a and c, 0 otherwise.
*/
int is_between(const double a,
	       const double b,
	       const double c)
{
  return ( ((b > a) && (b < c)) || ((b < a) && (b > c)) ) ? 1 : 0;
}




/* 
Approximate e using the power series expansion: 
1/0! + 1/1! + 1/2! + 1/3! + ... + 1/K!.
*/
double approximate_e(const int K)
{
  int i;

  /* Start at 1/0! = 1. */
  double E = 1;

  /* This will hold the value of i! for i=1, 2, ... */
  double F = 1;

  for (i=1; i<=K; ++i)
    {
      F *= i;
      E += 1 / F;
    }

  return E;
}


/* 
A crude way to determine if an integer is prime.
*/
int is_prime(const int P)
{
  int i;

  for (i=2; i<=sqrt(P); ++i)
    {
      if (P % i == 0) return 0;
    }

  return 1;
}


/* 
Returns 1 if P is an odd prime, 0 otherwise.
*/
int is_odd_prime(const int P)
{
  return ( (P != 2) && is_prime(P) );
}


/* 
Factor A as P*Q, where P and Q are any pair of nontrivial factors of
A.  The function returns one if such a factorization exists, and zero
otherwise.
*/
int factor_pair(const int A,
		int* P,
		int* Q)
{
  for (*P=2; *P<=sqrt(A); ++*P)
    {
      if (A % *P == 0)
	{
	  *Q = A / *P;
	  return 1;
	}
    }

  return 0;
}


/* 
A recursive function that returns the number of terms in the prime
factorization of A.
*/
int count_factors(const int A)
{
  int P, Q;
  int D = factor_pair(A, &P, &Q);

  if (D == 0) return 1;
  else return count_factors(P) + count_factors(Q);
}


/* 
Returns the number of terms of the harmonic series 1/1 + 1/2 + 1/3 +
 ... (which diverges to positive infinity) are needed to exceed Z.
*/
int harmonic_series(const double Z)
{
  int k = 0;
  double S = 0;

  while (S < Z)
    {
      ++k;
      S += 1 / (double)k;
    }

  return k;
}


/* 
Some people like to write "obfuscated" code like this...
*/
int harmonic_series1(const double Z)
{
  int k = 1;
  double S = 0;

  while ( (S += 1 / (double)k) < Z) ++k;

  return k;
}


/* 
or even like this.
*/
int harmonic_series2(const double Z)
{
  int k = 0;
  double S = 0;

  while ( (S += 1 / (double)(++k)) < Z);

  return k;
}



int main(int argc, 
	 char** argv)
{
  int k;
  
  printf("Check some prime numbers:\n");
  for (k=3; k<20; ++k)
    {
      printf("%d is%sprime\n", k, (is_prime(k) == 1) ? " " : " not ");
    }

  printf("\nPower series approximations to e:\n");
  for (k=1; k<18; ++k)
    {
      printf("approximate_E(%2d) = %16.13f\n", k, approximate_e(k));
    }

  printf("exp(1) = %16.13f\n\n", exp(1));

  printf("\nCount factors:\n");
  for (k=2; k<30; ++k)
    {
      int d = count_factors(k);
      printf("%d has %d factor%s\n", k, d, (d == 1) ? "" : "s");
    }

  printf("\nNumber of terms in a truncated harmonic series that"
	 " exceed a given value:\n");
  for (k=1; k<19; ++k)
    {
      printf("harmonic_series(%d) = %d\n", k, harmonic_series(k));
    }

  return 0;
}