Calculus


1   /*
2    * JScience - Java(TM) Tools and Libraries for the Advancement of Sciences.
3    * Copyright (C) 2007 - JScience (http://jscience.org/)
4    * All rights reserved.
5    * 
6    * Permission to use, copy, modify, and distribute this software is
7    * freely granted, provided that this notice is preserved.
8    */
9   package org.jscience.mathematics.numbers;
10  
11  /**
12   * <p> This class holds utilities upon arrays of positive <code>long</code>.</p>
13   *     
14   * @author <a HREF="mailto:jean-marie@dautelle.com">Jean-Marie Dautelle</a>
15   * @version 3.3, January 14, 2006
16   */
17  final class Calculus {
18  
19      /**
20       * Default constructor (private for utilities).
21       */
22      private Calculus() {
23      }
24  
25      static final long MASK_63 = 0x7FFFFFFFFFFFFFFFL;
26  
27      static final long MASK_32 = 0xFFFFFFFFL;
28  
29      static final long MASK_31 = 0x7FFFFFFFL;
30  
31      static final long MASK_8 = 0xFFL;
32          
33      /**
34       * Preconditions: xSize >= ySize
35       * @return z size
36       */
37      static int add(long[] x, int xSize, long[] y, int ySize, long[] z) {
38          long sum = 0;
39          int i = 0;
40          while (i < ySize) {
41              sum += x[i] + y[i];
42              z[i++] = sum & MASK_63;
43              sum >>>= 63;
44          }
45          while (true) {
46              if (sum == 0) {
47                  while (i < xSize) {
48                         z[i] = x[i++];
49                  }
50                  return xSize;
51              }
52              if (i == xSize) {
53                  z[xSize] = sum;
54                  return xSize + 1;
55              }
56              sum += x[i];
57              z[i++] = sum & MASK_63;
58              sum >>>= 63;
59          }
60      }
61  
62      /**
63       * Preconditions: x >= y
64       * @return z size
65       */
66      static int subtract(long[] x, int xSize, long[] y, int ySize,
67              long[] z) {
68          long diff = 0;
69          int i = 0;
70          while (i < ySize) {
71              diff += x[i] - y[i];
72              z[i++] = diff & MASK_63;
73              diff >>= 63; // Equals to -1 if borrow.
74          }
75          while (diff != 0) {
76              diff += x[i];
77              z[i++] = diff & MASK_63;
78              diff >>= 63; // Equals to -1 if borrow.
79          }        
80          // Copies rest of x to z.
81          while (i < xSize) {
82             z[i] = x[i++];
83          }
84          // Calculates size.
85          for (int j=xSize; j > 0;) { 
86              if (z[--j] != 0) return j + 1;
87          }
88          return 0;
89      }
90  
91      /**
92       * Preconditions: N/A
93       * @return 1, -1, 0 
94       */
95      static int compare(long[] left, long[] right, int size) {
96          for (int i = size; --i >= 0;) {
97              if (left[i] > right[i]) 
98                  return 1;
99              if (left[i] < right[i]) 
100                 return -1;
101         }
102         return 0;
103     }
104 
105     /**
106      * Preconditions: xSize != 0
107      * @return z size 
108      */
109     static int shiftLeft(int wordShift, int bitShift, long[] x,
110             int xSize, long[] z) {
111         final int shiftRight = 63 - bitShift;
112         int i = xSize;
113         int j = xSize + wordShift;
114         long tmp = x[--i];
115         long high = tmp >>> shiftRight;
116         if (high != 0) {
117             z[j] = high;
118         }
119         while (i > 0) {
120             z[--j] = ((tmp << bitShift) & MASK_63)
121                     | ((tmp = x[--i]) >>> shiftRight);
122         }
123         z[--j] = (tmp << bitShift) & MASK_63;
124         while (j > 0) {
125             z[--j] = 0;
126         }
127         return (high != 0) ? xSize + wordShift + 1 : xSize
128                 + wordShift;
129     }
130 
131     /**
132      * Preconditions: xSize > wordShift
133      * @return z size 
134      */
135     static int shiftRight(int wordShift, int bitShift, long[] x,
136             int xSize, long[] z) {
137         final int shiftLeft = 63 - bitShift;
138         int i = wordShift;
139         int j = 0;
140         long tmp = x[i];
141         while (i < xSize - 1) {
142             z[j++] = (tmp >>> bitShift) | ((tmp = x[++i]) << shiftLeft)
143                     & MASK_63;
144         }
145         tmp >>>= bitShift;
146         z[j] = tmp;
147         return (tmp != 0) ? j + 1 : j;
148     }
149 
150     /**
151      * Preconditions: y != 0, x != 0
152      * @return z size 
153      */
154     static int multiply(long[] x, int xSize, long y, long[] z) {
155         return multiply(x, xSize, y, z, 0);
156     }
157 
158     /**
159      * Preconditions: y != 0, xSize >= ySize
160      * @return z size 
161      */
162     static int multiply(long[] x, int xSize, long[] y, int ySize, long[] z) {
163         int zSize = 0;
164         for (int i = 0; i < ySize;) {
165             zSize = multiply(x, xSize, y[i], z, i++);
166         }
167         return zSize;
168     }
169 
170     // Multiplies by k, add to z if shift != 0
171     private static int multiply(long[] x, int xSize, long k, long[] z,
172             int shift) {
173         
174         final long kl = k & MASK_32; // 32 bits.
175         final long kh = k >> 32; // 31 bits
176         
177         long carry = 0; // 63 bits
178         for (int i = 0, j = shift; i < xSize;) {
179 
180             // Adds carry.
181             long zz = (shift == 0) ? carry : z[j] + carry; // 63 bits.
182             carry = zz >>> 63;
183             zz &= MASK_63; // 63 bits.
184 
185             // Splits words in [31 bits][32 bits]
186             final long w = x[i++];
187             final long wl = w & MASK_32; // 32 bits
188             final long wh = w >> 32; // 31 bits
189             
190             // Adds low.
191             long tmp = wl * kl; // 64 bits
192             carry += tmp >>> 63;
193             zz += tmp & MASK_63; // 64 bits.
194             carry += zz >>> 63;
195             zz &= MASK_63;
196             
197             // Adds middle.
198             tmp = wl * kh + wh * kl; // 64 bits.
199             carry += tmp >>> 31;
200             zz += (tmp << 32) & MASK_63; // 64 bits.
201             carry += zz >>> 63;
202             z[j++] = zz & MASK_63;
203             
204             // Adds high to carry.
205             carry += (wh * kh) << 1;
206             
207         }
208         int size = shift + xSize;
209         z[size] = carry;
210         if (carry == 0) return size;
211         return ++size;
212     }
213 
214     /**
215      * Preconditions: y is positive (31 bits).
216      * @return remainder 
217      */
218     static long divide(long[] x, int xSize, int y, long[] z) {
219         long r = 0;
220         for (int i = xSize; i > 0;) {
221             long w = x[--i];
222 
223             long wh = (r << 31) | (w >>> 32);
224             long qh = wh / y;
225             r = wh - qh * y;
226 
227             long wl = (r << 32) | (w & MASK_32);
228             long ql = wl / y;
229             r = wl - ql * y;
230 
231             z[i] = (qh << 32) | ql;
232         }
233         return r;
234     }
235     
236 }
A to Z: JavaDoc & Examples Daily Java News & Articles Open Source Projects Open Source Codes Free Computer Books Remove Frame
Popular Tags