View Javadoc
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 }