/ pow.s (emx+gcc) -- Copyright (c) 1991-1993 by Kolja Elsaesser / Modified 1992-1996 by Eberhard Mattes #include #define FUNC MATHSUFFIX1(pow) .globl LABEL(FUNC) .text ALIGN / double pow (double x, double y) #define cw_user -4(%ebp) #define cw_chop -2(%ebp) /define saved_ebp 0(%ebp) /define ret_addr 4(%ebp) #define x 8(%ebp) #if defined (LONG_DOUBLE) #define y 20(%ebp) #else #define y 16(%ebp) #endif LABEL(FUNC): pushl %ebp movl %esp, %ebp PROFILE_FRAME subl $4, %esp pushf /* Save flags */ fstcww cw_user movw cw_user, %ax orw $0x0c00, %ax /* Chop by truncating toward 0 */ movw %ax, cw_chop FLD y fxam fstsww %ax movb %ah, %al andb $0x47, %al cmpb $5, %al je Ly_pinf /* y = +INF */ cmpb $7, %al je Ly_minf /* y = -INF */ andb $0x45, %al cmpb $1, %al je Ly_nan /* y = NAN */ cmpb $0x40, %al je Ly_zero /* y = 0 */ movb %ah, %cl /* cl := fxam(y) */ / -INF < y < 0 or 0 < y < INF. Check x. FLD x fxam fstsww %ax movb %ah, %al andb $0x47, %al cmpb $5, %al je Lx_pinf /* x = +INF */ cmpb $7, %al je Lx_minf /* x = -INF */ andb $0x45, %al cmpb $1, %al je Lx_nan /* x = NAN */ cmpb $0x40, %al je Lx_zero /* x = 0 */ test $2, %ah jnz Lx_neg /* x < 0 */ / x > 0, y non-zero and finite Lx_pos: call Lpower jmp Ldone / x > 0, y non-zero and finite, overflow of y*log(x) ALIGN Lx_pos_ov: FLD x fcomps Lone fnstsww %ax sahf ja Ldone fchs jmp Ldone / x = 0, y non-zero and finite ALIGN Lx_zero:testb $2, %cl jnz Lx_zero_y_neg / x = 0, y > 0 fxch %st(1) /* st(0) := y, st(1) := x */ call Lis_odd_int fstp %st(0) /* Pop y */ cmpl $1, %eax /* Odd integer? */ je Ldone /* Yes: return x (+0 or -0) */ fstp %st(0) /* Pop x */ fldl __const_ZERO /* No: return +0 */ jmp Ldone / x = 0, y < 0 Lx_zero_y_neg: movb %ah, %ch /* ch := fxam(x) */ fxch %st(1) /* st(0) := y, st(1) := x */ call Lis_odd_int fstp %st(0) /* Pop y */ cmpl $1, %eax /* Odd integer? */ je L1 /* Yes: return 1/x (divide by zero) */ fstp %st(0) /* Pop x */ fldl __const_ZERO L1: fdivrs Lone jmp Ldone / x = +INF, y non-zero and finite ALIGN Lx_pinf:fstp %st(1) /* pop x, st(0):=x */ testb $2, %cl /* y > 0? */ jz Ldone /* Yes: return x = +INF */ fstp %st(0) /* Pop x */ fldl __const_ZERO /* Return +0 */ jmp Ldone / x = -INF, y non-zero and finite ALIGN Lx_minf:testb $2, %cl /* y < 0? */ jnz Lx_minf_y_neg /* Yes -> return +0 or -0 */ fxch %st(1) /* st(0) := y, st(1) := x */ call Lis_odd_int fstp %st(0) /* Pop y */ cmpl $1, %eax /* Is y an odd integer? */ je Ldone /* Yes: return x = -INF */ fchs /* No: return -x = +INF */ jmp Ldone / x = +INF, y < 0 and finite ALIGN Lx_minf_y_neg: fstp %st(0) /* Pop x */ call Lis_odd_int fstp %st(0) /* Pop y */ fldl __const_ZERO cmpl $1, %eax /* Is y an odd integer? */ jne Ldone /* No: return +0 */ fchs /* Yes: return -0 */ jmp Ldone / x = NAN, y non-zero and finite ALIGN Lx_nan: fstp %st(1) jmp Ldone / y = 0 ALIGN Ly_zero: fstp %st(0) flds Lone jmp Ldone / x < 0, y non-zero and finite ALIGN Lx_neg: fxch %st(1) /* st(0) := y, st(1) := x */ call Lis_odd_int cmp $0, %eax /* Is y an integer? */ je Lx_neg_invalid /* No -> return NAN */ fxch %st(1) /* st(0) := x, st(1) := y */ fchs /* st(0) := -x */ movl %eax, %ecx /* Remember whether y is odd or not */ call Lpower / Change sign if y is odd cmpl $1, %ecx jne Ldone fchs jmp Ldone / x < 0, y not an integer ALIGN Lx_neg_invalid: fstp %st(0) /* Pop y */ fstp %st(0) /* Pop x */ jmp Linvalid ALIGN Ly_nan: /* pow(x,NAN) = NAN */ Ldone: CONV(x) /* convert to double */ popf leave _xam j_inf Lerange j_nan Ledom Lreturn:EPILOGUE(FUNC) ALIGN Lerange:SET_ERRNO_CONST($ERANGE) jmp Lreturn ALIGN Ledom: SET_ERRNO_CONST($EDOM) jmp Lreturn / pow(x,+INF) ALIGN Ly_pinf:fstp %st(0) FLD x _xam j_nan Ldone fabs fcoms Lone fnstsww %ax fstp %st(0) sahf jz Linvalid ja Lret_pinf Lret_zero: fldl __const_ZERO jmp Ldone / pow(x,-INF) ALIGN Ly_minf:fstp %st(0) /* Pop y */ FLD x _xam j_nan Ldone fabs fcoms Lone fnstsww %ax fstp %st(0) /* Pop x */ sahf jz Linvalid ja Lret_zero Lret_pinf: fldl __HUGE_VAL jmp Ldone ALIGN Linvalid: fldl __const_ZERO fdivl __const_ZERO jmp Ldone / Check if x (in st(0)) is an odd integer / Returns: / 0 x is not an integer / 1 x is an odd integer / 2 x is an even integer (or too big to be odd) / x is not removed from the stack ALIGN Lis_odd_int: fld %st /* st(0)=x, st(1)=x */ fldcww cw_chop frndint /* st(0)=int(x), st(1)=x */ fldcww cw_user fcom %st(1) fnstsww %ax sahf jnz Lnot_int /* x not an integer */ fmuls Lhalf /* st(0)=x/2, st(1)=x */ fldcww cw_chop frndint fldcww cw_user fmuls Ltwo /* st(0)=2*(x/2), st(1)=x */ fcomp %st(1) fnstsww %ax sahf jnz Lodd movl $2, %eax ret ALIGN Lodd: movl $1, %eax ret ALIGN Lnot_int: fstp %st(0) /* pop int(x) */ movl $0, %eax ret / Compute st(0) ^ st(1) / Note: References x! ALIGN Lpower: fyl2x /* st(1) * log(st(0)) */ _xam j_inf Lpower_overflow fld %st fldcww cw_chop frndint fldcww cw_user fsubr %st, %st(1) fxch %st(1) f2xm1 fadds Lone /* TODO: to avoid overflow, scale 1.0, add later */ fscale fstp %st(1) /* Remove the scaling factor */ ret ALIGN Lpower_overflow: fabs FLD x fabs fcomps Lone fnstsww %ax sahf ja Lpower_ret fdivrs Lone Lpower_ret: ret Lhalf: .long 0x3f000000 Lone: .long 0x3f800000 Ltwo: .long 0x40000000