001: /*
002: * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
003: * Copyright (C) 2007 - JScience (http://jscience.org/)
004: * All rights reserved.
005: *
006: * Permission to use, copy, modify, and distribute this software is
007: * freely granted, provided that this notice is preserved.
008: */
009: package org.jscience.mathematics.number;
010:
011: import javolution.context.ObjectFactory;
012:
013: /**
014: * <p> This class holds utilities upon arrays of positive <code>long</code>.</p>
015: *
016: * @author <a href="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
017: * @version 3.3, January 14, 2006
018: */
019: final class Calculus {
020:
021: /**
022: * Default constructor (private for utilities).
023: */
024: private Calculus() {
025: }
026:
027: static final long MASK_63 = 0x7FFFFFFFFFFFFFFFL;
028:
029: static final long MASK_32 = 0xFFFFFFFFL;
030:
031: static final long MASK_31 = 0x7FFFFFFFL;
032:
033: static final long MASK_8 = 0xFFL;
034:
035: /**
036: * x += y
037: * @return x size
038: */
039: static int add(long[] x, int xSize, long y) {
040: long sum = x[0] + y;
041: x[0] = sum & MASK_63;
042: int i = 1;
043: sum >>>= 63;
044: while (sum != 0) {
045: if (i == xSize) {
046: x[xSize] = sum;
047: return xSize + 1;
048: }
049: sum += x[i];
050: x[i++] = sum & MASK_63;
051: sum >>>= 63;
052: }
053: return xSize;
054: }
055:
056: /**
057: * z = x + y
058: * Preconditions: xSize >= ySize
059: * @return z size
060: */
061: static int add(long[] x, int xSize, long[] y, int ySize, long[] z) {
062: long sum = 0;
063: int i = 0;
064: while (i < ySize) {
065: sum += x[i] + y[i];
066: z[i++] = sum & MASK_63;
067: sum >>>= 63;
068: }
069: while (true) {
070: if (sum == 0) {
071: while (i < xSize) {
072: z[i] = x[i++];
073: }
074: return xSize;
075: }
076: if (i == xSize) {
077: z[xSize] = sum;
078: return xSize + 1;
079: }
080: sum += x[i];
081: z[i++] = sum & MASK_63;
082: sum >>>= 63;
083: }
084: }
085:
086: /**
087: * z = x - y
088: * Preconditions: x >= y
089: * @return z size
090: */
091: static int subtract(long[] x, int xSize, long[] y, int ySize,
092: long[] z) {
093: long diff = 0;
094: int i = 0;
095: while (i < ySize) {
096: diff += x[i] - y[i];
097: z[i++] = diff & MASK_63;
098: diff >>= 63; // Equals to -1 if borrow.
099: }
100: while (diff != 0) {
101: diff += x[i];
102: z[i++] = diff & MASK_63;
103: diff >>= 63; // Equals to -1 if borrow.
104: }
105: // Copies rest of x to z.
106: while (i < xSize) {
107: z[i] = x[i++];
108: }
109: // Calculates size.
110: for (int j = xSize; j > 0;) {
111: if (z[--j] != 0)
112: return j + 1;
113: }
114: return 0;
115: }
116:
117: /**
118: * x.compare(y)
119: * Preconditions: xSize = ySize = size
120: * @return 1, -1, 0
121: */
122: static int compare(long[] x, long[] y, int size) {
123: for (int i = size; --i >= 0;) {
124: if (x[i] > y[i])
125: return 1;
126: if (x[i] < y[i])
127: return -1;
128: }
129: return 0;
130: }
131:
132: /**
133: * x << n
134: * Preconditions: xSize != 0
135: * @return z size
136: */
137: static int shiftLeft(int wordShift, int bitShift, long[] x,
138: int xSize, long[] z) {
139: final int shiftRight = 63 - bitShift;
140: int i = xSize;
141: int j = xSize + wordShift;
142: long tmp = x[--i];
143: long high = tmp >>> shiftRight;
144: if (high != 0) {
145: z[j] = high;
146: }
147: while (i > 0) {
148: z[--j] = ((tmp << bitShift) & MASK_63)
149: | ((tmp = x[--i]) >>> shiftRight);
150: }
151: z[--j] = (tmp << bitShift) & MASK_63;
152: while (j > 0) {
153: z[--j] = 0;
154: }
155: return (high != 0) ? xSize + wordShift + 1 : xSize + wordShift;
156: }
157:
158: /**
159: * x >> n
160: * Preconditions: xSize > wordShift
161: * @return z size
162: */
163: static int shiftRight(int wordShift, int bitShift, long[] x,
164: int xSize, long[] z) {
165: final int shiftLeft = 63 - bitShift;
166: int i = wordShift;
167: int j = 0;
168: long tmp = x[i];
169: while (i < xSize - 1) {
170: z[j++] = (tmp >>> bitShift) | ((tmp = x[++i]) << shiftLeft)
171: & MASK_63;
172: }
173: tmp >>>= bitShift;
174: z[j] = tmp;
175: return (tmp != 0) ? j + 1 : j;
176: }
177:
178: /**
179: * z = x * y
180: * Preconditions: y != 0, x != 0
181: * @return z size
182: */
183: static int multiply(long[] x, int xSize, long y, long[] z) {
184: return multiply(x, xSize, y, z, 0);
185: }
186:
187: /**
188: * z = x * y
189: * Preconditions: y != 0, xSize >= ySize
190: * @return z size
191: */
192: static int multiply(long[] x, int xSize, long[] y, int ySize,
193: long[] z) {
194: int zSize = 0;
195: for (int i = 0; i < ySize;) {
196: zSize = multiply(x, xSize, y[i], z, i++);
197: }
198: return zSize;
199: }
200:
201: // Multiplies by k, add to z if shift != 0
202: private static int multiply(long[] x, int xSize, long k, long[] z,
203: int shift) {
204:
205: final long kl = k & MASK_32; // 32 bits.
206: final long kh = k >> 32; // 31 bits
207:
208: long carry = 0; // 63 bits
209: for (int i = 0, j = shift; i < xSize;) {
210:
211: // Adds carry.
212: long zz = (shift == 0) ? carry : z[j] + carry; // 63 bits.
213: carry = zz >>> 63;
214: zz &= MASK_63; // 63 bits.
215:
216: // Splits words in [31 bits][32 bits]
217: final long w = x[i++];
218: final long wl = w & MASK_32; // 32 bits
219: final long wh = w >> 32; // 31 bits
220:
221: // Adds low.
222: long tmp = wl * kl; // 64 bits
223: carry += tmp >>> 63;
224: zz += tmp & MASK_63; // 64 bits.
225: carry += zz >>> 63;
226: zz &= MASK_63;
227:
228: // Adds middle.
229: tmp = wl * kh + wh * kl; // 64 bits.
230: carry += tmp >>> 31;
231: zz += (tmp << 32) & MASK_63; // 64 bits.
232: carry += zz >>> 63;
233: z[j++] = zz & MASK_63;
234:
235: // Adds high to carry.
236: carry += (wh * kh) << 1;
237:
238: }
239: int size = shift + xSize;
240: z[size] = carry;
241: if (carry == 0)
242: return size;
243: return ++size;
244: }
245:
246: /**
247: * z = x / y
248: * Preconditions: y is positive (31 bits).
249: * @return remainder
250: */
251: static long divide(long[] x, int xSize, int y, long[] z) {
252: long r = 0;
253: for (int i = xSize; i > 0;) {
254: long w = x[--i];
255:
256: long wh = (r << 31) | (w >>> 32);
257: long qh = wh / y;
258: r = wh - qh * y;
259:
260: long wl = (r << 32) | (w & MASK_32);
261: long ql = wl / y;
262: r = wl - ql * y;
263:
264: z[i] = (qh << 32) | ql;
265: }
266: return r;
267: }
268:
269: /**
270: * Multiplication logic (for concurrent context)
271: */
272: static final class MultiplyLogic implements Runnable {
273: private static final ObjectFactory<MultiplyLogic> FACTORY = new ObjectFactory<MultiplyLogic>() {
274: @Override
275: protected MultiplyLogic create() {
276: return new MultiplyLogic();
277: }
278: };
279: private LargeInteger _left, _right, _value;
280:
281: public static MultiplyLogic newInstance(LargeInteger left,
282: LargeInteger right) {
283: MultiplyLogic logic = FACTORY.object();
284: logic._left = left;
285: logic._right = right;
286: return logic;
287: }
288:
289: public void run() {
290: _value = _left.times(_right);// Recursive.
291: }
292:
293: public LargeInteger value() {
294: return _value;
295: }
296: };
297:
298: }
|