XProfan

'**************************************************************
WindowTitle "Calculate Function Gamma & Ln(Gamma) with complex arguments"
'(D)2017-04 Demo per Fortran>Basic>XProfan-11 by P.Specht, Vienna/Austria
'Q: https://jean-pierre.moreau.pagesperso-orange.fr/Basic/mcgama_bas.txt
WindowStyle 24:Window 0,0-%maxx,%maxy:showmax:set("Decimals",11)
'* ---------------------------------------------------------- *
'* Purpose: This program computes the gamma function G(z)     *
'*          or Complex Ln[G(z)] for a complex argument using  *
'*          subroutine CGAMA                                  *
'* Input :  x  --- Real part of z                             *
'*          y  --- Imaginary part of z                        *
'*          KF --- Function code:                             *
'*          KF=0 for Ln[G(z)]                                 *
'*          KF=1 for G(z)                                     *
'* Output:  GR --- Real part of Ln[G(z)] or G(z)              *
'*          GI --- Imaginary part of Ln[G(z)] or G(z)         *
'* Examples:                                                  *
'*    x         y           Re[G(z)]           Im[G(z)]       *
'*  --------------------------------------------------------  *
'*   2.50      5.00     .2267360319E-01    -.1172284404E-01   *
'*   5.00     10.00     .1327696517E-01     .3639011746E-02   *
'*   2.50     -5.00     .2267360319E-01     .1172284404E-01   *
'*   5.00    -10.00     .1327696517E-01    -.3639011746E-02   *
'*                                                            *
'*    x         y          Re[LnG(z)]         Im[LnG(z)]      *
'*  --------------------------------------------------------  *
'*   2.50      5.00    -.3668103262E+01     .5806009801E+01   *
'*   5.00     10.00    -.4285507444E+01     .1911707090E+02   *
'*   2.50     -5.00    -.3668103262E+01    -.5806009801E+01   *
'*   5.00    -10.00    -.4285507444E+01    -.1911707090E+02   *
'* ---------------------------------------------------------- *
'* SAMPLE RUNS:                                               *
'* Please enter KF, x and y: 1, 2.5,5                         *
'*    x    y          Re[G(z)]             Im[G(z)]           *
'*  --------------------------------------------------------  *
'*   2.5   5    2.267360318980015e-02  -1.172284404171513e-02 *
'*                                                            *
'* Please enter KF, x and y: 0, 2.5,5                          *
'*    x    y         Re[LnG(z)]           Im[LnG(z)]          *
'*  --------------------------------------------------------  *
'*   2.5   5   -3.66810326235451        5.806009800636287     *
'* ---------------------------------------------------------- *
'*    "Fortran Routines for Computation of Special Functions, *
'*     jin.ece.uiuc.edu/routines/routines.html".              *
'*                   QuickBasic Release By J-P Moreau, Paris. *
'*                                (www.jpmoreau.fr)           *
'**************************************************************
declare x!,y!,KF&,expx!,GR!,GI!, SINH!,COSH!
'... und per subroutine GAMMA:
declare A![10],G0!,GR1!,GI1!,SR!,SI!,T!,TH!,TH1!,TH2!,x0!,X1!,Y1!,Z1!,Z2!,PI!,J&,K&,NA&
CLS
MAIN:
PRINT "\n  Please enter Flag [1=Gamma(),0=Ln(Gamma())], Re(z), Im(z): "
locate %csrlin+1,5: INPUT KF&
locate %csrlin-1,25:input X!
locate %csrlin-1,50:input Y!
PRINT

IF KF& = 1

    PRINT "    x                   y              Re[G(z)]                  Im[G(z)]"

ELSE

    PRINT "    x                   y              Re[LnG(z)]               Im[LnG(z)]"

ENDIF

PRINT "   ---------------------------------------------------------------------------------"
CGAMA'(X!,Y!,KF&,GR!,GI!) '= GOSUB 1000 'call CGAMA(X,Y,KF,GR,GI)
PRINT tab(5);format$("%g",x!),tab(25);format$("%g",Y!),tab(40);format$("%g",GR!),tab(65);format$("%g",GI!)
PRINT
Goto "Main"
END'of Main Program
'Auxiliary functions
'500:
'Function SINH(xx)

proc sinh :parameters xx!

    expx! = EXP(xx!)
    SINH! = 0.5 * (expx! - 1 / expx!)
    'RETURN

endproc

'600:
'Function COSH(xx)

proc cosh :parameters xx!

    expx! = EXP(xx!)
    COSH! = .5 * (expx! + 1 / expx!)
    'RETURN

endproc

'1000:
'Subroutine CGAMA '(X,Y,KF,GR,GI)

PROC CGAMA

    ' ===========================================================
    '       Purpose: Compute the gamma function G(z) or Ln[G(z)]
    '                for a complex argument
    '       Input :  x  --- Real part of z
    '                y  --- Imaginary part of z
    '                KF --- Function code:
    '                       KF=0 for Ln[G(z)]
    '                       KF=1 for G(z)
    '       Output:  GR --- Real part of Ln[G(z)] or G(z)
    '                GI --- Imaginary part of Ln[G(z)] or G(z)
    ' ===========================================================
    ':::DECLARE A![10] 'DIM A(10)
    ':::DECLARE G0!,GR1!,GI1!,SR!,SI!,T!,TH!,TH1!,TH2!,x0!,X1!,Y1!,Z1!,Z2!,PI!
    ':::DECLARE J&,K&,NA&  '(alle nach Main gehoben per GLOBAL declare)
    PI! = 4 * ArcTaN(1)
    'Initialize table A
    A![1] = 8.333333333333333*10^-2: A![2] = -2.777777777777778*10^-3
    A![3] = 7.936507936507937*10^-4: A![4] = -5.952380952380952*10^-4
    A![5] = 8.417508417508418*10^-4: A![6] = -1.917526917526918*10^-3
    A![7] = 6.41025641025641*10^-3 : A![8] = -2.955065359477124*10^-2
    A![9] = 0.1796443723688307     : A![10]= -1.3924322169059

    IF (Y! = 0) AND (X! = INT(X!)) AND (X! <= 0)

        GR! = 1E30'arbitrary big number
        GI! = 0
        RETURN

    ELSEIF X! < 0

        X1! = X!
        Y1! = Y!
        X! = -X!
        Y! = -Y!

    ENDIF

    X0! = X!

    IF X! <= 7

        NA& = INT(7 - X!)
        X0! = X! + NA&

    ENDIF

    Z1! = SQRt(X0! * X0! + Y! * Y!)
    TH! = ArcTaN(Y! / X0!)
    GR! = (X0! - 0.5) * Ln(Z1!) - TH! * Y! - X0! + 0.5 * Ln(2 * PI!)
    GI! = TH! * (X0! - 0.5) + Y! * Ln(Z1!) - Y!
    'FOR K& = 1 TO 10

    whileloop 10:k&=&Loop

        T! = Z1! ^ (1 - 2 * K&)
        GR! = GR! + A![K&] * T! * COS((2 * K& - 1) * TH!)
        GI! = GI! - A![K&] * T! * SIN((2 * K& - 1) * TH!)
        'NEXT K

    endwhile

    IF X! <= 7'THEN

        GR1! = 0
        GI1! = 0
        'FOR J = 0 TO NA - 1

        whileloop 0,na&-1:J&=&Loop

            GR1! = GR1! + 0.5 * Ln((X! + J&) ^ 2 + Y! * Y!)
            GI1! = GI1! + ArcTaN(Y! / (X! + J&))
            'NEXT J

        endwhile

        GR! = GR! - GR1!
        GI! = GI! - GI1!

    ENDIF

    IF X1! < 0'THEN

        Z1! = SQRt(X! * X! + Y! * Y!)
        TH1! = ArcTaN(Y! / X!)
        xx! = PI! * Y!
        ::: sinh(xx!)'= GOSUB 500
        ::: cosh(xx!)'= GOSUB 600
        SR! = -SIN(PI! * X!) * COSH!
        SI! = -COS(PI! * X!) * SINH!
        Z2! = SQRt(SR! * SR! + SI! * SI!)
        TH2! = ArcTaN(SI! / SR!)
        case SR! < 0: TH2! = PI! + TH2!
        GR! = Ln(PI! / (Z1! * Z2!)) - GR!
        GI! = -TH1! - TH2! - GI!
        X! = X1!
        Y! = Y1!

    ENDIF

    IF KF& = 1'THEN

        G0! = EXP(GR!)
        GR! = G0! * COS(GI!)
        GI! = G0! * SIN(GI!)

    ENDIF

    'RETURN

EndProc

'end of file mcgama.prf