/* double_cover.cpp
 * Computes the Hodge diamond of a branched double cover of P^n.
 *
 * Copyright (C) 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
 *
 *    double_cover n d
 *
 * Computes the Hodge diamond of the double cover of P^n branched over a 
 * hypersurface of degree d. */


#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-(const hilbert_poly &);
  hilbert_poly operator()(int);
};


// the program
int main(int argc, char **argv) {
  const char *syntax_error = "usage: double_cover n d\n"
    "Computes the Hodge diamond of the double cover of P^n branched "
    "over a hypersurface of degree d.\n";

  if (argc != 3) { fputs(syntax_error, stderr); return 1; }

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

  int d = atoi(argv[2]); // degree of the branch divisor
  if (d < 2) { fputs(syntax_error, stderr); return 1; }
  if (d % 2) { fputs("Degree must be even.\n", stderr); return 1; }


  /* Let pi: X -> P^n be a doule cover branched over a hypersurface of 
   * degree d.  By a theorem of Lazarsfeld, pi^*: H^i(P^n) -> H^i(X) is 
   * an isomorphism for i < n, so it is enough to find the Euler 
   * characteristics of O_X, Omega^1_X, Omega^2_X, etc. */

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

  // the sheaves of holomorphic p-forms on P^n
  hilbert_poly *Omega_P = new hilbert_poly[n+1];
  Omega_P[0] = O;

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

  // the sheaves of holomorphic p-forms on X, pushed down
  hilbert_poly *pi_Omega_X = new hilbert_poly[n+1];
  pi_Omega_X[0] = O + O(-d/2);

  /* Note that X is a hypersurface in the total space of L := O(d/2)
   * cut out by a section of pi^* O(d).
   *
   * Start with the normal sequence 0 -> TX -> TL|_X -> p^* O(d) -> 0.
   * Dualize: 0 -> pi^* O(-d) -> Omega_L^1|_X -> Omega_X^1 -> 0.
   * Exterior power: 0 -> Omega_X^{p-1} ** pi^* O(-d) -> Omega_L^p|_X -> Omega_X^p -> 0.
   * Push down: 0 -> (pi_* Omega_X^{p-1})(-d) -> pi_* Omega_L^p|_X -> pi_* Omega_X^p -> 0.
   * 
   * To get the middle term, start with 0 -> pi^* L -> TL|_X -> pi^* TP^n -> 0.
   * Dualize: 0 -> pi^* Omega_P^1 -> Omega_L^1|_X -> pi^* O(-d/2) -> 0.
   * Exterior power: 0 -> pi^* Omega_P^p -> Omega_L^p|_X -> pi^* (Omega_P^{p-1}(-d/2)) -> 0.
   * Push down: 0 -> Omega_P^p + Omega_P^p(-d/2)
   *              -> pi_* Omega_L^p|_X 
   *              -> Omega_P^{p-1}(-d/2) + Omega_P^{p-1}(-d) -> 0. */

  for (int p = 1; p <= n; p++)
    pi_Omega_X[p] = Omega_P[p] + Omega_P[p](-d/2)
                   + Omega_P[p-1](-d/2) + Omega_P[p-1](-d)
                   - pi_Omega_X[p-1](-d);

  // use Lazarsfeld theorem
  int *middle_line = new int[n+1];
  for (int p = 0; p <= n; p++)
    middle_line[p] = minus_1_to_the(n-p) * pi_Omega_X[p].euler_char(n)
                     - (2*p == n ? 0 : minus_1_to_the(n));


  // print the result
  int *len = new int[n+1], cell_width = 2;
  for (int p = 0; p <= n; 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*n; i++) {
    if (i < n) {
      printf("%*s", cell_width*(n-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 == n)
      for (int j = 0; j <= n; j++) {
        print_centered(middle_line[j], len[j], cell_width);
        printf("%*s", cell_width, "");
      }
    else {
      printf("%*s", cell_width*(i-n), "");
      for (int j = 0; j <= 2*n-i; j++) {
        print_centered(2*j == 2*n-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]) 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;
}

// add
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;
  }
}

// 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;
}