Wrong Sign: even powers of negative numbers

[Charles Pegge]

I was working on a demo to solve the roots of functions (example: x^3-3*x^2+2*x=0)
and found that when negative numbers were powered by 2 4 6... they were coming out negative, when of course, they should be positive.

Even with the FPU, powers are not straight forward. A simple a^b requires a bucket full of code.

This is my candidate solution to the problem. It checks whether the exponent is odd or even, and whether the number was positive or negative and acts accordingly.

(There is still a problem with fractional powers of negative numbers - but only complex numbers can solve that one.)

If you can avoid using powers, and use multiplication wherever possible, your code will run a lot faster.

Code: OxygenBasic

1.  
2. deff pow
3. '=======
4.  
5. sub esp,16
6. fstcw [esp]
7. fstcw [esp+2]
8. or [esp],0xc00
9. fldcw [esp]
10. '
11. 'CHECK FOR ODD OR EVEN EXPONENT
12. '
13. (
14.   fist dword [esp+8]
15.   and byte [esp+8],1
16.   jnz exit
17.   fxch
18.   fabs
19.   fyl2x
20.   fld st0
21.   frndint
22.   fsub st1,st0
23.   fxch
24.   f2xm1
25.   fld1
26.   faddp st1,st0
27.   fscale
28.   fstp st1
29.   jmp fwd nng
30. )
31. fxch
32. fst dword [esp+8]
33. fabs
34. fyl2x
35. fld st0
36. frndint
37. fsub st1,st0
38. fxch
39. f2xm1
40. fld1
41. faddp st1,st0
42. fscale
43. fstp st1
44. (
45.   and byte [esp+11],128
46.   jz exit
47.   fchs
48. )
49. '
50. .nng 'DONE
51. '
52. fldcw [esp+2]
53. add esp,16
54.  
55. end deff
56.  
57.  
58. print pow(-2,3)
59.  
60.  

Charles

Rare and Valuable reference for FPU programming:

http://www.website.masmforum.com/tutorials/fptute/fpuchap11.htm

[Charles Pegge]

Could be useful.

I think those FPU instructions are fast though, it's mostly multiplication and you get the full 80bit precision. But x*x*x will always be faster then x^3.

Here is a version supporting negative integer powers:

Code: OxygenBasic

1. function IntPower(double x,sys n) as double
2.   sys i,inv
3.   double a=1
4.   if n<0 then n=-n : inv=1
5.   for i=1 to n
6.     a*=x
7.   next
8.   if inv then return 1/a
9.   return a
10. end function
11.  
12. '=====
13. 'TESTS
14. '=====
15.  
16. print IntPower 2,0
17. print IntPower 2,1
18. print IntPower 2,-1
19. print IntPower -2,2
20. print IntPower -2,3
21.  

Fractional powers of negative numbers are quite interesting. Are you familiar with complex numbers Peter?

Charles

[Charles Pegge]

sys a[ 100]
sys i=1

IndexBase 0

#show a[ i]=10

sys a[100]

IndexBase 2

#show a[ i]=10


You can see how arrays are accessed. When the index base is more than 0 then the index is pre-decremented by the index-base value.