/*

  Modular symbols tools

 */

#include "pari.h"

GEN ellcuspint(GEN,GEN,GEN);
GEN ellcuspintx(GEN,GEN);

GEN ellcuspint(GEN table, GEN a, GEN b) {
  long ltop = avma;
  GEN outval;

  outval = gsub(ellcuspintx(table,b),ellcuspintx(table,a));

  return gerepileupto(ltop,outval);
}

GEN ellcuspintx(GEN table, GEN x) {
  long N,kappa,k,qqkp1,qqkp2,qqkp3,msym1,msym2;
  GEN SUM,aa;

  N = lg(table) - 1;

  if (typ(x) == t_VEC) SUM = (GEN) ((GEN) table[1])[N];
  else if (typ(x) == t_INT) {
    SUM = cgetg(t_VEC,3); SUM[1] = zero; SUM[2] = zero;
  } else {

    aa = gcf(x);
    kappa = lg(aa) - 2;
    SUM = cgetg(t_VEC,3); SUM[1] = zero; SUM[2] = zero;
    k = 1;
    qqkp1 = 0; qqkp2 = 1;

    while (k <= kappa) {
      qqkp3 = itos((GEN) aa[k+1])*qqkp2 + qqkp1;
      msym1 = ((k % 2 == 0) ? -1 : 1)*qqkp3 % N;
      msym2 = qqkp2 % N;
      if (msym1 <= 0) msym1 = msym1 + N;
      if (msym2 <= 0) msym2 = msym2 + N;
      SUM = gadd(SUM,(GEN) ((GEN) table[msym2])[msym1]);
      k++; qqkp1 = qqkp2; qqkp2 = qqkp3;
    }

  }

  return SUM;
}