1 | ###
|
2 | Legendre function
|
3 |
|
4 | Example
|
5 |
|
6 | legendre(x,3,0)
|
7 |
|
8 | Result
|
9 |
|
10 | 5 3 3
|
11 | --- x - --- x
|
12 | 2 2
|
13 |
|
14 | The computation uses the following recurrence relation.
|
15 |
|
16 | P(x,0) = 1
|
17 |
|
18 | P(x,1) = x
|
19 |
|
20 | n*P(x,n) = (2*(n-1)+1)*x*P(x,n-1) - (n-1)*P(x,n-2)
|
21 |
|
22 | In the "for" loop we have i = n-1 so the recurrence relation becomes
|
23 |
|
24 | (i+1)*P(x,n) = (2*i+1)*x*P(x,n-1) - i*P(x,n-2)
|
25 |
|
26 | For m > 0
|
27 |
|
28 | P(x,n,m) = (-1)^m * (1-x^2)^(m/2) * d^m/dx^m P(x,n)
|
29 | ###
|
30 |
|
31 |
|
32 |
|
33 | Eval_legendre = ->
|
34 | # 1st arg
|
35 |
|
36 | push(cadr(p1))
|
37 | Eval()
|
38 |
|
39 | # 2nd arg
|
40 |
|
41 | push(caddr(p1))
|
42 | Eval()
|
43 |
|
44 | # 3rd arg (optional)
|
45 |
|
46 | push(cadddr(p1))
|
47 | Eval()
|
48 |
|
49 | p2 = pop()
|
50 | if (p2 == symbol(NIL))
|
51 | push_integer(0)
|
52 | else
|
53 | push(p2)
|
54 |
|
55 | legendre()
|
56 |
|
57 | #define X p1
|
58 | #define N p2
|
59 | #define M p3
|
60 | #define Y p4
|
61 | #define Y0 p5
|
62 | #define Y1 p6
|
63 |
|
64 |
|
65 | legendre = ->
|
66 | save()
|
67 | __legendre()
|
68 | restore()
|
69 |
|
70 | __legendre = ->
|
71 | m = 0
|
72 | n = 0
|
73 |
|
74 | p3 = pop()
|
75 | p2 = pop()
|
76 | p1 = pop()
|
77 |
|
78 | push(p2)
|
79 | n = pop_integer()
|
80 |
|
81 | push(p3)
|
82 | m = pop_integer()
|
83 |
|
84 | if (n < 0 || isNaN(n) || m < 0 || isNaN(m))
|
85 | push_symbol(LEGENDRE)
|
86 | push(p1)
|
87 | push(p2)
|
88 | push(p3)
|
89 | list(4)
|
90 | return
|
91 |
|
92 | if (issymbol(p1))
|
93 | __legendre2(n, m)
|
94 | else
|
95 | p4 = p1; # do this when X is an expr
|
96 | p1 = symbol(SECRETX)
|
97 | __legendre2(n, m)
|
98 | p1 = p4
|
99 | push(symbol(SECRETX))
|
100 | push(p1)
|
101 | subst()
|
102 | Eval()
|
103 |
|
104 | __legendre3(m)
|
105 |
|
106 | __legendre2 = (n, m) ->
|
107 | i = 0
|
108 |
|
109 | push_integer(1)
|
110 | push_integer(0)
|
111 |
|
112 | p6 = pop()
|
113 |
|
114 | # i=1 p5 = 0
|
115 | # p6 = 1
|
116 | # ((2*i+1)*x*p6 - i*p5) / i = x
|
117 | #
|
118 | # i=2 p5 = 1
|
119 | # p6 = x
|
120 | # ((2*i+1)*x*p6 - i*p5) / i = -1/2 + 3/2*x^2
|
121 | #
|
122 | # i=3 p5 = x
|
123 | # p6 = -1/2 + 3/2*x^2
|
124 | # ((2*i+1)*x*p6 - i*p5) / i = -3/2*x + 5/2*x^3
|
125 |
|
126 | for i in [0...n]
|
127 |
|
128 | p5 = p6
|
129 |
|
130 | p6 = pop()
|
131 |
|
132 | push_integer(2 * i + 1)
|
133 | push(p1)
|
134 | multiply()
|
135 | push(p6)
|
136 | multiply()
|
137 |
|
138 | push_integer(i)
|
139 | push(p5)
|
140 | multiply()
|
141 |
|
142 | subtract()
|
143 |
|
144 | push_integer(i + 1)
|
145 | divide()
|
146 |
|
147 | for i in [0...m]
|
148 | push(p1)
|
149 | derivative()
|
150 |
|
151 | # moveTos tos * (-1)^m * (1-x^2)^(m/2)
|
152 |
|
153 | __legendre3 = (m) ->
|
154 | if (m == 0)
|
155 | return
|
156 |
|
157 | if (car(p1) == symbol(COS))
|
158 | push(cadr(p1))
|
159 | sine()
|
160 | square()
|
161 | else if (car(p1) == symbol(SIN))
|
162 | push(cadr(p1))
|
163 | cosine()
|
164 | square()
|
165 | else
|
166 | push_integer(1)
|
167 | push(p1)
|
168 | square()
|
169 | subtract()
|
170 |
|
171 | push_integer(m)
|
172 | push_rational(1, 2)
|
173 | multiply()
|
174 | power()
|
175 | multiply()
|
176 |
|
177 | if (m % 2)
|
178 | negate()
|
179 |
|