| 1 | # power function for double precision floating point
|
| 2 |
|
| 3 |
|
| 4 |
|
| 5 | dpow = ->
|
| 6 | a = 0.0
|
| 7 | b = 0.0
|
| 8 | base = 0.0
|
| 9 | expo = 0.0
|
| 10 | result = 0.0
|
| 11 | theta = 0.0
|
| 12 |
|
| 13 | expo = pop_double()
|
| 14 | base = pop_double()
|
| 15 |
|
| 16 | # divide by zero?
|
| 17 |
|
| 18 | if (base == 0.0 && expo < 0.0)
|
| 19 | stop("divide by zero")
|
| 20 |
|
| 21 | # nonnegative base or integer power?
|
| 22 |
|
| 23 | if (base >= 0.0 || (expo % 1.0) == 0.0)
|
| 24 | result = Math.pow(base, expo)
|
| 25 | push_double(result)
|
| 26 | return
|
| 27 |
|
| 28 | result = Math.pow(Math.abs(base), expo)
|
| 29 |
|
| 30 | theta = Math.PI * expo
|
| 31 |
|
| 32 | # this ensures the real part is 0.0 instead of a tiny fraction
|
| 33 |
|
| 34 | if ((expo % 0.5) == 0.0)
|
| 35 | a = 0.0
|
| 36 | b = Math.sin(theta)
|
| 37 | else
|
| 38 | a = Math.cos(theta)
|
| 39 | b = Math.sin(theta)
|
| 40 |
|
| 41 | push_double(a * result)
|
| 42 | push_double(b * result)
|
| 43 | push(imaginaryunit)
|
| 44 | multiply()
|
| 45 | add()
|