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.statistics.descriptive;
18
19 import java.math.BigInteger;
20 import java.nio.ByteBuffer;
21
22 /**
23 * A mutable 192-bit unsigned integer.
24 *
25 * <p>This is a specialised class to implement an accumulator of squared {@code long} values.
26 *
27 * @since 1.1
28 */
29 final class UInt192 {
30 /** Mask for the lower 32-bits of a long. */
31 private static final long MASK32 = 0xffff_ffffL;
32
33 // Data is stored using integers to allow efficient sum-with-carry addition
34
35 /** bits 32-1 (low 32-bits). */
36 private int f;
37 /** bits 64-33. */
38 private int e;
39 /** bits 96-65. */
40 private int d;
41 /** bits 128-97. */
42 private int c;
43 /** bits 192-129 (high 64-bits). */
44 private long ab;
45
46 /**
47 * Create an instance.
48 */
49 private UInt192() {
50 // No-op
51 }
52
53 /**
54 * Create an instance using a direct binary representation.
55 * This is package-private for testing.
56 *
57 * @param hi High 64-bits.
58 * @param mid Middle 64-bits.
59 * @param lo Low 64-bits.
60 */
61 UInt192(long hi, long mid, long lo) {
62 this.f = (int) lo;
63 this.e = (int) (lo >>> Integer.SIZE);
64 this.d = (int) mid;
65 this.c = (int) (mid >>> Integer.SIZE);
66 this.ab = hi;
67 }
68
69 /**
70 * Create an instance using a direct binary representation.
71 *
72 * @param ab bits 192-129 (high 64-bits).
73 * @param c bits 128-97.
74 * @param d bits 96-65.
75 * @param e bits 64-33.
76 * @param f bits 32-1.
77 */
78 private UInt192(long ab, int c, int d, int e, int f) {
79 this.ab = ab;
80 this.c = c;
81 this.d = d;
82 this.e = e;
83 this.f = f;
84 }
85
86 /**
87 * Create an instance. The initial value is zero.
88 *
89 * @return the instance
90 */
91 static UInt192 create() {
92 return new UInt192();
93 }
94
95 /**
96 * Adds the squared value {@code x * x}.
97 *
98 * @param x Value.
99 */
100 void addSquare(long x) {
101 final long lo = x * x;
102 final long hi = IntMath.squareHigh(x);
103
104 // Sum with carry.
105 long s = (lo & MASK32) + (f & MASK32);
106 f = (int) s;
107 s = (s >>> Integer.SIZE) + (lo >>> Integer.SIZE) + (e & MASK32);
108 e = (int) s;
109 s = (s >>> Integer.SIZE) + (hi & MASK32) + (d & MASK32);
110 d = (int) s;
111 s = (s >>> Integer.SIZE) + (hi >>> Integer.SIZE) + (c & MASK32);
112 c = (int) s;
113 ab += s >>> Integer.SIZE;
114 }
115
116 /**
117 * Adds the value.
118 *
119 * @param x Value.
120 */
121 void add(UInt192 x) {
122 // Avoid issues adding to itself
123 final int ff = x.f;
124 final int ee = x.e;
125 final int dd = x.d;
126 final int cc = x.c;
127 final long aabb = x.ab;
128 // Sum with carry.
129 long s = (ff & MASK32) + (f & MASK32);
130 f = (int) s;
131 s = (s >>> Integer.SIZE) + (ee & MASK32) + (e & MASK32);
132 e = (int) s;
133 s = (s >>> Integer.SIZE) + (dd & MASK32) + (d & MASK32);
134 d = (int) s;
135 s = (s >>> Integer.SIZE) + (cc & MASK32) + (c & MASK32);
136 c = (int) s;
137 ab += (s >>> Integer.SIZE) + aabb;
138 }
139
140
141 /**
142 * Multiply by the unsigned value.
143 * Any overflow bits are lost.
144 *
145 * @param x Value.
146 * @return the product
147 */
148 UInt192 unsignedMultiply(int x) {
149 final long xx = x & MASK32;
150 // Multiply with carry.
151 long product = xx * (f & MASK32);
152 final int ff = (int) product;
153 product = (product >>> Integer.SIZE) + xx * (e & MASK32);
154 final int ee = (int) product;
155 product = (product >>> Integer.SIZE) + xx * (d & MASK32);
156 final int dd = (int) product;
157 product = (product >>> Integer.SIZE) + xx * (c & MASK32);
158 final int cc = (int) product;
159 // Possible overflow here and bits are lost
160 final long aabb = (product >>> Integer.SIZE) + xx * ab;
161 return new UInt192(aabb, cc, dd, ee, ff);
162 }
163
164 /**
165 * Subtracts the value.
166 * Any overflow bits (negative result) are lost.
167 *
168 * @param x Value.
169 * @return the difference
170 */
171 UInt192 subtract(UInt128 x) {
172 // Difference with carry.
173 // Subtract common part.
174 long diff = (f & MASK32) - (x.lo32() & MASK32);
175 final int ff = (int) diff;
176 diff = (diff >> Integer.SIZE) + (e & MASK32) - (x.mid32() & MASK32);
177 final int ee = (int) diff;
178 diff = (diff >> Integer.SIZE) + (d & MASK32) - (x.hi64() & MASK32);
179 final int dd = (int) diff;
180 diff = (diff >> Integer.SIZE) + (c & MASK32) - (x.hi64() >>> Integer.SIZE);
181 final int cc = (int) diff;
182 // Possible overflow here and bits are lost containing info on the
183 // magnitude of the true negative value
184 final long aabb = (diff >> Integer.SIZE) + ab;
185 return new UInt192(aabb, cc, dd, ee, ff);
186 }
187
188 /**
189 * Convert to a BigInteger.
190 *
191 * @return the value
192 */
193 BigInteger toBigInteger() {
194 final ByteBuffer bb = ByteBuffer.allocate(Integer.BYTES * 6)
195 .putLong(ab)
196 .putInt(c)
197 .putInt(d)
198 .putInt(e)
199 .putInt(f);
200 // Sign is always positive. This works for zero.
201 return new BigInteger(1, bb.array());
202 }
203
204 /**
205 * Convert to a double.
206 *
207 * @return the value
208 */
209 double toDouble() {
210 final long h = hi64();
211 final long m = mid64();
212 final long l = lo64();
213 if (h == 0) {
214 return IntMath.uint128ToDouble(m, l);
215 }
216 // For correct rounding we use a sticky bit to represent magnitude
217 // lost from the low 64-bits. The result is scaled by 2^64.
218 return IntMath.uint128ToDouble(h, m | ((l == 0) ? 0 : 1)) * 0x1.0p64;
219 }
220
221 /**
222 * Convert to an {@code int}; throwing an exception if the value overflows an {@code int}.
223 *
224 * @return the value
225 * @throws ArithmeticException if the value overflows an {@code int}.
226 * @see Math#toIntExact(long)
227 */
228 int toIntExact() {
229 return Math.toIntExact(toLongExact());
230 }
231
232 /**
233 * Convert to a {@code long}; throwing an exception if the value overflows a {@code long}.
234 *
235 * @return the value
236 * @throws ArithmeticException if the value overflows a {@code long}.
237 */
238 long toLongExact() {
239 // Test if we have more than 63-bits
240 if ((ab | c | d) != 0 || e < 0) {
241 throw new ArithmeticException("long integer overflow");
242 }
243 return lo64();
244 }
245
246 /**
247 * Return the lower 64-bits as a {@code long} value.
248 *
249 * @return the low 64-bits
250 */
251 long lo64() {
252 return (f & MASK32) | ((e & MASK32) << Integer.SIZE);
253 }
254
255 /**
256 * Return the middle 64-bits as a {@code long} value.
257 *
258 * @return bits 128-65
259 */
260 long mid64() {
261 return (d & MASK32) | ((c & MASK32) << Integer.SIZE);
262 }
263
264 /**
265 * Return the higher 64-bits as a {@code long} value.
266 *
267 * @return bits 192-129
268 */
269 long hi64() {
270 return ab;
271 }
272 }