サブ　eigen.js　　　Last modified: Mar 25, 2004
// 固有値(e)，固有ベクトル(v)を求める関数 eigen.js

function eigen(a, n, ne, eps, e, v)
{
  var sw, i, j, k, k1, m, nea, nm1, nm2, nv
  var nva, a1, ajk, d, ei, eps1 = 1e-8, eps2 = 1e-8, f, q
  var r, s, sr, t, wj1, w, lw
  var pass1 = 0, pass2 = 0
  lw = new Array(n)
  w = new Array(6)
  for (i = 0; i < 6; i++) {
    w[i] = new Array(n)
  }
  nv = nva = nea = Math.abs(ne)
  if (nea == 0 || nea > n) return 1
  nm1 = n-1; nm2 = n-2
  if (eps <= 0) {
    eps = 1e-16
  }
  if (n == 1) {
    e[0] = a[0][0]
    if (nv != 0) {
      v[0][0] = 1
    }
    return 0
  }
  if (n == 2) {
    w[0][0] = a[0][1]; t = 0.5*(a[0][0]+a[1][1])
    r = a[0][0]*a[1][1]-a[0][1]*a[0][1]
    d = t*t-r
    q = Math.abs(t)+Math.sqrt(d)
    if (t < 0) {
      q = -q
    }
    t *= ne
    if (t >= 0) {
      e[0] = q
      if (nea == 2) {
        e[1] = r/q
      }
    }
    else {
      e[0] = r/q
      if (nea == 2) {
        e[1] = q
      }
    }
  }
  else {
    for (k = 1; k <= nm2; k++) {
      k1 = k-1
      s = 0
      for (i = k; i < n; i++) {
        s += a[k1][i]*a[k1][i]
      }
      w[0][k1] = 0
      if (s != 0) {
        sr = Math.sqrt(s)
        a1 = a[k1][k]
        if (a1 < 0) {
          sr = -sr
        }
        w[0][k1] = -sr
        r = 1/(s+a1*sr)
        a[k1][k] = a1+sr
        for (i = k; i < n; i++) {
          s = 0
          for (j = k; j <= i; j++) {
            s += a[j][i]*a[k1][j]
          }
          if (i != nm1) {
            for (j = i+1; j < n; j++) {
              s += a[i][j]*a[k1][j]
            }
          }
          w[0][i] = s*r
        }
        s = 0
        for (i = k; i < n; i++) {
          s += a[k1][i]*w[0][i]
        }
        t = 0.5*s*r
        for (i = k; i < n; i++) {
          w[0][i] -= t*a[k1][i]
        }
        for (j = k; j < n; j++) {
          wj1 = w[0][j]
          ajk = a[k1][j]
          for (i = j; i < n; i++) {
            a[j][i] -= a[k1][i]*wj1+w[0][i]*ajk
          }
        }
      }
    }
    w[0][nm2] = a[nm2][nm1]
    for (i = 0; i < n; i++) {
      w[5][i] = a[i][i]
    }
    r = Math.abs(w[5][0])+Math.abs(w[0][0])
    f = Math.abs(w[0][nm2])+Math.abs(w[5][nm1])
    r = Math.max(r, f)
    for (i = 1; i < nm1; i++) {
      f = Math.abs(w[0][i-1])+Math.abs(w[5][i])+Math.abs(w[0][i])
      r = Math.max(r, f)
    }
    eps1 = r*1e-16
    eps2 = r*eps
    for (i = 0; i < nm1; i++) {
      w[1][i] = w[0][i]*w[0][i]
    }
    if (ne < 0) {
      r = -r
    }
    f = r
    for (i = 0; i < nea; i++) {
      e[i] = -r
    }
    pass2 = 0
    for (k = 1; k <= nea; k++) {
      if (!pass2) {
        k1 = k-1
        d = e[k1]
      }
      pass2 = 0
      t = 0.5*(d+f)
      if (Math.abs(d-f) > eps2 && t != d && t != f) {
        j = 0
        i = 1
        for (;;) {
          if (!pass1) {
            q = w[5][i-1]-t
          }
          pass1 = 0
          if (q >= 0) {
            j++
          }
          if (q != 0) {
            if (++i > n) {
              break
            }
            q = w[5][i-1]-t-w[1][i-2]/q
            pass1 = 1
            continue
          }
          else {
            i += 2
          }
          if (i > n) {
            break
          }
        }
        if (ne < 0) {
          j = n-j
        }
        if (j < k) {
          f = t
          pass2 = 1
          k--
          continue
        }
        d = t
        m = Math.min(j, nea)
        for (i = k-1; i < m; i++) {
          e[i] = t
        }
        pass2 = 1
        k--
        continue
      }
      e[k1] = t
    }
  }
  if (nv != 0) {
    if (n == 2) {
      w[5][0] = a[0][0]
      w[5][1] = a[1][1]
    }
    w[0][nm1] = 0
    for (i = 0; i < nva; i++) {
      ei = e[i]
      for (j = 0; j < n; j++) {
        w[1][j] = w[5][j]-ei
        w[2][j] = w[0][j]
        v[i][j] = 1
      }
      sw = 0
      for (j = 0; j < nm1; j++) {
        if (Math.abs(w[1][j]) >= Math.abs(w[0][j])) {
          if (w[1][j] == 0) {
            w[1][j] = 1e-30
          }
          w[4][j] = w[0][j]/w[1][j]
          lw[j] = 0
          w[1][j+1] -= w[4][j]*w[2][j]
          w[3][j] = 0
        }
        else {
          w[4][j] = w[1][j]/w[0][j]
          lw[j] = 1
          w[1][j] = w[0][j]
          t = w[2][j]
          w[2][j] = w[1][j+1]
          w[3][j] = w[2][j+1]
          w[1][j+1] = t-w[4][j]*w[2][j]
          w[2][j+1] = -w[4][j]*w[3][j]
        }
      }
      if (w[1][nm1] == 0) {
        w[1][nm1] = 1e-30
      }
      if (i > 0 && Math.abs(e[i]-e[i-1]) < eps1) {
        for (j = 0; j < n; j++) {
          v[i][j] = Math.random()
        }
      }
      for (;;) {
        t = v[i][nm1]
        r = v[i][nm2]
        v[i][nm1] = t/w[1][nm1]
        v[i][nm2] = (r-w[2][nm2]*v[i][nm1])/w[1][nm2]
        for (k = nm2-1; k >= 0; k--) {
          v[i][k] = (v[i][k]-w[2][k]*v[i][k+1]-w[3][k]*v[i][k+2])/w[1][k]
        }
        if (sw == 0) {
          sw = 1
          for (j = 0; j < nm1; j++) {
            if (lw[j] == 1) {
              t = v[i][j]
              v[i][j] = v[i][j+1]
              v[i][j+1] = t-w[4][j]*v[i][j+1]
            }
            else {
              v[i][j+1] -= w[4][j]*v[i][j]
            }
          }
          continue
        }
        else
        break
      }
    }
    if (n > 2) {
      for (i = 0; i < nm2; i++) {
        w[0][i] *= -a[i][i+1]
      }
      for (i = 0; i < nva; i++) {
        for (k = nm2-1; k >= 0; k--) {
          r = w[0][k]
          if (r != 0) {
            r = 1/r
            s = 0
            for (j = k+1; j < n; j++) {
              s += a[k][j]*v[i][j]
            }
            r *= s
            for (j = k+1; j < n; j++) {
              v[i][j] -= r*a[k][j]
            }
          }
        }
      }
    }
    for (i = 0; i < nva; i++) {
      if (i == 0) {
        m = 1
      }
      else {
        if (Math.abs(e[i]-e[i-1]) >= eps1) {
          m = i+1
        }
        else {
          for (j = m-1; j < i; j++) {
            s = 0
            for (k = 0; k < n; k++) {
              s += v[j][k]*v[i][k]
            }
            for (k = 0; k < n; k++) {
              v[i][k] -= s*v[j][k]
            }
          }
        }
      }
      t = s = 0
      for (j = 0; j < n; j++) {
        s += v[i][j]*v[i][j]
      }
      if (s != 0) {
        t = Math.sqrt(1/s)
      }
      for (j = 0; j < n; j++) {
        v[i][j] *= t
      }
    }
  }
  return 0
}