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