OSDN -- Develop and Download Open Source Software
  • R/O
  • HTTP
  • SSH
  • HTTPS
Fork

karinto: Commit

Karinto and KarintoTest


Commit MetaInfo

Revisiond01d6f45db3bcdebd99b5edf03422bb0537b6310 (tree)
Time2013-06-11 03:09:39
Authorkimikage <kimikage_ceo@hotm...>
Commiterkimikage

Log Message

MathExクラスの追加，誤差関数の実装(#23869)

Change Summary

Incremental Difference

--- a/Karinto/Karinto.csproj
+++ b/Karinto/Karinto.csproj
@@ -47,6 +47,7 @@
4747 <Compile Include="Fft.cs" />
4848 <Compile Include="HiokiHicorderData.cs" />
4949 <Compile Include="HiokiHicorderDataReader.cs" />
50+ <Compile Include="MathEx.cs" />
5051 <Compile Include="Numeric\Complex.cs" />
5152 <Compile Include="Numeric\Matrix.cs" />
5253 <Compile Include="Numeric\Matrix2x2.cs" />
--- /dev/null
+++ b/Karinto/MathEx.cs
@@ -0,0 +1,231 @@
1+/*
2+ * Karinto Library Project
3+ *
4+ * This software is distributed under a zlib-style license.
5+ * See license.txt for more information.
6+ */
7+
8+using System;
9+using System.Collections.Generic;
10+using System.Text;
11+
12+namespace Karinto
13+{
14+ public static class MathEx
15+ {
16+ #region fields
17+ private static readonly Polynomial[] erfcTable =
18+ {
19+ new Polynomial( //[0]: 1 <= x < 2
20+ 4.7019003296642450457550391704e-1,
21+ -5.4129916560282193120142282350e-4,
22+ 3.2466045488989207057247716494e-2,
23+ -1.5551574967900995569013129661e-3,
24+ 1.2991371652846980287617056420e-3,
25+ -1.2476738559919909068871115854e-4,
26+ 4.2164948133638407169423314936e-5,
27+ -5.5944555649744218134732620248e-6,
28+ 1.2391258283306184550290288637e-6,
29+ -1.8629332352031266041799875695e-7,
30+ 3.3508489330180889238527074083e-8,
31+ -5.1453759994172323811411458769e-9,
32+ 8.2778054007509566713316121080e-10,
33+ -1.2419584612081402943997565209e-10,
34+ 1.8628225410158766603724913483e-11,
35+ -2.6916642541339090011279256698e-12,
36+ 3.8321798025710167662264220379e-13,
37+ -5.2985573947392393836199793226e-14,
38+ 7.2214593891835651145391699616e-15,
39+ -1.0552316968012789176241975783e-15,
40+ 1.3848365749828243613160210837e-16
41+ ),
42+ new Polynomial( //[1]: 2 <= x < 4
43+ 0.28246919208993532,//2.8246919208993525100243945507e-1,
44+ 1.0295429295263581287408466969e-4,
45+ 2.6352174925955979906965931429e-2,
46+ -2.0840262584872572824603082922e-3,
47+ 1.6107317263933041407504226538e-3,
48+ -2.6413257091217055241174636048e-4,
49+ 9.3720096524725604112989062675e-5,
50+ -2.0394452777501772979594808178e-5,
51+ 5.4818843984699154429750547408e-6,
52+ -1.2769242863171188612896799307e-6,
53+ 3.0855108190100308040276886416e-7,
54+ -7.1110449954162265412516889674e-8,
55+ 1.6277854940656666077194365644e-8,
56+ -3.6388350444865150940102551512e-9,
57+ 8.0125451384130244729739875906e-10,
58+ -1.7333139256038089278659882792e-10,
59+ 3.6913251414732819339098560977e-11,
60+ -7.7391112861863152629855204763e-12,
61+ 1.5985771117554783386097643078e-12,
62+ -3.2558812235239520370430316160e-13,
63+ 6.5474397002885039401566856717e-14,
64+ -1.2908539064863451579569117003e-14,
65+ 2.4663418649675516649805424542e-15,
66+ -5.0102905606542145531985288605e-16,
67+ 1.1637629874200535139854357149e-16,
68+ -1.7702372520803693617903759643e-17
69+ ),
70+ new Polynomial( //[2]: 4 <= x < 8
71+ 2.6609916698112647641230096605e-1,
72+ 1.0052997692961233689565319255e-1,
73+ 0.063819586227166714,//6.3819586227166696717161903119e-2,
74+ 1.6510111747337704484097091109e-2,
75+ 6.7745775475928297665975539362e-3,
76+ 1.2613799574673550665257600345e-3,
77+ 4.6314050254419785813179921962e-4,
78+ 5.4816124629655721757719475954e-5,
79+ 2.5069740734847835425277961803e-5,
80+ 7.6494811710229433994592919775e-7,
81+ 1.3069891449561428965827894398e-6,
82+ -1.1211722505095090109921762128e-7,
83+ 7.7627072370605473720667209070e-8,
84+ -1.5067764445724103557487307330e-8,
85+ 5.3921138750029687859330065848e-9,
86+ -1.3591247804716958040748207366e-9,
87+ 4.0341965879492537779635899830e-10,
88+ -1.0884398300756496354533131958e-10,
89+ 3.0363217123340332343882922235e-11,
90+ -8.2439513886349055414081918841e-12,
91+ 2.2756305305902268859393373330e-12,
92+ -6.1246789459726400033402401708e-13,
93+ 1.4333003386843162197015125496e-13,
94+ -3.7815369812690037058307445308e-14,
95+ 1.7401564599128446023009328203e-14,
96+ -4.5875278594765951588791904401e-15
97+ ),
98+ new Polynomial( //[3]: 8 <= x < 16
99+ 4.6854221014893762619745618301e-2,
100+ -1.5511450952249084104775598329e-2,
101+ 5.1178905303441648577722059649e-3,
102+ -1.6829798529769531049332840624e-3,
103+ 5.5160777130644620308102757353e-4,
104+ -1.8020184996880086091090661904e-4,
105+ 5.8678514133519807272309470851e-5,
106+ -1.9045977453678276537004114523e-5,
107+ 6.1623270905278656151589976968e-6,
108+ -1.9875419915464967477768116933e-6,
109+ 6.3904356643246431611609250453e-7,
110+ -2.0483278667408340516487849749e-7,
111+ 6.5453904816333940501972558207e-8,
112+ -2.0852118199579889358624923805e-8,
113+ 6.6229050456699429842283832512e-9,
114+ -2.0972627311983894339808366077e-9,
115+ 6.6237821648144864340143708663e-10,
116+ -2.0853485442305791192334927397e-10,
117+ 6.5162153190849185531320932146e-11,
118+ -2.0380108335416826794297357698e-11,
119+ 6.6463403803777970165860627602e-12,
120+ -2.0777605497301746376639943675e-12,
121+ 4.6614566309518191650253075708e-13,
122+ -1.4074713019480904125485292343e-13,
123+ 1.0831110710204089408511450325e-13,
124+ -3.3906934201035073659517665140e-14
125+ ),
126+ new Polynomial( //[4]: 16 <= x < 28
127+ 2.8262089312268317148835811995e-2,
128+ -5.4152623609484073274229170763e-3,
129+ 1.6241399412435775154107409724e-3,
130+ -5.0744105591774430405597196162e-3,
131+ -1.7616157737971578130497741968e-2,
132+ 2.3083427166129394399962687197e-2,
133+ 1.1164742403248709989628167597e-1,
134+ -3.5273994341997794105095275190e-2,
135+ -3.0869428535741202544195312887e-1,
136+ -1.3020536214165453885942308553e-2,
137+ 4.3332671745834241123676787595e-1,
138+ 9.4385604975718760529994282754e-2,
139+ -3.0238469418692319645741714650e-1,
140+ -9.5995136295717931771565803096e-2,
141+ 8.3313274670729987166299100438e-2,
142+ 3.1593997744278794785645602547e-2
143+ )
144+ };
145+ #endregion
146+
147+ #region public methods
148+ public static double Erf(double x)
149+ {
150+ double a = Math.Abs(x);
151+ if (a >= 1.0)
152+ {
153+ return 1.0 - Erfc(x);
154+ }
155+ double y;
156+ double aa = a * a;
157+ if (a < 0.0625)
158+ {
159+ double c1 = a * aa;
160+ double c2 = c1 * aa;
161+ double c3 = c2 * aa;
162+ double c4 = c3 * aa;
163+
164+ double s1 = c1 * -315.0;// (3*5*7*9 / -3 / 1!);
165+ double s2 = s1 * 0 + c2 * 189.0 / 2.0;// (3*5*7*9 / +5 / 2!)
166+ double s3 = s2 * 0 + c3 * -45.0 / 2.0;// (3*5*7*9 / -7 / 3!)
167+ double s4 = s3 * 0 + c4 * 35.0 / 8.0;// (3*5*7*9 / 9 / 4!)
168+ y = a + (s1 + s2 + s3 + s4) / 945.0; //(3*5*7*9);
169+ }
170+ else
171+ {
172+ double c = a;
173+ double s = 0;
174+ double b = 1.0;
175+ for (double i = -1.0; i > -20.0; i-=1.0)
176+ {
177+ c *= aa / i;
178+ b += 2.0;
179+ s += c / b;
180+ }
181+ y = a + s;
182+ }
183+ y *= 1.1283791670955126; // 2 / pi^0.5
184+ return Math.Sign(x) > 0 ? y : -y;
185+ }
186+
187+ public static double Erfc(double x)
188+ {
189+ double a = Math.Abs(x);
190+ if (a < 1.0)
191+ {
192+ return 1.0 - Erf(x);
193+ }
194+ if (a > 27.226017111108501)
195+ {
196+ return x > 0 ? 0 : 2.0;
197+ }
198+ double y = 0.0;
199+
200+ if (a < 2.0)
201+ {
202+ y = Math.Exp(-1.1688327773406599 * a * a);
203+ y *= erfcTable[0].GetValueAt(a + a - 3.0);
204+ }
205+ else if (a < 4.0)
206+ {
207+ decimal b = -1.05068636145441m * (decimal)a * (decimal)a;
208+ y = Math.Exp((double)b);
209+ y *= erfcTable[1].GetValueAt(a - 3.0);
210+ }
211+ else if (a < 8.0)
212+ {
213+ decimal b = -1.0292687483818601m * (decimal)a * (decimal)a;
214+ y = Math.Exp((double)b);
215+ y *= erfcTable[2].GetValueAt(a * 0.5 - 3);
216+ }
217+ else if (a < 16.0)
218+ {
219+ y = Math.Exp(-a * a);
220+ y *= erfcTable[3].GetValueAt(a * 0.25 - 3);
221+ }
222+ else
223+ {
224+ y = Math.Exp(-a * a);
225+ y *= erfcTable[4].GetValueAt(a * 0.125 - 3);
226+ }
227+ return x > 0 ? y : 2.0 - y;
228+ }
229+ #endregion
230+ }
231+}
--- a/KarintoTest/KarintoTest.csproj
+++ b/KarintoTest/KarintoTest.csproj
@@ -50,6 +50,7 @@
5050 </ItemGroup>
5151 <ItemGroup>
5252 <Compile Include="CsvWriterTest.cs" />
53+ <Compile Include="MathExTest.cs" />
5354 <Compile Include="NumericTest\ComplexTest.cs" />
5455 <Compile Include="NumericTest\Matrix2x2Test.cs" />
5556 <Compile Include="NumericTest\OperatorTest.cs" />
--- /dev/null
+++ b/KarintoTest/MathExTest.cs
@@ -0,0 +1,174 @@
1+/*
2+ * Karinto Library Project
3+ *
4+ * This software is distributed under a zlib-style license.
5+ * See license.txt for more information.
6+ */
7+
8+using System;
9+using System.Collections.Generic;
10+using System.Text;
11+using NUnit.Framework;
12+using Karinto;
13+using System.Diagnostics;
14+
15+namespace KarintoTest
16+{
17+ class MathExTest
18+ {
19+ private string DoubleToRational(double d)
20+ {
21+ int e;
22+ ulong a;
23+ unsafe
24+ {
25+ ulong i64;
26+ *(&i64) = *(ulong*)((void*)&d);
27+ e = (int)(i64 >> 52) & 0x7FF;
28+ a = i64 & 0xFFFFFFFFFFFFFul | 0x10000000000000ul;
29+ }
30+ string s = d < 0 ? "-" : "";
31+ return String.Format("{0}{1} * 2b0^{2}", s, a, e -52 - 1023);
32+ }
33+
34+ [Test]
35+ public void Erf()
36+ {
37+ PointList p = new PointList();
38+ p.Add(-6.0, -1.0);
39+ p.Add(-1.0, -0.84270079294971489);
40+ p.Add(0.0, 0.0);
41+ p.Add(1e-300, 1.1283791670955126e-300);
42+ p.Add(1e-100, 1.1283791670955126e-100);
43+ p.Add(1e-50, 1.1283791670955126e-50);
44+ p.Add(1e-10, 1.1283791670955126e-10);
45+ p.Add(1e-9, 1.1283791670955126e-9);
46+ p.Add(1e-8, 1.1283791670955126e-8);
47+ p.Add(3e-8, 3.3851375012865365e-8);
48+ p.Add(1e-7, 1.1283791670955088e-7);
49+ p.Add(3e-7, 3.3851375012864358e-7);
50+ p.Add(1e-6, 1.1283791670951364e-6);
51+ p.Add(3e-6, 3.3851375012763826e-6);
52+ p.Add(1e-5, 1.1283791670579000e-5);
53+ p.Add(3e-5, 3.3851375002709968e-5);
54+ p.Add(5e-5, 5.6418958307759834e-5);
55+ p.Add(1e-4, 1.1283791633342487e-4);
56+ p.Add(3e-4, 3.3851373997324151e-4);
57+ p.Add(5e-4, 5.6418953653196121e-4);
58+ p.Add(1e-3, 1.1283787909692363e-3);
59+ p.Add(3e-3, 3.3851273459014537e-3);
60+ p.Add(5e-3, 5.6418488200315501e-3);
61+ p.Add(1e-2, 1.1283415555849616e-2);
62+ p.Add(3e-2, 3.3841222341735429e-2);
63+ p.Add(5e-2, 5.6371977797016630e-2);
64+ p.Add(1e-1, 1.1246291601828490e-1);
65+ p.Add(0.2, 0.22270258921047847);
66+ p.Add(0.3, 0.32862675945912739);
67+ p.Add(0.4, 0.42839235504666845);
68+ p.Add(0.5, 0.52049987781304652);
69+ p.Add(0.6, 0.60385609084792591);
70+ p.Add(0.7, 0.67780119383741844);
71+ p.Add(0.8, 0.74210096470766052);
72+ p.Add(0.9, 0.79690821242283216);
73+ p.Add(0.99999999999999989, 0.84270079294971478);
74+ p.Add(1.0, 0.84270079294971489);
75+ p.Add(1.0000000000000002, 0.84270079294971501);
76+ p.Add(2.0, 0.99532226501895271);
77+ p.Add(3.0, 0.99997790950300136);
78+ p.Add(4.0, 0.99999998458274209);
79+ p.Add(5.0, 0.99999999999846256);
80+ p.Add(6.0, 1.0);
81+ foreach (Point t in p)
82+ {
83+ double erf = MathEx.Erf(t.X);
84+ string m = String.Format("x = {0}, y = {1}",
85+ t.X, DoubleToRational(erf));
86+ Assert.AreEqual(t.Y, erf, Math.Abs(t.Y) * 1.7763568394002505e-15, m);
87+ }
88+ }
89+
90+ [Test]
91+ public void Erfc()
92+ {
93+ PointList p = new PointList();
94+ p.Add(-28.0, 2.0);
95+ p.Add(-1.0, 1.8427007929497148);
96+ p.Add(0.0, 1.0);
97+ p.Add(0.1, 0.88753708398171505);
98+ p.Add(0.5, 0.47950012218695348);
99+ p.Add(1.0, 0.15729920705028513);
100+ p.Add(1.0000000000000002, 0.15729920705028505);
101+ p.Add(1.1, 0.11979493042591827);
102+ p.Add(1.2, 0.089686021770364638);
103+ p.Add(1.3, 0.065992055059347549);
104+ p.Add(1.4, 0.047714880237351202);
105+ p.Add(1.5, 0.033894853524689274);
106+ p.Add(1.6, 0.023651616655355985);
107+ p.Add(1.7, 0.016209541409225439);
108+ p.Add(1.8, 0.010909498364269283);
109+ p.Add(1.9, 0.0072095707647425325);
110+ p.Add(1.9999999999999998, 0.0046777349810472706);
111+ p.Add(2.0, 0.0046777349810472662);
112+ p.Add(2.1, 0.0029794666563329841);
113+ p.Add(2.2, 0.0018628462979818898);
114+ p.Add(2.3, 0.0011431765973566525);
115+ p.Add(2.4, 0.0006885138966450789);
116+ p.Add(2.5, 0.00040695201744495892);
117+ p.Add(2.6, 0.00023603441652934908);
118+ p.Add(2.7, 0.00013433273994052419);
119+ p.Add(2.8, 7.5013194665459108e-05);
120+ p.Add(2.9, 4.1097878099458858e-05);
121+ p.Add(3.0, 2.2090496998585441e-05);
122+ p.Add(3.1, 1.1648657367199589e-05);
123+ p.Add(3.2, 6.0257611517620876e-06);
124+ p.Add(3.3, 3.0577097964381650e-06);
125+ p.Add(3.4, 1.5219933628622864e-06);
126+ p.Add(3.5, 7.4309837234141278e-07);
127+ p.Add(4.0, 1.5417257900280020e-08);
128+
129+ foreach (Point t in p)
130+ {
131+ double erfc = MathEx.Erfc(t.X);
132+ string m = String.Format("x = {0}, y = {1}",
133+ t.X, DoubleToRational(erfc));
134+ Assert.AreEqual(t.Y, erfc, Math.Abs(t.Y) * 1.7763568394002505e-15, m);
135+ }
136+
137+ p.Clear();
138+ p.Add(4.5, 1.9661604415428876e-10);
139+ p.Add(5.0, 1.5374597944280349e-12);
140+ p.Add(5.5, 7.3578479179743983e-15);
141+ p.Add(6.0, 2.1519736712498913e-17);
142+ p.Add(6.5, 3.8421483271206475e-20);
143+ p.Add(7.0, 4.1838256077794142e-23);
144+ p.Add(7.5, 2.7766493860305689e-26);
145+ p.Add(8.0, 1.1224297172982926e-29);
146+ p.Add(10.0, 2.0884875837625449e-45);
147+ p.Add(12.0, 1.3562611692059042e-64);
148+ p.Add(14.0, 3.0372298477503115e-87);
149+
150+ foreach (Point t in p)
151+ {
152+ double erfc = MathEx.Erfc(t.X);
153+ string m = String.Format("x = {0}, y = {1}",
154+ t.X, DoubleToRational(erfc));
155+ Assert.AreEqual(t.Y, erfc, Math.Abs(t.Y) * 1e-14, m);
156+ }
157+
158+ p.Clear();
159+ p.Add(16.0, 2.3284857515715299e-113);
160+ p.Add(20.0, 5.3958656116079012e-176);
161+ p.Add(24.0, 1.6489825831519335e-252);
162+ p.Add(28.0, 0.0);
163+ p.Add(100.0, 0.0);
164+
165+ foreach (Point t in p)
166+ {
167+ double erfc = MathEx.Erfc(t.X);
168+ string m = String.Format("x = {0}, y = {1}",
169+ t.X, DoubleToRational(erfc));
170+ Assert.AreEqual(t.Y, erfc, Math.Abs(t.Y) * 0.5, m);
171+ }
172+ }
173+ }
174+}
Show on old repository browser