/*----------------------------------------------------------------------------
Programm sgnsum:                                   Duesseldorf, den 12.01.1995
Berechnung des exakten Vorzeichens einer Summe a[1]+...+a[m]-b[1]-...-b[n].
Eingabe der Summanden von Typ float.
----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <unistd.h>
#include <sys/time.h>
#include <sys/resource.h>

/*----------------------------------------------------------------------------
typedef tStat:
Typ fuer Statistikdaten.
m..Anzahl der positiven Summanden (Liste a),
n..Anzahl der negativen Summanden (Liste b),
CaseSteps[i]..Anzahl der in Schritt 3 ausgefuehrten 3 Unterfaelle,
iteratione..Anzahl der Iterationen,
nTerm..Anzahl der Elemente in Liste b zur Zeit des Abbruchs,
mTerm..Anzahl der Elemente in Liste a zur Zeit des Abbruchs,
runtime..Laufzeit zur Berechnung des Vorzeichens.
----------------------------------------------------------------------------*/
typedef struct
{
 int n,m,CasesStep3[3],iterations,nTerm,mTerm,runtime;
} tStat;

/*----------------------------------------------------------------------------
Funktion GetCPUUsage:
Rueckgabe: Die vom Beginn des Programms verbrauchte CPU-Zeit in Mikrosekunden.
----------------------------------------------------------------------------*/
int GetCPUUsage(void)
{
 struct rusage usage;
 getrusage(RUSAGE_SELF,&usage);
 return usage.ru_utime.tv_sec*100+usage.ru_utime.tv_usec/10000;
}

/*----------------------------------------------------------------------------
Funktion sort:
Beschreibung: Heapsort aus dem Buch "Numerical Recipes in C" von Press, 
   Flannery, Teukolsky und Vetterling, Cambridge University Press Seite 247.
Eingabe: ra..Feld von doubles, n..Laenge von ra.
Ausgabe: ra wird absteigend sortiert.
----------------------------------------------------------------------------*/
void sort(int n, double *ra)
{
 int l,j,ir,i;
 double rra;
 
 l=(n >> 1)+1;
 ir=n;
 for(;;)
 {
  if (l>1)
   rra=ra[--l];
  else
  {
   rra=ra[ir];
   ra[ir]=ra[1];
   if (--ir==1)
   {
    ra[1]=rra;
    return;
   }
  } /* (l<=1) */
  i=l;
  j=l<<1;
  while (j<=ir)
  {
   if ( j<ir && ra[j]>ra[j+1] ) ++j;
   if ( rra>ra[j] )
   {
    ra[i]=ra[j];
    j+=(i=j);
   }
   else
    j=ir+1;
  }
  ra[i]=rra;
 } /* for(;;) */
}

/*----------------------------------------------------------------------------
Funktion swap:
Beschreibung: swap(a,b) vertauscht a und b.
----------------------------------------------------------------------------*/
void swap(double *a, double *b)
{
 double tmp=*a;
 *a=*b;
 *b=tmp;
}

void BuildHeapFromTop(int n, double *ra)
{
 int i=1,m;
 double top=ra[1];

 while (2*i<=n)
 {
  m= 2*i;
  if ( ra[m]<ra[m+1] ) if (m<n) m++;
  if (top<ra[m]) { ra[i]=ra[m]; i=m; }
  else break; 
 }
 ra[i]=top;
}

void BuildHeapFromBelow(int n, double *ra)
{
 int i=n,m;
 double last=ra[n];
 
 while (i/2>0)
 {
  m= i/2;
  if ( ra[m]<last ) { ra[i]=ra[m]; i=m;}
  else break;
 }
 ra[i]= last;
}

/*----------------------------------------------------------------------------
Funktion RandomDouble:
Rueckgabe: Einen zufaellig erzeugten Wert vom Typ double.
----------------------------------------------------------------------------*/
double RandomDouble(void)
{
 int sgn,EXP;
 if (random()%2) sgn=1;
 else sgn=-1;
 EXP=random()%60;
 return sgn*ldexp(drand48(),EXP);
}

/*----------------------------------------------------------------------------
Funktion :sgnsum
Beschreibung: sgnsum berechnet ohne Rundungsfehler das Vorzeichen der Summe 
   der S[i]. 
Eingabe: S..Ein Feld von Summanden des Typs double, nS..Laenge von S.
Ausgabe: st..Statistik (Typbeschreibung siehe oben).
Rueckgabe: Das Vorzeichen der Summe.
Lokale Variablen: n..Laegne Liste b; m..Laenge Liste a; E..Exponent von a[1];
   F..Exponent von b[1]; sg..-1 -> Summe negativ, 0 -> Summe 0, 1 -> Summe 
   positiv, -2 -> es sind weitere Berechnungen noetig; die restlichen
   Variablen sind Hilfsvariablen.
----------------------------------------------------------------------------*/
int sgnsum(double *S, int nS)
{
 int n,m,E,F,i,sg=-2;
 double *a,*b,as,ass,bs,bss,uu,u,v;
 tStat * st;
 
 st= new tStat;
 
/*----------------------------------------------------------------------------
Initialisierung von st und den Listen a und b.
----------------------------------------------------------------------------*/
 st->runtime=-GetCPUUsage();
 if ( (a=(double*)calloc(nS+3,sizeof(double)))==NULL )
 {
  printf("No mem. Program terminated.\n");
  exit(1);
 }
 if ( (b=(double*)calloc(nS+3,sizeof(double)))==NULL )
 {
  printf("No mem. Program terminated.\n");
  exit(1);
 }
 
/*----------------------------------------------------------------------------
Splitten von S in positive Summanden a und negative Summanden b.
----------------------------------------------------------------------------*/
 n=m=0;
 for (i=1;i<=nS;i++)
  if (S[i]>0) a[++m]=S[i];
  else if (S[i]<0) b[++n]=fabs(S[i]);
 
 st->CasesStep3[0]=st->CasesStep3[1]=st->CasesStep3[2]=0;
 st->iterations=0;
 st->n=n;
 st->m=m;
 
/*----------------------------------------------------------------------------
Faelle abfangen, in denen es nichts zu berechnen gibt.
----------------------------------------------------------------------------*/
 if ( n==0 && m==0) sg=0;
 if ( n==0 ) sg=1;
 if ( m==0 ) sg=-1;

/*----------------------------------------------------------------------------
Sortieren der Listen a und b in absteigender Reihenfolge.
----------------------------------------------------------------------------*/
 if ( m>1 && sg==-2 ) sort(m,a);
 if ( n>1 && sg==-2 ) sort(n,b);
 
/*----------------------------------------------------------------------------
Hauptschleife.
----------------------------------------------------------------------------*/
 while (sg==-2)
 {

/*
  for (i=1;i<=m;i++) printf("a[%i]: %e\n",i,a[i]);
  for (i=1;i<=n;i++) printf("b[%i]: %e\n",i,b[i]);
*/
  st->iterations++;

/*----------------------------------------------------------------------------
Wenn eine Liste zu lang wird, wird das Programm abgebrochen.
----------------------------------------------------------------------------*/
  if (m>nS||n>nS)
  {
   printf("List full. Program aborted\n");
   exit(1);
  }

/*----------------------------------------------------------------------------
 Step 1: (Termination Criteria)
----------------------------------------------------------------------------*/
  if ( n==0 && m==0 ) { sg=0; break; }
  else if ( n==0 || a[1]>n*b[1] && m>0 ) { sg=1; break;}
  else if ( m==0 || b[1]>m*a[1] && n>0 ) { sg=-1; break;}
  
/*----------------------------------------------------------------------------
 Step 2: (Auxiliary variables)
----------------------------------------------------------------------------*/
  as=ass=bs=bss=0;
  
/*----------------------------------------------------------------------------
 Step 3: (Comparision and processing of the leading summands of the lists)
----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------
Berechnug der Exponenten E von a[1] und F von b[1] zur Basis 2.
a[1]-2^E (b[1]-2^F) nullt dann das fuehrende Bit von a[1] (b[1]).
----------------------------------------------------------------------------*/
  frexp(a[1],&E);
  E--;
  frexp(b[1],&F);
  F--;

/*----------------------------------------------------------------------------
Step 3, case (i):
----------------------------------------------------------------------------*/
  if (E==F) 
  {
   st->CasesStep3[0]++;
   if (a[1]>=b[1]) as=a[1]-b[1];
   else bs= b[1]-a[1];
  }

/*----------------------------------------------------------------------------
Step 3, case (ii):
----------------------------------------------------------------------------*/
  else if (E>F)
  {
   st->CasesStep3[1]++;
   uu= ldexp(1,F);
   if (b[1]==uu) u=uu;
   else u=uu*2;
   as=a[1]-u;
   ass=u-b[1];
  }
  
/*----------------------------------------------------------------------------
Step 3, case (iii):
----------------------------------------------------------------------------*/
  else if (F>E)
  {
   st->CasesStep3[2]++;
   uu=ldexp(1,E);
   if (a[1]==uu) v=uu;
   else v=uu*2;
   bs= b[1]-v;
   bss= v-a[1];
  }

/*----------------------------------------------------------------------------
 Step 4: (Rearrangement of the lists, keeping S constant.)
----------------------------------------------------------------------------*/

  if (as==0 && ass==0)
  {
   a[1]=a[m];
   m--;
   BuildHeapFromTop(m,a);
  }
  if (as==0 && ass!=0)
  {
   a[1]=ass;
   BuildHeapFromTop(m,a);
  }
  if (as!=0 && ass==0)
  {
   a[1]=as;
   BuildHeapFromTop(m,a);
  }
  if (as!=0 && ass!=0)
  {
   a[1]=as;
   BuildHeapFromTop(m,a);
   a[++m]=ass;
   BuildHeapFromBelow(m,a);
  }
  
  if (bs==0 && bss==0)
  {
   b[1]=b[n];
   n--;
   BuildHeapFromTop(n,b);
  }
  if (bs==0 && bss!=0)
  {
   b[1]=bss;
   BuildHeapFromTop(n,b);
  }
  if (bs!=0 && bss==0)
  {
   b[1]=bs;
   BuildHeapFromTop(n,b);
  }
  if (bs!=0 && bss!=0)
  {
   b[1]=bs;
   BuildHeapFromTop(n,b);
   b[++n]=bss;
   BuildHeapFromBelow(n,b);
  }
  for (i=2;i<=m;i++)
   if (a[1]<a[i]||a[i]<=0) 
   {
    printf("Fehler in Liste a: %e\n",a[i]);
    getchar();
   }
 } /* while(sg==-2) */
 free(b);
 free(a);
 delete st;

/*----------------------------------------------------------------------------
Update der Statistikvariablen.
----------------------------------------------------------------------------*/
 st->nTerm=n;
 st->mTerm=m;
 st->runtime+=GetCPUUsage();
 return sg;
}