public MatrixValue Function(ScalarValue j1, ScalarValue j2)
        {
            if (isNotHalf(j1))
                throw new ArgumentException("0, +-0.5, +-1, +-1.5, ...", j1.Re.ToString());

            if (isNotHalf(j2))
                throw new ArgumentException("0, +-0.5, +-1, +-1.5, ...", j2.Re.ToString());

            var l = 1;
            var M = new MatrixValue();

            for (var m1 = -j1.Re; m1 <= j1.Re; m1 += 1.0)
            {
                for (var m2 = -j2.Re; m2 <= j2.Re; m2 += 1.0)
                {
                    var m = m1 + m2;
                    var ja = j1.Re + j2.Re;

                    for (var j = Math.Abs(m); j <= ja; j += 1.0)
                    {
                        var v = CGCoefficients(j1.Re, j2.Re, j, m1, m2);
                        M[l, 1] = new ScalarValue(m);
                        M[l, 2] = new ScalarValue(m1);
                        M[l, 3] = new ScalarValue(m2);
                        M[l, 4] = new ScalarValue(j);
                        M[l, 5] = new ScalarValue(v);
                        l++;
                    }
                }
            }

            return M;
        }

        public MatrixValue Function(MatrixValue M)
        {
            if (M.DimensionX != 2 && M.DimensionY != 2)
                throw new YAMPMatrixDimensionException(2, M.DimensionX, M.DimensionY, M.DimensionX);

            var isTransposed = M.DimensionY != 2;

            if (isTransposed)
            {
                M = M.Transpose();
            }

            var m = new MatrixValue(2, M.DimensionX);

            for (var i = 1; i <= M.DimensionX; i++)
            {
                var x = M[1, i].Re;
                var y = M[2, i].Re;
                var phi = Math.Atan2(y, x);
                var r = x == 0.0 ? y : (y == 0.0 ? x : x / Math.Cos(phi));
                m[1, i] = new ScalarValue(r * Math.Cos(phi));
                m[2, i] = new ScalarValue(r * Math.Sin(phi));
            }

            return isTransposed ? m.Transpose() : m;
        }

        public MatrixValue Function(MatrixValue Y, MatrixValue x)
        {
            var X = new Double[x.Length];

            for (var i = 0; i < x.Length; i++)
            {
                X[i] = x[i + 1].Re;
            }

            if (!Y.IsVector)
            {
                var M = new MatrixValue();

                for (var i = 1; i <= Y.DimensionX; i++)
                {
                    var N = YMath.Histogram(Y.GetColumnVector(i), X);

                    for (var j = 1; j <= N.Length; j++)
                    {
                        M[j, i] = N[j];
                    }
                }

                return M;
            }

            return YMath.Histogram(Y, X);
        }

        public MatrixValue Function(MatrixValue M)
        {
            if (M.DimensionX != 3 && M.DimensionY != 3)
                throw new YAMPMatrixDimensionException(3, M.DimensionX, M.DimensionY, M.DimensionX);

            var isTransposed = M.DimensionY != 3;

            if (isTransposed)
            {
                M = M.Transpose();
            }

            var m = new MatrixValue(3, M.DimensionX);

            for (var i = 1; i <= M.DimensionX; i++)
            {
                var x = M[1, i].Re;
                var y = M[2, i].Re;
                var z = M[3, i].Re;
                var r = Math.Sqrt(x * x + y * y + z * z);
                var phi = Math.Atan2(y, x);
                var theta = Math.Acos(z / r);
                m[1, i] = new ScalarValue(r);
                m[2, i] = new ScalarValue(theta);
                m[3, i] = new ScalarValue(phi);
            }

            return isTransposed ? m.Transpose() : m;
        }

        /// <summary>
        /// Creates a new instance.
        /// </summary>
        /// <param name="samples">The Nx2 matrix containing the sample data.</param>
        public SplineInterpolation(MatrixValue samples)
            : base(samples)
        {
            a = new double[Np];
            h = new double[Np];

            for (int i = 2; i <= Np; i++)
                h[i - 1] = samples[i, 1].Re - samples[i - 1, 1].Re;

            if (Np > 2)
            {
                double[] sub = new double[Np - 1];
                double[] diag = new double[Np - 1];
                double[] sup = new double[Np - 1];

                for (int i = 2; i < Np; i++)
                {
                    int j = i - 1;
                    diag[j] = (h[j] + h[j + 1]) / 3;
                    sup[j] = h[j + 1] / 6;
                    sub[j] = h[j] / 6;
                    a[j] = (samples[i + 1, 2].Re - samples[i, 2].Re) / h[j + 1] - (samples[i, 2].Re - samples[i - 1, 2].Re) / h[j];

                }

                SolveTridiag(sub, diag, sup, ref a, Np - 2);
            }
        }

        public BarPlotValue Function(MatrixValue Y, MatrixValue x)
        {
            var bp = new BarPlotValue();
            var X = new Double[x.Length];

            for (var i = 0; i < x.Length; i++)
            {
                X[i] = x[i + 1].Re;
            }

            if (Y.IsVector)
            {
                bp.AddPoints(YMath.Histogram(Y, X));
            }
            else
            {
                var M = new MatrixValue();

                for (var i = 1; i <= Y.DimensionX; i++)
                {
                    var N = YMath.Histogram(Y.GetColumnVector(i), X);

                    for (var j = 1; j <= N.Length; j++)
                    {
                        M[j, i] = N[j];
                    }
                }

                bp.AddPoints(M);
            }

            return bp;
        }

        /// <summary>
        /// Creates the right QR decomposition (Givens or Householder) depending on the given matrix.
        /// </summary>
        /// <param name="A">The matrix to decompose.</param>
        /// <returns>The right QR decomposition implementation.</returns>
        public static QRDecomposition Create(MatrixValue A)
        {
            if (A.IsComplex)
                return new GivensDecomposition(A);

            return new HouseholderDecomposition(A);
        }

        public MatrixValue Function(ArgumentsValue values)
        {
            var m = new MatrixValue();

            for (var i = 1; i <= values.Length; i++)
            {
                var sy = m.DimensionY;
                var sx = m.DimensionX;

                if (values[i] is ScalarValue)
                {
                    var s = (ScalarValue)values[i];
                    m[sy + 1, sx + 1] = s.Clone();
                }
                else if (values[i] is MatrixValue)
                {
                    var n = (MatrixValue)values[i];

                    for (var l = 1; l <= n.DimensionX; l++)
                    {
                        for (var k = 1; k <= n.DimensionY; k++)
                        {
                            m[sy + k, sx + l] = n[k, l].Clone();
                        }
                    }
                }
                else
                {
                    throw new YAMPArgumentInvalidException(Name, values[i].Header, i);
                }
            }

            return m;
        }

        public ArgumentsValue Function(ScalarValue n)
        {
            var nn = n.GetIntegerOrThrowException("n", Name);
            int dim = nn + 1;

            if (dim < 2)
                throw new YAMPArgumentRangeException("n", 1.0);

            var X = new MatrixValue(dim, dim); // x = sin(phi) * cos(theta)
            var Y = new MatrixValue(dim, dim); // y = sin(phi) * sin(theta)
            var Z = new MatrixValue(dim, dim); // z = cos(phi)

            var stheta = Table(0.0, 2.0 * Math.PI, dim, Math.Sin);
            var ctheta = Table(0.0, 2.0 * Math.PI, dim, Math.Cos);
            var sphi = Table(0.0, Math.PI, dim, Math.Sin);
            var cphi = Table(0.0, Math.PI, dim, Math.Cos);

            for (var j = 0; j < dim; j++)
            {
                var col = j + 1;

                for (var i = 0; i < dim; i++)
                {
                    var row = i + 1;

                    X[row, col] = new ScalarValue(sphi[j] * ctheta[i]);
                    Y[row, col] = new ScalarValue(sphi[j] * stheta[i]);
                    Z[row, col] = new ScalarValue(cphi[j]);
                }
            }

            return new ArgumentsValue(X, Y, Z);
        }

        public override Value Perform(Value argument)
        {
            if (argument is ScalarValue)
            {
                return argument;
            }

            if (argument is MatrixValue)
            {
                var m = argument as MatrixValue;

                if (m.DimensionX == 1)
                {
                    return GetVectorMin(m.GetColumnVector(1));
                }
                else if (m.DimensionY == 1)
                {
                    return GetVectorMin(m.GetRowVector(1));
                }
                else
                {
                    var M = new MatrixValue(1, m.DimensionX);

                    for (var i = 1; i <= m.DimensionX; i++)
                    {
                        M[1, i] = GetVectorMin(m.GetColumnVector(i));
                    }

                    return M;
                }
            }

            throw new YAMPOperationInvalidException("min", argument);
        }

        public MatrixValue Function(MatrixValue M)
        {
            if (M.DimensionX != 3 && M.DimensionY != 3)
                throw new YAMPMatrixDimensionException(3, M.DimensionX, M.DimensionY, M.DimensionX);

            var isTransposed = M.DimensionY != 3;

            if (isTransposed)
            {
                M = M.Transpose();
            }

            var m = new MatrixValue(3, M.DimensionX);

            for (var i = 1; i <= M.DimensionX; i++)
            {
                var r = M[1, i].Re;
                var theta = M[2, i].Re;
                var phi = M[3, i].Re;
                var rt = r * Math.Sin(theta);
                m[1, i] = new ScalarValue(rt * Math.Cos(phi));
                m[2, i] = new ScalarValue(rt * Math.Sin(phi));
                m[3, i] = new ScalarValue(r * Math.Cos(theta));
            }

            return isTransposed ? m.Transpose() : m;
        }

 /// <summary>
 /// QR Decomposition, computed by Householder reflections.
 /// </summary>
 /// <param name="A">Rectangular matrix</param>
 /// <returns>Structure to access R and the Householder vectors and compute Q.</returns>
 protected QRDecomposition(MatrixValue A)
 {
     // Initialize.
     FullRank = true;
     m = A.DimensionY;
     n = A.DimensionX;
 }

        public ArgumentsValue Function(MatrixValue X, MatrixValue Y)
        {
            if (X.Length != Y.Length)
                throw new YAMPDifferentLengthsException(X.Length, Y.Length);

            var x1 = new ScalarValue();
            var y1 = new ScalarValue();
            var xy = new ScalarValue();
            var x2 = new ScalarValue();
            var slope = new ScalarValue();
            var offset = new ScalarValue();

            for (var i = 1; i <= X.Length; i++)
            {
                x1 = x1 + X[i];
                y1 = y1 + Y[i];
                xy = xy + X[i] * Y[i];
                x2 = x2 + X[i] * X[i];
            }

            var J = ((Double)X.Length * x2) - (x1 * x1);

            if (J.Re != 0.0)
            {
                slope = (((Double)X.Length * xy) - (x1 * y1)) / J.Re;
                offset = ((y1 * x2) - (x1 * xy)) / J.Re;
            }

            return new ArgumentsValue(slope, offset);
        }

        Plot2DValue Plot(IFunction f, Double minx, Double maxx, Double precision)
        {
            var cp = new Plot2DValue();
            var N = (Int32)((maxx - minx) / precision) + 1;
            var M = new MatrixValue(N, 2);
            var x = new ScalarValue(minx);

            for (var i = 0; i < N; i++)
            {
                var row = i + 1;
                var y = f.Perform(Context, x);
                M[row, 1] = x.Clone();

                if (y is ScalarValue)
                {
                    M[row, 2] = (ScalarValue)y;
                }
                else if (y is MatrixValue)
                {
                    var Y = (MatrixValue)y;

                    for (var j = 1; j <= Y.Length; j++)
                    {
                        M[row, j + 1] = Y[j];
                    }
                }

                x.Re += precision;
            }

            cp.AddPoints(M);
            return cp;
        }

        static MatrixValue Generate(Int32 n, Int32 m)
        {
            var distribution = new DiscreteUniformDistribution(Rng);
            var l = n * m;
            var M = new MatrixValue(n, m);
            var numbers = new List<Int32>();
            distribution.Alpha = 0;

            for (var i = 1; i <= l; i++)
            {
                numbers.Add(i);
            }

            for (var j = 1; j <= n; j++)
            {
                for (var i = 1; i <= m; i++)
                {
                    distribution.Beta = numbers.Count - 1;
                    var index = distribution.Next();
                    index = Math.Max(Math.Min(0, index), numbers.Count - 1);
                    M[j, i] = new ScalarValue(numbers[index]);
                    numbers.RemoveAt(index);
                }
            }

            return M;
        }

        public MatrixValue Function(MatrixValue x, MatrixValue y, ScalarValue n)
        {
            if (x.Length != y.Length)
                throw new YAMPDifferentLengthsException(x.Length, y.Length);

            var nn = n.GetIntegerOrThrowException("n", Name);
            var m = nn + 1;

            if (m < 2)
                throw new YAMPArgumentRangeException("n", 0.0);

            var M = new MatrixValue(x.Length, m);
            var b = new MatrixValue(x.Length, 1);

            for (var j = 1; j <= M.Rows; j++)
            {
                var el = ScalarValue.One;
                var z = x[j];

                for (var i = 1; i <= M.Columns; i++)
                {
                    M[j, i] = el;
                    el *= z;
                }

                b[j, 1] = y[j];
            }

            var qr = QRDecomposition.Create(M);
            return qr.Solve(b);
        }

        /// <summary>
        /// Performs the function - maps each entry of a matrix to a matrix.
        /// </summary>
        /// <param name="argument">Either a Scalar or a Matrix.</param>
        /// <returns>The scalar or matrix.</returns>
        public override Value Perform(Value argument)
        {
            if (argument is ScalarValue)
            {
                return GetValue((ScalarValue)argument);
            }
            else if (argument is MatrixValue)
            {
                var A = (MatrixValue)argument;
                var M = new MatrixValue(A.DimensionY, A.DimensionX);

                for (var j = 1; j <= A.DimensionY; j++)
                {
                    for (var i = 1; i <= A.DimensionX; i++)
                    {
                        M[j, i] = GetValue(A[j, i]);
                    }
                }

                return M;
            }
            else if (argument is ArgumentsValue)
            {
                throw new YAMPArgumentNumberException(Name, ((ArgumentsValue)argument).Length, 1);
            }

            throw new YAMPOperationInvalidException(Name, argument);
        }

 public MatrixValue Function(MatrixValue original, ScalarValue x)
 {
     var spline = new SplineInterpolation(original);
     var M = new MatrixValue(1, 2);
     M[1, 1].Re = x.Re;
     M[1, 2].Re = spline.ComputeValue(x.Re);
     return M;
 }

 public SurfacePlotValue Function(MatrixValue X, MatrixValue Y, MatrixValue Z)
 {
     var splot = new SurfacePlotValue();
     splot.IsSurf = true;
     splot.IsMesh = false;
     splot.AddPoints(X, Y, Z);
     return splot;
 }

        /// <summary>
        /// Check for symmetry, then construct the eigenvalue decomposition
        /// </summary>
        /// <param name="Arg">Square matrix</param>
        /// <returns>Structure to access D and V.</returns>
        public Eigenvalues(MatrixValue Arg)
        {
            var A = Arg.GetRealMatrix();
            n = Arg.DimensionX;
            V = new double[n][];

            for (int i = 0; i < n; i++)
                V[i] = new double[n];

            d = new double[n];
            e = new double[n];

            issymmetric = true;

            for (int j = 0; (j < n) && issymmetric; j++)
            {
                for (int i = 0; (i < n) && issymmetric; i++)
                    issymmetric = (A[i][j] == A[j][i]);
            }

            if (issymmetric)
            {
                for (int i = 0; i < n; i++)
                {
                    for (int j = 0; j < n; j++)
                        V[i][j] = A[i][j];
                }

                // Tridiagonalize.
                tred2();

                // Diagonalize.
                tql2();
            }
            else
            {
                H = new double[n][];

                for (int i2 = 0; i2 < n; i2++)
                    H[i2] = new double[n];

                ort = new double[n];

                for (int j = 0; j < n; j++)
                {
                    for (int i = 0; i < n; i++)
                        H[i][j] = A[i][j];
                }

                // Reduce to Hessenberg form.
                orthes();

                // Reduce Hessenberg to real Schur form.
                hqr2();
            }
        }