/* Beregning af koefficienter til q-ekspansionen for cuspformen i Knapps
eksempel side 304. Beregning af det tilsvarende gitter (udspÊndt af
omega1,omega2).

For at benytte typen long_double_complex er man n¯dt til at bruge c++ 
compileren. Kompileres med "g++ qexp.cpp". Filen "a.out" kan nu k¯res.

DesvÊrre kan c++'s standard biblioteker kun regne med den prÊcision der
underst¯ttes af hardwaren. (dvs omkring 52-64 bit). SÂ der er ingen
garanti for at resultaterne passer bare nogenlunde, og ingen kontrol med hvor
mange bits der smides vÊk ved afrundinger.

Beregningen af koefficineterne for q-expansionen tager kvadratisk tid.
Her viser C++ sin styrke. Man har fuld kontrol over hvad der sker og er sikker
p at der ikke spildes tid med allokere og deallokere hukommelse.

Den bedste kombination af brug af C og et symbolsk matematikprogram vil
vÊre at lade c-koden beregne koefficienterne og lade det symbolske
matematikprogram evaluere H.

Skrevet af Anders Nedergaard Jensen. December 1999.
*/
#define nterms 300
#include <stdio.h>
#include <complex.h>
#include <iostream.h>
#include <math.h>

typedef long_double_complex compl;

compl i(0.0,1.0);
#define pi 3.14159265358979323846L
#define V4 8,1,-33,-4
#define V6 9,1,-55,-6

// Hack, for at f ulovlige indices til at give 0
int dummy[nterms*30];
#define ctable ((&(dummy[0]))+nterms*29)

void calckoeff()
{
  /* Beregning af c_n'erne */
  int m,n;

  ctable[0]=0;
  ctable[1]=1;
  for(m=2;m<nterms;m++)ctable[m]=0;
  for(n=1;n<nterms;n++)
    {
      for(m=nterms-1;m>=0;m--)
      ctable[m]=ctable[m]
    -2*ctable[m-n]
    +1*ctable[m-2*n]
    -2*ctable[m-11*n]
    +4*ctable[m-12*n]
    -2*ctable[m-13*n]
    +1*ctable[m-22*n]
    -2*ctable[m-23*n]
    +1*ctable[m-24*n];
    }
}

compl H(compl q)
{
  /* evaluer H i punktet q */
  compl h(0.0,0.0);
  // tabel til at opl¯fte q til nte potens med kun et logaritmisk antal 
  // operationer (for ikke at miste n¯jagtighed)
  // (Noget tilsvarende burde man mÂske overveje at g¯re for summen af n led)
  compl qpow[20];

  int n;

  qpow[1]=q;
  for(n=2;n<20;n++)qpow[n]=qpow[n-1]*qpow[n-1];

  for(n=1;n<nterms;n++)
    {
      compl a(1.0,0.0);
      int j=1,N=n;
      while(N)
	{
	  if(N&1)a*=qpow[j];
	  N>>=1;
	  j++;
	}
      h+=(ctable[n]*a)/n;
    }
  return h;
}

compl transform(compl tau, int a, int b, int c, int d)
{
  /* Lad matricen virke p tau..... */
  return (a*tau+b)/(c*tau+d);
}

void print(compl a)
{
  /* Udskrivning af a */
  printf("%28.18G+%28.18Gi\n",(double)real(a),(double)imag(a));
}

void calculate(compl tau)
{
  /* Udf¯r beregningerne fra side 305 p det givne tau */
  compl pi2i(0.0,2*pi),h1,h2,h3,omega1,omega2,tmp;

  printf("\ntau:");print(tau);
  print(h1=H(exp(pi2i*tau)));
  tmp=transform(tau,V4);
  printf("V4otau:");print(tmp);
  print(h2=H(exp(pi2i*tmp)));
  tmp=transform(tau,V6);
  printf("V6otau:");print(tmp);
  print(h3=H(exp(pi2i*tmp)));

  omega1=(h2-h1)/pi2i;
  omega2=(h3-h1)/pi2i;

  printf("omega1: ");print(omega1);
  printf("omega2: ");print(omega2);
}

void main()
{
  /* Hovedprogram */

  printf("NTerms %i\n",nterms);

  calckoeff();
  for(int n=0;n<nterms;n++)printf("c[%i]=%i \n",n,ctable[n]);

  calculate(compl(0.125,0.025));
  calculate(compl(0,1));
  calculate(compl(0.12,0.02));
  calculate(compl(0.225,0.075));
  calculate(compl(0.525,0.075));
  calculate(compl(0.625,0.05));
  calculate(compl(-0.225,0.15));
  calculate(compl(-0.425,0.25)); // et godt tau
  calculate(compl(-0.125,0.15));
}