/* See dgpulse.h */
/* Last edited on 2004-11-02 21:51:28 by stolfi */

#include "dgtree.h"
#include "dgspline.h"
#include "dgpulse.h"
#include "dgbezier.h"
#include "dglocus.h"
#include "js.h"
#include <math.h>

/* INTERNAL PROTOTYPES */

double dg_mother_H_pulse_poly
  ( dg_Cont c, 
    dg_PulseIndex pix, 
    double z
  ); 
  /* Computes the mother H-pulse of continuity {c}, degree {g=2*c+1},
    and index {pix} (which must lie in {0..c}) for the argument {z}
    which must lie in {[0_1]}. */

double dg_mother_B_pulse_poly
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix, 
    double z
  ); 
  /* Computes the mother B-pulse polynomial of continuity {c}, degree
    {g}, and index {pix} (which must lie in {0..c}) for the argument
    {z} which must lie in {[0_1]}. */

dg_GridSize dg_mother_H_pulse_supp_count
  ( dg_Cont c, 
    dg_PulseIndex pix
  ); 
  /* Returns the number of unit-grid cells that comprise the support 
    of the mother H-pulse. */

dg_Interval dg_mother_H_pulse_sum_range
  ( dg_Cont c, 
    dg_PulseIndex pix,
    dg_GridSize sz
  );
  /* An interval containing {P(z) = dg_mother_H_pulse_sum(c,pix,sz,z)}
    in [-1_+1]. Takes wrap-around into account. */

double dg_mother_H_pulse_sum_max
  ( dg_Cont c, 
    dg_PulseIndex pix,
    dg_GridSize sz
  );
  /* The maximum absolute value of {P(z) = dg_mother_H_pulse_sum(c,pix,sz,z)} in 
    [-1_+1]. Takes wrap-around into account. */

dg_GridSize dg_mother_B_pulse_supp_count
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix
  );
  /* Returns the number of unit-grid cells that comprise the support 
    of the mother H-pulse. */
    
dg_Interval dg_mother_B_pulse_sum_range
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz
  );
  /* The range of {P(z) = dg_mother_B_pulse_sum(c,g,pix,sz,z)} in 
    [-1_+1]. Takes wrap-around into account. */
    
double dg_mother_B_pulse_sum_max
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz
  );
  /* The maximum absolute value of {P(z) = dg_mother_B_pulse_sum(c,g,pix,sz,z)} in 
    [-1_+1]. Takes wrap-around into account. */

/* IMPLEMENTATIONS */

double dg_pulse_eval
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g,
    dg_PulseIndex pix, 
    dg_GridSize sz,
    dg_GridPos j, 
    double x
  )
  { /* Scale argument {[0 _ 1)} to {[0 _ sz)} and shift from {j} to 0: */
    double y = x*sz - j;
    return dg_mother_pulse_sum(kind, c, g, pix, sz, y);
  }

dg_GridSize dg_pulse_supp_count
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g,
    dg_PulseIndex pix, 
    dg_GridSize sz
  )
  { dg_GridSize m = dg_mother_pulse_supp_count(kind, c, g, pix);
    return (m <= sz ? m : sz);
  }
  
/* GENERIC MOTHER PULSES */

double dg_mother_pulse_sum
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz, 
    double z
  )
  {
    /* Reduce {z} to {[-sz+1 _ 1)} modulo sz: */
    z -= floor(z/sz)*sz;
    if (z >= 1) { z -= sz; }
    affirm(z >= -(double)(sz+1), "bad frac - too small");
    affirm(z < 1.0, "bad frac - too big");
    
    /* Discover low endpoint {lo} of mother pulse's support: */
    dg_GridPos lo = -(dg_mother_pulse_supp_count(kind, c, g, pix) - 1);
    /* Add {w(z)} for all {z} in {[lo _ 1)} congruent to {y} mod {sz}: */
    double f = 0.0;
    while (z >= lo) { f += dg_mother_pulse_eval(kind, c, g, pix, z); z -= sz; }
    return f;
  }

double dg_mother_pulse_eval
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix, 
    double z
  )
  { affirm(g >= 0, "negative degree");
    switch(kind)
      { case DG_PK_H: 
          affirm(g == 2*c + 1, "bad degree for H-pulse");
          return dg_mother_H_pulse_eval(c, pix, z);
        case DG_PK_B: 
          return dg_mother_B_pulse_eval(c, g, pix, z);
        default:
          affirm(FALSE, "unknown pulse kind");
          return 0.0; /* To keep the compiler happy. */
      }
  }

dg_Interval dg_mother_pulse_sum_range
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz
  )
  { affirm(g >= 0, "negative degree");
    switch(kind)
      { case DG_PK_H: 
          affirm(g == 2*c + 1, "bad degree for H-pulse");
          return dg_mother_H_pulse_sum_range(c, pix, sz);
        case DG_PK_B: 
          return dg_mother_B_pulse_sum_range(c, g, pix, sz);
        default:
          affirm(FALSE, "unknown pulse kind");
          /* Compiler pacifier: return an empty interval: */
          return (dg_Interval){{ +1.0, -1.0 }}; 
      }
  }

dg_GridSize dg_mother_pulse_supp_count
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix
  )
  { affirm(g >= 0, "negative degree");
    switch(kind)
      { case DG_PK_H: 
          affirm(g == 2*c + 1, "bad degree for H-pulse");
          return dg_mother_H_pulse_supp_count(c, pix);
        case DG_PK_B: 
          return dg_mother_B_pulse_supp_count(c, g, pix);
        default:
          affirm(FALSE, "unknown pulse kind");
          return 0.0; /* To keep the compiler happy. */
      }
  }
  
void dg_mother_pulse_to_bezier
  ( dg_PulseKind kind,
    dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    double **bp
  )
  {
    affirm(g >= 0, "negative degree");

    int wd = dg_mother_pulse_supp_count(kind, c, g, pix);
    int k;
    
    /* Allocate Bézier coeff vectors as needed: */
    for (k = 0; k < wd; k++)
      { if (bp[k] == NULL) 
          { bp[k] = notnull(malloc((g+1)*sizeof(double))); }
      }

    switch(kind)
      { case DG_PK_H: 
          affirm(g == 2*c + 1, "bad degree for H-pulse");
          dg_mother_H_pulse_to_bezier(c, pix, bp);
        case DG_PK_B: 
          dg_mother_B_pulse_to_bezier(c, g, pix, bp);
        default:
          affirm(FALSE, "unknown pulse kind");
          return; /* To keep the compiler happy. */
      }
  }

/* H-PULSES */

double dg_mother_H_pulse_eval
  ( dg_Cont c, 
    dg_PulseIndex pix, 
    double z
  )
  { affirm(c >= 0, "bad contin for H-pulse");
    affirm((pix >= 0) && (pix <= c), "invalid H-pulse index");
    /* Evaluate pulse: */
    if ((z <= -1.0) || (z >= +1.0))
      { return 0.0; }
    else if (z >=  0.0) 
      { return dg_mother_H_pulse_poly(c, pix, z); }
    else
      { double s = (pix % 2 == 0 ? 1.0 : -1.0);
        return s * dg_mother_H_pulse_poly(c, pix, -z);
      }
  }

double dg_mother_H_pulse_poly
  ( dg_Cont c, 
    dg_PulseIndex pix, 
    double z
  )
  { /* Assumes {c} and {pix} are OK. */
    /* !!! Needs a general implementation !!! */

    /* 
      HR[0,0] = 1-x

      HR[1,0] = (2*x+1)*u^2
      HR[1,1] = 3*x*u^2

      HR[2,0] = (6*x^2+3*x+1)*u^3
      HR[2,1] = 5*x*(3*x+1)*u^3
      HR[2,2] = 5*x^2*u^3
    */
  
    double u = 1.0 - z;
    switch (c) 
      { case 0: 
          switch (pix)
            { case 0: return u;
            }
        case 1:
          switch (pix)
            { case 0: return (2*z + 1)*u*u;
              case 1: return 3*z*u*u;
            }
        case 2:
          switch (pix)
            { case 0: return ((6*z + 3)*z + 1)*u*u*u;
              case 1: return 5*z*(3*z + 1)*u*u*u;
              case 2: return 5*(z*z*u*u*u);
            }
      }
    affirm(FALSE, "not implemented yet");
    return 0.0;
  }

dg_Interval dg_mother_H_pulse_sum_range
  ( dg_Cont c, 
    dg_PulseIndex pix,
    dg_GridSize sz
  )
  { affirm(c >= 0, "bad contin for H-pulse");
    affirm((pix >= 0) && (pix <= c), "invalid H-pulse index");
    double pmax = dg_mother_H_pulse_sum_max(c, pix, sz);
    double pmin = (pix % 2 == 0 ? 0.0 : -pmax);
    return (dg_Interval){{ pmin, pmax }};
  }
  
double dg_mother_H_pulse_sum_max
  ( dg_Cont c, 
    dg_PulseIndex pix,
    dg_GridSize sz
  )
  { /* Assumes {c} and {pix} are OK. */
    /* !!! Needs a general implementation !!! */
  
    if (sz > 1)
      { /* No wrap-around: */
        switch (c) 
          { case 0: 
              switch (pix)
                { case 0: return 1.0;
                }
            case 1:
              switch (pix)
                { case 0: return 1.0;
                  case 1: return 12.0/27.0;
                }
            case 2:
              switch (pix)
                { case 0: return 1.0;
                  case 1: return 80.0/81.0;
                  case 2: return 1080.0/3125.0;
                }
          }
      }
    else
      { /* Wrap-around pulse: */
        switch (c) 
          { case 0: 
              switch (pix)
                { case 0: return 1.0;
                }
            case 1:
              switch (pix)
                { case 0: return 1.0;
                  case 1: return 0.28867513465254790917;
                    /* sqrt(3.0)/6.0 */
                }
            case 2:
              switch (pix)
                { case 0: return 1.0;
                  case 1: return 0.73374572609092830121;
                    /* 1/90*sqrt(15)*sqrt(15-2*sqrt(30))*(3+sqrt(30)) */
                  case 2: return 10.0/16.0;
                }
          }
      }
    affirm(FALSE, "not implemented yet");
    return 1.0;
  }

dg_GridSize dg_mother_H_pulse_supp_count
  ( dg_Cont c, 
    dg_PulseIndex pix
  )
  { affirm(c >= 0, "bad contin for H-pulse");
    affirm((pix >= 0) && (pix <= c), "invalid H-pulse index");
    return 2;
  }

void dg_mother_H_pulse_to_bezier
  ( dg_Cont c, 
    dg_PulseIndex pix,
    double **bp
  )
  {
    double *bp0 = bp[0];
    double *bp1 = bp[1];
    
    switch (c) 
      { case 0: /* {g == 1} */
          switch (pix)
            { case 0: 
                bp0[0] = 0; bp0[1] = 1;
                bp1[0] = 1; bp1[1] = 0;
                return;
            }
        case 1: /* {g == 3} */
          switch (pix)
            { case 0: 
                bp0[0] = 00; bp0[1] = 00; bp0[2] = +1; bp0[3] = +1;
                bp1[0] = +1; bp1[1] = +1; bp1[2] = 00; bp1[3] = 00;
                return;
              case 1: 
                bp0[0] = 00; bp0[1] = 00; bp0[2] = -1; bp0[3] = 00;
                bp1[0] = 00; bp1[1] = +1; bp1[2] = 00; bp1[3] = 00;
                return;
            }
        case 2: /* {g == 5} */
          switch (pix)
            { case 0:
                /* ((6*z + 3)*z + 1)*u*u*u; */
                bp0[0] = 00; bp0[1] = 00; bp0[2] = 00; bp0[3] = +1; bp0[4] = +1; bp0[5] = +1;
                bp1[0] = +1; bp1[1] = +1; bp1[2] = +1; bp1[3] = 00; bp1[4] = 00; bp1[5] = 00;
                return;
              case 1: 
                /*  5*z*(3*z + 1)*u*u*u; */
                bp0[0] = 00; bp0[1] = 00; bp0[2] = 00; bp0[3] = -2; bp0[4] = -1; bp0[5] = 00;
                bp1[0] = 00; bp1[1] = +1; bp1[2] = +2; bp1[3] = 00; bp1[4] = 00; bp1[5] = 00;
                return;
              case 2: 
                /*  5*(z*z*u*u*u); */
                bp0[0] = 00; bp0[1] = 00; bp0[2] = 00; bp0[3] = +1; bp0[4] = 00; bp0[5] = 00;
                bp1[0] = 00; bp1[1] = 00; bp1[2] = +1; bp1[3] = 00; bp1[4] = 00; bp1[5] = 00;
                return;
            }
      }
    affirm(FALSE, "not implemented yet");
    return 1.0;
  }

/* B-PULSES */

double dg_mother_B_pulse_eval
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix, 
    double z
  )
  { affirm(c >= -1, "bad contin for B-pulse");
    affirm(g >= 0, "negative degree for B-pulse");
    affirm(g >= 2*c + 1, "bad degree for B-pulse");
    affirm(g <= DG_MAX_B_PULSE_DEGREE, "degree too high for B-pulse");
    affirm((pix >= 0) && (pix <= g-c-1), "invalid index");
    if (pix > c)
      { /* A Bernstein-Bézier element with single-cell support: */
        if ((z <= 0.0) || (z >= 1.0)) 
          { return 0.0; }
        else
          { return dg_bezier_poly(g, pix, z); }
      }
    else
      { /* A two-interval pulse: */
        if ((z <= -1.0) || (z >= +1.0))
          { return 0.0; }
        else
          { if (z >=  0.0) 
              { return dg_mother_B_pulse_poly(c, g, pix, z); }
            else
              { return dg_mother_B_pulse_poly(c, g, c-pix, -z); }
          }
      }
  }

double dg_mother_B_pulse_poly
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix, 
    double z
  )
  { /* Assumes {g}, {c} and {pix} are OK. */
    double u = 1.0 - z;
    /* Compute the nonzero Bézier coeffs, namely {a[0..pix]}. */
    double a[DG_MAX_B_PULSE_DEGREE + 1];
    int r, s;
    a[pix] = 1.0;
    for (r = pix-1; r >= 0; r--) { a[r] = a[r+1]/2; }
    for (s = pix+1; s <= c; s++)
      { a[pix] = a[pix]/2;
        for (r = pix-1; r >= 0; r--) { a[r] = (a[r] + a[r+1])/2; }
      }
    /* Apply DeCasteljau on coeffs {a[0..pix]}: */
    for (s = 1; s <= g; s++)
      { for (r = 0; (r < pix) && (r <= g-s); r++)
          { a[r] = u*a[r] + z*a[r+1]; }
        if (s <= g - pix) { a[pix] = u*a[pix]; }
      }
    return a[0];
  }

dg_Interval dg_mother_B_pulse_sum_range
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz
  )
  { affirm(c >= -1, "bad contin for B-pulse");
    affirm(g >= 0, "negative degree for B-pulse");
    affirm(g >= 2*c + 1, "bad degree for B-pulse");
    affirm(g <= DG_MAX_B_PULSE_DEGREE, "degree too high for B-pulse");
    affirm((pix >= 0) && (pix <= g-c-1), "invalid index");
    double pmax = dg_mother_B_pulse_sum_max(c, g, pix, sz);
    return (dg_Interval){{ 0.0, pmax }};
  }

double dg_mother_B_pulse_sum_max
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    dg_GridSize sz
  )
  { /* Assumes {g,c,pix} are OK. */
    /* !!! Needs a general implementation !!! */
    if (pix > c)
      { /* Single-cell Bézier pulse: */
        return dg_bezier_poly_max(g, pix);
      }
    else if (sz > 1)
      { /* No wrap-around: */
        switch (c) 
          { case 0: 
              return 1.0;
            case 1:
              return 16.0/25.0;
            case 2:
              switch (pix)
                { case 0: return 0.53079492556952189365;
                  case 1: return 0.50;
                  case 2: return 0.53079492556952189365;
                }
          }
      }
    else
      { /* Wrap-around pulse: */
        switch (c) 
          { case 0: 
              return 1.0;
            case 1:
              return 0.64433756729740644112;
                /* 1/2+sqrt(3)/12 */
            case 2:
              switch (pix)
                { case 0: return 0.5314173849806585844;
                  case 1: return 0.50; 
                  case 2: return 0.5314173849806585844;
                }
          }
      }
    affirm(FALSE, "not implemented yet");
    return 1.0;
  }

dg_GridSize dg_mother_B_pulse_supp_count
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix
  )
  { affirm(c >= -1, "bad contin for B-pulse");
    affirm(g >= 0, "negative degree for B-pulse");
    affirm(g >= 2*c + 1, "bad degree for B-pulse");
    affirm(g <= DG_MAX_B_PULSE_DEGREE, "degree too high for B-pulse");
    affirm((pix >= 0) && (pix <= g-c-1), "invalid index");
    return (pix > c ? 1 : 2);
  }

void dg_mother_B_pulse_to_bezier
  ( dg_Cont c, 
    dg_Degree g, 
    dg_PulseIndex pix,
    double **bp
  )
  {
    affirm(g >= 0, "negative degree");
    /* To be implemented */
    affirm(FALSE, "not implemented yet");
    return 1.0;
  }