Post Reply 
41 MCODE - Hermite and Lagrange Polynomials
12-25-2014, 06:00 PM (This post was last modified: 12-27-2014 06:25 AM by Ángel Martin.)
Post: #1
41 MCODE - Hermite and Lagrange Polynomials
As usual when you thought a project was done new ideas pop up that just can't be ignored... like this sweet & short MCODE implementation for the Hermite and Lagrange orthogonal polynomials using the recurrence approach - as an afterthought to be included in the SandMatrix module.


0D4    "T"    Hermite Polynomials
04D    "M"    n,x in {Y,X}
048    "H"    Ángel Martin
284    CLRF 7    
02B    JNC +05    
0C7    "G"    Legendre Polynomials
045    "E"    n,x in {Y,X}
04C    "L"    Ángel Martin
288    SETF 7    
1A5    ?NC XQ    Check for valid entries
100    ->4069    [CHKST2] - sets DEC mode
04E    C=0  ALL    
35C    PT=12         C= 1
050    LD@PT- 1    
168    WRIT 5(M)    
0F8    READ 3(X)    x
10E    A=C ALL    
28C    ?FSET 7    Legendre?
01D    ?NC XQ    no, ergo Hermite
060    ->1807    [AD2_10]
1A8    WRIT 6(N)    2x or x
0B8    READ 2(Y)    n
088    SETF 5    do integer portion
0ED    ?NC XQ    leaves result in 13-digit form
064    ->193B    [INTFRC] - doesn't need DEC
070    N=C ALL    update counter
2FA    ?C#0 M    n = 0?
18B    JNC +49d    yes, result in (5)M
009    ?NC XQ    {A,B} = {A,B}-1
060    ->1802    [SUBONE]
2FA    ?C#0 M    was n=1 ?
19B    JNC +51d    yes, result in (6)N
28C    ?FSET 7    Legendre?
025    ?NC XQ          No, 2 (n-1)
060    ->1809    [AD1-10 ]
178    READ 5(M)    H(n-2)
13D    ?NC XQ      
060    ->184F    [MP1_10]
2BE    C=-C-1 MS    sign change
11E    A=C MS    ditto in 13-digit form
089    ?NC XQ    
064    ->1922    [STSCR]
28C    ?FSET 7    Legendre?
05B    JNC +11d    
0B0    C=N ALL    n
10E    A=C ALL    
01D    ?NC XQ    2n
060    ->1807    [AD2_10]
009    ?NC XQ    2n-1
060    ->1802    [SUBONE]
0F8    READ 3(X)    x
13D    ?NC XQ      (2n-1)*x
060    ->184F    [MP1_10]
02B    JNC +05    
0F8    READ 3(X)    x
10E    A=C ALL    
01D    ?NC XQ    2x
060    ->1807    [AD2_10]
1B8    READ 6(N)    H(n-1) or: P(n-1)
168    WRIT 5(M)    the new H(n-2) / P(n-2)
13D    ?NC XQ      2x* H(n-1) or: (2n-1)*x*P(n-1)
060    ->184F    [MP1_10]
0D1    ?NC XQ    recall first term
064    ->1934    [RCSCR]
031    ?NC XQ    add to complete it
060    ->180C    [AD2-13]
28C    ?FSET 7    Legendre?
023    JNC +04    no, skip
0B0    C=N ALL    n
269    ?NC XQ    
060    ->189A    [DV1-10]
1A8    WRIT 6(N)    the new H(n-1)
0B0    C=N ALL    
1FD    ?NC XQ    {A,B} = C-1
100    ->407F    [DECC10] - sets DEC
273    JNC - 50d    
178    READ 5(M)    
070    N=C ALL    parameter passing
260    SETHEX    
3AD    PORT DEP:    Abandon ship
08C    GO    in orderly fashion
13E    ->AD3E    [NFRX2]
1B8    READ 6(N)    
3D3    JNC -06

Only valid for integer orders. More complex formulas exist for non-integer cases, see JM Baillard's pages for details:

Note that the listing above is for a sub-function sitting in a bank-switched page.
Comments and suggestions for improvement are welcome.

Enjoy, as usual ;-)
Find all posts by this user
Quote this message in a reply
Post Reply 

User(s) browsing this thread: 1 Guest(s)