/* hodge.cpp -- computes the Hodge diamond of a complete intersection in P^n
 *
 * Copyright (C) 2007 Nicolas Addington (adding@math.wisc.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
 *
 *    hodge n d_1 d_2 ...
 *
 * Computes the Hodge diamond of the complete intersection of hypersurfaces of
 * degrees d_1, d_2, ... in P^n. */

/* CAVEAT OVERFLOW
 *
 * I'm too lazy to check for integer overflows.  But if the numbers are that big,
 * what were you going to do with them anyway? */


#include <iostream>

using namespace std;

// global variables
int n; // we're in P^n
int *choose_cache = NULL;

// function prototypes
int choose(int, int);
inline int minus_1_to_the(int);
void print_centered(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)
   * etc. */

  hilbert_poly(int);
  hilbert_poly() {}

  int euler_char();
  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: hodge_no n d_1 d_2 ...\n"
    "Finds the Hodge numbers of the complete intersection of hypersurfaces of "
    "degrees d_1, d_2, ... in P^n.\n";

  if (argc < 2) { cerr << syntax_error; return 1; }

  n = atoi(argv[1]); // we're in P^n
  if (n < 1) { cerr << syntax_error; return 1; }

  int m = argc - 2; // number of hypersurfaces
  int *d = new int[m] - 1; // degrees of hypersurfaces.  1-based array
  for (int i = 1; i <= m; i++) {
    d[i] = atoi(argv[i+1]);
    if (d[i] < 1) { cerr << syntax_error; return 1; }
  }

  int dim = n - m; // dimension of intersection
  if (dim < 0) { cerr << "too many hypersurfaces" << endl; 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.
   * rinse and repeat. */

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

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

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

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

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

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

  /* we will need a bunch of binomial coefficients.
   * (max_len-2 choose n) is the biggest we will need,
   * where max_len is the max of Omega[p].len.
   * choose_cache[j] will be (j choose n) ultimately. */
  int max_len = 0;
  for (int p = 0; p <= dim/2; p++)
    if (Omega[p].len > max_len)
      max_len = Omega[p].len;
  if (max_len-2 > n) {
    int choose_cache_len = max_len - n - 1; // why? you'll see...
    choose_cache = new int[choose_cache_len];
    for (int i = 0; i < choose_cache_len; i++)
      choose_cache[i] = 1;
    for (int i = 1; i <= n; i++)
      for (int j = 1; j < choose_cache_len; j++)
        choose_cache[j] += choose_cache[j-1];
    choose_cache -= n; // ...here.
  }

  // Lefschetz hyperplane theorem
  int *middle_line = new int[dim + 1];
  for (int p = 0; p <= dim/2; p++) {
    middle_line[p] = minus_1_to_the(dim - p) * Omega[p].euler_char();
    if (p < dim/2 || dim % 2 == 1)
      middle_line[p] -= minus_1_to_the(dim);

    // it's Kaehler
    middle_line[dim - p] = middle_line[p];
  }


  // print the result.
  // this is almost as much work as the computation.

  int width = 0;
  for (int p = 0; p <= dim/2; p++) {
    int w = snprintf(NULL, 0, "%d", middle_line[p]) + 1;
    if (w > width) width = w;
  }

  char *spaces = new char[width+1];
  for (int i = 0; i < width; i++)
    spaces[i] = ' ';
  spaces[width] = 0;

  for (int i = 0; i <= 2*dim; i++) {
    if (i < dim) {
      for (int j = 0; j < dim-i; j++)
        cout << spaces;
      for (int j = 0; j <= i; j++) {
        print_centered(i % 2 == 0 && j == i/2 ? 1 : 0, width);
        cout << spaces;
      }
    } else if (i == dim)
      for (int j = 0; j <= dim; j++) {
        print_centered(middle_line[j], width);
        cout << spaces;
      }
    else {
      for (int j = 0; j < i-dim; j++)
        cout << spaces;
      for (int j = 0; j <= 2*dim-i; j++) {
        print_centered(i % 2 == 0 && j == dim-i/2 ? 1 : 0, width);
        cout << spaces;
      }
    }
    cout << endl;
  }


  // let the garbage collect itself when the process terminates

  return 0;
}


// binomial coefficient
// a < 0 is permitted
int choose(int a, int b) {
  if (a == b)
    return 1;
  if (0 <= a && a < b)
    return 0;
  if (a < 0 && b == n && choose_cache)
    return minus_1_to_the(n) * choose_cache[n-1-a];

  long long ret = 1;
  for (int i = 0; i < b; i++)
    ret *= a-i;
  for (int i = 1; i <= b; i++)
    ret /= i;
  return 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 width) {
  int len = 0 <= k && k < 10 ? 1 : snprintf(NULL, 0, "%d", k);
  for (int i = 0; i <= width; i++)
    if (i == (width-len)/2) {
      cout << k;
      i += len;
    } else
      cout << ' ';
}


// hilbert_poly constructor
hilbert_poly::hilbert_poly(int _len) {
  len = _len;
  coeffs = new int[len];
}

// plug in t = 0
int hilbert_poly::euler_char() {
  int ret = 0;
  for (int i = 0; i < len; i++)
    if (coeffs[i]) 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) { cerr << "can't twist that way" << endl; 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;
}