; ; 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 fact label tbyte include fact.inc 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, Factorial; .code ; ----------------------------------------------- ; 'C' calling conventions: ; double y = Factorial(int n) ; ----------------------------------------------- Factorial proc arg n:dword uses ebx xor eax,eax mov eax,n shl eax,1 mov ebx,eax shl eax,2 add ebx,eax fld tbyte ptr fact[ebx] ret Factorial endp ; ----------------------------------------------- ; '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