#! /usr/bin/gawk -f
# Last edited on 2012-09-25 21:34:14 by stolfilocal

# Tries to generate a northogonal polynomial basis for R^n
# Under the norm <f,g> = (1/n)*(SUM [f[i]**g[i] : i in 0..n-1])

BEGIN{
  if (n == "") { arg_error("must define {n}"); }
  m = n; # For now.
  if (op == "") { arg_error("must define {op}"); }
  split("", U);    # The first basis element.
  split("", D);    # The diagonal operator, equivalent to multiplying by {x}.
  pi = 3.1415926;
  for (j = 0; j < n; j++) 
    {
      if (op == "hcos") 
        { U[j] = 1.00; D[j] = -cos(pi*j/n); }
      else if (op == "hqua")
        { U[j] = 1.00; D[j] = -cos(pi*(j+0.5)/n); }
      else if (op == "poly")
        { U[j] = 1.00; D[j] = (j + 0.5)/n - 0.5; }
      else if (op == "hran")
        { U[j] = 1.00; D[j] = 0.8*prand(j); }
      else 
        { arg_error(("invalid {op} = \"" op "\"")); }
    }
  if (op == "hran") { dsort(D,n); }
  
  split("", phi);  # The first {m} elements of the basis.
  
  for (i = 0; i < m; i++)
    { # Generate element {phi[i,0..n-1]} {i}:
      if (i == 0)
        { # Basis element 0 is the unit vector: 
          for (j = 0; j < n; j++) { phi[i,j] = U[j]; }
        }
      else
        { # Basis element i is {b*D(phi[i-1]) + a*phi[i-1] + c*phi[i-2]}: 
          cDD1 = dot_DD(phi, D, n, i-1, i-1);
          cD1 = dot_D(phi, D, n, i-1, i-1);
          cD2 = (i == 1 ? 0 : dot_D(phi, D, n, i-1, i-2));
          h = cDD1 - cD1*cD1 - cD2*cD2;
          printf "%4d cDD1 = %+12.6f  cD1 = %+12.6f  cD2 = %+12.6f  h = %+12.6f\n", i, cDD1, cD1, cD2, h > "/dev/stderr";
          if (abs(h) < 1.0e-6) { num_error(i, "division by small {h}"); }
          b = 1/sqrt(h);
          a = -b*cD1;
          c = -b*cD2;
          for (j = 0; j < n; j++) { phi[i,j] = (a + b*D[j])*phi[i-1,j] + c*phi[i-2,j]; }
        }
    }

  # Write out:
  for (j = 0; j < n; j++) 
    { printf "%4d", j;
      for (i = 0; i < m; i++)
        { printf " %+9.6f", phi[i,j];
        }
      printf "\n";
    }

  check_basis(phi,m,n);
}

function check_basis(phi,m,n,  M,i,j,bug,tmp)
  {
    # Checks whether the basis {phi[0..m-1]} with {n} columns is orthonormal.
    bug = 0;
    split("",M);  # Scalar product matrix:
    printf "Checking orthogonality...\n" > "/dev/stderr";
    for (i = 0; i < n; i++)
      { for (j = 0; j < n; j++) 
          { M[i,j] = dot(phi,n,i,j);
            tmp = (i == j);
            if (abs(M[i,j] - tmp) > 1.0e-6) 
              { printf "** <phi[%d],phi[%d]> = %+9.6f\n", i, j, M[i,j] > "/dev/stderr";
                bug=1;
              }
          }
      }

    if (bug)
      { 
        printf "Scalar product matrix:\n" > "/dev/stderr";
        for (i = 0; i < n; i++)
          { for (j = 0; j < n; j++) 
              { M[i,j] = dot(phi,n,i,j);
                printf " %+9.6f", M[i,j] > "/dev/stderr";
              }
            printf "\n" > "/dev/stderr";
          }
      }
  }

function prand(j,  z)
  { 
    # A pseudo-random function of the nonnegative integer {j},
    # with values in {[-1 _ +1]}:
    z = 3.1415927*log(371*j + 417);
    z = z - int(z);
    z = 21.71819*abs(sin(49*z));
    z = z - int(z);
    return 2*z-1;
  }

function dsort(D,n,   m,i)
  {
    # Sorts {D[0..n-1]} in order of increasing value.}
    m = asort(D); # Sorts and renumbers {1..n}.
    if (n != m) { prog_error(("inconsistent sort")); }
    for (i = 0; i < n; i++) { D[i] = D[i+1]; }
    delete D[n];
  }
  
function dot(phi,n,i,j,   sum,k)
  { 
    # Dot product of rows {i,j} of the matrix {phi} with {n} columns.
    sum = 0;
    for (k = 0; k < n; k++) { sum += phi[i,k]*phi[j,k]; }
    return sum/n;
  }

function dot_D(phi,D,n,i,j,   sum,k)
  { 
    # Dot product of rows {i,j} of the matrix {phi} with {n} columns, with 
    # row {i} transformed by the diagonal operator {D}.
    sum = 0;
    for (k = 0; k < n; k++) { sum += D[k]*phi[i,k]*phi[j,k]; }
    return sum/n;
  }

function dot_DD(phi,D,n,i,j,   sum,k)
  { 
    # Dot product of rows {i,j} of the matrix {phi} with {n} columns, with 
    # both rows transformed by the diagonal operator {D}.
    sum = 0;
    for (k = 0; k < n; k++) { sum += D[k]*D[k]*phi[i,k]*phi[j,k]; }
    return sum/n;
  }

function abs(x)
  { 
    if (x < 0) { return -x; } else { return x; }
  }

function num_error(i,msg)
  { 
    printf "** element %d: %s\n", i, msg > "/dev/stderr";
    exit 1;
  }

function arg_error(msg)
  { 
    printf "** argument error: %s\n", msg > "/dev/stderr";
    exit 1;
  }

function prog_error(msg)
  { 
    printf "** program error: %s\n", msg > "/dev/stderr";
    exit 1;
  }