/* svc: A fast C implementation of the membership function for the
   Smith-Volterra-Cantor set (also known as the "fat Cantor set"), a
   nowhere-dense fractal subset of [0,1] that has Lebesgue measure 1/2.
   See http://en.wikipedia.org/wiki/Smith-Volterra-Cantor_set .

   Of course, we can't really work with real numbers, but this code
   provides a discrete equivalent, e.g., a subset of the same shape
   that contains half of the integers in an interval [0,max).  The
   runtime of a membership test is O(log max).  

   On my Thinkpad T61 laptop, unsigned longs are 32 bits, so this code
   works up to max=2**31, which can be expressed as ((ulong)1)<<31.
   The function can individually check all 2 billion points in this
   range in less than a minute -- about 24 ns per point.

   (I wrote this because I had a possible construction that involved
   SVC membership tests on random real numbers, and I wanted to do a
   sanity check on its desired properties before trying to prove
   anything.)

   Compile as: gcc -O3 svc.c -o svc

   Author: Jason Eisner <jason@cs.jhu.edu>, October 2007
   I hereby place this code snippet in the public domain; 
   acknowledgment is welcome but not required.
*/

#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

typedef unsigned long int ulong;
int svc1(ulong i, ulong max, ulong del);

/* This is the function to call: a membership test for the 
   discrete Smith-Volterra-Cantor set over [0,max).  
   Returns true iff i is in this set.

   The caller must ensure that 0 <= i < max, 
   and that max is a power of 2 and at least 8. 
   We return true for exactly half the points i in this range. 

   Note: A more careful implementation that paid attention to rounding
   could probably relax these restrictions on max, requiring only that
   max is even.  I think it may already turn out to work for max=4 and
   max=2.
*/

int svc(ulong i, ulong max) {
  return svc1(i, max, max/4);
}


/* recursive helper function that does the real work.
   If we're in the middle del points of [0,max), we return false 
   since these are to be deleted.  Otherwise we have to recurse
   with del <- del/4.  The tricky part is to get the base case right
   to ensure that we delete a total of half the points overall.
*/
int svc1(ulong i, ulong max, ulong del) {
  if (del==0)  /* Really del=0.5, but we rounded down.  
		  This is the base case if the original max had the form 2^(2k+1).
		  We won't be able to recurse further, but if we were doing real
		  numbers, we would delete a total of 1 from here on: 
		  0.5 + 2*(0.5/4) + 4*(0.5/16) + ...
		  So delete the middle 1 (max will be odd in this case) and stop. */
    return (i!=max/2);
  else if (del==1)  /* This is the base case if the original max had the form 2^(2k). 
		       We won't be able to recurse further, but if we were doing
		       real numbers, we would delete a total of 2 from here on.
		       So delete the middle 2 (max will be even in this case) and stop. */
    return (i!=max/2-1 && i!=max/2);
  else {   /* Delete the middle del points and then recurse. */
    ulong midlo = (max-del)/2;   /* smallest point to delete on this round */
    if (i >= midlo)
      if (i < midlo+del)
	return 0;   /* we're in the middle del points that we're deleting on this round */
      else 
	i=(max-1)-i;    /* we're in the upper half of [0,max); exploit symmetry and flip to the lower half */
    return svc1(i, midlo, del/4);   /* recurse */
  }
}


/* A little code to test the svc() function. */

int main(void) {
  ulong i, max, count;

  /* Graphically show our discrete SVC set for some small ranges,
     and just summary statistics for larger ranges. */
  for (max=8; max > 0; max*=2) {
    count = 0;
    for (i=0; i < max; i++) {
      int result = svc(i,max);
      count += result;
      if (max < 80)     
	putchar(result ? '*' : ' '); 
    }
    assert(count == max/2);   /* check that we're hitting exactly half */
    printf("| (%lu/%lu)\n", count, max);
  }

  return 0;
}