/***********************************************************
ci.c -- 余弦積分
************************************************************
${\rm Ci}(x) = \gamma + \log x
+ \int_{0}^{x} \frac{\cos t - 1}{t} \, dt$
x Ci(x)
1 0.33740392290096813466
2 0.42298082877486499570
3 0.11962978600800032763
4 -0.14098169788693041164
5 -0.19002974965664387862
10 -0.045456433004455372635
20 0.044419820845353316540
30 -0.033032417282071143779
40 0.019020007896208766962
50 -0.0056283863241163054402
***********************************************************/
Const EULER = 0.577215664901532860606512090082 /* Eulerの定数 $\gamma$ */
Function Ci_series( x As Double) As Double /* 級数展開 */
Dim k As Long
Dim s, t, u
s = EULER + Log(x)
x = - x * x: t = 1
k=2
While k < 1000
t = t * (x / ((k - 1) * k))
u = s: s = s + (t / k)
If (s = u) Then
Ci_series = s
Exit Function
End If
k = k + 2
Wend
Print "Si_series(): 収束しません."
Ci_series = s
End Function
Function Ci_asympt(x As Double) As Double /* 漸近展開 */
Dim k As Long, flag As Long
Dim t, f, g, fmax, fmin, gmax, gmin
fmax = 2: gmax = 2: fmin = 0:gmin = 0
f = 0:g = 0: t = 1 / x
k = 0:flag = 0
while (flag <> 15)
f = f+t: k=k+1:t = t * (k / x)
if (f < fmax) Then fmax = f else flag = flag or 1
g= g+ t: k=k+1:t = t * (k / x)
if (g < gmax) Then gmax = g else flag = flag or 2
f = f-t: k=k+1:t = t * (k / x)
if (f > fmin) Then fmin = f else flag = flag or 4
g = g-t: k=k+1:t = t * (k / x)
if (g > gmin) Then gmin = g else flag = flag or 8
Wend
Ci_asympt = 0.5 * ((fmax + fmin) * Sin(x) - (gmax + gmin) * Cos(x))
End Function
Function Ci(x As Double) As Double
if (x < 0) Then
Ci = -Ci(-x)
Exit Function
End If
if (x < 18) Then
Ci = Ci_series(x)
Exit Function
End If
Ci = Ci_asympt(x)
End Function
#N88BASIC
Dim x As Double
Print (" x Ci(x)")
For x = 1 To 50
Print x, Ci(x)
Next
