1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.apache.commons.math4.legacy.linear;
18
19 import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
20 import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
21
22 /**
23 * This class implements inverses of Hilbert Matrices as
24 * {@link RealLinearOperator}.
25 */
26 public class InverseHilbertMatrix
27 extends RealLinearOperator {
28
29 /** The size of the matrix. */
30 private final int n;
31
32 /**
33 * Creates a new instance of this class.
34 *
35 * @param n Size of the matrix to be created.
36 */
37 public InverseHilbertMatrix(final int n) {
38 this.n = n;
39 }
40
41 /** {@inheritDoc} */
42 @Override
43 public int getColumnDimension() {
44 return n;
45 }
46
47 /**
48 * Returns the {@code (i, j)} entry of the inverse Hilbert matrix. Exact
49 * arithmetic is used; in case of overflow, an exception is thrown.
50 *
51 * @param i Row index (starts at 0).
52 * @param j Column index (starts at 0).
53 * @return The coefficient of the inverse Hilbert matrix.
54 */
55 public long getEntry(final int i, final int j) {
56 long val = i + j + 1;
57 long aux = BinomialCoefficient.value(n + i, n - j - 1);
58 val = Math.multiplyExact(val, aux);
59 aux = BinomialCoefficient.value(n + j, n - i - 1);
60 val = Math.multiplyExact(val, aux);
61 aux = BinomialCoefficient.value(i + j, i);
62 val = Math.multiplyExact(val, aux);
63 val = Math.multiplyExact(val, aux);
64 return ((i + j) & 1) == 0 ? val : -val;
65 }
66
67 /** {@inheritDoc} */
68 @Override
69 public int getRowDimension() {
70 return n;
71 }
72
73 /** {@inheritDoc} */
74 @Override
75 public RealVector operate(final RealVector x) {
76 if (x.getDimension() != n) {
77 throw new DimensionMismatchException(x.getDimension(), n);
78 }
79 final double[] y = new double[n];
80 for (int i = 0; i < n; i++) {
81 double pos = 0.;
82 double neg = 0.;
83 for (int j = 0; j < n; j++) {
84 final double xj = x.getEntry(j);
85 final long coeff = getEntry(i, j);
86 final double daux = coeff * xj;
87 // Positive and negative values are sorted out in order to limit
88 // catastrophic cancellations (do not forget that Hilbert
89 // matrices are *very* ill-conditioned!
90 if (daux > 0.) {
91 pos += daux;
92 } else {
93 neg += daux;
94 }
95 }
96 y[i] = pos + neg;
97 }
98 return new ArrayRealVector(y, false);
99 }
100 }