UNPKG

11.9 kBtext/coffeescriptView Raw
1# Symbolic multiplication
2
3# multiplication is commutative, so it can't be used
4# e.g. on two matrices.
5# But it can be used, say, on a scalar and a matrix.,
6# so the output of a multiplication is not
7# always a scalar.
8
9#extern void append(void)
10#static void parse_p1(void)
11#static void parse_p2(void)
12#static void __normalize_radical_factors(int)
13
14Eval_multiply = ->
15 push(cadr(p1))
16 Eval()
17 p1 = cddr(p1)
18 while (iscons(p1))
19 push(car(p1))
20 Eval()
21 multiply()
22 p1 = cdr(p1)
23
24# this one doesn't eval the factors,
25# so you pass i*(-1)^(1/2), it wouldnt't
26# give -1, because i is not evalled
27multiply = ->
28 if (esc_flag)
29 stop("escape key stop")
30 if (isnum(stack[tos - 2]) && isnum(stack[tos - 1]))
31 multiply_numbers()
32 else
33 save()
34 yymultiply()
35 restore()
36
37yymultiply = ->
38 h = 0
39 i = 0
40 n = 0
41
42 # pop operands
43
44 p2 = pop()
45 p1 = pop()
46
47 h = tos
48
49 # is either operand zero?
50
51 if (iszero(p1) || iszero(p2))
52 if evaluatingAsFloats then push_double(0.0) else push(zero)
53 return
54
55 # is either operand a sum?
56
57 #console.log("yymultiply: expanding: " + expanding)
58 if (expanding && isadd(p1))
59 p1 = cdr(p1)
60 if evaluatingAsFloats then push_double(0.0) else push(zero)
61 while (iscons(p1))
62 push(car(p1))
63 push(p2)
64 multiply()
65 add()
66 p1 = cdr(p1)
67 return
68
69 if (expanding && isadd(p2))
70 p2 = cdr(p2)
71 if evaluatingAsFloats then push_double(0.0) else push(zero)
72 while (iscons(p2))
73 push(p1)
74 push(car(p2))
75 multiply()
76 add()
77 p2 = cdr(p2)
78 return
79
80 # scalar times tensor?
81
82 if (!istensor(p1) && istensor(p2))
83 push(p1)
84 push(p2)
85 scalar_times_tensor()
86 return
87
88 # tensor times scalar?
89
90 if (istensor(p1) && !istensor(p2))
91 push(p1)
92 push(p2)
93 tensor_times_scalar()
94 return
95
96 # adjust operands
97
98 if (car(p1) == symbol(MULTIPLY))
99 p1 = cdr(p1)
100 else
101 push(p1)
102 list(1)
103 p1 = pop()
104
105 if (car(p2) == symbol(MULTIPLY))
106 p2 = cdr(p2)
107 else
108 push(p2)
109 list(1)
110 p2 = pop()
111
112 # handle numerical coefficients
113
114 if (isnum(car(p1)) && isnum(car(p2)))
115 push(car(p1))
116 push(car(p2))
117 multiply_numbers()
118 p1 = cdr(p1)
119 p2 = cdr(p2)
120 else if (isnum(car(p1)))
121 push(car(p1))
122 p1 = cdr(p1)
123 else if (isnum(car(p2)))
124 push(car(p2))
125 p2 = cdr(p2)
126 else
127 if evaluatingAsFloats then push_double(1.0) else push(one)
128
129 parse_p1()
130 parse_p2()
131
132 while (iscons(p1) && iscons(p2))
133
134 # if (car(p1)->gamma && car(p2)->gamma) {
135 # combine_gammas(h)
136 # p1 = cdr(p1)
137 # p2 = cdr(p2)
138 # parse_p1()
139 # parse_p2()
140 # continue
141 # }
142
143 if (caar(p1) == symbol(OPERATOR) && caar(p2) == symbol(OPERATOR))
144 push_symbol(OPERATOR)
145 push(cdar(p1))
146 push(cdar(p2))
147 append()
148 cons()
149 p1 = cdr(p1)
150 p2 = cdr(p2)
151 parse_p1()
152 parse_p2()
153 continue
154
155 switch (cmp_expr(p3, p4))
156 when -1
157 push(car(p1))
158 p1 = cdr(p1)
159 parse_p1()
160 when 1
161 push(car(p2))
162 p2 = cdr(p2)
163 parse_p2()
164 when 0
165 combine_factors(h)
166 p1 = cdr(p1)
167 p2 = cdr(p2)
168 parse_p1()
169 parse_p2()
170 else
171 stop("internal error 2")
172
173 # push remaining factors, if any
174
175 while (iscons(p1))
176 push(car(p1))
177 p1 = cdr(p1)
178
179 while (iscons(p2))
180 push(car(p2))
181 p2 = cdr(p2)
182
183 # normalize radical factors
184
185 # example: 2*2(-1/2) -> 2^(1/2)
186
187 # must be done after merge because merge may produce radical
188
189 # example: 2^(1/2-a)*2^a -> 2^(1/2)
190
191 __normalize_radical_factors(h)
192
193 # this hack should not be necessary, unless power returns a multiply
194
195 #for (i = h; i < tos; i++) {
196 # if (car(stack[i]) == symbol(MULTIPLY)) {
197 # multiply_all(tos - h)
198 # return
199 # }
200 #}
201
202 if (expanding)
203 for i in [h...tos]
204 if (isadd(stack[i]))
205 multiply_all(tos - h)
206 return
207
208 # n is the number of result factors on the stack
209
210 n = tos - h
211
212 if (n == 1)
213 return
214
215 # discard integer 1
216
217 if (isrational(stack[h]) && equaln(stack[h], 1))
218 if (n == 2)
219 p7 = pop()
220 pop()
221 push(p7)
222 else
223 stack[h] = symbol(MULTIPLY)
224 list(n)
225 return
226
227 list(n)
228 p7 = pop()
229 push_symbol(MULTIPLY)
230 push(p7)
231 cons()
232
233# Decompose a factor into base and power.
234#
235# input: car(p1) factor
236#
237# output: p3 factor's base
238#
239# p5 factor's power (possibly 1)
240
241parse_p1 = ->
242 p3 = car(p1)
243 p5 = if evaluatingAsFloats then one_as_double else one
244 if (car(p3) == symbol(POWER))
245 p5 = caddr(p3)
246 p3 = cadr(p3)
247
248# Decompose a factor into base and power.
249#
250# input: car(p2) factor
251#
252# output: p4 factor's base
253#
254# p6 factor's power (possibly 1)
255
256parse_p2 = ->
257 p4 = car(p2)
258 p6 = if evaluatingAsFloats then one_as_double else one
259 if (car(p4) == symbol(POWER))
260 p6 = caddr(p4)
261 p4 = cadr(p4)
262
263# h an integer
264combine_factors = (h) ->
265 push(p4)
266 push(p5)
267 push(p6)
268 add()
269 power()
270 p7 = pop()
271 if (isnum(p7))
272 push(stack[h])
273 push(p7)
274 multiply_numbers()
275 stack[h] = pop()
276 else if (car(p7) == symbol(MULTIPLY))
277 # power can return number * factor (i.e. -1 * i)
278 if (isnum(cadr(p7)) && cdddr(p7) == symbol(NIL))
279 push(stack[h])
280 push(cadr(p7))
281 multiply_numbers()
282 stack[h] = pop()
283 push(caddr(p7))
284 else
285 push(p7)
286 else
287 push(p7)
288
289gp = [
290 [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
291 [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
292 [0,0,1,-6,-7,-8,-3,-4,-5,13,14,15,-16,9,10,11,-12],
293 [0,0,6,-1,-11,10,-2,-15,14,12,-5,4,-9,16,-8,7,-13],
294 [0,0,7,11,-1,-9,15,-2,-13,5,12,-3,-10,8,16,-6,-14],
295 [0,0,8,-10,9,-1,-14,13,-2,-4,3,12,-11,-7,6,16,-15],
296 [0,0,3,2,15,-14,1,11,-10,16,-8,7,13,12,-5,4,9],
297 [0,0,4,-15,2,13,-11,1,9,8,16,-6,14,5,12,-3,10],
298 [0,0,5,14,-13,2,10,-9,1,-7,6,16,15,-4,3,12,11],
299 [0,0,13,12,-5,4,16,-8,7,-1,-11,10,-3,-2,-15,14,-6],
300 [0,0,14,5,12,-3,8,16,-6,11,-1,-9,-4,15,-2,-13,-7],
301 [0,0,15,-4,3,12,-7,6,16,-10,9,-1,-5,-14,13,-2,-8],
302 [0,0,16,-9,-10,-11,-13,-14,-15,-3,-4,-5,1,-6,-7,-8,2],
303 [0,0,9,-16,8,-7,-12,5,-4,-2,-15,14,6,-1,-11,10,3],
304 [0,0,10,-8,-16,6,-5,-12,3,15,-2,-13,7,11,-1,-9,4],
305 [0,0,11,7,-6,-16,4,-3,-12,-14,13,-2,8,-10,9,-1,5],
306 [0,0,12,13,14,15,9,10,11,-6,-7,-8,-2,-3,-4,-5,-1]
307]
308
309#if 0
310
311# h an int
312combine_gammas = (h) ->
313 n = gp[Math.floor(p1.gamma)][Math.floor(p2.gamma)]
314 if (n < 0)
315 n = -n
316 push(stack[h])
317 negate()
318 stack[h] = pop()
319 if (n > 1)
320 push(_gamma[n])
321
322#endif
323
324multiply_noexpand = ->
325 prev_expanding = expanding
326 expanding = 0
327 multiply()
328 expanding = prev_expanding
329
330# multiply n factors on stack
331
332# n an integer
333multiply_all = (n) ->
334 i = 0
335 if (n == 1)
336 return
337 if (n == 0)
338 push if evaluatingAsFloats then one_as_double else one
339 return
340 h = tos - n
341 push(stack[h])
342 for i in [1...n]
343 push(stack[h + i])
344 multiply()
345 stack[h] = pop()
346 moveTos h + 1
347
348# n an integer
349multiply_all_noexpand = (n) ->
350 prev_expanding = expanding
351 expanding = 0
352 multiply_all(n)
353 expanding = prev_expanding
354
355#-----------------------------------------------------------------------------
356#
357# Symbolic division, or numeric division if doubles are found.
358#
359# Input: Dividend and divisor on stack
360#
361# Output: Quotient on stack
362#
363#-----------------------------------------------------------------------------
364
365divide = ->
366 if (isnum(stack[tos - 2]) && isnum(stack[tos - 1]))
367 divide_numbers()
368 else
369 inverse()
370 multiply()
371
372# this is different from inverse of a matrix (inv)!
373inverse = ->
374 if (isnum(stack[tos - 1]))
375 invert_number()
376 else
377 push_integer(-1)
378 power()
379
380reciprocate = ->
381 inverse()
382
383negate = ->
384 if (isnum(stack[tos - 1]))
385 negate_number()
386 else
387 if evaluatingAsFloats then push_double(-1.0) else push_integer(-1)
388 multiply()
389
390negate_expand = ->
391 prev_expanding = expanding
392 expanding = 1
393 negate()
394 expanding = prev_expanding
395
396negate_noexpand = ->
397 prev_expanding = expanding
398 expanding = 0
399 negate()
400 expanding = prev_expanding
401
402#-----------------------------------------------------------------------------
403#
404# Normalize radical factors
405#
406# Input: stack[h] Coefficient factor, possibly 1
407#
408# stack[h + 1] Second factor
409#
410# stack[tos - 1] Last factor
411#
412# Output: Reduced coefficent and normalized radicals (maybe)
413#
414# Example: 2*2^(-1/2) -> 2^(1/2)
415#
416# (power number number) is guaranteed to have the following properties:
417#
418# 1. Base is an integer
419#
420# 2. Absolute value of exponent < 1
421#
422# These properties are assured by the power function.
423#
424#-----------------------------------------------------------------------------
425
426#define A p1
427#define B p2
428
429#define BASE p3
430#define EXPO p4
431
432#define TMP p5
433
434
435# h is an int
436__normalize_radical_factors = (h) ->
437
438 i = 0
439 # if coeff is 1 or floating then don't bother
440
441 if (isplusone(stack[h]) || isminusone(stack[h]) || isdouble(stack[h]))
442 return
443
444 # if no radicals then don't bother
445
446 for i in [(h + 1)...tos]
447 if (__is_radical_number(stack[i]))
448 break
449
450 if (i == tos)
451 return
452
453 # ok, try to simplify
454
455 save()
456
457 # numerator
458
459 push(stack[h])
460 mp_numerator()
461 #console.log("__normalize_radical_factors numerator: " + stack[tos-1])
462 p1 = pop(); # p1 is A
463
464 for i in [(h + 1)...tos]
465
466 if (isplusone(p1) || isminusone(p1)) # p1 is A
467 break
468
469 if (!__is_radical_number(stack[i]))
470 continue
471
472 p3 = cadr(stack[i]); #p3 is BASE
473 p4 = caddr(stack[i]); #p4 is EXPO
474
475 # exponent must be negative
476
477 if (!isnegativenumber(p4)) #p4 is EXPO
478 continue
479
480 # numerator divisible by p3 (base)?
481
482 push(p1); # p1 is A
483 push(p3); #p3 is BASE
484 divide()
485
486 p5 = pop(); #p5 is TMP
487
488 if (!isinteger(p5)) #p5 is TMP
489 continue
490
491 # reduce numerator
492
493 p1 = p5; # p1 is A #p5 is TMP
494
495 # invert radical
496
497 push_symbol(POWER)
498 push(p3); #p3 is BASE
499 push if evaluatingAsFloats then one_as_double else one
500 push(p4); #p4 is EXPO
501 add()
502 list(3)
503 stack[i] = pop()
504
505 # denominator
506
507 push(stack[h])
508 mp_denominator()
509 #console.log("__normalize_radical_factors denominator: " + stack[tos-1])
510 p2 = pop(); # p2 is B
511
512 for i in [(h + 1)...tos]
513
514 if (isplusone(p2)) # p2 is B
515 break
516
517 if (!__is_radical_number(stack[i]))
518 continue
519
520 p3 = cadr(stack[i]); #p3 is BASE
521 p4 = caddr(stack[i]); #p4 is EXPO
522
523 # exponent must be positive
524
525 if (isnegativenumber(p4)) #p4 is EXPO
526 continue
527
528 # denominator divisible by p3? #p3 is BASE
529
530 push(p2); # p2 is B
531 push(p3); #p3 is BASE
532 divide()
533
534 p5 = pop(); #p5 is TMP
535
536 if (!isinteger(p5)) #p5 is TMP
537 continue
538 #console.log("__new radical p5: " + p5.toString())
539 #console.log("__new radical top stack: " + stack[tos-1])
540
541 # reduce denominator
542
543 p2 = p5; # p2 is B #p5 is TMP
544
545 # invert radical
546
547 push_symbol(POWER)
548 push(p3); #p3 is BASE
549 push(p4); #p4 is EXPO
550 #console.log("__new radical p3: " + p3.toString())
551 #console.log("__new radical p4: " + p4.toString())
552 push(one)
553 subtract()
554
555 if dontCreateNewRadicalsInDenominatorWhenEvalingMultiplication
556 if (isinteger(p3) and !isinteger(stack[tos-1]) and isnegativenumber(stack[tos - 1]))
557 # bail out,
558 # we want to avoid going ahead with the subtraction of
559 # the exponents, because that would turn a perfectly good
560 # integer exponent in the denominator into a fractional one
561 # i.e. a radical.
562 # Note that this only prevents new radicals ending up
563 # in the denominator, it doesn't fix existing ones.
564 pop()
565 pop()
566 pop()
567
568 push(p1); # p1 is A
569 push(p3); #p3 is BASE
570 divide()
571 p1 = pop()
572
573 break
574 #console.log("__new radical exponent: " + stack[tos-1])
575
576 list(3)
577 stack[i] = pop()
578
579 # reconstitute the coefficient
580
581 push(p1); # p1 is A
582 push(p2); # p2 is B
583 divide()
584 stack[h] = pop()
585
586 restore()
587
588# don't include i
589
590# p is a U
591__is_radical_number = (p) ->
592 # don't use i
593
594 if (car(p) == symbol(POWER) && isnum(cadr(p)) && isnum(caddr(p)) && !isminusone(cadr(p)))
595 return 1
596 else
597 return 0
598
599#-----------------------------------------------------------------------------
600#
601# > a*hilbert(2)
602# ((a,1/2*a),(1/2*a,1/3*a))
603#
604# Note that "a" is presumed to be a scalar. Is this correct?
605#
606# Yes, because "*" has no meaning if "a" is a tensor.
607# To multiply tensors, "dot" or "outer" should be used.
608#
609# > dot(a,hilbert(2))
610# dot(a,((1,1/2),(1/2,1/3)))
611#
612# In this case "a" could be a scalar or tensor so the result is not
613# expanded.
614#
615#-----------------------------------------------------------------------------
616
617