/* complete_intersection.cpp
 * Computes the Hodge diamond of a complete intersection in P^n.
 *
 * Copyright (C) 2007, 2010 Nicolas Addington (adding@math.duke.edu)
 *
 * This program is free software; you can redistribute it and/or modify 
 * it under the terms of the GNU General Public License as published by 
 * the Free Software Foundation; either version 2 of the License, or (at 
 * your option) any later version. */

/* SYNTAX
 *
 *    complete_intersection n d_1 d_2 ...
 *
 * Computes the Hodge diamond of the complete intersection of 
 * hypersurfaces of degrees d_1, d_2, ... in P^n. */


#include <cstdio>
#include <cstdlib>
#include <climits>
#include <cmath>

// function prototypes
int choose(int, int);
inline int minus_1_to_the(int);
inline void print_centered(int, int, int);

// our class
class hilbert_poly {
public:
  int len;
  int *coeffs;
  /* coeffs[0] is the coefficient of (t+n choose n)
   * coeffs[1] is the coefficient of (t+n-1 choose n)
   * coeffs[2] is the coefficient of (t+n-2 choose n)
   * etc. */

  hilbert_poly(int _len) : len(_len), coeffs(new int[len]) {}
  hilbert_poly() : len(0) {}

  int euler_char(int);
  hilbert_poly operator*(int);
  hilbert_poly operator-(const hilbert_poly &);
  hilbert_poly operator()(int);
};


// the program
int main(int argc, char **argv) {
  const char *syntax_error = "usage: complete_intersection n d_1 d_2 ...\n"
    "Computes the Hodge diamond of the complete intersection of " 
    "hypersurfaces of degrees d_1, d_2, ... in P^n.\n";

  if (argc < 2) { fputs(syntax_error, stderr); return 1; }

  int n = atoi(argv[1]); // we're in P^n
  if (n < 0) { fputs(syntax_error, stderr); return 1; }

  int m = argc - 2; // number of hypersurfaces
  int *d = new int[m]; // degrees of hypersurfaces
  for (int i = 0; i < m; i++) {
    d[i] = atoi(argv[i+2]);
    if (d[i] < 1) { fputs(syntax_error, stderr); return 1; }
  }

  int dim = n - m; // dimension of intersection
  if (dim < 0) { fputs("Too many hypersurfaces.\n", stderr); return 1; }

  /* We start with O_{P^n} and Omega^p_{P^n} (the structure sheaf and sheaves of holomorphic p-forms).
   * Next we restrict them to V_1 := the first hypersurface.
   * Next we compute O_{V_1} and Omega^p_{V_1}.
   * Next we restrict these to V_2 := the intersection of V_1 and the second hypersurface.
   * And so on. */

  hilbert_poly *Omega = new hilbert_poly[dim + 1],
               *Restr = new hilbert_poly[dim + 1];

  // the structure sheaf of P^n
  hilbert_poly O(1);
  O.coeffs[0] = 1;

  Omega[0] = O;
  // the Euler sequence and its exterior powers
  for (int p = 1; p <= dim; p++)
    Omega[p] = O(-p) * choose(n+1, p) - Omega[p-1];

  // for each hypersurface,
  for (int i = 0; i < m; i++) {
    // restrict Omega^p
    for (int p = 0; p <= dim; p++)
      Restr[p] = Omega[p] - Omega[p](-d[i]);

    // the structure sheaf is the restriction of the structure sheaf
    Omega[0] = Restr[0];

    // the adjunction formula and its exterior powers
    for (int p = 1; p <= dim; p++)
      Omega[p] = Restr[p] - Omega[p-1](-d[i]);
  }

  // Lefschetz hyperplane theorem
  int *middle_line = new int[dim+1];
  for (int p = 0; p <= dim; p++)
    middle_line[p] = minus_1_to_the(dim-p) * Omega[p].euler_char(n)
                     - (2*p == dim ? 0 : minus_1_to_the(dim));


  // print the result
  int *len = new int[dim+1], cell_width = 2;
  for (int p = 0; p <= dim; p++) {
    len[p] = snprintf(NULL, 0, "%d", middle_line[p]);
    if (len[p] + 1 > cell_width) cell_width = len[p] + 1;
  }
  for (int i = 0; i <= 2*dim; i++) {
    if (i < dim) {
      printf("%*s", cell_width*(dim-i), "");
      for (int j = 0; j <= i; j++) {
        print_centered(2*j == i ? 1 : 0, 1, cell_width);
        printf("%*s", cell_width, "");
      }
    } else if (i == dim)
      for (int j = 0; j <= dim; j++) {
        print_centered(middle_line[j], len[j], cell_width);
        printf("%*s", cell_width, "");
      }
    else {
      printf("%*s", cell_width*(i-dim), "");
      for (int j = 0; j <= 2*dim-i; j++) {
        print_centered(2*j == 2*dim-i ? 1 : 0, 1, cell_width);
        printf("%*s", cell_width, "");
      }
    }
    putchar('\n');
  }


  // let the garbage collect itself when the process terminates

  return 0;
}


/* The binomial coefficient a choose b; we permit a < 0.
 * The use of floating-point is a little quick and dirty,
 * but it avoids unnecessary overflows. */
int choose(int a, int b) {
  if (a < 0) return minus_1_to_the(b) * choose(b-a-1, b);
  if (a < b) return 0;
  if (a-b < b) b = a-b;
  double ret = 1;
  for (int i = 0; i < b; i++) {
    ret /= b-i;
    if (INT_MAX / (a-i) < ret) {
      fputs("Overflow.\n", stderr);
      exit(1);
      // I'm too lazy to check for overflows anywhere else.
    }
    ret *= a-i;
  }
  return (int)rint(ret);
}

// self-explanatory
int minus_1_to_the(int p) {
  return p % 2 == 0 ? 1 : -1;
}

// print an integer k, padded with spaces on either side
void print_centered(int k, int len, int cell_width) {
  int left_space = (cell_width-len)/2;
  printf("%*s%-*d", left_space, "", cell_width - left_space, k);
}


// plug in t = 0
int hilbert_poly::euler_char(int n) {
  int ret = 0;
  for (int i = 0; i < len; i++)
    if (coeffs[i] != 0) ret += coeffs[i] * choose(n-i, n);
  return ret;
}

// multiply by an integer
hilbert_poly hilbert_poly::operator*(int k) {
  hilbert_poly ret(len);
  for (int i = 0; i < len; i++)
    ret.coeffs[i] = coeffs[i] * k;
  return ret;
}

// subtract
hilbert_poly hilbert_poly::operator-(const hilbert_poly &rhs) {
  if (len < rhs.len) {
    hilbert_poly ret(rhs.len);
    for (int i = 0; i < len; i++)
      ret.coeffs[i] = coeffs[i] - rhs.coeffs[i];
    for (int i = len; i < rhs.len; i++)
      ret.coeffs[i] = -rhs.coeffs[i];
    return ret;
  } else {
    hilbert_poly ret(len);
    for (int i = 0; i < rhs.len; i++)
      ret.coeffs[i] = coeffs[i] - rhs.coeffs[i];
    for (int i = rhs.len; i < len; i++)
      ret.coeffs[i] = coeffs[i];
    return ret;
  }
}

// twist
hilbert_poly hilbert_poly::operator()(int k) {
  if (k > 0) { fputs("Can't twist that way.\n", stderr); exit(1); }
  hilbert_poly ret(len - k);
  for (int i = 0; i < -k; i++)
    ret.coeffs[i] = 0;
  for (int i = 0; i < len; i++)
    ret.coeffs[i-k] = coeffs[i];
  return ret;
}