UNPKG

9.62 kBtext/coffeescriptView Raw
1# derivative
2
3
4
5#define F p3
6#define X p4
7#define N p5
8
9Eval_derivative = ->
10
11 # evaluate 1st arg to get function F
12
13 i = 0
14 p1 = cdr(p1)
15 push(car(p1))
16 Eval()
17
18 # evaluate 2nd arg and then...
19
20 # example result of 2nd arg what to do
21 #
22 # d(f) nil guess X, N = nil
23 # d(f,2) 2 guess X, N = 2
24 # d(f,x) x X = x, N = nil
25 # d(f,x,2) x X = x, N = 2
26 # d(f,x,y) x X = x, N = y
27
28 p1 = cdr(p1)
29 push(car(p1))
30 Eval()
31
32 p2 = pop()
33 if (p2 == symbol(NIL))
34 guess()
35 push(symbol(NIL))
36 else if (isnum(p2))
37 guess()
38 push(p2)
39 else
40 push(p2)
41 p1 = cdr(p1)
42 push(car(p1))
43 Eval()
44
45 p5 = pop(); # p5 is N
46 p4 = pop(); # p4 is X
47 p3 = pop(); # p3 is F
48
49 while (1)
50
51 # p5 (N) might be a symbol instead of a number
52
53 if (isnum(p5)) # p5 is N
54 push(p5); # p5 is N
55 n = pop_integer()
56 if (isNaN(n))
57 stop("nth derivative: check n")
58 else
59 n = 1
60
61 push(p3); # p3 is F
62
63 if (n >= 0)
64 for i in [0...n]
65 push(p4); # p4 is X
66 derivative()
67 else
68 n = -n
69 for i in [0...n]
70 push(p4); # p4 is X
71 integral()
72
73 p3 = pop(); # p3 is F
74
75 # if p5 (N) is nil then arglist is exhausted
76
77 if (p5 == symbol(NIL)) # p5 is N
78 break
79
80 # otherwise...
81
82 # N arg1 what to do
83 #
84 # number nil break
85 # number number N = arg1, continue
86 # number symbol X = arg1, N = arg2, continue
87 #
88 # symbol nil X = N, N = nil, continue
89 # symbol number X = N, N = arg1, continue
90 # symbol symbol X = N, N = arg1, continue
91
92 if (isnum(p5)) # p5 is N
93 p1 = cdr(p1)
94 push(car(p1))
95 Eval()
96 p5 = pop(); # p5 is N
97 if (p5 == symbol(NIL)) # p5 is N
98 break; # arglist exhausted
99 if (isnum(p5)) # p5 is N
100 doNothing = 1 # N = arg1
101 else
102 p4 = p5; # p5 is N # p4 is X # X = arg1
103 p1 = cdr(p1)
104 push(car(p1))
105 Eval()
106 p5 = pop(); # p5 is N # N = arg2
107 else
108 p4 = p5; # p5 is N # p4 is X # X = N
109 p1 = cdr(p1)
110 push(car(p1))
111 Eval()
112 p5 = pop(); # p5 is N # N = arg1
113
114 push(p3); # p3 is F # final result
115
116derivative = ->
117 save()
118 p2 = pop()
119 p1 = pop()
120 if (isnum(p2))
121 stop("undefined function")
122 if (istensor(p1))
123 if (istensor(p2))
124 d_tensor_tensor()
125 else
126 d_tensor_scalar()
127 else
128 if (istensor(p2))
129 d_scalar_tensor()
130 else
131 d_scalar_scalar()
132 restore()
133
134d_scalar_scalar = ->
135 if (issymbol(p2))
136 d_scalar_scalar_1()
137 else
138 # Example: d(sin(cos(x)),cos(x))
139 # Replace cos(x) <- X, find derivative, then do X <- cos(x)
140 push(p1); # sin(cos(x))
141 push(p2); # cos(x)
142 push(symbol(SECRETX)); # X
143 subst(); # sin(cos(x)) -> sin(X)
144 push(symbol(SECRETX)); # X
145 derivative()
146 push(symbol(SECRETX)); # X
147 push(p2); # cos(x)
148 subst(); # cos(X) -> cos(cos(x))
149
150d_scalar_scalar_1 = ->
151 # d(x,x)?
152
153 if (equal(p1, p2))
154 push(one)
155 return
156
157 # d(a,x)?
158
159 if (!iscons(p1))
160 push(zero)
161 return
162
163 if (isadd(p1))
164 dsum()
165 return
166
167 if (car(p1) == symbol(MULTIPLY))
168 dproduct()
169 return
170
171 if (car(p1) == symbol(POWER))
172 dpower()
173 return
174
175 if (car(p1) == symbol(DERIVATIVE))
176 dd()
177 return
178
179 if (car(p1) == symbol(LOG))
180 dlog()
181 return
182
183 if (car(p1) == symbol(SIN))
184 dsin()
185 return
186
187 if (car(p1) == symbol(COS))
188 dcos()
189 return
190
191 if (car(p1) == symbol(TAN))
192 dtan()
193 return
194
195 if (car(p1) == symbol(ARCSIN))
196 darcsin()
197 return
198
199 if (car(p1) == symbol(ARCCOS))
200 darccos()
201 return
202
203 if (car(p1) == symbol(ARCTAN))
204 darctan()
205 return
206
207 if (car(p1) == symbol(SINH))
208 dsinh()
209 return
210
211 if (car(p1) == symbol(COSH))
212 dcosh()
213 return
214
215 if (car(p1) == symbol(TANH))
216 dtanh()
217 return
218
219 if (car(p1) == symbol(ARCSINH))
220 darcsinh()
221 return
222
223 if (car(p1) == symbol(ARCCOSH))
224 darccosh()
225 return
226
227 if (car(p1) == symbol(ARCTANH))
228 darctanh()
229 return
230
231 if (car(p1) == symbol(ABS))
232 dabs()
233 return
234
235 if (car(p1) == symbol(SGN))
236 dsgn()
237 return
238
239 if (car(p1) == symbol(HERMITE))
240 dhermite()
241 return
242
243 if (car(p1) == symbol(ERF))
244 derf()
245 return
246
247 if (car(p1) == symbol(ERFC))
248 derfc()
249 return
250
251 if (car(p1) == symbol(BESSELJ))
252 if (iszero(caddr(p1)))
253 dbesselj0()
254 else
255 dbesseljn()
256 return
257
258 if (car(p1) == symbol(BESSELY))
259 if (iszero(caddr(p1)))
260 dbessely0()
261 else
262 dbesselyn()
263 return
264
265 if (car(p1) == symbol(INTEGRAL) && caddr(p1) == p2)
266 derivative_of_integral()
267 return
268
269 dfunction()
270
271dsum = ->
272 h = tos
273 p1 = cdr(p1)
274 while (iscons(p1))
275 push(car(p1))
276 push(p2)
277 derivative()
278 p1 = cdr(p1)
279 add_all(tos - h)
280
281dproduct = ->
282 i = 0
283 j = 0
284 n = 0
285 n = length(p1) - 1
286 for i in [0...n]
287 p3 = cdr(p1)
288 for j in [0...n]
289 push(car(p3))
290 if (i == j)
291 push(p2)
292 derivative()
293 p3 = cdr(p3)
294 multiply_all(n)
295 add_all(n)
296
297#-----------------------------------------------------------------------------
298#
299# v
300# y = u
301#
302# log y = v log u
303#
304# 1 dy v du dv
305# - -- = - -- + (log u) --
306# y dx u dx dx
307#
308# dy v v du dv
309# -- = u (- -- + (log u) --)
310# dx u dx dx
311#
312#-----------------------------------------------------------------------------
313
314dpower = ->
315 push(caddr(p1)); # v/u
316 push(cadr(p1))
317 divide()
318
319 push(cadr(p1)); # du/dx
320 push(p2)
321 derivative()
322
323 multiply()
324
325 push(cadr(p1)); # log u
326 logarithm()
327
328 push(caddr(p1)); # dv/dx
329 push(p2)
330 derivative()
331
332 multiply()
333
334 add()
335
336 push(p1); # u^v
337
338 multiply()
339
340dlog = ->
341 push(cadr(p1))
342 push(p2)
343 derivative()
344 push(cadr(p1))
345 divide()
346
347# derivative of derivative
348#
349# example: d(d(f(x,y),y),x)
350#
351# p1 = d(f(x,y),y)
352#
353# p2 = x
354#
355# cadr(p1) = f(x,y)
356#
357# caddr(p1) = y
358
359dd = ->
360 # d(f(x,y),x)
361
362 push(cadr(p1))
363 push(p2)
364 derivative()
365
366 p3 = pop()
367
368 if (car(p3) == symbol(DERIVATIVE))
369
370 # sort dx terms
371
372 push_symbol(DERIVATIVE)
373 push_symbol(DERIVATIVE)
374 push(cadr(p3))
375
376 if (lessp(caddr(p3), caddr(p1)))
377 push(caddr(p3))
378 list(3)
379 push(caddr(p1))
380 else
381 push(caddr(p1))
382 list(3)
383 push(caddr(p3))
384
385 list(3)
386
387 else
388 push(p3)
389 push(caddr(p1))
390 derivative()
391
392# derivative of a generic function
393
394dfunction = ->
395 p3 = cdr(p1); # p3 is the argument list for the function
396
397 if (p3 == symbol(NIL) || Find(p3, p2))
398 push_symbol(DERIVATIVE)
399 push(p1)
400 push(p2)
401 list(3)
402 else
403 push(zero)
404
405dsin = ->
406 push(cadr(p1))
407 push(p2)
408 derivative()
409 push(cadr(p1))
410 cosine()
411 multiply()
412
413dcos = ->
414 push(cadr(p1))
415 push(p2)
416 derivative()
417 push(cadr(p1))
418 sine()
419 multiply()
420 negate()
421
422dtan = ->
423 push(cadr(p1))
424 push(p2)
425 derivative()
426 push(cadr(p1))
427 cosine()
428 push_integer(-2)
429 power()
430 multiply()
431
432darcsin = ->
433 push(cadr(p1))
434 push(p2)
435 derivative()
436 push_integer(1)
437 push(cadr(p1))
438 push_integer(2)
439 power()
440 subtract()
441 push_rational(-1, 2)
442 power()
443 multiply()
444
445darccos = ->
446 push(cadr(p1))
447 push(p2)
448 derivative()
449 push_integer(1)
450 push(cadr(p1))
451 push_integer(2)
452 power()
453 subtract()
454 push_rational(-1, 2)
455 power()
456 multiply()
457 negate()
458
459# Without simplify With simplify
460#
461# d(arctan(y/x),x) -y/(x^2*(y^2/x^2+1)) -y/(x^2+y^2)
462#
463# d(arctan(y/x),y) 1/(x*(y^2/x^2+1)) x/(x^2+y^2)
464
465darctan = ->
466 push(cadr(p1))
467 push(p2)
468 derivative()
469 push_integer(1)
470 push(cadr(p1))
471 push_integer(2)
472 power()
473 add()
474 inverse()
475 multiply()
476 simplify()
477
478dsinh = ->
479 push(cadr(p1))
480 push(p2)
481 derivative()
482 push(cadr(p1))
483 ycosh()
484 multiply()
485
486dcosh = ->
487 push(cadr(p1))
488 push(p2)
489 derivative()
490 push(cadr(p1))
491 ysinh()
492 multiply()
493
494dtanh = ->
495 push(cadr(p1))
496 push(p2)
497 derivative()
498 push(cadr(p1))
499 ycosh()
500 push_integer(-2)
501 power()
502 multiply()
503
504darcsinh = ->
505 push(cadr(p1))
506 push(p2)
507 derivative()
508 push(cadr(p1))
509 push_integer(2)
510 power()
511 push_integer(1)
512 add()
513 push_rational(-1, 2)
514 power()
515 multiply()
516
517darccosh = ->
518 push(cadr(p1))
519 push(p2)
520 derivative()
521 push(cadr(p1))
522 push_integer(2)
523 power()
524 push_integer(-1)
525 add()
526 push_rational(-1, 2)
527 power()
528 multiply()
529
530darctanh = ->
531 push(cadr(p1))
532 push(p2)
533 derivative()
534 push_integer(1)
535 push(cadr(p1))
536 push_integer(2)
537 power()
538 subtract()
539 inverse()
540 multiply()
541
542dabs = ->
543 push(cadr(p1))
544 push(p2)
545 derivative()
546 push(cadr(p1))
547 sgn()
548 multiply()
549
550dsgn = ->
551 push(cadr(p1))
552 push(p2)
553 derivative()
554 push(cadr(p1))
555 dirac()
556 multiply()
557 push_integer(2)
558 multiply()
559
560dhermite = ->
561 push(cadr(p1))
562 push(p2)
563 derivative()
564 push_integer(2)
565 push(caddr(p1))
566 multiply()
567 multiply()
568 push(cadr(p1))
569 push(caddr(p1))
570 push_integer(-1)
571 add()
572 hermite()
573 multiply()
574
575derf = ->
576 push(cadr(p1))
577 push_integer(2)
578 power()
579 push_integer(-1)
580 multiply()
581 exponential()
582 if evaluatingAsFloats
583 push_double(Math.PI)
584 else
585 push_symbol(PI)
586 push_rational(-1,2)
587 power()
588 multiply()
589 push_integer(2)
590 multiply()
591 push(cadr(p1))
592 push(p2)
593 derivative()
594 multiply()
595
596derfc = ->
597 push(cadr(p1))
598 push_integer(2)
599 power()
600 push_integer(-1)
601 multiply()
602 exponential()
603 if evaluatingAsFloats
604 push_double(Math.PI)
605 else
606 push_symbol(PI)
607 push_rational(-1,2)
608 power()
609 multiply()
610 push_integer(-2)
611 multiply()
612 push(cadr(p1))
613 push(p2)
614 derivative()
615 multiply()
616
617dbesselj0 = ->
618 push(cadr(p1))
619 push(p2)
620 derivative()
621 push(cadr(p1))
622 push_integer(1)
623 besselj()
624 multiply()
625 push_integer(-1)
626 multiply()
627
628dbesseljn = ->
629 push(cadr(p1))
630 push(p2)
631 derivative()
632 push(cadr(p1))
633 push(caddr(p1))
634 push_integer(-1)
635 add()
636 besselj()
637 push(caddr(p1))
638 push_integer(-1)
639 multiply()
640 push(cadr(p1))
641 divide()
642 push(cadr(p1))
643 push(caddr(p1))
644 besselj()
645 multiply()
646 add()
647 multiply()
648
649dbessely0 = ->
650 push(cadr(p1))
651 push(p2)
652 derivative()
653 push(cadr(p1))
654 push_integer(1)
655 besselj()
656 multiply()
657 push_integer(-1)
658 multiply()
659
660dbesselyn = ->
661 push(cadr(p1))
662 push(p2)
663 derivative()
664 push(cadr(p1))
665 push(caddr(p1))
666 push_integer(-1)
667 add()
668 bessely()
669 push(caddr(p1))
670 push_integer(-1)
671 multiply()
672 push(cadr(p1))
673 divide()
674 push(cadr(p1))
675 push(caddr(p1))
676 bessely()
677 multiply()
678 add()
679 multiply()
680
681derivative_of_integral = ->
682 push(cadr(p1))
683
684