330 lines
4.1 KiB
ArmAsm
330 lines
4.1 KiB
ArmAsm
.globl gamma, _gamma, signgam, _signgam
|
|
.globl log, sin
|
|
half = 040000
|
|
one = 40200
|
|
two = 40400
|
|
eight = 41000
|
|
large = 77777
|
|
ldfps = 170100^tst
|
|
stfps = 170200^tst
|
|
/
|
|
/ gamma computes the log of the abs of the gamma function.
|
|
/ gamma accepts its argument and returns its result
|
|
/ in fr0. The carry bit is set if the result is
|
|
/ too large to represent.
|
|
/ The sign of the gamma function is
|
|
/ returned in the globl cell signgam.
|
|
/
|
|
/ The coefficients for expansion around zero
|
|
/ are #5243 from Hart & Cheney; for expansion
|
|
/ around infinity they are #5404.
|
|
/
|
|
/ movf arg,fr0
|
|
/ jsr pc,gamma
|
|
/ movf fr0,...
|
|
/
|
|
|
|
_gamma:
|
|
mov r5,-(sp)
|
|
mov sp,r5
|
|
movf 4(r5),fr0
|
|
jsr pc,gamma
|
|
mov (sp)+,r5
|
|
rts pc
|
|
gamma:
|
|
stfps -(sp)
|
|
ldfps $200
|
|
clr signgam
|
|
movf fr1,-(sp)
|
|
tstf fr0
|
|
cfcc
|
|
ble negative
|
|
cmpf $eight,fr0
|
|
cfcc
|
|
blt asymptotic
|
|
jsr pc,regular
|
|
/
|
|
lret:
|
|
jsr pc,log
|
|
bec ret
|
|
4
|
|
ret:
|
|
movf (sp)+,fr1
|
|
ldfps (sp)+
|
|
clc
|
|
rts pc
|
|
/
|
|
erret:
|
|
movf $large,fr0
|
|
movf (sp)+,fr1
|
|
ldfps (sp)+
|
|
sec
|
|
rts pc
|
|
|
|
/
|
|
/ here for positive x > 8
|
|
/ the log of the gamma function is
|
|
/ approximated directly by the asymptotic series.
|
|
/
|
|
asymptotic:
|
|
movf fr0,-(sp)
|
|
movf fr0,fr1
|
|
jsr pc,log
|
|
subf $half,fr1
|
|
mulf fr1,fr0
|
|
subf (sp),fr0
|
|
addf goobie,fr0
|
|
/
|
|
movf $one,fr1
|
|
divf (sp)+,fr1
|
|
movf fr0,-(sp)
|
|
movf fr1,-(sp)
|
|
mulf fr1,fr1
|
|
/
|
|
mov r0,-(sp)
|
|
mov $p5p,r0
|
|
mov $5,-(sp)
|
|
movf *(r0)+,fr0
|
|
1:
|
|
mulf fr1,fr0
|
|
addf *(r0)+,fr0
|
|
dec (sp)
|
|
bne 1b
|
|
tst (sp)+
|
|
mov (sp)+,r0
|
|
mulf (sp)+,fr0
|
|
addf (sp)+,fr0
|
|
br ret
|
|
|
|
/
|
|
/ here on negative
|
|
/ the negative gamma is computed
|
|
/ in terms of the pos gamma.
|
|
/
|
|
negative:
|
|
absf fr0
|
|
movf fr0,fr1
|
|
mulf pi,fr0
|
|
jsr pc,sin
|
|
negf fr0
|
|
cfcc
|
|
beq erret
|
|
bgt 1f
|
|
inc signgam
|
|
absf fr0
|
|
1:
|
|
mov signgam,-(sp)
|
|
mulf fr1,fr0
|
|
divf pi,fr0
|
|
jsr pc,log
|
|
movf fr0,-(sp)
|
|
movf fr1,fr0
|
|
jsr pc,gamma
|
|
addf (sp)+,fr0
|
|
negf fr0
|
|
mov (sp)+,signgam
|
|
br ret
|
|
|
|
/
|
|
/ control comes here for arguments less than 8.
|
|
/ if the argument is 2<x<3 then compute by
|
|
/ a rational approximation.
|
|
/ if x<2 or x>3 then the argument
|
|
/ is reduced in range by the formula
|
|
/ gamma(x+1) = x*gamma(x)
|
|
/
|
|
regular:
|
|
subf $two,fr0
|
|
cfcc
|
|
bge 1f
|
|
addf $two,fr0
|
|
movf fr0,-(sp)
|
|
addf $one,fr0
|
|
movf fr0,-(sp)
|
|
addf $one,fr0
|
|
jsr pc,regular
|
|
divf (sp)+,fr0
|
|
divf (sp)+,fr0
|
|
rts pc
|
|
1:
|
|
cmpf $one,fr0
|
|
cfcc
|
|
bgt 1f
|
|
addf $one,fr0
|
|
movf fr0,-(sp)
|
|
subf $two,fr0
|
|
jsr pc,1b
|
|
mulf (sp)+,fr0
|
|
rts pc
|
|
1:
|
|
movf fr2,-(sp)
|
|
mov r0,-(sp)
|
|
mov $p4p,r0
|
|
mov $6,-(sp)
|
|
movf fr0,fr2
|
|
movf *(r0)+,fr0
|
|
1:
|
|
mulf fr2,fr0
|
|
addf *(r0)+,fr0
|
|
dec (sp)
|
|
bne 1b
|
|
mov $7,(sp)
|
|
movf fr2,fr1
|
|
br 2f
|
|
1:
|
|
mulf fr2,fr1
|
|
2:
|
|
addf *(r0)+,fr1
|
|
dec (sp)
|
|
bne 1b
|
|
tst (sp)+
|
|
mov (sp)+,r0
|
|
divf fr1,fr0
|
|
movf (sp)+,fr2
|
|
rts pc
|
|
/
|
|
.data
|
|
p4p:
|
|
p6;p5;p4;p3;p2;p1;p0
|
|
q6;q5;q4;q3;q2;q1;q0
|
|
|
|
/ p6 = -.67449 50724 59252 89918 d1
|
|
/ p5 = -.50108 69375 29709 53015 d2
|
|
/ p4 = -.43933 04440 60025 67613 d3
|
|
/ p3 = -.20085 27401 30727 91214 d4
|
|
/ p2 = -.87627 10297 85214 89560 d4
|
|
/ p1 = -.20886 86178 92698 87364 d5
|
|
/ p0 = -.42353 68950 97440 89647 d5
|
|
/ q7 = 1.0 d0
|
|
/ q6 = -.23081 55152 45801 24562 d2
|
|
/ q5 = +.18949 82341 57028 01641 d3
|
|
/ q4 = -.49902 85266 21439 04834 d3
|
|
/ q3 = -.15286 07273 77952 20248 d4
|
|
/ q2 = +.99403 07415 08277 09015 d4
|
|
/ q1 = -.29803 85330 92566 49932 d4
|
|
/ q0 = -.42353 68950 97440 90010 d5
|
|
p1:
|
|
143643
|
|
26671
|
|
36161
|
|
72154
|
|
p2:
|
|
143410
|
|
165327
|
|
54121
|
|
172630
|
|
p3:
|
|
142773
|
|
10340
|
|
74264
|
|
152066
|
|
p4:
|
|
142333
|
|
125113
|
|
176657
|
|
75740
|
|
p5:
|
|
141510
|
|
67515
|
|
65111
|
|
24263
|
|
p6:
|
|
140727
|
|
153242
|
|
163350
|
|
32217
|
|
p0:
|
|
144045
|
|
70660
|
|
101665
|
|
164444
|
|
q1:
|
|
143072
|
|
43052
|
|
50302
|
|
136745
|
|
q2:
|
|
43433
|
|
50472
|
|
145404
|
|
175462
|
|
q3:
|
|
142677
|
|
11556
|
|
144553
|
|
154177
|
|
q4:
|
|
142371
|
|
101646
|
|
141245
|
|
11264
|
|
q5:
|
|
42075
|
|
77614
|
|
43022
|
|
27573
|
|
q6:
|
|
141270
|
|
123404
|
|
76130
|
|
12650
|
|
q0:
|
|
144045
|
|
70660
|
|
101665
|
|
164444
|
|
|
|
p5p:
|
|
s5;s4;s3;s2;s1;s0
|
|
/
|
|
/ s5 = -.16334 36431 d-2
|
|
/ s4 = +.83645 87892 2 d-3
|
|
/ s3 = -.59518 96861 197 d-3
|
|
/ s2 = +.79365 05764 93454 d-3
|
|
/ s1 = -.27777 77777 35865 004 d-2
|
|
/ s0 = +.83333 33333 33331 01837 d-1
|
|
/ goobie = 0.91893 85332 04672 74178 d0
|
|
s5:
|
|
135726
|
|
14410
|
|
15074
|
|
17706
|
|
s4:
|
|
35533
|
|
42714
|
|
111634
|
|
76770
|
|
s3:
|
|
135434
|
|
3200
|
|
171173
|
|
156141
|
|
s2:
|
|
35520
|
|
6375
|
|
12373
|
|
111437
|
|
s1:
|
|
136066
|
|
5540
|
|
132625
|
|
63540
|
|
s0:
|
|
37252
|
|
125252
|
|
125252
|
|
125047
|
|
goobie:
|
|
40153
|
|
37616
|
|
41445
|
|
172645
|
|
pi:
|
|
40511
|
|
7732
|
|
121041
|
|
64302
|
|
.bss
|
|
_signgam:
|
|
signgam:.=.+2
|