Package cern.jet.math

Examples of cern.jet.math.Functions


      if (norm > 1)
         s = (int) Math.ceil(Num.log2(norm));
      else
         s = 0;

      Functions F = Functions.functions;    // alias F
      // B <-- B/2^s
      double v = 1.0 / Math.pow(2.0, s);
      if (v <= 0)
          throw new IllegalArgumentException ("   v <= 0");
      B.assign (F.mult(v));

      DoubleFactory2D fac = DoubleFactory2D.dense;
      final DoubleMatrix2D B0 = fac.identity(n);    // B^0 = I
      final DoubleMatrix2D B2 = alge.mult(B, B);    // B^2
      final DoubleMatrix2D B4 = alge.mult(B2, B2)// B^4

      DoubleMatrix2D T = B2.copy();          // T = work matrix
      DoubleMatrix2D W = B4.copy();          // W = work matrix
      W.assign (F.mult(cPade[9]));           // W <-- W*cPade[9]
      W.assign (T, F.plusMult(cPade[7]));    // W <-- W + T*cPade[7]
      DoubleMatrix2D U = alge.mult(B4, W);   // U <-- B4*W

      // T = B2.copy();
      W = B4.copy();
      W.assign (F.mult(cPade[5]));           // W <-- W*cPade[5]
      W.assign (T, F.plusMult(cPade[3]));    // W <-- W + T*cPade[3]
      W.assign (B0, F.plusMult(cPade[1]));   // W <-- W + B0*cPade[1]
      U.assign (W, F.plus);                  // U <-- U + W
      U = alge.mult(B, U);                   // U <-- B*U

      // T = B2.copy();
      W = B4.copy();
      W.assign (F.mult(cPade[8]));           // W <-- W*cPade[8]
      W.assign (T, F. plusMult(cPade[6]));   // W <-- W + T*cPade[6]
      DoubleMatrix2D V = alge.mult(B4, W);   // V <-- B4*W

      // T = B2.copy();
      W = B4.copy();
      W.assign (F.mult(cPade[4]));           // W <-- W*cPade[4]
      W.assign (T, F.plusMult(cPade[2]));    // W <-- W + T*cPade[2]
      W.assign (B0, F.plusMult(cPade[0]));   // W <-- W + B0*cPade[0]
      V.assign (W, F.plus);                  // V <-- V + W

      W = V.copy();
      W.assign(U, F.plus);                   // W = V + U, Padé numerator
      T = V.copy();
      T.assign(U, F.minus);                  // T = V - U, Padé denominator

      // Compute Padé approximant for exponential = W / T
      LUDecomposition lu = new LUDecomposition(T);
      B = lu.solve(W);

      if (false) {
         // This overflows for large |mu|
         // B <-- B^(2^s)
         for(int i = 0; i < s; i++)
            B = alge.mult(B, B);
         /*
         if (bal > 0) {
            throw new UnsupportedOperationException ("   balancing");
         } */
         v = Math.exp(mu);
         B.assign (F.mult(v));               // B <-- B*e^mu

      } else {
         // equivalent to B^(2^s) * e^mu, but only if no balancing
         double r = mu * v;
         r = Math.exp(r);
         B.assign (F.mult(r));
         for (int i = 0; i < s; i++)
            B = alge.mult(B, B);
      }

      return B;
View Full Code Here

TOP

Related Classes of cern.jet.math.Functions

Copyright © 2018 www.massapicom. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.