;
;	Math function library
;
;	(C) 1993, C. Bond.  All rights reserved.
;
    .486
    .model flat, C

	.data
;   Note that the following coefficients are sqrt(2*pi) times
;	the coefficients of the term in the Lanczos expansion for
;	(small) gamma = 8, and order = 8.
gcoeff	dt	 2.50662827463085895669
	dt 	 4951.52807661845316286
	dt	-11022.6030401376228508
	dt	 8679.53339641626425256
	dt	-2900.13167318763062839
	dt	 387.369697577684282966
	dt	-15.6756300917512932238
	dt	 8.683645856906762224733e-2
	dt	-1.7105384786443111054e-6

gm      dt  7.5 ; gamma - 0.5
half	dt	0.5
fpf     dq  5.5

BigNum	dq	1e308

lgcoeff dq	1.000000000190015
	dq	76.18009172947146
	dq	-86.5053203291677
	dq	24.01409824083091
	dq	-1.231739572450155
	dq	0.1208650973866179e-2
	dq	-0.5395239384953e-5
lgconst	dq	2.5066282746310008

    public Gamma, LnGamma

	.code
;	-----------------------------------------------
;	'C' calling conventions:
;		double y = Gamma(double x)
;	-----------------------------------------------

;	This is an implementation of the precision gamma
;	function approximation by C. Lanczos. It is
;	valid for all real arguments x, where
;	? < x < 171+, and x <> 0, -1, -2, . . .
;
;	The approximation can also be used for complex
;	arguments, 'z'. (See CGamma.)
;
;	Coefficients are computed for the 8th order
;	approximation and are stored in extended real
;	format.
;

Gamma   proc
	arg x:qword
    uses ebx
	local status_old:word,status_new:word,sign_flg:word
	
;
;	Check sign of argument, and adjust 'x' if necessary.
;	Use 386/486 instead of co-processor for speed.
;
    mov ax,word ptr [x+6]
	and ax,08000h
	mov sign_flg,ax
	xor word ptr [x+6],ax
;
;	Sign flag is zero if argument is non-negative
;	Argument now equals |argument|
;
;   Initialize coefficient array pointer [ebx]
;	assuming 'ds' points to data segment.


    mov ebx,offset tbyte ptr gcoeff
                            ;-----------regs---------------
                            ; 0   1   2   3   4   5   6   7
                            ;------------------------------
    fld tbyte ptr [ebx]     ; C0
    fld1                    ; 1   C0
    fld x                   ; x   1   C0
	mov cx,8
;
;	compute the series part of the solution
;
gloop:
    add ebx,10
    fld tbyte ptr [ebx]     ; C1  x   1   C0
    fdiv st,st(1)           ; k1  x   1   C0
	faddp st(3),st			; x   1   C
	fadd st,st(1)			; x++ 1   C
	dec cx
	jne gloop

    fstp st(0)              ; 1   C
    fstp st(0)              ; C
    fld x                   ; x   C
	fld tbyte ptr gm		; gm  x   C
    fadd                    ; K   C
    fldln2                  ; l2  K   C

    fld st(1)               ; K   l2  K   C
    fyl2x                   ; le  K   C
    fld x                   ; x   le  K   C
	fld tbyte ptr half		; .5  x   le  C
    fsub                    ; x-  le  K   C
    fmul                    ; K2  K   C
    fsubr                   ; K3  C

;
;                                      A   B
;       Multiply argument by log e so 2 = e
;                               2
;
    fldl2e                  ; l2  K3  C
    fmul                    ; K4  C

;
;	Set rounding control to round down to handle negative numbers properly
;	
	fstcw status_old
	fstcw status_new
	and status_new, 0f3ffh
	or status_new, 0400h
	fldcw status_new

    fld st(0)               ; K4  K4  C 
    frndint                 ; KI  K4  C
	fldcw status_old	
	fsub st(1),st			; KI  KF  C
    fxch                    ; KF  KI  C
    f2xm1                   ; 2^- KI  C
    fld1                    ; 1   2^- KI  C
    fadd                    ; 2^  KI  C
    fscale                  ; y'  KI  C
    fxch                    ; KI  y'  C
    fstp st(0)              ; y'  C
    fmul                    ; y
	test sign_flg,08000h
	jz gx
;
;	Analytic continuation of gamma function for
;	negative arguments.
;	Gamma(-x) = pi * Csc(pi * x)/(-x * Gamma(x))
;
    fld x                   ; x   y'
    fmul                    ; d
    fchs                    ; -d
    fldpi                   ; pi  -d
    fmul x                  ; K6  -d
    fsin                    ; K7  -d
    fmul                    ; d'
    fldpi                   ; pi  d'
    fdivr                   ; y
gx:
	ret
Gamma	endp

;	-----------------------------------------------
;	'C' calling conventions:
;		double y = LnGamma(double x)
;	-----------------------------------------------
LnGamma proc
	arg x:qword
    uses ebx

                            ;-----------regs---------------
                            ; 0   1   2   3   4   5   6   7
                            ;------------------------------
    fld x                   ; x
    fld st(0)               ; x   x
    fadd qword ptr fpf      ; K1  x
    fldln2                  ; l2  K1  x
    fld st(1)               ; K1  l2  K1  x
    fyl2x                   ; K2  K1  x
    fld half                ; .5  K2  K1  x
    fadd st,st(3)           ; K3  K2  K1  x
    fmul                    ; K4  K1  x
    fsub                    ; K5  x
    mov ebx, offset qword ptr lgcoeff
    fld qword ptr [ebx]     ; C0  K5  x
    fld1                    ; 1   C0  K5  x
    fld st(3)               ; x   1   C0  K5  x 
	mov cx,6
lg1:
    add ebx,8
	fadd st(0),st(1)		; x++ 1   C0  K5  x
    fld qword ptr [ebx]     ; C1  x++ 1   C0  K5  x
	fdiv st(0),st(1)		; k1  x++ 1   C0  K5  x
	faddp st(3),st			; x++ 1   C   K5  x
    dec cx                  ;
    jne lg1                 ;
                            ;
    fstp st(0)              ; 1   C   K5  x 
    fstp st(0)              ; C   K5  x
	fmul lgconst			; C1  K5  x
	fdiv st(0),st(2)		; C2  K5  x
    ffree st(2)             ; C2  K5
    fldln2                  ; l2  C2  K5 
    fxch                    ; C2  l2  K5 
    fyl2x                   ; C3  K5
    fsubr                   ; y
	ret 
LnGamma endp
	end