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 * https://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
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.linear;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathIllegalArgumentException;
26 import org.hipparchus.util.ArithmeticUtils;
27 import org.hipparchus.util.CombinatoricsUtils;
28
29 /**
30 * This class implements inverses of Hilbert Matrices as
31 * {@link RealLinearOperator}.
32 */
33 public class InverseHilbertMatrix
34 implements RealLinearOperator {
35
36 /** The size of the matrix. */
37 private final int n;
38
39 /**
40 * Creates a new instance of this class.
41 *
42 * @param n Size of the matrix to be created.
43 */
44 public InverseHilbertMatrix(final int n) {
45 this.n = n;
46 }
47
48 /** {@inheritDoc} */
49 @Override
50 public int getColumnDimension() {
51 return n;
52 }
53
54 /**
55 * Returns the {@code (i, j)} entry of the inverse Hilbert matrix. Exact
56 * arithmetic is used; in case of overflow, an exception is thrown.
57 *
58 * @param i Row index (starts at 0).
59 * @param j Column index (starts at 0).
60 * @return The coefficient of the inverse Hilbert matrix.
61 */
62 public long getEntry(final int i, final int j) {
63 long val = i + j + 1;
64 long aux = CombinatoricsUtils.binomialCoefficient(n + i, n - j - 1);
65 val = ArithmeticUtils.mulAndCheck(val, aux);
66 aux = CombinatoricsUtils.binomialCoefficient(n + j, n - i - 1);
67 val = ArithmeticUtils.mulAndCheck(val, aux);
68 aux = CombinatoricsUtils.binomialCoefficient(i + j, i);
69 val = ArithmeticUtils.mulAndCheck(val, aux);
70 val = ArithmeticUtils.mulAndCheck(val, aux);
71 return ((i + j) & 1) == 0 ? val : -val;
72 }
73
74 /** {@inheritDoc} */
75 @Override
76 public int getRowDimension() {
77 return n;
78 }
79
80 /** {@inheritDoc} */
81 @Override
82 public RealVector operate(final RealVector x) {
83 if (x.getDimension() != n) {
84 throw new MathIllegalArgumentException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
85 x.getDimension(), n);
86 }
87 final double[] y = new double[n];
88 for (int i = 0; i < n; i++) {
89 double pos = 0.;
90 double neg = 0.;
91 for (int j = 0; j < n; j++) {
92 final double xj = x.getEntry(j);
93 final long coeff = getEntry(i, j);
94 final double daux = coeff * xj;
95 // Positive and negative values are sorted out in order to limit
96 // catastrophic cancellations (do not forget that Hilbert
97 // matrices are *very* ill-conditioned!
98 if (daux > 0.) {
99 pos += daux;
100 } else {
101 neg += daux;
102 }
103 }
104 y[i] = pos + neg;
105 }
106 return new ArrayRealVector(y, false);
107 }
108 }