UNPKG

4.45 kBtext/coffeescriptView Raw
1# Rational power function
2
3
4
5qpow = ->
6 save()
7 qpowf()
8 restore()
9
10#define BASE p1
11#define EXPO p2
12
13qpowf = ->
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
188normalize_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
235is_small_integer = (p) ->
236 return p.q.a.isSmall