1 | 'Classes/ActiveBasic/Math/Math.ab
|
---|
2 |
|
---|
3 | Namespace ActiveBasic
|
---|
4 | Namespace Math
|
---|
5 | '----
|
---|
6 | '浮動小数点数補助
|
---|
7 | Function 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)
|
---|
16 | End Function
|
---|
17 |
|
---|
18 | Function 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
|
---|
29 | End Function
|
---|
30 |
|
---|
31 | Function 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
|
---|
43 | End Function
|
---|
44 |
|
---|
45 | '----
|
---|
46 | '冪乗
|
---|
47 | Function 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
|
---|
63 | End Function
|
---|
64 |
|
---|
65 | Function 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
|
---|
77 | End Function
|
---|
78 |
|
---|
79 | Function 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
|
---|
103 | End Function
|
---|
104 |
|
---|
105 | 'xの符号だけをyのものにした値を返す。
|
---|
106 | '引数 x 元となる絶対値
|
---|
107 | '引数 y 元となる符号
|
---|
108 | Function CopySign(x As Double, y As Double) As Double
|
---|
109 | SetQWord(VarPtr(CopySign), (GetQWord(VarPtr(x)) And &h7fffffffffffffff) Or (GetQWord(VarPtr(y)) And &h8000000000000000))
|
---|
110 | End Function
|
---|
111 |
|
---|
112 | Function CopySign(x As Single, y As Single) As Single
|
---|
113 | SetDWord(VarPtr(CopySign), (GetDWord(VarPtr(x)) And &h7fffffff) Or (GetDWord(VarPtr(y)) And &h80000000))
|
---|
114 | End Function
|
---|
115 |
|
---|
116 | '----
|
---|
117 | '絶対値
|
---|
118 | Function Abs(n As Double) As Double
|
---|
119 | If n < 0 Then
|
---|
120 | Abs = -n
|
---|
121 | Else
|
---|
122 | Abs = n
|
---|
123 | End If
|
---|
124 | End Function
|
---|
125 |
|
---|
126 | Function Abs(n As Single) As Single
|
---|
127 | If n < 0 Then
|
---|
128 | Abs = -n
|
---|
129 | Else
|
---|
130 | Abs = n
|
---|
131 | End If
|
---|
132 | End Function
|
---|
133 |
|
---|
134 | Function Abs(n As Int64) As Int64
|
---|
135 | If n < 0 Then
|
---|
136 | Abs = -n
|
---|
137 | Else
|
---|
138 | Abs = n
|
---|
139 | End If
|
---|
140 | End Function
|
---|
141 |
|
---|
142 | Function Abs(n As Long) As Long
|
---|
143 | If n < 0 Then
|
---|
144 | Abs = -n
|
---|
145 | Else
|
---|
146 | Abs = n
|
---|
147 | End If
|
---|
148 | End Function
|
---|
149 |
|
---|
150 | Function Abs(n As Integer) As Integer
|
---|
151 | If n < 0 Then
|
---|
152 | Abs = -n
|
---|
153 | Else
|
---|
154 | Abs = n
|
---|
155 | End If
|
---|
156 | End Function
|
---|
157 |
|
---|
158 | Function Abs(n As SByte) As SByte
|
---|
159 | If n < 0 Then
|
---|
160 | Abs = -n
|
---|
161 | Else
|
---|
162 | Abs = n
|
---|
163 | End If
|
---|
164 | End Function
|
---|
165 |
|
---|
166 | '----
|
---|
167 | '指数・対数
|
---|
168 |
|
---|
169 | Function 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)
|
---|
198 | End Function
|
---|
199 |
|
---|
200 | Function 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
|
---|
210 | End Function
|
---|
211 |
|
---|
212 | Function 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
|
---|
228 | End Function
|
---|
229 |
|
---|
230 | Function Log10(x As Double) As Double
|
---|
231 | Return Log(x) * Detail._System_InverseLn10
|
---|
232 | End Function
|
---|
233 |
|
---|
234 | '----
|
---|
235 | '三角関数
|
---|
236 | Function 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
|
---|
254 | End Function
|
---|
255 |
|
---|
256 | Function 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)
|
---|
264 | End Function
|
---|
265 |
|
---|
266 | Function 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
|
---|
285 | End Function
|
---|
286 |
|
---|
287 | '--
|
---|
288 | '三角関数の逆関数
|
---|
289 | Function 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
|
---|
295 | End Function
|
---|
296 |
|
---|
297 | Function 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
|
---|
303 | End Function
|
---|
304 |
|
---|
305 | Function 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
|
---|
340 | End Function
|
---|
341 |
|
---|
342 | Function 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
|
---|
351 | End Function
|
---|
352 |
|
---|
353 | '----
|
---|
354 | '双曲線関数
|
---|
355 | Function 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
|
---|
363 | End Function
|
---|
364 |
|
---|
365 | Function 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
|
---|
373 | End Function
|
---|
374 |
|
---|
375 |
|
---|
376 | '----
|
---|
377 | '浮動小数点数判定
|
---|
378 | Function 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
|
---|
386 | End Function
|
---|
387 |
|
---|
388 | Function 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)
|
---|
393 | End Function
|
---|
394 |
|
---|
395 | Function 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
|
---|
399 | End Function
|
---|
400 |
|
---|
401 | '----
|
---|
402 | 'その他
|
---|
403 | Function 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
|
---|
419 | End Function
|
---|
420 |
|
---|
421 | Namespace Detail
|
---|
422 |
|
---|
423 | Function 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
|
---|
430 | End Function
|
---|
431 |
|
---|
432 | #ifdef _WIN64
|
---|
433 |
|
---|
434 | Function GetNaN() As Double
|
---|
435 | SetQWord(VarPtr(GetNaN) As *QWord, &H7FF8000000000000)
|
---|
436 | End Function
|
---|
437 |
|
---|
438 | Function 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)
|
---|
442 | End Function
|
---|
443 |
|
---|
444 | #else
|
---|
445 |
|
---|
446 | Function GetNaN() As Double
|
---|
447 | Dim p = VarPtr(GetNaN) As *DWord
|
---|
448 | p[0] = 0
|
---|
449 | p[1] = &H7FF80000
|
---|
450 | End Function
|
---|
451 |
|
---|
452 | Function 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
|
---|
458 | End Function
|
---|
459 |
|
---|
460 | #endif
|
---|
461 |
|
---|
462 | Function 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)
|
---|
471 | End Function
|
---|
472 |
|
---|
473 | Function _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)
|
---|
491 | End Function
|
---|
492 |
|
---|
493 | Const _System_D = 4.4544551033807686783083602485579e-6 As Double
|
---|
494 | Const _System_UrTan_N = 19 As Long
|
---|
495 | Const _System_EPS5 = 0.001 As Double
|
---|
496 | Const _System_Atan_N = 20 As Long
|
---|
497 | Const _System_HalfPI = (_System_PI * 0.5)
|
---|
498 | Const _System_InverseHalfPI = (2 / _System_PI) '1 / (PI / 2)
|
---|
499 | Const _System_InverseLn10 = 0.43429448190325182765112891891661 '1 / (ln 10)
|
---|
500 | Const _System_InverseSqrt2 = 0.70710678118654752440084436210485 '1 / (√2)
|
---|
501 | Const _System_LOG2 = 0.6931471805599453094172321214581765680755
|
---|
502 | Const _System_Log_N = 7 As Long
|
---|
503 |
|
---|
504 | End Namespace 'Detail
|
---|
505 | End Namespace 'Math
|
---|
506 | End Namespace 'ActiveBasic
|
---|