/* Verify inequality associated with quasi-random graphs
   and complete minors.  */
/* Copyright 2002 Joseph Samuel Myers.
   All rights reserved.
  
   Redistribution and use in source and binary forms,
   with or without modification, are permitted provided that
   the following conditions are met:
   1. Redistributions of source code must retain the above
      copyright notice, this list of conditions and the
      following disclaimer.
   2. Redistributions in binary form must reproduce the above
      copyright notice, this list of conditions and the
      following disclaimer in the documentation and/or
      other materials provided with the distribution.
   3. The name of the author may not be used to endorse
      or promote products derived from this software without
      specific prior written permission.
  
   THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND
   ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING,
   BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
   MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
   ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE
   FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,
   OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
   DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
   AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
   STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
   OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

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

static long double q0; /* q_0 = 0.543689012692...,
                          real root of
                          q^3 + q^2 + q - 1 = 0.  */
static long double p0; /* p_0 = 1 - q_0,
                          real root of
                          p^3 - 4p^2 + 6p - 2 = 0.  */
static long double em1; /* e^-1.  */
static long double em2; /* e^-2.  */
static long double em2x4; /* 4e^-2.  */

typedef struct {
  long double min;
  long double max;
} BOUNDS;

/* Print an error and exit.  */
static void
die(const char *format, ...)
{
  va_list args;
  fprintf(stderr, "minors-quasi-ineq: ");
  va_start(args, format);
  vfprintf(stderr, format, args);
  va_end(args);
  exit(EXIT_FAILURE);
}










/* Return -r log r.  */
static long double
mrlogr (long double r)
{
  if (r < 0.0L || r > 1.0L)
    die("mrlogr out of range");
  if (r == 0.0L)
    return 0.0L;
  return -r * logl(r);
}

/* Return r (log r)^2.  */
static long double
rlogr2 (long double r)
{
  long double l;
  if (r < 0.0L || r > 1.0L)
    die("rlogr2 out of range");
  if (r == 0.0L)
    return 0.0L;
  l = logl(r);
  return r * l * l;
}

/* Return x^(1/y^2).  */
static long double
powov2 (long double x, long double y)
{
  if (x < 0.0L || x >= 1.0L || y < 0.0L || y > 1.0L)
    die("powov2 out of range");
  if (y == 0.0L || x == 0.0L)
    return 0.0L;
  return powl(x, 1.0L / y / y);
}







/* Return the bounds of -r log r.  */
static BOUNDS
mlog_bounds (BOUNDS r)
{
  BOUNDS ret;
  long double left;
  long double right;
  left = mrlogr(r.min);
  right = mrlogr(r.max);
  ret.min = (left < right ? left : right);
  if (r.min <= em1 && em1 <= r.max)
    ret.max = em1;
  else
    ret.max = (left > right ? left : right);
  return ret;
}

/* Return the bounds of r (log r)^2.  */
static BOUNDS
log2_bounds (BOUNDS r)
{
  BOUNDS ret;
  long double left;
  long double right;
  left = rlogr2(r.min);
  right = rlogr2(r.max);
  ret.min = (left < right ? left : right);
  if (r.min <= em2 && em2 <= r.max)
    ret.max = em2x4;
  else
    ret.max = (left > right ? left : right);
  return ret;
}








/* Return the bounds of d^2f/d\alpha^2.  */
static BOUNDS
d2f_bounds (BOUNDS qb, BOUNDS ab)
{
  BOUNDS qoa; /* Bounds of q^(1/\alpha^2).  */
  BOUNDS qoal1; /* ... * - log self.  */
  BOUNDS qoal2; /* ... * log^2 self.  */
  BOUNDS qo1a; /* Bounds of q^(1/(1-\alpha)^2).  */
  BOUNDS qo1al1; /* ... * - log self.  */
  BOUNDS qo1al2; /* ... * log^2 self.  */
  BOUNDS ret;
  qoa.min = powov2(qb.min, ab.min);
  qoa.max = powov2(qb.max, ab.max);
  qo1a.min = powov2(qb.min, 1.0L - ab.max);
  qo1a.max = powov2(qb.max, 1.0L - ab.min);
  qoal1 = mlog_bounds(qoa);
  qoal2 = log2_bounds(qoa);
  qo1al1 = mlog_bounds(qo1a);
  qo1al2 = log2_bounds(qo1a);
  ret.min = (-4.0L
             + 2.0L * qoa.min + 2.0L * qoal1.min
             + 4.0L * qoal2.min
             + 2.0L * qo1a.min + 2.0L * qo1al1.min
             + 4.0L * qo1al2.min);
  ret.max = (-4.0L
             + 2.0L * qoa.max + 2.0L * qoal1.max
             + 4.0L * qoal2.max
             + 2.0L * qo1a.max + 2.0L * qo1al1.max
             + 4.0L * qo1al2.max);
  return ret;
}










/* Print a fraction in reduced form.
   The only factors to remove are powers of 2.  */
static void
print_reduced_fraction(int n, int d)
{
  while (n % 2 == 0 && d % 2 == 0) {
    n /= 2;
    d /= 2;
  }
  if (n == 0)
    printf("$0$");
  else
    printf("$\\frac{%d}{%d}$", n, d);
}

/* Attempt to prove second derivative always negative
   for given q.  */
static void
prove_d2f_neg (BOUNDS qb)
{
  long double amin, astep;
  amin = 0.0L;
  astep = 0.5L;
  while (amin < 0.5L) {
    BOUNDS ab;
    BOUNDS tb;
    ab.min = amin;
    ab.max = amin + astep;
    tb = d2f_bounds(qb, ab);
    if (tb.max < 0.0L) {
      printf("$%Lg$ & $%Lg$ & ",
             qb.min, qb.max);
      print_reduced_fraction(ab.min / astep, 1.0L / astep);
      printf(" & ");
      print_reduced_fraction(ab.max / astep, 1.0L / astep);
      printf(" & $%.6Lf$ & $%.6Lf$ \\\\\n",
             tb.min, tb.max);
      amin += astep;
    } else {


      astep /= 2;
      if (astep < 1.0L / 4096.0L)
        break;
    }
  }
  if (amin >= 0.5L) {
    printf("Succeeded in proving second derivative "
           "negative for all alpha,\n"
           "%Lf <= q <= %Lf\n", qb.min, qb.max);
  } else {
    die("FAILED to prove second derivative "
        "negative for all alpha,\n"
        "%Lf <= q <= %Lf", qb.min, qb.max);
  }
}


/* Attempt to prove f always positive for given q.  */
static void
prove_f_pos (BOUNDS qb)
{
  /* Divide [ 0, 1/2 ] into some number of parts.  Bound
     second derivative on each part.  We know first
     derivative is zero at centre; bound it on each part.
     We know f is zero at alpha = 0; bound it at end of
     each part.  Want: f > 0 in centre.  Use this from say
     q = 0.4 up, but not too far up (too near q0).  */
  int parts = 2;
  BOUNDS d2b[32];
  BOUNDS d1b[32]; /* 1st derivative on a region.  */
  BOUNDS d1bleft[32]; /* 1st derivative to
                         the left of a region.  */
  BOUNDS fb[32], fbright[32];
  while (parts <= 32) {
    int i;
    long double d = 0.5L / (long double)parts;
    BOUNDS rb, lb;
    for (i = 0; i < parts; i++) {
      BOUNDS ab;
      ab.min = d * i;
      ab.max = d * (i + 1);
      d2b[i] = d2f_bounds(qb, ab);
    }
    rb.min = 0.0L;
    rb.max = 0.0L;
    for (i = parts - 1; i >= 0; i--) {
      d1bleft[i].min = rb.min - d * d2b[i].max;
      d1bleft[i].max = rb.max - d * d2b[i].min;
      d1b[i].min = (rb.min < d1bleft[i].min
                    ? rb.min
                    : d1bleft[i].min);
      d1b[i].max = (rb.max > d1bleft[i].max
                    ? rb.max
                    : d1bleft[i].max);
      rb = d1bleft[i];
    }
    lb.min = 0.0L;
    lb.max = 0.0L;
    for (i = 0; i < parts; i++) {
      fbright[i].min = lb.min + d * d1b[i].min;
      fbright[i].max = lb.max + d * d1b[i].max;
      fb[i].min = (lb.min < fbright[i].min
                   ? lb.min
                   : fbright[i].min);
      fb[i].max = (lb.max > fbright[i].max
                   ? lb.max
                   : fbright[i].max);
      if (fbright[i].min <= 0.0L)
        break;
      lb = fbright[i];
    }
    if (i == parts)
      break;
    parts *= 2;
  }
  if (parts <= 32) {
    int i;
    for (i = 0; i < parts; i++) {
      printf("$%Lg$ & $%Lg$ "
             "& $\\frac{%d}{%d}$ & $\\frac{%d}{%d}$ "
             "& $%.6Lf$ & $%.6Lf$ & $%.6Lf$ & $%.6Lf$ "
             "& $%.6Lf$ & $%.6Lf$ \\\\\n",
             qb.min, qb.max, i, 2 * parts, i + 1, 2 * parts,
             d2b[i].min, d2b[i].max,
             d1bleft[i].min, d1bleft[i].max,
             fbright[i].min, fbright[i].max);
    }
    printf("Proved f always positive for "
           "%Lf <= q <= %Lf, %d steps.\n",
           qb.min, qb.max, parts);
  } else {
    die("FAILED to prove f always positive for "
        "%Lf <= q <= %Lf.",
        qb.min, qb.max);
  }
}

/* Attempt to prove f of a certain shape for given q.  */
static void
prove_f_shape (BOUNDS qb)
{
  /* Divide [ 0, 1/2 ] into some number of parts.  Bound
     second derivative on each part.  We know first
     derivative is zero at centre; bound it on each part.
     Want: first a region with second derivative negative,
     then the first derivative negative until the centre
     (zero at the centre - so then allow a region with
     second derivative positive).  */
  int parts = 2;
  BOUNDS d2b[32];
  BOUNDS d1b[32]; /* 1st derivative on a region.  */
  BOUNDS d1bleft[32]; /* 1st derivative to
                         the left of a region.  */
  while (parts <= 32) {
    int i;
    long double d = 0.5L / (long double)parts;
    BOUNDS rb;
    for (i = 0; i < parts; i++) {
      BOUNDS ab;
      ab.min = d * i;
      ab.max = d * (i + 1);
      d2b[i] = d2f_bounds(qb, ab);
    }
    rb.min = 0.0L;
    rb.max = 0.0L;
    for (i = parts - 1; i >= 0; i--) {
      d1bleft[i].min = rb.min - d * d2b[i].max;
      d1bleft[i].max = rb.max - d * d2b[i].min;
      d1b[i].min = (rb.min < d1bleft[i].min
                    ? rb.min
                    : d1bleft[i].min);
      d1b[i].max = (rb.max > d1bleft[i].max
                    ? rb.max
                    : d1bleft[i].max);
      rb = d1bleft[i];
    }
    for (i = 0; i < parts; i++) {
      if (d2b[i].max >= 0)
        break;
    }
    for (; i < parts; i++) {
      if (d1b[i].max >= 0)
        break;
    }
    for (; i < parts; i++) {
      if (d2b[i].min <= 0)
        break;
    }
    if (i == parts)
      break;
    parts *= 2;
  }
  if (parts <= 32) {
    int i;
    for (i = 0; i < parts; i++) {
      printf("$%Lg$ & $%Lg$ "
             "& $\\frac{%d}{%d}$ & $\\frac{%d}{%d}$ "
             "& $%.6Lf$ & $%.6Lf$ & $%.6Lf$ & $%.6Lf$ \\\\\n",
             qb.min, qb.max, i, 2 * parts, i + 1, 2 * parts,
             d2b[i].min, d2b[i].max,
             d1bleft[i].min, d1bleft[i].max);
    }
    printf("Proved f correct shape for "
           "%Lf <= q <= %Lf, %d steps.\n",
           qb.min, qb.max, parts);
  } else {
    die("FAILED to prove f correct shape for "
        "%Lf <= q <= %Lf.",
        qb.min, qb.max);
  }
}


int
main (void)
{
  /* Compute q0 by cubic formula.  */
  q0 = (cbrtl(3.0L * sqrtl(33.0L) + 17.0L)
        - cbrtl(3.0L * sqrtl(33.0L) - 17.0L)
        - 1.0L) / 3.0L;
  p0 = 1.0L - q0;
  em1 = expl(-1.0L);
  em2 = expl(-2.0L);
  em2x4 = 4 * em2;
  prove_d2f_neg((BOUNDS) { 0.0L, 0.4L });
  prove_f_pos((BOUNDS) { 0.4L, 0.48L });
  prove_f_shape((BOUNDS) { 0.48L, 0.5L });
  prove_f_shape((BOUNDS) { 0.5L, 0.55L });
  exit(EXIT_SUCCESS);
}