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

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

Circle, Sqrtを壊していたので修正。プロンプト画面の描画ステートメントの座標の引数をDoubleへ変更。

File size: 10.3 KB
Line 
1'Classes/ActiveBasic/Math/Math.ab
2
3Namespace ActiveBasic
4Namespace Math
5
6Const PI = _System_PI
7
8'----
9'浮動小数点数補助
10Function 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)
19End Function
20
21Function 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
32End Function
33
34Function 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
46End Function
47
48'----
49'冪乗
50Function 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
66End Function
67
68Function 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
80End Function
81
82Function 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
106End Function
107
108'xの符号だけをyのものにした値を返す。
109'引数 x 元となる絶対値
110'引数 y 元となる符号
111Function CopySign(x As Double, y As Double) As Double
112 SetQWord(VarPtr(CopySign), (GetQWord(VarPtr(x)) And &h7fffffffffffffff) Or (GetQWord(VarPtr(y)) And &h8000000000000000))
113End Function
114
115Function CopySign(x As Single, y As Single) As Single
116 SetDWord(VarPtr(CopySign), (GetDWord(VarPtr(x)) And &h7fffffff) Or (GetDWord(VarPtr(y)) And &h80000000))
117End Function
118
119'----
120'絶対値
121Function Abs(n As Double) As Double
122 If n < 0 Then
123 Abs = -n
124 Else
125 Abs = n
126 End If
127End Function
128
129Function Abs(n As Single) As Single
130 If n < 0 Then
131 Abs = -n
132 Else
133 Abs = n
134 End If
135End Function
136
137Function Abs(n As Int64) As Int64
138 If n < 0 Then
139 Abs = -n
140 Else
141 Abs = n
142 End If
143End Function
144
145Function Abs(n As Long) As Long
146 If n < 0 Then
147 Abs = -n
148 Else
149 Abs = n
150 End If
151End Function
152
153Function Abs(n As Integer) As Integer
154 If n < 0 Then
155 Abs = -n
156 Else
157 Abs = n
158 End If
159End Function
160
161Function Abs(n As SByte) As SByte
162 If n < 0 Then
163 Abs = -n
164 Else
165 Abs = n
166 End If
167End Function
168
169'----
170'指数・対数
171
172Function 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)
201End Function
202
203Function 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
213End Function
214
215Function 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
231End Function
232
233Function Log10(x As Double) As Double
234 Return Log(x) * Detail._System_InverseLn10
235End Function
236
237'----
238'三角関数
239Function 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
257End Function
258
259Function 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)
267End Function
268
269Function 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
288End Function
289
290'--
291'三角関数の逆関数
292Function 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
298End Function
299
300Function 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
306End Function
307
308Function 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
343End Function
344
345Function 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
354End Function
355
356'----
357'双曲線関数
358Function 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
366End Function
367
368Function 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
376End Function
377
378
379'----
380'浮動小数点数判定
381Function 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
389End Function
390
391Function 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)
396End Function
397
398Function 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
402End Function
403
404'----
405'その他
406Function 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
422End Function
423
424Namespace Detail
425
426Function 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
433End Function
434
435#ifdef _WIN64
436
437Function GetNaN() As Double
438 SetQWord(VarPtr(GetNaN) As *QWord, &H7FF8000000000000)
439End Function
440
441Function 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)
445End Function
446
447#else
448
449Function GetNaN() As Double
450 Dim p = VarPtr(GetNaN) As *DWord
451 p[0] = 0
452 p[1] = &H7FF80000
453End Function
454
455Function 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
461End Function
462
463#endif
464
465Function 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)
474End Function
475
476Function _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)
494End Function
495
496Const _System_D = 4.4544551033807686783083602485579e-6 As Double
497Const _System_UrTan_N = 19 As Long
498Const _System_EPS5 = 0.001 As Double
499Const _System_Atan_N = 20 As Long
500Const _System_HalfPI = (_System_PI * 0.5)
501Const _System_InverseHalfPI = (2 / _System_PI) '1 / (PI / 2)
502Const _System_InverseLn10 = 0.43429448190325182765112891891661 '1 / (ln 10)
503Const _System_InverseSqrt2 = 0.70710678118654752440084436210485 '1 / (√2)
504Const _System_LOG2 = 0.6931471805599453094172321214581765680755
505Const _System_Log_N = 7 As Long
506
507End Namespace 'Detail
508End Namespace 'Math
509End Namespace 'ActiveBasic
Note: See TracBrowser for help on using the repository browser.