「余弦積分」の編集履歴（バックアップ）一覧はこちら

「余弦積分」(2010/07/16 (金) 15:18:04) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

#asciiart(blockquote){
/***********************************************************
	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
}