/* discontinuous.c  two functions to test finding minimum code
 *                  f(x)= (x^2-1) if bottom bit of x=0
 *                      = 1.0     if bottom bit of x=1   -1 < x < 1
 *                        64 bit IEEE floating point
 *                        this is numerically everywhere discontimuous
 *
 *                  a recursive function that is continuous,
 *                  yet nowhere analytically differentiable 0 < x < 1
 *                  f(x) =
 *                  if(x<=1.0/3.0) return (3.0/4.0)*f(3.0*x)
 *             else if(x>=2.0/3.0) return (1.0/4.0)+(3.0/4.0)*f(3.0*x-2.0)
 *                else             return (1.0/4.0)+(1.0/2.0)*f(2.0-3.0*x)
 */

#include <stdio.h>
#include <stdlib.h>
#include "udrnrt.h"

double fd(double x) /* for 0 < x < 1 */
{
  int debug = 0;
  static double y; /*  retains value between calls */
  static int n = 30;
  if(n<=0) {n=30; return y;}
  n--;
  if(debug) printf("fd x=%20.12f \n", x);
  if(x<=1.0/3.0) {y=(3.0/4.0)*fd(3.0*x);
                  if(debug) printf("fd1/3 x=%20.12f, y=%20.12f \n", x, y);
                  return y;}
  else if(x>=2.0/3.0) {y=(1.0/4.0)+(3.0/4.0)*fd(3.0*x-2.0);
                       if(debug) printf("fd2/3 x=%20.12f, y=%20.12f \n", x, y);
                       return y;}
  else /* 1/3 < x < 2/3 */ {y=(1.0/4.0)+(1.0/2.0)*fd(2.0-3.0*x);
                       if(debug) printf("fdmid x=%20.12f, y=%20.12f \n", x, y);
                            return y;}
}

double f(double x) /* for -1 < x < 1 */
{
  union xix{double x; long ix;} overlay;
  overlay.x = x;
  if(overlay.ix & 1L) return 1.0;
  return (x*x-1.0);
}

int main(int argc, char *argv[])
{
  double x, y, dx;
  FILE * fout, * fdout;
  int n = 100;
  int i;

  printf("discontinuous.c running \n");
  fout = fopen("fdiscontinuous.out","w");
  fdout = fopen("fddiscontinuous.out","w");

  printf("a few points of y = f(x) \n");
  x = 0.0000001;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.000001;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.1;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.3;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.5;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.7;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.9;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.999999;
  y = f(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);

  printf("a few points of y = fd(x) \n");
  x = 0.0000001;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.000001;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.01;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (1.0/3.0)- 0.01;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (1.0/3.0)+ 0.01;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (1.0/3.0)+ 0.02;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (1.0/3.0)+ 0.03;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.5;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.50000001;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (2.0/3.0)- 0.001;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = (2.0/3.0)+ 0.001;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 0.99999999;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);
  x = 1.0;
  y = fd(x);
  printf("x=%20.12f, = %LX, y=%20.12f \n", x, x, y);

  printf("write fdiscontinuous.out, %d points \n", 2*n+1);
  for(i=0; i<=2*n; i++)
  {
    x = 2.0*udrnrt()-1; /*random -1 to 1 */
    y = f(x);
    fprintf(fout, "%f %f\n", x, y);
  }
  fclose(fout);

  printf("write fddiscontinuous.out, %d points \n", n+1);
  dx = 0.01;
  x = 0.0;
  for(i=0; i<=n; i++)
  {
    y = fd(x);
    fprintf(fdout, "%f %f\n", x, y);
    x = x+dx;
  }
  fclose(fdout);

  printf("discontinuous.c finished \n");
  return 0;
} /* end discontinuous.c */