home *** CD-ROM | disk | FTP | other *** search
/ minnie.tuhs.org / unixen.tar / unixen / PDP-11 / Trees / V6 / usr / source / s3 / gamma.s < prev    next >
Encoding:
Text File  |  1975-05-13  |  4.1 KB  |  330 lines
  1. .globl gamma, _gamma, signgam, _signgam
  2. .globl    log, sin
  3. half = 040000
  4. one = 40200
  5. two = 40400
  6. eight = 41000
  7. large = 77777
  8. ldfps = 170100^tst
  9. stfps = 170200^tst
  10. /
  11. /    gamma computes the log of the abs of the gamma function.
  12. /    gamma accepts its argument and returns its result
  13. /    in fr0.  The carry bit is set if the result is
  14. /    too large to represent.
  15. /    The sign of the gamma function is
  16. /    returned in the globl cell signgam.
  17. /
  18. /    The coefficients for expansion around zero
  19. /    are #5243 from Hart & Cheney; for expansion
  20. /    around infinity they are #5404.
  21. /
  22. /    movf    arg,fr0
  23. /    jsr    pc,gamma
  24. /    movf    fr0,...
  25. /
  26.  
  27. _gamma:
  28.     mov    r5,-(sp)
  29.     mov    sp,r5
  30.     movf    4(r5),fr0
  31.     jsr    pc,gamma
  32.     mov    (sp)+,r5
  33.     rts    pc
  34. gamma:
  35.     stfps    -(sp)
  36.     ldfps    $200
  37.     clr    signgam
  38.     movf    fr1,-(sp)
  39.     tstf    fr0
  40.     cfcc
  41.     ble    negative
  42.     cmpf    $eight,fr0
  43.     cfcc
  44.     blt    asymptotic
  45.     jsr    pc,regular
  46. /
  47. lret:
  48.     jsr    pc,log
  49.     bec    ret
  50.     4
  51. ret:
  52.     movf    (sp)+,fr1
  53.     ldfps    (sp)+
  54.     clc
  55.     rts    pc
  56. /
  57. erret:
  58.     movf    $large,fr0
  59.     movf    (sp)+,fr1
  60.     ldfps    (sp)+
  61.     sec
  62.     rts    pc
  63.  
  64. /
  65. /    here for positive x > 8
  66. /    the log of the gamma function is
  67. /    approximated directly by the asymptotic series.
  68. /
  69. asymptotic:
  70.     movf    fr0,-(sp)
  71.     movf    fr0,fr1
  72.     jsr    pc,log
  73.     subf    $half,fr1
  74.     mulf    fr1,fr0
  75.     subf    (sp),fr0
  76.     addf    goobie,fr0
  77. /
  78.     movf    $one,fr1
  79.     divf    (sp)+,fr1
  80.     movf    fr0,-(sp)
  81.     movf    fr1,-(sp)
  82.     mulf    fr1,fr1
  83. /
  84.     mov    r0,-(sp)
  85.     mov    $p5p,r0
  86.     mov    $5,-(sp)
  87.     movf    *(r0)+,fr0
  88. 1:
  89.     mulf    fr1,fr0
  90.     addf    *(r0)+,fr0
  91.     dec    (sp)
  92.     bne    1b
  93.     tst    (sp)+
  94.     mov    (sp)+,r0
  95.     mulf    (sp)+,fr0
  96.     addf    (sp)+,fr0
  97.     br    ret
  98.  
  99. /
  100. /    here on negative
  101. /    the negative gamma is computed
  102. /    in terms of the pos gamma.
  103. /
  104. negative:
  105.     absf    fr0
  106.     movf    fr0,fr1
  107.     mulf    pi,fr0
  108.     jsr    pc,sin
  109.     negf    fr0
  110.     cfcc
  111.     beq    erret
  112.     bgt    1f
  113.     inc    signgam
  114.     absf    fr0
  115. 1:
  116.     mov    signgam,-(sp)
  117.     mulf    fr1,fr0
  118.     divf    pi,fr0
  119.     jsr    pc,log
  120.     movf    fr0,-(sp)
  121.     movf    fr1,fr0
  122.     jsr    pc,gamma
  123.     addf    (sp)+,fr0
  124.     negf    fr0
  125.     mov    (sp)+,signgam
  126.     br    ret
  127.  
  128. /
  129. /    control comes here for arguments less than 8.
  130. /    if the argument is 2<x<3 then compute by
  131. /    a rational approximation.
  132. /    if x<2 or x>3 then the argument
  133. /    is reduced in range by the formula
  134. /    gamma(x+1) = x*gamma(x)
  135. /
  136. regular:
  137.     subf    $two,fr0
  138.     cfcc
  139.     bge    1f
  140.     addf    $two,fr0
  141.     movf    fr0,-(sp)
  142.     addf    $one,fr0
  143.     movf    fr0,-(sp)
  144.     addf    $one,fr0
  145.     jsr    pc,regular
  146.     divf    (sp)+,fr0
  147.     divf    (sp)+,fr0
  148.     rts    pc
  149. 1:
  150.     cmpf    $one,fr0
  151.     cfcc
  152.     bgt    1f
  153.     addf    $one,fr0
  154.     movf    fr0,-(sp)
  155.     subf    $two,fr0
  156.     jsr    pc,1b
  157.     mulf    (sp)+,fr0
  158.     rts    pc
  159. 1:
  160.     movf    fr2,-(sp)
  161.     mov    r0,-(sp)
  162.     mov    $p4p,r0
  163.     mov    $6,-(sp)
  164.     movf    fr0,fr2
  165.     movf    *(r0)+,fr0
  166. 1:
  167.     mulf    fr2,fr0
  168.     addf    *(r0)+,fr0
  169.     dec    (sp)
  170.     bne    1b
  171.     mov    $7,(sp)
  172.     movf    fr2,fr1
  173.     br    2f
  174. 1:
  175.     mulf    fr2,fr1
  176. 2:
  177.     addf    *(r0)+,fr1
  178.     dec    (sp)
  179.     bne    1b
  180.     tst    (sp)+
  181.     mov    (sp)+,r0
  182.     divf    fr1,fr0
  183.     movf    (sp)+,fr2
  184.     rts    pc
  185. /
  186. .data
  187. p4p:
  188.     p6;p5;p4;p3;p2;p1;p0
  189.     q6;q5;q4;q3;q2;q1;q0
  190.  
  191. /    p6 = -.67449 50724 59252 89918 d1
  192. /    p5 = -.50108 69375 29709 53015 d2
  193. /    p4 = -.43933 04440 60025 67613 d3
  194. /    p3 = -.20085 27401 30727 91214 d4
  195. /    p2 = -.87627 10297 85214 89560 d4
  196. /    p1 = -.20886 86178 92698 87364 d5
  197. /    p0 = -.42353 68950 97440 89647 d5
  198. /    q7 = 1.0 d0
  199. /    q6 = -.23081 55152 45801 24562 d2
  200. /    q5 = +.18949 82341 57028 01641 d3
  201. /    q4 = -.49902 85266 21439 04834 d3
  202. /    q3 = -.15286 07273 77952 20248 d4
  203. /    q2 = +.99403 07415 08277 09015 d4
  204. /    q1 = -.29803 85330 92566 49932 d4
  205. /    q0 = -.42353 68950 97440 90010 d5
  206. p1:
  207.     143643
  208.     26671
  209.     36161
  210.     72154
  211. p2:
  212.     143410
  213.     165327
  214.     54121
  215.     172630
  216. p3:
  217.     142773
  218.     10340
  219.     74264
  220.     152066
  221. p4:
  222.     142333
  223.     125113
  224.     176657
  225.     75740
  226. p5:
  227.     141510
  228.     67515
  229.     65111
  230.     24263
  231. p6:
  232.     140727
  233.     153242
  234.     163350
  235.     32217
  236. p0:
  237.     144045
  238.     70660
  239.     101665
  240.     164444
  241. q1:
  242.     143072
  243.     43052
  244.     50302
  245.     136745
  246. q2:
  247.     43433
  248.     50472
  249.     145404
  250.     175462
  251. q3:
  252.     142677
  253.     11556
  254.     144553
  255.     154177
  256. q4:
  257.     142371
  258.     101646
  259.     141245
  260.     11264
  261. q5:
  262.     42075
  263.     77614
  264.     43022
  265.     27573
  266. q6:
  267.     141270
  268.     123404
  269.     76130
  270.     12650
  271. q0:
  272.     144045
  273.     70660
  274.     101665
  275.     164444
  276.  
  277. p5p:
  278.     s5;s4;s3;s2;s1;s0
  279. /
  280. /    s5 = -.16334 36431 d-2
  281. /    s4 = +.83645 87892 2 d-3
  282. /    s3 = -.59518 96861 197 d-3
  283. /    s2 = +.79365 05764 93454 d-3
  284. /    s1 = -.27777 77777 35865 004 d-2
  285. /    s0 = +.83333 33333 33331 01837 d-1
  286. /    goobie = 0.91893 85332 04672 74178 d0
  287. s5:
  288.     135726
  289.     14410
  290.     15074
  291.     17706
  292. s4:
  293.     35533
  294.     42714
  295.     111634
  296.     76770
  297. s3:
  298.     135434
  299.     3200
  300.     171173
  301.     156141
  302. s2:
  303.     35520
  304.     6375
  305.     12373
  306.     111437
  307. s1:
  308.     136066
  309.     5540
  310.     132625
  311.     63540
  312. s0:
  313.     37252
  314.     125252
  315.     125252
  316.     125047
  317. goobie:
  318.     40153
  319.     37616
  320.     41445
  321.     172645
  322. pi:
  323.     40511
  324.     7732
  325.     121041
  326.     64302
  327. .bss
  328. _signgam:
  329. signgam:.=.+2
  330.