UNIX 6th Ed. - /usr/source/s3/gamma.s

   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

Defined functions

Defined variables