
Begin["System`"]

Unprotect[SinIntegral, CosIntegral, SinhIntegral, CoshIntegral,
	FresnelC, FresnelS];

SinIntegral::usage = "SinIntegral[x] gives the sine integral
Integrate[Sin[t]/t, {t, 0, x}]."

CosIntegral::usage = "CosIntegral[x] gives the cosine integral
EulerGamma + Log[x] + Integrate[(Cos[t] - 1)/t, {t, 0, x}]."

SinhIntegral::usage = "SinhIntegral[x] gives the hyperbolic sine integral
Integrate[Sinh[t]/t, {t, 0, x}]."

CoshIntegral::usage = "CoshIntegral[x] gives the hyperbolic cosine integral
EulerGamma + Log[x] + Integrate[(Cosh[t] - 1)/t, {t, 0, x}]."

FresnelS::usage = "FresnelS[x] gives the Fresnel integral
S[x] = Integrate[Sin[Pi t^2/2], {t, 0, x}]."

FresnelC::usage = "FresnelC[x] gives the Fresnel integral
C[x] = Integrate[Cos[Pi t^2/2], {t, 0, x}]."

Begin["`Private`"]
SinIntegral[z_]:=	(* evaluated for inexact z. *)
    Module[{x, y},
	x = Re[z];
	y = Im[z];
	Which[
	    z == 0,	z,
	    x + y < 0,	-SinIntegral[-z],
	    y == 0,	Im[ExpIntegralE[1,I x]] + N[Pi/2,Precision[z]] + I y/x,
	    x >= y,	-I/2 (ExpIntegralE[1, I z] - ExpIntegralE[1, -I z]) +
				N[Pi/2, Precision[z]],
	    True,	I/2 (ExpIntegralEi[-I z] + ExpIntegralE[1,-I z])
	]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

CosIntegral[z_]:=	(* evaluated for inexact z. *)
    Module[{x, y},
	x = Re[z];
	y = Im[z];
	Which[
	    z == 0,	-Infinity,
	    x + y < 0,	CosIntegral[-z] + If[y < 0, -I, I] N[Pi, Precision[z]],
	    y == 0,	-Re[ExpIntegralE[1, I x]] + I y/x,
	    x >= y,	-(ExpIntegralE[1, I z] + ExpIntegralE[1, -I z])/2,
	    True,	(ExpIntegralEi[-I z] - ExpIntegralE[1, -I z])/2 +
				I N[Pi/2, Precision[z]]
	]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

SinhIntegral[z_]:=	(* evaluated for inexact z. *)
    Module[{shi = -I SinIntegral[I z]},
	If[Head[z] === Real, Re[shi], shi]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

CoshIntegral[z_]:=	(* evaluated for inexact z. *)
    Module[{chi},
	If[z == 0, Return[-Infinity]];
	chi = CosIntegral[I z] + If[Re[z] < 0 && Im[z] >= 0, 3I, -I] *
		N[Pi/2, Precision[z]];
	If[Head[z] === Real && z > 0, Re[chi], chi]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

FresnelS[z_] :=
    Module[{ipiz2 = N[Pi, Precision[z]] I z^2/2},
	If[Head[z] === Real,
	    z Im[Hypergeometric1F1[1/2, 3/2, ipiz2]],
	    -I z (Hypergeometric1F1[1/2, 3/2, ipiz2] -
		Hypergeometric1F1[1/2, 3/2, -ipiz2])/2
	]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

FresnelC[z_] :=
    Module[{ipiz2 = N[Pi, Precision[z]] I z^2/2},
	If[Head[z] === Real,
	    z Re[Hypergeometric1F1[1/2, 3/2, ipiz2]],
	    z (Hypergeometric1F1[1/2, 3/2, ipiz2] +
		Hypergeometric1F1[1/2, 3/2, -ipiz2])/2
	]
    ] /; (NumberQ[z] && Precision[z] < Infinity)

SinIntegral[0]:=0
CosIntegral[0]:=-Infinity
SinhIntegral[0]:=0
CoshIntegral[0]:=-Infinity
FresnelS[0]:=0
FresnelC[0]:=0

SinIntegral[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	Pi/2,
	    I,	DirectedInfinity[I],
	    -1,	-Pi/2,
	    -I,	DirectedInfinity[-I]
	] /; MemberQ[{1, I, -1, -I}, d]

CosIntegral[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	0,
	    I,	DirectedInfinity[1] + I Pi/2,
	    -1,	-I Pi,
	    -I,	DirectedInfinity[1] - I Pi/2
	] /; MemberQ[{1, I, -1, -I}, d]

SinhIntegral[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	DirectedInfinity[1],
	    I,	I Pi/2,
	    -1,	DirectedInfinity[-1],
	    -I,	-I Pi/2
	] /; MemberQ[{1, I, -1, -I}, d]

CoshIntegral[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	DirectedInfinity[1],
	    I,	I Pi/2,
	    -1,	DirectedInfinity[1] + I Pi,
	    -I,	-I Pi/2
	] /; MemberQ[{1, I, -1, -I}, d]

FresnelS[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	1/2,
	    I,	-I/2,
	    -1,	-1/2,
	    -I,	I/2
	] /; MemberQ[{1, I, -1, -I}, d]

FresnelC[DirectedInfinity[d_]] :=
	Switch[d,
	    1, 	1/2,
	    I,	I/2,
	    -1,	-1/2,
	    -I,	-I/2
	] /; MemberQ[{1, I, -1, -I}, d]

SinIntegral/: SinIntegral' := Sin[#]/# &
CosIntegral/: CosIntegral' := Cos[#]/# &
SinhIntegral/: SinhIntegral' := Sinh[#]/# &
CoshIntegral/: CoshIntegral' := Cosh[#]/# &
FresnelS/: FresnelS' := Sin[Pi/2 #^2]&
FresnelC/: FresnelC' := Cos[Pi/2 #^2]&

End[]

Attributes[SinIntegral] = {Listable, ReadProtected};
Attributes[CosIntegral] = {Listable, ReadProtected};
Attributes[SinhIntegral] = {Listable, ReadProtected};
Attributes[CoshIntegral] = {Listable, ReadProtected};
Attributes[FresnelS] = {Listable, ReadProtected};
Attributes[FresnelC] = {Listable, ReadProtected};

Protect[SinIntegral, CosIntegral, SinhIntegral, CoshIntegral,
	FresnelS, FresnelC]

End[]

(* :Tests:

d1[z_] := Chop[1 - SinIntegral[z] / NIntegrate[Sin[t]/t, {t, 0, z},
		PrecisionGoal -> 13, Method -> DoubleExponential]];
d2[z_] := Chop[1 - SinhIntegral[z] / NIntegrate[Sinh[t]/t, {t, 0, z},
		PrecisionGoal -> 13, Method -> DoubleExponential]];
d3[z_] := Chop[1 - CosIntegral[z] / (N[EulerGamma] + Log[z] +
		NIntegrate[(Cos[t] - 1)/t, {t, 0, z},
			PrecisionGoal -> 13, Method -> DoubleExponential])];
d4[z_] := Chop[1 - CoshIntegral[z] / (N[EulerGamma] + Log[z] +
		NIntegrate[(Cosh[t] - 1)/t, {t, 0, z},
			PrecisionGoal -> 13, Method -> DoubleExponential])];
d5[z_] := Chop[1 - FresnelS[z] / NIntegrate[Sin[Pi/2 t^2], {t, 0, z},
		PrecisionGoal -> 13, Method -> DoubleExponential]];
d6[z_] := Chop[1 - FresnelC[z] / NIntegrate[Cos[Pi/2 t^2], {t, 0, z},
		PrecisionGoal -> 13, Method -> DoubleExponential]];

Do[zz = Random[Complex, {-20-20I, 20+20I}];
	Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];
Do[zz = Random[Real, {-20, 20}];
	Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];
Do[zz = Random[Real, {-20, 20}]I;
	Print[{zz, d1[zz], d2[zz], d3[zz], d4[zz], d5[zz], d6[zz]}], {20}];

*)
