|
Replies:
11
-
Last Post:
Jul 6, 2009 11:25 PM
by: mtj
|
Threads:
[
Previous
|
Next
]
|
|
|
|
|
|
IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Apr 11, 2009 2:55 PM
|
|
|
Forgive the large post of source code, but maybe this will help some of you that, like me, are fighting to find the missing Math functions from Java SE in ME. I ported the ones I needed from Sun's own FDLIBM C-library -- Sun should really provide the missing Math functionality as a standard extension to J2ME themselves for developers if it's not part of the device JVM. Google's Andriod is using the equivalent of Java5 -- Java ME needs to come out of the dark ages to compete.
Sun, or others, let me know if you'd like me to add more functions to this -- the exp, log, and pow were the most glaring omissions for implementing geographic projections like Mercator and Lambert.
I've battered this with millions of JUnit random number tests and it gives _exactly_ the same results as the Java SE JVM Math.exp, log, and pow methods since it is an exact port from the C-library underlying Java SE.
This is also available in SVN at http://sourceforge.net/projects/akmemobilemaps/ SVN.
package akme.mobile.util;
/* Portions copyright Sun Microsystems as below.
* This was ported from the FDLIBM C-library.
* ====================================================
* Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
*
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/**
* Additional Math functions missing from Java ME / J2ME.
*
* @see http://www.netlib.org/fdlibm/readme
*/
public class MathUtil {
/* Common constants. */
private static final double
zero = 0.0,
one = 1.0,
huge = 1.0e+300,
two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
private static final long
HI_MASK = 0xffffffff00000000L,
LO_MASK = 0x00000000ffffffffL;
private static final int
HI_SHIFT = 32;
/**
* Return Math.E to the exponent a.
* This in turn uses ieee7854_exp(double).
*/
public static final double exp(double a) {
return ieee754_exp(a);
}
/**
* Return the natural logarithm, ln(a), as it relates to Math.E.
* This in turn uses ieee7854_log(double).
*/
public static final double log(double a) {
return ieee754_log(a);
}
/**
* Return a to the power of b, sometime written as a ** b
* but not to be confused with the bitwise ^ operator.
* This in turn uses ieee7854_log(double).
*/
public static final double pow(double a, double b) {
return ieee754_pow(a, b);
}
/* __ieee754_exp(x)
* Returns the exponential of x.
*
* Method
* 1. Argument reduction:
* Reduce x to an r so that |r| ]] 7.09782712893383973096e+02 then exp(x) overflow
* if x < -7.45133219101941108420e+02 then exp(x) underflow
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
]]]]> HI_SHIFT); /* high word of x */
xsb = (hx>>31)&1; /* sign bit of x */
hx &= 0x7fffffff; /* high word of |x| */
/* filter out non-finite argument */
if(hx >= 0x40862E42) { /* if |x|>=709.78... */
if(hx>=0x7ff00000) {
lx = (int)((long)xl & LO_MASK); /* low word of x */
if(((hx&0xfffff)|lx)!=0)
return x+x; /* NaN */
else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
}
if(x > o_threshold) return huge*huge; /* overflow */
if(x < u_threshold) return twom1000*twom1000; /* 0x3fd62e42) { /* if |x| > 0.5 ln2 */
if(hx < 0x3FF0A2B2) { /*]]]]]one) return one+x;/* trigger inexact */
}
//else k = 0; // handled at declaration
/* x is now in primary range */
t = x*x;
c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
if(k==0) return one-((x*c)/(c-2.0)-x);
else y = one-((lo-(x*c)/(2.0-c))-hi);
yl = Double.doubleToLongBits(y);
if(k >= -1021) {
yl += ((long)k<<(20+HI_SHIFT)); /*]> HI_SHIFT); /* high word of x */
lx = (int)(xl & LO_MASK); /* low word of x */
k=0;
if (hx < 0x00100000) { /*> HI_SHIFT); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
//__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
x = Double.longBitsToDouble(((long)(hx|(i^0x3ff00000)) << HI_SHIFT) | (Double.doubleToLongBits(x) & LO_MASK));
k += (i>>20);
f = x-1.0;
if((0x000fffff&(2+hx))<3) { /*0) {
hfsq=0.5*f*f;
if(k==0) return f-(hfsq-s*(hfsq+R)); else
return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
} else {
if(k==0) return f-s*(f-R); else
return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
}
}
/* __ieee754_pow(x,y) return x**y
*
* n
* Method: Let x = 2 * (1+f)
* 1. Compute and return log2(x) in two pieces:
* log2(x) = w1 + w2,
* where w1 has 53-24 = 29 bit trailing zeros.
* 2. Perform y*log2(x) = n+y' by simulating muti-precision
* arithmetic, where |y'| 1) ** +INF is +INF
* 6. +-(|x| > 1) ** -INF is +0
* 7. +-(|x| < 1) ** +INF is +0
* 8. +-(|x| < 1) ** -INF is +INF
* 9. +-1 ** +-INF is NAN
* 10. +0 ** (+anything except 0, NAN) is +0
* 11. -0 ** (+anything except 0, NAN, odd integer) is +0
* 12. +0 ** (-anything except 0, NAN) is +INF
* 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
* 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
* 15. +INF ** (+anything except 0,NAN) is +INF
* 16. +INF ** (-anything except 0,NAN) is +0
* 17. -INF ** (anything) = -0 ** (-anything)
* 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
* 19. (-anything except 0 and inf) ** (non-integer) is NAN
*
* Accuracy:
* pow(x,y) returns x**y nearly rounded. In particular
* pow(integer,integer)
* always returns the correct integer provided it is
* representable.
*
* Constants :
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
]]]> (29+HI_SHIFT))^1;
//i1 = 1-i0;
hx = (int)(Double.doubleToLongBits(x) >> HI_SHIFT);
lx = (int)(Double.doubleToLongBits(x) & LO_MASK);
hy = (int)(Double.doubleToLongBits(y) >> HI_SHIFT);
ly = (int)(Double.doubleToLongBits(y) & LO_MASK);
ix = hx&0x7fffffff;
iy = hy&0x7fffffff;
/* y==zero: x**0 = 1 */
if((iy|ly)==0) return one;
/* +-NaN return x+y */
if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
return x+y;
/* determine if y is an odd int when x < 0
* yisint = 0 ... y is not an integer
* yisint = 1 ... y is an odd int
* yisint = 2 ... y is an even int
*/
=0x43400000) yisint = 2; /* even integer y */
else if(iy>=0x3ff00000) {
k = (iy>>20)-0x3ff; /* exponent */
if(k>20) {
j = ly>>(52-k);
if((j<<(52-k))=>(20-k);
if((j<<(20-k))== 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
return (hy>=0)? y: zero;
else /* (|x|<1)**-,+inf ==0) /* x >= +0 */
return Math.sqrt(x);
}
}
ax = Math.abs(x);
/* special value of x */
if(lx==0) {
if(ix==0x7ff00000||ix==0||ix==0x3ff00000){
z = ax; /*x is +-0,+-inf,+-1*/
if(hy<0) z = one/z; /*>31)+1;
/* (x<0)**(non-int) is NaN */
0x41e00000) { /* if |y| > 2**31 */
if(iy>0x43f00000){ /* if |y| > 2**64, must o/uflow */
if(ix=0x3ff00000) return (hy>0)? huge*huge:tiny*tiny;
}
/* over/underflow if x is not close to one */
if(ix<0x3fefffff) return (hy<0)? s*huge*huge:s*tiny*tiny;
if(ix>0x3ff00000) return (hy>0)? s*huge*huge:s*tiny*tiny;
/* now |1-x| is tiny > HI_SHIFT); }
n += ((ix)>>20)-0x3ff;
j = ix&0x000fffff;
/* determine interval */
ix = j|0x3ff00000; /* normalize ix */
if(j]]]]]>1)|0x20000000)+0x00080000+(k<<18);
t_h = Double.longBitsToDouble(((long)((int)((ix>>1)|0x20000000)+0x00080000+(k<<18)) << HI_SHIFT) | (Double.doubleToLongBits(t_h) & LO_MASK));
t_l = ax - (t_h-bp[k]);
s_l = v*((u-s_h*t_h)-s_h*t_l);
/*]]]> HI_SHIFT);
i = (int)(Double.doubleToLongBits(z) & LO_MASK);
if (j>=0x40900000) { /* z >= 1024 */
if(((j-0x40900000)|i)!=0) /* if z > 1024 */
return s*huge*huge; /* overflow */
else {
if(p_l+ovt>z-p_h) return s*huge*huge; /* overflow */
}
} else if((j&0x7fffffff)>=0x4090cc00 ) { /* z >20)-0x3ff;
n = 0;
if(i>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
n = j+(0x00100000>>(k+1));
k = ((n&0x7fffffff)>>20)-0x3ff; /* new k for n */
t = zero;
//__HI(t) = (n&~(0x000fffff>>k));
t = Double.longBitsToDouble(((long)((int)n&~(0x000fffff>>k)) << HI_SHIFT) | (Double.doubleToLongBits(t) & LO_MASK));
n = ((n&0x000fffff)|0x00100000)>>(20-k);
if(j<0) n = -n;
p_h -= t;
}
t = p_l+p_h;
//> HI_SHIFT);
j += (n<<20);
if((j>>20)> HI_SHIFT);
lx = (int)(Double.doubleToLongBits(x) & LO_MASK);
k = (hx&0x7ff00000)>>20; /* extract exponent */
if (k==0) { /* 0 or subnormal x */
if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
x *= two54;
hx = (int)(Double.doubleToLongBits(x) >> HI_SHIFT);
k = ((hx&0x7ff00000)>>20) - 54;
if (n< -50000) return tiny*x; /* 0x7fe) return huge*copysign(huge,x); /* overflow */
if (k > 0) /* normal result */
{
//__HI(x) = (hx&0x800fffff)|(k<<20);
x = Double.longBitsToDouble(((long)((int)(hx&0x800fffff)|(k<<20)) << HI_SHIFT) | (Double.doubleToLongBits(x) & LO_MASK));
return x;
}
if (k <= -54)
if (n > 50000) /* in case integer overflow in n+k */
return huge*copysign(huge,x); /*overflow*/
else return tiny*copysign(tiny,x); /*underflow*/
k += 54; /* subnormal result */
//__HI(x) = (hx&0x800fffff)|(k<<20);
x = Double.longBitsToDouble(((long)((int)(hx&0x800fffff)|(k<<20)) << HI_SHIFT) | (Double.doubleToLongBits(x) & LO_MASK));
return x*twom54;
}
/* 0.0D) ? -x : ((y > 0.0D && x < 0.0D) ? -x : x);
}
/*= 0.0D ? x : -x;
}
*/
}
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Apr 17, 2009 6:45 AM
in response to: kmashint
|
|
|
Hi,
I've downloaded the file from the sourceforge repository.
I have tested and the library function pow(x,y) does not work.
Test it with:
pow(10,3)
pow(10,4)
It fails instead of returning the righ value it returns:
pow(10,3)-> 1.1392378155556871E-305
pow(10,4)-> 9.113902524445497E-305
Instead pow(10,2) works right....
Any idea?
Thank you
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Apr 21, 2009 7:22 PM
in response to: kmashint
|
|
|
Thanks, I see that as well -- I literally ran millions of test cases but on random value and not simple integer exponents. I'll track down the issue and re-post.
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Apr 28, 2009 8:10 PM
in response to: kmashint
|
|
|
I'm splitting the common Java ME libraries from akmemobilemaps at sourceforge.net to this project at dev.java.net to be closer to Java ME developers looking for sample code. I found it far too painful to get started with Java ME development so hopefully this helps other along the way.
https://akmemobile.dev.java.net/
On a side note, it looks like the LWUIT project Sun is finally helping developers with a higher-level GUI library. It's a much need effort considering competition iPhone and Andriod.
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Apr 29, 2009 9:52 AM
in response to: kmashint
|
|
|
BTW this project has been approved and is in the M&E Community now.
Binky
mefeedback@mobileandembedded.org wrote:
> I'm splitting the common Java ME libraries from akmemobilemaps at sourceforge.net to this project at dev.java.net to be closer to Java ME developers looking for sample code. I found it far too painful to get started with Java ME development so hopefully this helps other along the way.
>
> https://akmemobile.dev.java.net/
>
> On a side note, it looks like the LWUIT project Sun is finally helping developers with a higher-level GUI library. It's a much need effort considering competition iPhone and Andriod.
> [Message sent by forum member 'kmashint' (kmashint)]
>
> http://forums.java.net/jive/thread.jspa?messageID=344245
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: feedback-unsubscribe@mobileandembedded.dev.java.net
> For additional commands, e-mail: feedback-help@mobileandembedded.dev.java.net
>
>
---------------------------------------------------------------------
To unsubscribe, e-mail: feedback-unsubscribe@mobileandembedded.dev.java.net
For additional commands, e-mail: feedback-help@mobileandembedded.dev.java.net
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Jun 8, 2009 11:05 PM
in response to: kmashint
|
|
|
Thanks for this information.
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Jun 18, 2009 9:33 PM
in response to: kmashint
|
|
|
this is great stuff, just what i needed!
any plans by the way to implement acos/asin?
thanks!
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Jun 19, 2009 7:53 PM
in response to: brenzosa
|
|
|
Glad to hear this is helpful for others. I can't believe these weren't eventually included in the Java ME already. I'm adding methods as people need them and doing test comparison against the full Java SE Math method where they should match exactly.
I'll put acos/asin next on the list and will post an Announcement to the akmemobile.dev.java.net project when they're ready. Please use the code there not what's posted here.
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Jul 6, 2009 7:21 AM
in response to: kmashint
|
|
|
Hi!
Thanks for your hard work providing the mathematical capabilities of J2SE to the J2ME users.
Do you have any eta on the acos/asin functions. It would be just what I need, and I'm not enough of a mathematician to figure it out myself how to calculate those.
Thanks once more.
|
|
|
|
|
|
|
|
|
Re: IEEE Math on Java ME / J2ME - the real exp, log, and pow methods
Posted:
Jul 6, 2009 11:25 PM
in response to: kmashint
|
|
|
Thank you!
I don't know how I missed it yesterday  I went through the svn and could not see them :P Maybe just tiredness in the afternoon, too little coffee and too much programming. Thanks again.
|
|
|
|
|
java.net RSS
