DECLARE SUB CISIA (X#, CI#, SI#)
'******************************************************************
'*      Purpose: This program computes the cosine and sine        *
'*               integrals using subroutine CISIA                 *
'*      Input :  x  --- Argument of Ci(x) and Si(x)               *
'*      Output:  CI --- Ci(x)                                     *
'*               SI --- Si(x)                                     *
'*      Example:                                                  *
'*                   x         Ci(x)          Si(x)               *
'*                ------------------------------------            *
'*                  0.0     - i             .00000000             *
'*                  5.0     -.19002975     1.54993124             *
'*                 10.0     -.04545643     1.65834759             *
'*                 20.0      .04441982     1.54824170             *
'*                 30.0     -.03303242     1.56675654             *
'*                 40.0      .01902001     1.58698512             *
'* -------------------------------------------------------------- *
'* REFERENCE: "Fortran Routines for Computation of Special        *
'*             Functions jin.ece.uiuc.edu/routines/routines.html" *
'*                                                                *
'*                            Basic Release By J-P Moreau, Paris. *
'*                                     (www.jpmoreau.fr)          *
'******************************************************************
DEFDBL A-H, O-Z
DEFINT I-N
'   PROGRAM MCISIA
	PRINT
	INPUT " Please enter x: ", X
	PRINT "   x         Ci(x)          Si(x)   "
	PRINT "------------------------------------"
	
	CALL CISIA(X, CI, SI)
	
	IF X <> 0! THEN PRINT " "; X; " "; CI; " "; SI
	IF X = 0! THEN PRINT " "; X; " i"; " O"
	END

'end of file mcisia.bas

	SUB CISIA (X, CI, SI)

'       =============================================
'       Purpose: Compute cosine and sine integrals
'                Si(x) and Ci(x)  ( x ò 0 )
'       Input :  x  --- Argument of Ci(x) and Si(x)
'       Output:  CI --- Ci(x)
'                SI --- Si(x)
'       =============================================

	DIM BJ(100)
	P2 = 1.570796326794897#
	EL = .5772156649015329#
	EPS = .000000000000001#
	X2 = X * X
	IF (X = 0!) THEN
	   CI = -1000000000000000#
	   SI = 0!
	ELSEIF (X <= 16!) THEN
	   XR = -.25 * X2
	   CI = EL + LOG(X) + XR
	   FOR K = 2 TO 40
	      XR = -.5 * XR * (K - 1) / (K * K * (2 * K - 1)) * X2
	      CI = CI + XR
	      IF (ABS(XR) < ABS(CI) * EPS) THEN GOTO 15
	   NEXT K
15         XR = X
	   SI = X
	   FOR K = 1 TO 40
	      XR = -.5# * XR * (2 * K - 1) / K / (4 * K * K + 4 * K + 1) * X2
	      SI = SI + XR
	  IF (ABS(XR) < ABS(SI) * EPS) THEN GOTO 20'return
	   NEXT K
	ELSEIF (X <= 32!) THEN
	   M = INT(47.2 + .82 * X)
	   XA1 = 0!
	   XA0 = .000000000000001#
	   FOR K = M TO 1 STEP -1
	      XA = 4! * K * XA0 / X - XA1
	      BJ(K) = XA
	      XA1 = XA0
	      XA0 = XA
	   NEXT K
	   XS = BJ(1)
	   FOR K = 3 TO M STEP 2
	      XS = XS + 2! * BJ(K)
	   NEXT K
	   BJ(1) = BJ(1) / XS
	   FOR K = 2 TO M
	      BJ(K) = BJ(K) / XS
	   NEXT K
	   XR = 1!
	   XG1 = BJ(1)
	   FOR K = 2 TO M
	      XR = .25 * XR * (2! * K - 3!) ^ 2 / ((K - 1!) * (2! * K - 1!) ^ 2) * X
	      XG1 = XG1 + BJ(K) * XR
	   NEXT K
	   XR = 1!
	   XG2 = BJ(1)
	   FOR K = 2 TO M
	      XR = .25 * XR * (2! * K - 5!) ^ 2 / ((K - 1!) * (2! * K - 3!) ^ 2) * X
	      XG2 = XG2 + BJ(K) * XR
	   NEXT K
	   XCS = COS(X / 2!)
	   XSS = SIN(X / 2!)
	   CI = EL + LOG(X) - X * XSS * XG1 + 2 * XCS * XG2 - 2 * XCS * XCS
	   SI = X * XCS * XG1 + 2 * XSS * XG2 - SIN(X)
	ELSE
	   XR = 1!
	   XF = 1!
		   FOR K = 1 TO 9
	      XR = -2! * XR * K * (2 * K - 1) / X2
	      XF = XF + XR
	   NEXT K
	   XR = 1! / X
	   XG = XR
	   FOR K = 1 TO 8
	      XR = -2! * XR * (2 * K + 1) * K / X2
	  XG = XG + XR
	   NEXT K
	   CI = XF * SIN(X) / X - XG * COS(X) / X
	   SI = P2 - XF * COS(X) / X - XG * SIN(X) / X
	END IF
20      END SUB