2 * Copyright (C) 1999-2000 Harri Porten (porten@kde.org)
3 * Copyright (C) 2007, 2008 Apple Inc. All Rights Reserved.
5 * This library is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU Lesser General Public
7 * License as published by the Free Software Foundation; either
8 * version 2 of the License, or (at your option) any later version.
10 * This library is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * Lesser General Public License for more details.
15 * You should have received a copy of the GNU Lesser General Public
16 * License along with this library; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 #include "MathObject.h"
25 #include "ObjectPrototype.h"
26 #include "Operations.h"
28 #include <wtf/Assertions.h>
29 #include <wtf/MathExtras.h>
30 #include <wtf/RandomNumber.h>
31 #include <wtf/RandomNumberSeed.h>
35 ASSERT_HAS_TRIVIAL_DESTRUCTOR(MathObject
);
37 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncAbs(ExecState
*);
38 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncACos(ExecState
*);
39 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncASin(ExecState
*);
40 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan(ExecState
*);
41 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan2(ExecState
*);
42 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCeil(ExecState
*);
43 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncCos(ExecState
*);
44 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncExp(ExecState
*);
45 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncFloor(ExecState
*);
46 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog(ExecState
*);
47 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncMax(ExecState
*);
48 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncMin(ExecState
*);
49 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncPow(ExecState
*);
50 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncRandom(ExecState
*);
51 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncRound(ExecState
*);
52 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncSin(ExecState
*);
53 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncSqrt(ExecState
*);
54 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncTan(ExecState
*);
55 static EncodedJSValue JSC_HOST_CALL
mathProtoFuncIMul(ExecState
*);
59 #include "MathObject.lut.h"
63 const ClassInfo
MathObject::s_info
= { "Math", &Base::s_info
, 0, ExecState::mathTable
, CREATE_METHOD_TABLE(MathObject
) };
65 /* Source for MathObject.lut.h
67 abs mathProtoFuncAbs DontEnum|Function 1
68 acos mathProtoFuncACos DontEnum|Function 1
69 asin mathProtoFuncASin DontEnum|Function 1
70 atan mathProtoFuncATan DontEnum|Function 1
71 atan2 mathProtoFuncATan2 DontEnum|Function 2
72 ceil mathProtoFuncCeil DontEnum|Function 1
73 cos mathProtoFuncCos DontEnum|Function 1
74 exp mathProtoFuncExp DontEnum|Function 1
75 floor mathProtoFuncFloor DontEnum|Function 1
76 log mathProtoFuncLog DontEnum|Function 1
77 max mathProtoFuncMax DontEnum|Function 2
78 min mathProtoFuncMin DontEnum|Function 2
79 pow mathProtoFuncPow DontEnum|Function 2
80 random mathProtoFuncRandom DontEnum|Function 0
81 round mathProtoFuncRound DontEnum|Function 1
82 sin mathProtoFuncSin DontEnum|Function 1
83 sqrt mathProtoFuncSqrt DontEnum|Function 1
84 tan mathProtoFuncTan DontEnum|Function 1
85 imul mathProtoFuncIMul DontEnum|Function 2
89 MathObject::MathObject(JSGlobalObject
* globalObject
, Structure
* structure
)
90 : JSNonFinalObject(globalObject
->vm(), structure
)
94 void MathObject::finishCreation(ExecState
* exec
, JSGlobalObject
* globalObject
)
96 Base::finishCreation(globalObject
->vm());
97 ASSERT(inherits(&s_info
));
99 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "E"), jsNumber(exp(1.0)), DontDelete
| DontEnum
| ReadOnly
);
100 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "LN2"), jsNumber(log(2.0)), DontDelete
| DontEnum
| ReadOnly
);
101 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "LN10"), jsNumber(log(10.0)), DontDelete
| DontEnum
| ReadOnly
);
102 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "LOG2E"), jsNumber(1.0 / log(2.0)), DontDelete
| DontEnum
| ReadOnly
);
103 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "LOG10E"), jsNumber(0.4342944819032518), DontDelete
| DontEnum
| ReadOnly
); // See ECMA-262 15.8.1.5
104 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "PI"), jsNumber(piDouble
), DontDelete
| DontEnum
| ReadOnly
);
105 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "SQRT1_2"), jsNumber(sqrt(0.5)), DontDelete
| DontEnum
| ReadOnly
);
106 putDirectWithoutTransition(exec
->vm(), Identifier(exec
, "SQRT2"), jsNumber(sqrt(2.0)), DontDelete
| DontEnum
| ReadOnly
);
109 bool MathObject::getOwnPropertySlot(JSCell
* cell
, ExecState
* exec
, PropertyName propertyName
, PropertySlot
&slot
)
111 return getStaticFunctionSlot
<JSObject
>(exec
, ExecState::mathTable(exec
), jsCast
<MathObject
*>(cell
), propertyName
, slot
);
114 bool MathObject::getOwnPropertyDescriptor(JSObject
* object
, ExecState
* exec
, PropertyName propertyName
, PropertyDescriptor
& descriptor
)
116 return getStaticFunctionDescriptor
<JSObject
>(exec
, ExecState::mathTable(exec
), jsCast
<MathObject
*>(object
), propertyName
, descriptor
);
119 // ------------------------------ Functions --------------------------------
121 EncodedJSValue JSC_HOST_CALL
mathProtoFuncAbs(ExecState
* exec
)
123 return JSValue::encode(jsNumber(fabs(exec
->argument(0).toNumber(exec
))));
126 EncodedJSValue JSC_HOST_CALL
mathProtoFuncACos(ExecState
* exec
)
128 return JSValue::encode(jsDoubleNumber(acos(exec
->argument(0).toNumber(exec
))));
131 EncodedJSValue JSC_HOST_CALL
mathProtoFuncASin(ExecState
* exec
)
133 return JSValue::encode(jsDoubleNumber(asin(exec
->argument(0).toNumber(exec
))));
136 EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan(ExecState
* exec
)
138 return JSValue::encode(jsDoubleNumber(atan(exec
->argument(0).toNumber(exec
))));
141 EncodedJSValue JSC_HOST_CALL
mathProtoFuncATan2(ExecState
* exec
)
143 double arg0
= exec
->argument(0).toNumber(exec
);
144 double arg1
= exec
->argument(1).toNumber(exec
);
145 return JSValue::encode(jsDoubleNumber(atan2(arg0
, arg1
)));
148 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCeil(ExecState
* exec
)
150 return JSValue::encode(jsNumber(ceil(exec
->argument(0).toNumber(exec
))));
153 EncodedJSValue JSC_HOST_CALL
mathProtoFuncCos(ExecState
* exec
)
155 return JSValue::encode(jsDoubleNumber(cos(exec
->argument(0).toNumber(exec
))));
158 EncodedJSValue JSC_HOST_CALL
mathProtoFuncExp(ExecState
* exec
)
160 return JSValue::encode(jsDoubleNumber(exp(exec
->argument(0).toNumber(exec
))));
163 EncodedJSValue JSC_HOST_CALL
mathProtoFuncFloor(ExecState
* exec
)
165 return JSValue::encode(jsNumber(floor(exec
->argument(0).toNumber(exec
))));
168 EncodedJSValue JSC_HOST_CALL
mathProtoFuncLog(ExecState
* exec
)
170 return JSValue::encode(jsDoubleNumber(log(exec
->argument(0).toNumber(exec
))));
173 EncodedJSValue JSC_HOST_CALL
mathProtoFuncMax(ExecState
* exec
)
175 unsigned argsCount
= exec
->argumentCount();
176 double result
= -std::numeric_limits
<double>::infinity();
177 for (unsigned k
= 0; k
< argsCount
; ++k
) {
178 double val
= exec
->argument(k
).toNumber(exec
);
179 if (std::isnan(val
)) {
183 if (val
> result
|| (!val
&& !result
&& !std::signbit(val
)))
186 return JSValue::encode(jsNumber(result
));
189 EncodedJSValue JSC_HOST_CALL
mathProtoFuncMin(ExecState
* exec
)
191 unsigned argsCount
= exec
->argumentCount();
192 double result
= +std::numeric_limits
<double>::infinity();
193 for (unsigned k
= 0; k
< argsCount
; ++k
) {
194 double val
= exec
->argument(k
).toNumber(exec
);
195 if (std::isnan(val
)) {
199 if (val
< result
|| (!val
&& !result
&& std::signbit(val
)))
202 return JSValue::encode(jsNumber(result
));
205 #if PLATFORM(IOS) && CPU(ARM_THUMB2)
207 static double fdlibmPow(double x
, double y
);
209 static ALWAYS_INLINE
bool isDenormal(double x
)
211 static const uint64_t signbit
= 0x8000000000000000ULL
;
212 static const uint64_t minNormal
= 0x0001000000000000ULL
;
213 return (bitwise_cast
<uint64_t>(x
) & ~signbit
) - 1 < minNormal
- 1;
216 static ALWAYS_INLINE
bool isEdgeCase(double x
)
218 static const uint64_t signbit
= 0x8000000000000000ULL
;
219 static const uint64_t infinity
= 0x7fffffffffffffffULL
;
220 return (bitwise_cast
<uint64_t>(x
) & ~signbit
) - 1 >= infinity
- 1;
223 static ALWAYS_INLINE
double mathPow(double x
, double y
)
225 if (!isDenormal(x
) && !isDenormal(y
)) {
226 double libmResult
= pow(x
,y
);
227 if (libmResult
|| isEdgeCase(x
) || isEdgeCase(y
))
230 return fdlibmPow(x
,y
);
235 ALWAYS_INLINE
double mathPow(double x
, double y
)
242 EncodedJSValue JSC_HOST_CALL
mathProtoFuncPow(ExecState
* exec
)
246 double arg
= exec
->argument(0).toNumber(exec
);
247 double arg2
= exec
->argument(1).toNumber(exec
);
249 if (std::isnan(arg2
))
250 return JSValue::encode(jsNaN());
251 if (std::isinf(arg2
) && fabs(arg
) == 1)
252 return JSValue::encode(jsNaN());
253 return JSValue::encode(jsNumber(mathPow(arg
, arg2
)));
256 EncodedJSValue JSC_HOST_CALL
mathProtoFuncRandom(ExecState
* exec
)
258 return JSValue::encode(jsDoubleNumber(exec
->lexicalGlobalObject()->weakRandomNumber()));
261 EncodedJSValue JSC_HOST_CALL
mathProtoFuncRound(ExecState
* exec
)
263 double arg
= exec
->argument(0).toNumber(exec
);
264 double integer
= ceil(arg
);
265 return JSValue::encode(jsNumber(integer
- (integer
- arg
> 0.5)));
268 EncodedJSValue JSC_HOST_CALL
mathProtoFuncSin(ExecState
* exec
)
270 return JSValue::encode(exec
->vm().cachedSin(exec
->argument(0).toNumber(exec
)));
273 EncodedJSValue JSC_HOST_CALL
mathProtoFuncSqrt(ExecState
* exec
)
275 return JSValue::encode(jsDoubleNumber(sqrt(exec
->argument(0).toNumber(exec
))));
278 EncodedJSValue JSC_HOST_CALL
mathProtoFuncTan(ExecState
* exec
)
280 return JSValue::encode(jsDoubleNumber(tan(exec
->argument(0).toNumber(exec
))));
283 EncodedJSValue JSC_HOST_CALL
mathProtoFuncIMul(ExecState
* exec
)
285 int32_t left
= exec
->argument(0).toInt32(exec
);
286 if (exec
->hadException())
287 return JSValue::encode(jsNull());
288 int32_t right
= exec
->argument(1).toInt32(exec
);
289 return JSValue::encode(jsNumber(left
* right
));
292 #if PLATFORM(IOS) && CPU(ARM_THUMB2)
294 // The following code is taken from netlib.org:
295 // http://www.netlib.org/fdlibm/fdlibm.h
296 // http://www.netlib.org/fdlibm/e_pow.c
297 // http://www.netlib.org/fdlibm/s_scalbn.c
299 // And was originally distributed under the following license:
302 * ====================================================
303 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
305 * Developed at SunSoft, a Sun Microsystems, Inc. business.
306 * Permission to use, copy, modify, and distribute this
307 * software is freely granted, provided that this notice
309 * ====================================================
312 * ====================================================
313 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
315 * Permission to use, copy, modify, and distribute this
316 * software is freely granted, provided that this notice
318 * ====================================================
321 /* __ieee754_pow(x,y) return x**y
324 * Method: Let x = 2 * (1+f)
325 * 1. Compute and return log2(x) in two pieces:
327 * where w1 has 53-24 = 29 bit trailing zeros.
328 * 2. Perform y*log2(x) = n+y' by simulating muti-precision
329 * arithmetic, where |y'|<=0.5.
330 * 3. Return x**y = 2**n*exp(y'*log2)
333 * 1. (anything) ** 0 is 1
334 * 2. (anything) ** 1 is itself
335 * 3. (anything) ** NAN is NAN
336 * 4. NAN ** (anything except 0) is NAN
337 * 5. +-(|x| > 1) ** +INF is +INF
338 * 6. +-(|x| > 1) ** -INF is +0
339 * 7. +-(|x| < 1) ** +INF is +0
340 * 8. +-(|x| < 1) ** -INF is +INF
341 * 9. +-1 ** +-INF is NAN
342 * 10. +0 ** (+anything except 0, NAN) is +0
343 * 11. -0 ** (+anything except 0, NAN, odd integer) is +0
344 * 12. +0 ** (-anything except 0, NAN) is +INF
345 * 13. -0 ** (-anything except 0, NAN, odd integer) is +INF
346 * 14. -0 ** (odd integer) = -( +0 ** (odd integer) )
347 * 15. +INF ** (+anything except 0,NAN) is +INF
348 * 16. +INF ** (-anything except 0,NAN) is +0
349 * 17. -INF ** (anything) = -0 ** (-anything)
350 * 18. (-anything) ** (integer) is (-1)**(integer)*(+anything**integer)
351 * 19. (-anything except 0 and inf) ** (non-integer) is NAN
354 * pow(x,y) returns x**y nearly rounded. In particular
355 * pow(integer,integer)
356 * always returns the correct integer provided it is
360 * The hexadecimal values are the intended ones for the following
361 * constants. The decimal values may be used, provided that the
362 * compiler will convert from decimal to binary accurately enough
363 * to produce the hexadecimal values shown.
366 #define __HI(x) *(1+(int*)&x)
367 #define __LO(x) *(int*)&x
371 dp_h
[] = { 0.0, 5.84962487220764160156e-01,}, /* 0x3FE2B803, 0x40000000 */
372 dp_l
[] = { 0.0, 1.35003920212974897128e-08,}, /* 0x3E4CFDEB, 0x43CFD006 */
376 two53
= 9007199254740992.0, /* 0x43400000, 0x00000000 */
380 two54
= 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
381 twom54
= 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
382 /* poly coefs for (3/2)*(log(x)-2s-2/3*s**3 */
383 L1
= 5.99999999999994648725e-01, /* 0x3FE33333, 0x33333303 */
384 L2
= 4.28571428578550184252e-01, /* 0x3FDB6DB6, 0xDB6FABFF */
385 L3
= 3.33333329818377432918e-01, /* 0x3FD55555, 0x518F264D */
386 L4
= 2.72728123808534006489e-01, /* 0x3FD17460, 0xA91D4101 */
387 L5
= 2.30660745775561754067e-01, /* 0x3FCD864A, 0x93C9DB65 */
388 L6
= 2.06975017800338417784e-01, /* 0x3FCA7E28, 0x4A454EEF */
389 P1
= 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
390 P2
= -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
391 P3
= 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
392 P4
= -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
393 P5
= 4.13813679705723846039e-08, /* 0x3E663769, 0x72BEA4D0 */
394 lg2
= 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
395 lg2_h
= 6.93147182464599609375e-01, /* 0x3FE62E43, 0x00000000 */
396 lg2_l
= -1.90465429995776804525e-09, /* 0xBE205C61, 0x0CA86C39 */
397 ovt
= 8.0085662595372944372e-0017, /* -(1024-log2(ovfl+.5ulp)) */
398 cp
= 9.61796693925975554329e-01, /* 0x3FEEC709, 0xDC3A03FD =2/(3ln2) */
399 cp_h
= 9.61796700954437255859e-01, /* 0x3FEEC709, 0xE0000000 =(float)cp */
400 cp_l
= -7.02846165095275826516e-09, /* 0xBE3E2FE0, 0x145B01F5 =tail of cp_h*/
401 ivln2
= 1.44269504088896338700e+00, /* 0x3FF71547, 0x652B82FE =1/ln2 */
402 ivln2_h
= 1.44269502162933349609e+00, /* 0x3FF71547, 0x60000000 =24b 1/ln2*/
403 ivln2_l
= 1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
405 inline double fdlibmScalbn (double x
, int n
)
410 k
= (hx
&0x7ff00000)>>20; /* extract exponent */
411 if (k
==0) { /* 0 or subnormal x */
412 if ((lx
|(hx
&0x7fffffff))==0) return x
; /* +-0 */
415 k
= ((hx
&0x7ff00000)>>20) - 54;
416 if (n
< -50000) return tiny
*x
; /*underflow*/
418 if (k
==0x7ff) return x
+x
; /* NaN or Inf */
420 if (k
> 0x7fe) return huge
*copysign(huge
,x
); /* overflow */
421 if (k
> 0) /* normal result */
422 {__HI(x
) = (hx
&0x800fffff)|(k
<<20); return x
;}
424 if (n
> 50000) /* in case integer overflow in n+k */
425 return huge
*copysign(huge
,x
); /*overflow*/
426 else return tiny
*copysign(tiny
,x
); /*underflow*/
428 k
+= 54; /* subnormal result */
429 __HI(x
) = (hx
&0x800fffff)|(k
<<20);
433 double fdlibmPow(double x
, double y
)
435 double z
,ax
,z_h
,z_l
,p_h
,p_l
;
436 double y1
,t1
,t2
,r
,s
,t
,u
,v
,w
;
437 int i0
,i1
,i
,j
,k
,yisint
,n
;
441 i0
= ((*(int*)&one
)>>29)^1; i1
=1-i0
;
442 hx
= __HI(x
); lx
= __LO(x
);
443 hy
= __HI(y
); ly
= __LO(y
);
444 ix
= hx
&0x7fffffff; iy
= hy
&0x7fffffff;
446 /* y==zero: x**0 = 1 */
447 if((iy
|ly
)==0) return one
;
449 /* +-NaN return x+y */
450 if(ix
> 0x7ff00000 || ((ix
==0x7ff00000)&&(lx
!=0)) ||
451 iy
> 0x7ff00000 || ((iy
==0x7ff00000)&&(ly
!=0)))
454 /* determine if y is an odd int when x < 0
455 * yisint = 0 ... y is not an integer
456 * yisint = 1 ... y is an odd int
457 * yisint = 2 ... y is an even int
461 if(iy
>=0x43400000) yisint
= 2; /* even integer y */
462 else if(iy
>=0x3ff00000) {
463 k
= (iy
>>20)-0x3ff; /* exponent */
466 if(static_cast<unsigned>(j
<<(52-k
))==ly
) yisint
= 2-(j
&1);
469 if((j
<<(20-k
))==iy
) yisint
= 2-(j
&1);
474 /* special value of y */
476 if (iy
==0x7ff00000) { /* y is +-inf */
477 if(((ix
-0x3ff00000)|lx
)==0)
478 return y
- y
; /* inf**+-1 is NaN */
479 else if (ix
>= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
480 return (hy
>=0)? y
: zero
;
481 else /* (|x|<1)**-,+inf = inf,0 */
482 return (hy
<0)?-y
: zero
;
484 if(iy
==0x3ff00000) { /* y is +-1 */
485 if(hy
<0) return one
/x
; else return x
;
487 if(hy
==0x40000000) return x
*x
; /* y is 2 */
488 if(hy
==0x3fe00000) { /* y is 0.5 */
489 if(hx
>=0) /* x >= +0 */
495 /* special value of x */
497 if(ix
==0x7ff00000||ix
==0||ix
==0x3ff00000){
498 z
= ax
; /*x is +-0,+-inf,+-1*/
499 if(hy
<0) z
= one
/z
; /* z = (1/|x|) */
501 if(((ix
-0x3ff00000)|yisint
)==0) {
502 z
= (z
-z
)/(z
-z
); /* (-1)**non-int is NaN */
504 z
= -z
; /* (x<0)**odd = -(|x|**odd) */
512 /* (x<0)**(non-int) is NaN */
513 if((n
|yisint
)==0) return (x
-x
)/(x
-x
);
515 s
= one
; /* s (sign of result -ve**odd) = -1 else = 1 */
516 if((n
|(yisint
-1))==0) s
= -one
;/* (-ve)**(odd int) */
519 if(iy
>0x41e00000) { /* if |y| > 2**31 */
520 if(iy
>0x43f00000){ /* if |y| > 2**64, must o/uflow */
521 if(ix
<=0x3fefffff) return (hy
<0)? huge
*huge
:tiny
*tiny
;
522 if(ix
>=0x3ff00000) return (hy
>0)? huge
*huge
:tiny
*tiny
;
524 /* over/underflow if x is not close to one */
525 if(ix
<0x3fefffff) return (hy
<0)? s
*huge
*huge
:s
*tiny
*tiny
;
526 if(ix
>0x3ff00000) return (hy
>0)? s
*huge
*huge
:s
*tiny
*tiny
;
527 /* now |1-x| is tiny <= 2**-20, suffice to compute
528 log(x) by x-x^2/2+x^3/3-x^4/4 */
529 t
= ax
-one
; /* t has 20 trailing zeros */
530 w
= (t
*t
)*(0.5-t
*(0.3333333333333333333333-t
*0.25));
531 u
= ivln2_h
*t
; /* ivln2_h has 21 sig. bits */
532 v
= t
*ivln2_l
-w
*ivln2
;
537 double ss
,s2
,s_h
,s_l
,t_h
,t_l
;
539 /* take care subnormal number */
541 {ax
*= two53
; n
-= 53; ix
= __HI(ax
); }
542 n
+= ((ix
)>>20)-0x3ff;
544 /* determine interval */
545 ix
= j
|0x3ff00000; /* normalize ix */
546 if(j
<=0x3988E) k
=0; /* |x|<sqrt(3/2) */
547 else if(j
<0xBB67A) k
=1; /* |x|<sqrt(3) */
548 else {k
=0;n
+=1;ix
-= 0x00100000;}
551 /* compute ss = s_h+s_l = (x-1)/(x+1) or (x-1.5)/(x+1.5) */
552 u
= ax
-bp
[k
]; /* bp[0]=1.0, bp[1]=1.5 */
557 /* t_h=ax+bp[k] High */
559 __HI(t_h
)=((ix
>>1)|0x20000000)+0x00080000+(k
<<18);
560 t_l
= ax
- (t_h
-bp
[k
]);
561 s_l
= v
*((u
-s_h
*t_h
)-s_h
*t_l
);
562 /* compute log(ax) */
564 r
= s2
*s2
*(L1
+s2
*(L2
+s2
*(L3
+s2
*(L4
+s2
*(L5
+s2
*L6
)))));
569 t_l
= r
-((t_h
-3.0)-s2
);
570 /* u+v = ss*(1+...) */
573 /* 2/(3log2)*(ss+...) */
577 z_h
= cp_h
*p_h
; /* cp_h+cp_l = 2/(3*log2) */
578 z_l
= cp_l
*p_h
+p_l
*cp
+dp_l
[k
];
579 /* log2(ax) = (ss+..)*2/(3*log2) = n + dp_h + z_h + z_l */
581 t1
= (((z_h
+z_l
)+dp_h
[k
])+t
);
583 t2
= z_l
-(((t1
-t
)-dp_h
[k
])-z_h
);
586 /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
589 p_l
= (y
-y1
)*t1
+y
*t2
;
594 if (j
>=0x40900000) { /* z >= 1024 */
595 if(((j
-0x40900000)|i
)!=0) /* if z > 1024 */
596 return s
*huge
*huge
; /* overflow */
598 if(p_l
+ovt
>z
-p_h
) return s
*huge
*huge
; /* overflow */
600 } else if((j
&0x7fffffff)>=0x4090cc00 ) { /* z <= -1075 */
601 if(((j
-0xc090cc00)|i
)!=0) /* z < -1075 */
602 return s
*tiny
*tiny
; /* underflow */
604 if(p_l
<=z
-p_h
) return s
*tiny
*tiny
; /* underflow */
608 * compute 2**(p_h+p_l)
613 if(i
>0x3fe00000) { /* if |z| > 0.5, set n = [z+0.5] */
614 n
= j
+(0x00100000>>(k
+1));
615 k
= ((n
&0x7fffffff)>>20)-0x3ff; /* new k for n */
617 __HI(t
) = (n
&~(0x000fffff>>k
));
618 n
= ((n
&0x000fffff)|0x00100000)>>(20-k
);
625 v
= (p_l
-(t
-p_h
))*lg2
+t
*lg2_l
;
629 t1
= z
- t
*(P1
+t
*(P2
+t
*(P3
+t
*(P4
+t
*P5
))));
630 r
= (z
*t1
)/(t1
-two
)-(w
+z
*w
);
634 if((j
>>20)<=0) z
= fdlibmScalbn(z
,n
); /* subnormal output */
635 else __HI(z
) += (n
<<20);