/*
      [web:reg] Hodrick Prescott Filter - DLL
    
      © 2006, Kurt Annen
      email:  annen@web-reg.de
      web:    www.web-reg.de
      
      'Cardware'
      This DLL was written by Kurt Annen (www.web-reg.de/contact.html)
      If you like the program, please send me a postcard.
*/


#include "hp_filter.h"


DLLEXPORT int hp_filter(double *data, UINT nobs,double lambda, double *output, BOOL ret)
{
   /* 
         *  Schneller Algorithmus zur Lösung des LGS BxX=Y
         * unter Berücksichtigung, dass B eine pentadiagonale Matrix ist
         * siehe z.B. Helmut Späth "Numerik, Eine Einführung für
         * Mathematiker und Informatiker", S. 110 ff, Braunschweig & 
         * Wiesbaden, 1994
         *  
   */
          
   double* a;
   double* b;
   double *c; 
   register UINT  i=0; 
   double H1=0;
        double H2=0;
        double H3=0;
        double H4=0;
        double H5=0;
        double HH1=0;
        double HH2=0;
        double HH3=0;
        double HH4=0;
        double HH5=0;
        double Z;
        double HB;
        double HC;
   
   /* Testen, ob nobs >5 */
   
   if (nobs <= 5) return 1;
   
   /* Speicher in den HEAP */
   
   a = (double *) calloc(nobs, sizeof(double));
   b = (double *) calloc(nobs, sizeof(double)) ;
   c = (double *) calloc(nobs, sizeof(double));
   
   /* Testen auf Speicherproblem */
   if (a == NULL || b == NULL || c == NULL) return 2; 
   
   /* Da bei diesem Algorithmus, die Daten überschrieben wrden, müssen wir
      ein neuen Array füllem.
   */
      
   for(i=0;i<nobs;i++) output[i]=data[i];
   
   /* Koeffizientenmatrix erstellen */
   
   
   a[0]=1.0+lambda;
   b[0]=-2.0*lambda;
   c[0]=lambda;
   
   for(i=1;i<nobs-2;i++)
   {
   a[i]=6.0*lambda+1;
   b[i]=-4.0*lambda;
   c[i]=lambda;                                             
   }
   a[1]=5.0*lambda+1;
   a[nobs-1]=1.0+lambda;
   a[nobs-2]=5.0*lambda+1.0;

   b[nobs-2]=-2.0*lambda;
   b[nobs-1]=0.0;
   c[nobs-2]=0.0;
   c[nobs-1]=0.0;
   
  /* Vorwärts du junge Garde */
   for (i=0;i<nobs;i++)
{
Z=a[i]-H4*H1-HH5*HH2;
if (Z==0) return 3;
HB=b[i];
HH1=H1;
H1=(HB-H4*H2)/Z;
b[i]=H1;
HC=c[i];
HH2=H2;
H2=HC/Z;
c[i]=H2;
a[i]=(output[i]-HH3*HH5-H3*H4)/Z;
HH3=H3;
H3=a[i];
H4=HB-H5*HH1;
HH5=H5;
H5=HC;
}

        
 /* Rückwärts muss die junge Garde auch */         
H2=0;
H1=a[nobs-1];
output[nobs-1]=H1;


for (i=nobs-1;i>0;i--)
{

output[i-1]=a[i-1]-b[i-1]*H1-c[i-1]*H2;
H2=H1;
H1=output[i-1];
}

 /* Trend oder zyklische Komponente? */
   if (ret==FALSE) for(i=0;i<nobs;i++) output[i]=data[i]-output[i];
   
 /* Wer will schon Speicherlücken? */
   
   free(a);
   free(b);
   free(c);
   
   return 0;
}


BOOL APIENTRY DllMain (HINSTANCE hInst     /* Library instance handle. */ ,
                       DWORD reason        /* Reason this function is being called. */ ,
                       LPVOID reserved     /* Not used. */ )
{
    switch (reason)
    {
      case DLL_PROCESS_ATTACH:
        break;

      case DLL_PROCESS_DETACH:
        break;

      case DLL_THREAD_ATTACH:
        break;

      case DLL_THREAD_DETACH:
        break;
    }

    /* Returns TRUE on success, FALSE on failure */
    return TRUE;
}