1 | # Rational power function
|
2 |
|
3 |
|
4 |
|
5 | qpow = ->
|
6 | save()
|
7 | qpowf()
|
8 | restore()
|
9 |
|
10 | #define BASE p1
|
11 | #define EXPO p2
|
12 |
|
13 | qpowf = ->
|
14 | expo = 0
|
15 | #unsigned int a, b, *t, *x, *y
|
16 |
|
17 | p2 = pop(); # p2 is EXPO
|
18 | p1 = pop(); # p1 is BASE
|
19 |
|
20 | # if base is 1 or exponent is 0 then return 1
|
21 | if (isplusone(p1) || iszero(p2)) # p1 is BASE # p2 is EXPO
|
22 | push_integer(1)
|
23 | return
|
24 |
|
25 | # if (-1)^(1/2) -> leave it as is
|
26 | if (isminusone(p1) and isoneovertwo(p2)) # p1 is BASE # p2 is EXPO
|
27 | push(imaginaryunit)
|
28 | return
|
29 |
|
30 | # if base is zero then return 0
|
31 | if (iszero(p1)) # p1 is BASE
|
32 | if (isnegativenumber(p2)) # p2 is EXPO
|
33 | stop("divide by zero")
|
34 | push(zero)
|
35 | return
|
36 |
|
37 | # if exponent is 1 then return base
|
38 | if (isplusone(p2)) # p2 is EXPO
|
39 | push(p1); # p1 is BASE
|
40 | return
|
41 |
|
42 | # if exponent is integer then power
|
43 | if (isinteger(p2)) # p2 is EXPO
|
44 | push(p2); # p2 is EXPO
|
45 | expo = pop_integer()
|
46 | if isNaN(expo)
|
47 | # expo greater than 32 bits
|
48 | push_symbol(POWER)
|
49 | push(p1); # p1 is BASE
|
50 | push(p2); # p2 is EXPO
|
51 | list(3)
|
52 | return
|
53 |
|
54 | x = mpow(p1.q.a, Math.abs(expo)); # p1 is BASE
|
55 | y = mpow(p1.q.b, Math.abs(expo)); # p1 is BASE
|
56 | if (expo < 0)
|
57 | t = x
|
58 | x = y
|
59 | y = t
|
60 | x = makeSignSameAs(x,y)
|
61 | y = makePositive(y)
|
62 |
|
63 | p3 = new U()
|
64 | p3.k = NUM
|
65 | p3.q.a = x
|
66 | p3.q.b = y
|
67 | push(p3)
|
68 | return
|
69 |
|
70 | # from here on out the exponent is NOT an integer
|
71 |
|
72 | # if base is -1 then normalize polar angle
|
73 | if (isminusone(p1)) # p1 is BASE
|
74 | push(p2); # p2 is EXPO
|
75 | normalize_angle()
|
76 | return
|
77 |
|
78 | # if base is negative then (-N)^M -> N^M * (-1)^M
|
79 | if (isnegativenumber(p1)) # p1 is BASE
|
80 | push(p1); # p1 is BASE
|
81 | negate()
|
82 | push(p2); # p2 is EXPO
|
83 | qpow()
|
84 |
|
85 | push_integer(-1)
|
86 | push(p2); # p2 is EXPO
|
87 | qpow()
|
88 |
|
89 | multiply()
|
90 | return
|
91 |
|
92 | # if p1 (BASE) is not an integer then power numerator and denominator
|
93 | if (!isinteger(p1)) # p1 is BASE
|
94 | push(p1); # p1 is BASE
|
95 | mp_numerator()
|
96 | push(p2); # p2 is EXPO
|
97 | qpow()
|
98 | push(p1); # p1 is BASE
|
99 | mp_denominator()
|
100 | push(p2); # p2 is EXPO
|
101 | negate()
|
102 | qpow()
|
103 | multiply()
|
104 | return
|
105 |
|
106 | # At this point p1 (BASE) is a positive integer.
|
107 |
|
108 | # If p1 (BASE) is small then factor it.
|
109 | if (is_small_integer(p1)) # p1 is BASE
|
110 | push(p1); # p1 is BASE
|
111 | push(p2); # p2 is EXPO
|
112 | quickfactor()
|
113 | return
|
114 |
|
115 | # At this point p1 (BASE) is a positive integer and p2 (EXPO) is not an integer.
|
116 |
|
117 | if ( !p2.q.a.isSmall || !p2.q.b.isSmall ) # p2 is EXPO
|
118 | push_symbol(POWER)
|
119 | push(p1) # p1 is BASE
|
120 | push(p2); # p2 is EXPO
|
121 | list(3)
|
122 | return
|
123 |
|
124 | a = p2.q.a; # p2 is EXPO
|
125 | b = p2.q.b; # p2 is EXPO
|
126 |
|
127 | x = mroot(p1.q.a, b); # p1 is BASE
|
128 |
|
129 | if (x == 0)
|
130 | push_symbol(POWER)
|
131 | push(p1); # p1 is BASE
|
132 | push(p2); # p2 is EXPO
|
133 | list(3)
|
134 | return
|
135 |
|
136 | y = mpow(x, a)
|
137 |
|
138 | #mfree(x)
|
139 |
|
140 | p3 = new U()
|
141 |
|
142 | p3.k = NUM
|
143 |
|
144 | if (p2.q.a.isNegative()) # p2 is EXPO
|
145 | p3.q.a = bigInt(1)
|
146 | p3.q.b = y
|
147 | else
|
148 | p3.q.a = y
|
149 | p3.q.b = bigInt(1)
|
150 |
|
151 | push(p3)
|
152 |
|
153 | #-----------------------------------------------------------------------------
|
154 | #
|
155 | # Normalize the angle of unit imaginary, i.e. (-1) ^ N
|
156 | #
|
157 | # Input: N on stack (must be rational, not float)
|
158 | #
|
159 | # Output: Result on stack
|
160 | #
|
161 | # Note:
|
162 | #
|
163 | # n = q * d + r
|
164 | #
|
165 | # Example:
|
166 | # n d q r
|
167 | #
|
168 | # (-1)^(8/3) -> (-1)^(2/3) 8 3 2 2
|
169 | # (-1)^(7/3) -> (-1)^(1/3) 7 3 2 1
|
170 | # (-1)^(5/3) -> -(-1)^(2/3) 5 3 1 2
|
171 | # (-1)^(4/3) -> -(-1)^(1/3) 4 3 1 1
|
172 | # (-1)^(2/3) -> (-1)^(2/3) 2 3 0 2
|
173 | # (-1)^(1/3) -> (-1)^(1/3) 1 3 0 1
|
174 | #
|
175 | # (-1)^(-1/3) -> -(-1)^(2/3) -1 3 -1 2
|
176 | # (-1)^(-2/3) -> -(-1)^(1/3) -2 3 -1 1
|
177 | # (-1)^(-4/3) -> (-1)^(2/3) -4 3 -2 2
|
178 | # (-1)^(-5/3) -> (-1)^(1/3) -5 3 -2 1
|
179 | # (-1)^(-7/3) -> -(-1)^(2/3) -7 3 -3 2
|
180 | # (-1)^(-8/3) -> -(-1)^(1/3) -8 3 -3 1
|
181 | #
|
182 | #-----------------------------------------------------------------------------
|
183 |
|
184 | #define A p1
|
185 | #define Q p2
|
186 | #define R p3
|
187 |
|
188 | normalize_angle = ->
|
189 | save()
|
190 |
|
191 | p1 = pop(); # p1 is A
|
192 |
|
193 | # integer exponent?
|
194 |
|
195 | if (isinteger(p1)) # p1 is A
|
196 | if (p1.q.a.isOdd()) # p1 is A
|
197 | push_integer(-1); # odd exponent
|
198 | else
|
199 | push_integer(1); # even exponent
|
200 | restore()
|
201 | return
|
202 |
|
203 | # floor
|
204 |
|
205 | push(p1); # p1 is A
|
206 | bignum_truncate()
|
207 | p2 = pop(); # p2 is Q
|
208 |
|
209 | if (isnegativenumber(p1)) # p1 is A
|
210 | push(p2) # p2 is Q
|
211 | push_integer(-1)
|
212 | add()
|
213 | p2 = pop(); # p2 is Q
|
214 |
|
215 | # remainder (always positive)
|
216 |
|
217 | push(p1); # p1 is A
|
218 | push(p2); # p2 is Q
|
219 | subtract()
|
220 | p3 = pop(); # p3 is R
|
221 |
|
222 | # remainder becomes new angle
|
223 | push_symbol(POWER)
|
224 | push_integer(-1)
|
225 | push(p3) # p3 is R
|
226 | list(3)
|
227 |
|
228 | # negate if quotient is odd
|
229 |
|
230 | if (p2.q.a.isOdd()) # p2 is Q
|
231 | negate()
|
232 |
|
233 | restore()
|
234 |
|
235 | is_small_integer = (p) ->
|
236 | return p.q.a.isSmall
|