View Javadoc

1   // MathSupport.java, created Mon Feb  5 23:23:21 2001 by joewhaley
2   // Copyright (C) 2001-3 John Whaley <jwhaley@alum.mit.edu>
3   // Licensed under the terms of the GNU LGPL; see COPYING for details.
4   package joeq.Runtime;
5   
6   import joeq.Class.PrimordialClassLoader;
7   import joeq.Class.jq_Class;
8   import joeq.Class.jq_StaticField;
9   import joeq.Class.jq_StaticMethod;
10  
11  /*
12   * @author  John Whaley <jwhaley@alum.mit.edu>
13   * @version $Id: MathSupport.java 1456 2004-03-09 22:01:46Z jwhaley $
14   */
15  public abstract class MathSupport {
16  
17      public static boolean ucmp(/*unsigned*/int a, /*unsigned*/int b) {
18          if ((b < 0) && (a >= 0)) return true;
19          if ((b >= 0) && (a < 0)) return false;
20          return a<b;
21      }
22      public static boolean ulcmp(/*unsigned*/long a, /*unsigned*/long b) {
23          if ((b < 0L) && (a >= 0L)) return true;
24          if ((b >= 0L) && (a < 0L)) return false;
25          return a<b;
26      }
27      public static /*unsigned*/int udiv(/*unsigned*/int a, /*unsigned*/int b) {
28          if (b < 0) {
29              if (a < 0) {
30                  if (a >= b) return 1;
31              }
32              return 0;
33          }
34          if (a >= 0) return a/b;
35          // hard case
36          int a2 = a >>> 1;
37          int d = a2/b;
38          int m = ((a2%b) << 1) + (a&1);
39          return (d<<1) + m/b;
40      }
41      
42      public static /*unsigned*/int urem(/*unsigned*/int a, /*unsigned*/int b) {
43          if (b < 0) {
44              if (a < 0) {
45                  if (a >= b) return a-b;
46              }
47              return a;
48          }
49          if (a >= 0) return a%b;
50          // hard case
51          int a2 = a >>> 1;
52          int m = ((a2%b) << 1) + (a&1);
53          return m%b;
54      }
55      
56      public static /*unsigned*/long umul(/*unsigned*/int a, /*unsigned*/int b) {
57          char a_1 = (char)a;
58          char a_2 = (char)(a>>16);
59          char b_1 = (char)b;
60          char b_2 = (char)(b>>16);
61          int r1;
62          long result = 0L;
63          r1 = a_1*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)    ; r1 -= Integer.MAX_VALUE; }
64          result += ((long)r1)   ;
65          r1 = a_2*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<16; r1 -= Integer.MAX_VALUE; }
66          result += ((long)r1)<<16;
67          r1 = a_1*b_2; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<16; r1 -= Integer.MAX_VALUE; }
68          result += ((long)r1)<<16;
69          return result;
70      }
71      
72      public static long imul(int u, int v) {
73          /*unsigned*/int u1, u0, v1, v0, udiff, vdiff, high, mid, low;
74          /*unsigned*/int prodh, prodl, was;
75          boolean neg;
76  
77          u1 = HHALF(u);
78          u0 = LHALF(u);
79          v1 = HHALF(v);
80          v0 = LHALF(v);
81  
82          low = u0 * v0;
83  
84          /* This is the same small-number optimization as before. */
85          if (u1 == 0 && v1 == 0) return ((long)low)&0xFFFFFFFF;
86  
87          if (!ucmp(u1, u0)) {
88              udiff = u1 - u0; neg = false;
89          } else {
90              udiff = u0 - u1; neg = true;
91          }
92          if (!ucmp(v0, v1)) {
93              vdiff = v0 - v1;
94          } else {
95              vdiff = v1 - v0; neg = !neg;
96          }
97          mid = udiff * vdiff;
98  
99          high = u1 * v1;
100 
101         /* prod = (high << 2N) + (high << N); */
102         prodh = high + HHALF(high);
103         prodl = LHUP(high);
104 
105         /* if (neg) prod -= mid << N; else prod += mid << N; */
106         if (neg) {
107             was = prodl;
108             prodl -= LHUP(mid);
109             prodh -= HHALF(mid) + (ucmp(was, prodl)?1:0);
110         } else {
111             was = prodl;
112             prodl += LHUP(mid);
113             prodh += HHALF(mid) + (ucmp(prodl, was)?1:0);
114         }
115 
116         /* prod += low << N */
117         was = prodl;
118         prodl += LHUP(low);
119         prodh += HHALF(low) + (ucmp(prodl, was)?1:0);
120         /* ... + low; */
121         if (ucmp((prodl += low), low))
122             prodh++;
123 
124         /* return 4N-bit product */
125         return COMBINEQ(prodh, prodl);
126     }
127     
128     public static long lmul(long a, long b) {
129         long u, v, low, prod;
130         /*unsigned*/int high, mid, udiff, vdiff;
131         boolean negall, negmid;
132         if (a >= 0) { u = a; negall = false; }
133         else { u = -a; negall = true; }
134         if (b >= 0) { v = b; }
135         else { v = -b; negall = !negall; }
136 
137         /*unsigned*/int u1, u0, v1, v0;
138         u1 = HHALFQ(u);
139         u0 = LHALFQ(u);
140         v1 = HHALFQ(v);
141         v0 = LHALFQ(v);
142         if (u1 == 0 && v1 == 0) {
143             prod = imul(u0, v0);
144         } else {
145             /*
146              * Compute the three intermediate products, remembering
147              * whether the middle term is negative.  We can discard
148              * any upper bits in high and mid, so we can use native
149              * u_long * u_long => u_long arithmetic.
150              */
151             low = imul(u0, v0);
152             if (!ucmp(u1, u0)) {
153                 negmid = false; udiff = u1 - u0;
154             } else {
155                 negmid = true; udiff = u0 - u1;
156             }
157             if (!ucmp(v0, v1)) {
158                 vdiff = v0 - v1;
159             } else {
160                 vdiff = v1 - v0; negmid = !negmid;
161             }
162             mid = udiff * vdiff;
163             
164             high = (/*unsigned*/ int)umul(u1, v1);
165 
166             /*
167              * Assemble the final product.
168              */
169             prod = COMBINEQ(high + (negmid ? -mid : mid) + LHALFQ(low) + HHALFQ(low), HHALFQ(low));
170         }
171         return (negall ? -prod : prod);
172     }
173     
174     public static long ulmul(long a, long b) {
175         char a_1 = (char)a;
176         char a_2 = (char)(a>>16);
177         char a_3 = (char)(a>>32);
178         char a_4 = (char)(a>>48);
179         char b_1 = (char)b;
180         char b_2 = (char)(b>>16);
181         char b_3 = (char)(b>>32);
182         char b_4 = (char)(b>>48);
183         int r1;
184         long result = 0L;
185         r1 = a_1*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)    ; r1 -= Integer.MAX_VALUE; }
186         result += ((long)r1)   ;
187         r1 = a_2*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<16; r1 -= Integer.MAX_VALUE; }
188         result += ((long)r1)<<16;
189         r1 = a_1*b_2; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<16; r1 -= Integer.MAX_VALUE; }
190         result += ((long)r1)<<16;
191         r1 = a_3*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<32; r1 -= Integer.MAX_VALUE; }
192         result += ((long)r1)<<32;
193         r1 = a_2*b_2; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<32; r1 -= Integer.MAX_VALUE; }
194         result += ((long)r1)<<32;
195         r1 = a_1*b_3; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<32; r1 -= Integer.MAX_VALUE; }
196         result += ((long)r1)<<32;
197         r1 = a_4*b_1; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<48; r1 -= Integer.MAX_VALUE; }
198         result += ((long)r1)<<48;
199         r1 = a_3*b_2; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<48; r1 -= Integer.MAX_VALUE; }
200         result += ((long)r1)<<48;
201         r1 = a_2*b_3; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<48; r1 -= Integer.MAX_VALUE; }
202         result += ((long)r1)<<48;
203         r1 = a_1*b_4; if (r1 < 0) { result += ((long)Integer.MAX_VALUE)<<48; r1 -= Integer.MAX_VALUE; }
204         result += ((long)r1)<<48;
205         return result;
206     }
207     
208     public static long ldiv(long uq, long vq) {
209         boolean neg = false;
210         if (uq < 0L) { uq = -uq; neg = !neg; }
211         if (vq < 0L) { vq = -vq; neg = !neg; }
212         uq = uldivrem(uq, vq, false);
213         if (neg) uq = -uq;
214         return uq;
215     }
216 
217     public static long lrem(long uq, long vq) {
218         boolean neg = false;
219         if (uq < 0L) { uq = -uq; neg = !neg; }
220         if (vq < 0L) { vq = -vq; /*neg = !neg;*/ }
221         uq = uldivrem(uq, vq, true);
222         if (neg) uq = -uq;
223         return uq;
224     }
225 
226     private static int HHALFQ(long v) { return (int)(v>>32); }
227     private static int LHALFQ(long v) { return (int)v; }
228     private static char HHALF(int v) { return (char)(v>>16); }
229     private static char LHALF(int v) { return (char)v; }
230     private static int LHUP(int v) { return v<<16; }
231     private static /*unsigned*/int COMBINE(char hi, char lo) { return (hi<<16) | lo; }
232     private static /*unsigned*/long COMBINEQ(/*unsigned*/ int hi, /*unsigned*/ int lo) {
233         return (((long)hi)<<32) | (((long)lo) & 0xFFFFFFFFL);
234     }
235 
236     /*
237      * Multiprecision divide.  This algorithm is from Knuth vol. 2 (2nd ed),
238      * section 4.3.1, pp. 257--259.
239      *
240      * NOTE: the version here is adapted from NetBSD C source code (author unknown).
241      */
242     public static /*unsigned*/long uldivrem(/*unsigned*/long uq, /*unsigned*/long vq, boolean rem) {
243         final int HALF_BITS = 16;
244         final int B = 1 << HALF_BITS;
245         
246         char u[], v[], q[];
247         char v1, v2;
248         /*unsigned*/int qhat, rhat, t;
249         int m, n, d, j, i;
250         int ui=0, vi=0, qi=0;
251         if (vq == 0) throw new ArithmeticException();
252         if (ulcmp(uq, vq)) {
253             if (rem) return uq;
254             return 0L;
255         }
256         u = new char[5];
257         v = new char[5];
258         q = new char[5];
259         /*
260          * Break dividend and divisor into digits in base B, then
261          * count leading zeros to determine m and n.  When done, we
262          * will have:
263          *      u = (u[1]u[2]...u[m+n]) sub B
264          *      v = (v[1]v[2]...v[n]) sub B
265          *      v[1] != 0
266          *      1 < n <= 4 (if n = 1, we use a different division algorithm)
267          *      m >= 0 (otherwise u < v, which we already checked)
268          *      m + n = 4
269          * and thus
270          *      m = 4 - n <= 2
271          */
272         u[0] = 0;
273         u[1] = HHALF(HHALFQ(uq));
274         u[2] = LHALF(HHALFQ(uq));
275         u[3] = HHALF(LHALFQ(uq));
276         u[4] = LHALF(LHALFQ(uq));
277         v[1] = HHALF(HHALFQ(vq));
278         v[2] = LHALF(HHALFQ(vq));
279         v[3] = HHALF(LHALFQ(vq));
280         v[4] = LHALF(LHALFQ(vq));
281         for (n = 4; v[vi+1] == 0; vi++) {
282             if (--n == 1) {
283                 /*unsigned*/int rbj;    /* r*B+u[j] (not root boy jim) */
284                 char q1, q2, q3, q4;
285 
286                 /*
287                  * Change of plan, per exercise 16.
288                  *      r = 0;
289                  *      for j = 1..4:
290                  *              q[j] = floor((r*B + u[j]) / v),
291                  *              r = (r*B + u[j]) % v;
292                  * We unroll this completely here.
293                  */
294                 t = v[vi+2];    /* nonzero, by definition */
295                 q1 = (char)udiv(u[ui+1], t);
296                 rbj = COMBINE((char)urem(u[ui+1], t), u[ui+2]);
297                 q2 = (char)udiv(rbj, t);
298                 rbj = COMBINE((char)urem(rbj, t), u[ui+3]);
299                 q3 = (char)udiv(rbj, t);
300                 rbj = COMBINE((char)urem(rbj, t), u[ui+4]);
301                 q4 = (char)udiv(rbj, t);
302                 if (rem) return (long)urem(rbj, t);
303                 return COMBINEQ(COMBINE(q1, q2), COMBINE(q3, q4));
304             }
305         }
306 
307         /*
308          * By adjusting q once we determine m, we can guarantee that
309          * there is a complete four-digit quotient at &qspace[1] when
310          * we finally stop.
311          */
312         for (m = 4 - n; u[ui+1] == 0; ++ui)
313             m--;
314         for (i = 4 - m; --i >= 0;)
315             q[qi+i] = 0;
316         qi += 4 - m;
317 
318         /*
319          * Here we run Program D, translated from MIX to C and acquiring
320          * a few minor changes.
321          *
322          * D1: choose multiplier 1 << d to ensure v[1] >= B/2.
323          */
324         d = 0;
325         for (t = v[vi+1]; ucmp(t, B/2); t <<= 1)
326             d++;
327         if (d > 0) {
328             shl(u, ui, m + n, d);               /* u <<= d */
329             shl(v, vi+1, n - 1, d);             /* v <<= d */
330         }
331         /*
332          * D2: j = 0.
333          */
334         j = 0;
335         v1 = v[vi+1];   /* for D3 -- note that v[1..n] are constant */
336         v2 = v[vi+2];   /* for D3 */
337         do {
338             char uj0, uj1, uj2;
339                 
340             /*
341              * D3: Calculate qhat (\^q, in TeX notation).
342              * Let qhat = min((u[j]*B + u[j+1])/v[1], B-1), and
343              * let rhat = (u[j]*B + u[j+1]) mod v[1].
344              * While rhat < B and v[2]*qhat > rhat*B+u[j+2],
345              * decrement qhat and increase rhat correspondingly.
346              * Note that if rhat >= B, v[2]*qhat < rhat*B.
347              */
348             uj0 = u[ui+j + 0];  /* for D3 only -- note that u[j+...] change */
349             uj1 = u[ui+j + 1];  /* for D3 only */
350             uj2 = u[ui+j + 2];  /* for D3 only */
351             boolean sim_goto = false;
352             if (uj0 == v1) {
353                 qhat = B;
354                 rhat = uj1;
355                 //goto qhat_too_big;
356                 sim_goto = true;
357             } else {
358                 /*unsigned*/int n2 = COMBINE(uj0, uj1);
359                 qhat = udiv(n2, v1);
360                 rhat = urem(n2, v1);
361             }
362             while (sim_goto || (ucmp(COMBINE((char)rhat, uj2), (int)umul(v2, qhat)))) {
363         //qhat_too_big:
364                 sim_goto = false;
365                 qhat--;
366                 if ((rhat += v1) >= B)
367                     break;
368             }
369             /*
370              * D4: Multiply and subtract.
371              * The variable `t' holds any borrows across the loop.
372              * We split this up so that we do not require v[0] = 0,
373              * and to eliminate a final special case.
374              */
375             for (t = 0, i = n; i > 0; i--) {
376                 t = u[ui+i + j] - (int)umul(v[vi+i], qhat) - t;
377                 u[ui+i + j] = LHALF(t);
378                 t = (B - HHALF(t)) & (B - 1);
379             }
380             t = u[ui+j] - t;
381             u[ui+j] = LHALF(t);
382             /*
383              * D5: test remainder.
384              * There is a borrow if and only if HHALF(t) is nonzero;
385              * in that (rare) case, qhat was too large (by exactly 1).
386              * Fix it by adding v[1..n] to u[j..j+n].
387              */
388             if (HHALF(t) != 0) {
389                 qhat--;
390                 for (t = 0, i = n; i > 0; i--) { /* D6: add back. */
391                     t += u[ui+i + j] + v[vi+i];
392                     u[ui+i + j] = LHALF(t);
393                     t = HHALF(t);
394                 }
395                 u[ui+j] = LHALF(u[ui+j] + t);
396             }
397             q[qi+j] = (char)qhat;
398         } while (++j <= m);             /* D7: loop on j. */
399 
400         /*
401          * If caller wants the remainder, we have to calculate it as
402          * u[m..m+n] >>> d (this is at most n digits and thus fits in
403          * u[m+1..m+n], but we may need more source digits).
404          */
405         if (rem) {
406             if (d != 0) {
407                 for (i = m + n; i > m; --i)
408                     u[ui+i] = (char)((u[ui+i] >>> d) |
409                         LHALF(u[ui+i - 1] << (HALF_BITS - d)));
410                 u[ui+i] = 0;
411             }
412             return COMBINEQ(COMBINE(u[1], u[2]), COMBINE(u[3], u[4]));
413         }
414         return COMBINEQ(COMBINE(q[1], q[2]), COMBINE(q[3], q[4]));
415     }
416     private static void shl(char[] p, int off, int len, int sh)
417     {
418         int i;
419         final int HALF_BITS = 16;
420         for (i = 0; i < len; i++)
421                 p[off+i] = (char) (LHALF(p[off+i] << sh) | (p[off+i + 1] >>> (HALF_BITS - sh)));
422         p[off+i] = LHALF(p[off+i] << sh);
423     }
424 
425     // greatest double that can be rounded to an int
426     public static double maxint = (double)Integer.MAX_VALUE;
427     // least double that can be rounded to an int
428     public static double minint = (double)Integer.MIN_VALUE;
429   
430     // greatest double that can be rounded to a long
431     public static double maxlong = (double)Long.MAX_VALUE;
432     // least double that can be rounded to a long
433     public static double minlong = (double)Long.MIN_VALUE;
434     
435     public static final jq_Class _class;
436     public static final jq_StaticMethod _lmul;
437     public static final jq_StaticMethod _ldiv;
438     public static final jq_StaticMethod _lrem;
439     public static final jq_StaticField _maxint;
440     public static final jq_StaticField _minint;
441     public static final jq_StaticField _maxlong;
442     public static final jq_StaticField _minlong;
443     static {
444         _class = (jq_Class)PrimordialClassLoader.loader.getOrCreateBSType("Ljoeq/Runtime/MathSupport;");
445         _lmul = _class.getOrCreateStaticMethod("lmul", "(JJ)J");
446         _ldiv = _class.getOrCreateStaticMethod("ldiv", "(JJ)J");
447         _lrem = _class.getOrCreateStaticMethod("lrem", "(JJ)J");
448         _maxint = _class.getOrCreateStaticField("maxint", "D");
449         _minint = _class.getOrCreateStaticField("minint", "D");
450         _maxlong = _class.getOrCreateStaticField("maxlong", "D");
451         _minlong = _class.getOrCreateStaticField("minlong", "D");
452     }
453     
454 }