source: trunk/ab5.0/ablib/src/Classes/ActiveBasic/Math/Math.ab@ 589

Last change on this file since 589 was 589, checked in by イグトランス (egtra), 16 years ago

数学関数をActiveBasic.Mathへ統合

File size: 10.3 KB
Line 
1'Classes/ActiveBasic/Math/Math.ab
2
3Namespace ActiveBasic
4Namespace Math
5'----
6'浮動小数点数補助
7Function ldexp(x As Double, n As Long) As Double
8 If x = 0 Then
9 ldexp = 0
10 Exit Function
11 End If
12 Dim pSrc = VarPtr(x) As *QWord
13 Dim pDest = VarPtr(ldexp) As *QWord
14 n += (pSrc[0] >> 52) As DWord And &h7FF
15 pDest[0] = n << 52 Or (pSrc[0] And &h800FFFFFFFFFFFFF)
16End Function
17
18Function frexp(x As Double, ByRef n As Long) As Double
19 If x = 0 Then
20 n = 0
21 frexp = 0
22 Exit Function
23 End If
24
25 Dim pSrc = VarPtr(x) As *QWord
26 Dim pDest = VarPtr(frexp) As *QWord
27 n = ((pSrc[0] >> 52) As DWord And &h7FF) - 1022
28 pDest[0] = (pSrc[0] And &h800FFFFFFFFFFFFF) Or &h3FE0000000000000
29End Function
30
31Function frexp(x As Single, ByRef n As Long) As Single
32 If x = 0 Then
33 n = 0
34 frexp = 0
35 Exit Function
36 End If
37
38 Dim pSrc As *DWord, pDest As *DWord
39 pSrc = VarPtr(x) As *DWord
40 pDest = VarPtr(frexp) As *DWord
41 n = ((pSrc[0] >> 23) And &hFF) - 126
42 pDest[0] = (pSrc[0] And &h807FFFFF) Or &h7E000000
43End Function
44
45'----
46'冪乗
47Function pow(x As Double, n As Long) As Double
48 Dim abs_n As Long
49 Dim r = 1 As Double
50
51 abs_n=Abs(n) As Long
52 While abs_n<>0
53 If abs_n and 1 Then r *= x
54 x = x * x
55 abs_n >>= 1 ' abs_n \= 2
56 Wend
57
58 If n>=0 Then
59 pow=r
60 Else
61 pow=1/r
62 End If
63End Function
64
65Function pow(x As Double, y As Double) As Double
66 Dim yi = y As Long
67 If y = yi Then
68 pow = pow(x, yi)
69 ElseIf x>0 Then
70 pow = Exp(y * Log(x))
71 Exit Function
72 ElseIf x<>0 or y<=0 Then
73 pow = Detail.GetNaN()
74 Else
75 pow = 0
76 End If
77End Function
78
79Function Sqrt(x As Double) As Double
80 If x > 0 Then
81 If IsInf(x) Then
82 Sqrt = x
83 Else
84 Sqrt = x
85 Dim i = (VarPtr(Sqrt) + 6) As *Word
86 Dim jj = GetWord(i) As Long
87 Dim j = jj >> 5 As Long
88 Dim k = (jj And &h0000001f) As Long
89 j = (j + 511) << 4 + k
90 SetWord(i, j)
91 Dim last As Double
92 Do
93 last = Sqrt
94 Sqrt = (x / Sqrt + Sqrt) * 0.5
95 Loop While Sqrt <> last
96 End If
97 ElseIf x < 0 Then
98 Sqrt = Detail.GetNaN()
99 Else
100 'x = 0 Or NaN
101 Sqrt = x
102 End If
103End Function
104
105'xの符号だけをyのものにした値を返す。
106'引数 x 元となる絶対値
107'引数 y 元となる符号
108Function CopySign(x As Double, y As Double) As Double
109 SetQWord(VarPtr(CopySign), (GetQWord(VarPtr(x)) And &h7fffffffffffffff) Or (GetQWord(VarPtr(y)) And &h8000000000000000))
110End Function
111
112Function CopySign(x As Single, y As Single) As Single
113 SetDWord(VarPtr(CopySign), (GetDWord(VarPtr(x)) And &h7fffffff) Or (GetDWord(VarPtr(y)) And &h80000000))
114End Function
115
116'----
117'絶対値
118Function Abs(n As Double) As Double
119 If n < 0 Then
120 Abs = -n
121 Else
122 Abs = n
123 End If
124End Function
125
126Function Abs(n As Single) As Single
127 If n < 0 Then
128 Abs = -n
129 Else
130 Abs = n
131 End If
132End Function
133
134Function Abs(n As Int64) As Int64
135 If n < 0 Then
136 Abs = -n
137 Else
138 Abs = n
139 End If
140End Function
141
142Function Abs(n As Long) As Long
143 If n < 0 Then
144 Abs = -n
145 Else
146 Abs = n
147 End If
148End Function
149
150Function Abs(n As Integer) As Integer
151 If n < 0 Then
152 Abs = -n
153 Else
154 Abs = n
155 End If
156End Function
157
158Function Abs(n As SByte) As SByte
159 If n < 0 Then
160 Abs = -n
161 Else
162 Abs = n
163 End If
164End Function
165
166'----
167'指数・対数
168
169Function Exp(x As Double) As Double
170 If IsNaN(x) Then
171 Return x
172 Else If IsInf(x) Then
173 If 0 > x Then
174 Return 0
175 Else
176 Return x
177 End If
178 End If
179 Dim k As Long
180 If x >= 0 Then
181 k = Fix(x / Detail._System_LOG2 + 0.5)
182 Else
183 k = Fix(x / Detail._System_LOG2 - 0.5)
184 End If
185
186 x -= k * Detail._System_LOG2
187
188 Dim x2 = x * x
189 Dim w = x2 / 22
190
191 Dim i = 18
192 While i >= 6
193 w = x2 / (w + i)
194 i -= 4
195 Wend
196
197 Return ldexp((2 + w + x) / (2 + w - x), k)
198End Function
199
200Function Log1p(x As Double) As Double
201 If x < -1 Or IsNaN(x) Then
202 Log1p = Detail.GetNaN()
203 ElseIf x = 0 Then
204 x = 0
205 ElseIf IsInf(x) Then
206 Log1p = x
207 Else
208 Log1p = Detail.Log1p(x)
209 End If
210End Function
211
212Function Log(x As Double) As Double
213 If x = 0 Then
214 Log = Detail.GetInf(True)
215 ElseIf x < 0 Or IsNaN(x) Then
216 Log = Detail.GetNaN()
217 ElseIf IsInf(x) Then
218 Log = x
219 Else
220 Dim tmp = x * Detail._System_InverseSqrt2
221 Dim p = VarPtr(tmp) As *QWord
222 Dim m = GetQWord(p) And &h7FF0000000000000
223 Dim k = ((m >> 52) As DWord) As Long - 1022
224 SetQWord(p, m + &h0010000000000000)
225 x /= tmp
226 Log = Detail._System_LOG2 * k + Detail.Log1p(x - 1)
227 End If
228End Function
229
230Function Log10(x As Double) As Double
231 Return Log(x) * Detail._System_InverseLn10
232End Function
233
234'----
235'三角関数
236Function Sin(x As Double) As Double
237 If IsNaN(x) Then
238 Return x
239 ElseIf IsInf(x) Then
240 Return Detail.GetNaN()
241 Exit Function
242 End If
243
244 Dim k As Long
245 Dim t As Double
246
247 t = Detail._Support_tan((x * 0.5) As Double, k)
248 t = 2 * t / (1 + t * t)
249 If (k And 1) = 0 Then 'k mod 2 = 0 Then
250 Return t
251 Else
252 Return -t
253 End If
254End Function
255
256Function Cos(x As Double) As Double
257 If IsNaN(x) Then
258 Return x
259 ElseIf IsInf(x) Then
260 Return Detail.GetNaN()
261 End If
262
263 Return Sin((Detail._System_HalfPI - Abs(x)) As Double)
264End Function
265
266Function Tan(x As Double) As Double
267 If IsNaN(x) Then
268 Tan = x
269 Exit Function
270 ElseIf IsInf(x) Then
271 Tan = Detail.GetNaN()
272 Exit Function
273 End If
274
275 Dim k As Long
276 Dim t As Double
277 t = Detail._Support_tan(x, k)
278 If (k And 1) = 0 Then 'k mod 2 = 0 Then
279 Return t
280 ElseIf t <> 0 Then
281 Return -1 / t
282 Else
283 Return CopySign(Detail.GetInf(False), -t)
284 End If
285End Function
286
287'--
288'三角関数の逆関数
289Function Asin(x As Double) As Double
290 If x < -1 Or x > 1 Then
291 Asin = Detail.GetNaN()
292 Else
293 Asin = Atan(x / Sqrt(1 - x * x))
294 End If
295End Function
296
297Function Acos(x As Double) As Double
298 If x < -1 Or x > 1 Then
299 Acos = Detail.GetNaN()
300 Else
301 Acos = Detail._System_HalfPI - Asin(x)
302 End If
303End Function
304
305Function Atan(x As Double) As Double
306 If IsNaN(x) Then
307 Atan = x
308 Exit Function
309 ElseIf IsInf(x) Then
310 Atan = CopySign(_System_PI, x)
311 Exit Function
312 End If
313 Dim i As Long
314 Dim sgn As Long
315 Dim dbl = 0 As Double
316
317 If x > 1 Then
318 sgn = 1
319 x = 1 / x
320 ElseIf x < -1 Then
321 sgn = -1
322 x = 1 / x
323 Else
324 sgn = 0
325 End If
326
327 For i = Detail._System_Atan_N To 1 Step -1
328 Dim t As Double
329 t = i * x
330 dbl = (t * t) / (2 * i + 1 + dbl)
331 Next
332
333 If sgn > 0 Then
334 Atan = Detail._System_HalfPI - x / (1 + dbl)
335 ElseIf sgn < 0 Then
336 Atan = -Detail._System_HalfPI - x / (1 + dbl)
337 Else
338 Atan = x / (1 + dbl)
339 End If
340End Function
341
342Function Atan2(y As Double, x As Double) As Double
343 If x = 0 Then
344 Atan2 = Sgn(y) * Detail._System_HalfPI
345 Else
346 Atan2 = Atn(y / x)
347 If x < 0 Then
348 Atan2 += CopySign(_System_PI, y)
349 End If
350 End If
351End Function
352
353'----
354'双曲線関数
355Function Sinh(x As Double) As Double
356 If Abs(x) > Detail._System_EPS5 Then
357 Dim t As Double
358 t = Exp(x)
359 Return (t - 1 / t) * 0.5
360 Else
361 Return x * (1 + x * x * 0.166666666666667) ' x * (1 + x * x / 6)
362 End If
363End Function
364
365Function Tanh(x As Double) As Double
366 If x > Detail._System_EPS5 Then
367 Return 2 / (1 + Exp(-2 * x)) - 1
368 ElseIf x < -Detail._System_EPS5 Then
369 Return 1 - 2 / (Exp(2 * x) + 1)
370 Else
371 Return x * (1 - x * x * 0.333333333333333) 'x * (1 - x * x / 3)
372 End If
373End Function
374
375
376'----
377'浮動小数点数判定
378Function IsNaN(ByVal x As Double) As Boolean
379 Dim p = VarPtr(x) As *DWord
380 IsNaN = False
381 If (p[1] And &H7FF00000) = &H7FF00000 Then
382 If (p[0] <> 0) Or ((p[1] And &HFFFFF) <> 0) Then
383 IsNaN = True
384 End If
385 End If
386End Function
387
388Function IsInf(x As Double) As Boolean
389 Dim p = VarPtr(x) As *DWord
390 p[1] And= &h7fffffff
391 Dim inf = Detail.GetInf(False)
392 IsInf = (memcmp(p As *Byte, VarPtr(inf), SizeOf (Double)) = 0)
393End Function
394
395Function IsFinite(x As Double) As Boolean
396 Dim p = VarPtr(x) As *DWord
397 p[1] And= &H7FF00000
398 IsFinite = ( p[1] And &H7FF00000 ) = &H7FF00000
399End Function
400
401'----
402'その他
403Function Hypot(x As Double, y As Double) As Double
404 If x = 0 Then
405 Hypot = Abs(y)
406 Else If y = 0 Then
407 Hypot = Abs(x)
408 Else
409 Dim ax = Abs(x)
410 Dim ay = Abs(y)
411 If ay > ax Then
412 Dim t = x / y
413 Hypot = ay * Sqr(1 + t * t)
414 Else
415 Dim t = y / x
416 Hypot = ax * Sqr(1 + t * t)
417 End If
418 End If
419End Function
420
421Namespace Detail
422
423Function SetSign(x As Double, isNegative As Long) As Double
424#ifdef _WIN64
425 SetQWord(AddressOf(CopySign), (GetQWord(VarPtr(x)) And &h7fffffffffffffff) Or (isNegative << 63))
426#else
427 SetDWord(AddressOf(CopySign), GetDWord(VarPtr(x)))
428 SetDWord(AddressOf(CopySign) + SizeOf (DWord), GetQWord(VarPtr(x) + SizeOf (DWord)) And &h7fffffff Or (isNegative << 31))
429#endif
430End Function
431
432#ifdef _WIN64
433
434Function GetNaN() As Double
435 SetQWord(VarPtr(GetNaN) As *QWord, &H7FF8000000000000)
436End Function
437
438Function GetInf(sign As Boolean) As Double
439 Dim s = 0 As QWord
440 If sign Then s = 1 << 63
441 SetQWord(VarPtr(GetInf) As *QWord, &h7FF0000000000000 Or s)
442End Function
443
444#else
445
446Function GetNaN() As Double
447 Dim p = VarPtr(GetNaN) As *DWord
448 p[0] = 0
449 p[1] = &H7FF80000
450End Function
451
452Function GetInf(sign As Boolean) As Double
453 Dim s = 0 As DWord
454 If sign Then s = (1 As DWord) << 31
455 Dim p = VarPtr(GetInf) As *DWord
456 p[0] = 0
457 p[1] = &h7FF00000 Or s
458End Function
459
460#endif
461
462Function Log1p(x As Double) As Double
463 Dim s = 0 As Double
464 Dim i = _System_Log_N As Long
465 While i >= 1
466 Dim t = (i * x) As Double
467 s = t * (2 * i + 1 + s) / (2 + t)
468 i--
469 Wend
470 Return x / (1 + s)
471End Function
472
473Function _Support_tan(x As Double, ByRef k As Long) As Double
474 If x >= 0 Then
475 k = ( Fix(x * _System_InverseHalfPI) + 0.5 ) As Long
476 Else
477 k = ( Fix(x * _System_InverseHalfPI) - 0.5 ) As Long
478 End If
479
480 x = (x - (3217.0 / 2048.0) * k) + _System_D * k
481
482 Dim x2 = x * x
483 Dim t = 0 As Double
484
485 Dim i As Long
486 For i = _System_UrTan_N To 3 Step -2
487 t = x2 / (i - t)
488 Next
489
490 _Support_tan = x / (1 - t)
491End Function
492
493Const _System_D = 4.4544551033807686783083602485579e-6 As Double
494Const _System_UrTan_N = 19 As Long
495Const _System_EPS5 = 0.001 As Double
496Const _System_Atan_N = 20 As Long
497Const _System_HalfPI = (_System_PI * 0.5)
498Const _System_InverseHalfPI = (2 / _System_PI) '1 / (PI / 2)
499Const _System_InverseLn10 = 0.43429448190325182765112891891661 '1 / (ln 10)
500Const _System_InverseSqrt2 = 0.70710678118654752440084436210485 '1 / (√2)
501Const _System_LOG2 = 0.6931471805599453094172321214581765680755
502Const _System_Log_N = 7 As Long
503
504End Namespace 'Detail
505End Namespace 'Math
506End Namespace 'ActiveBasic
Note: See TracBrowser for help on using the repository browser.