Karinto/MathEx.cs

   1 /*\r
   2  *      Karinto Library Project\r
   3  *\r
   4  *      This software is distributed under a zlib-style license.\r
   5  *      See license.txt for more information.\r
   6  */\r
   7 \r
   8 using System;\r
   9 using System.Collections.Generic;\r
  10 using System.Text;\r
  11 \r
  12 namespace Karinto\r
  13 {\r
  14     public static class MathEx\r
  15     {\r
  16         #region fields\r
  17         private static readonly Polynomial[] erfcTable = \r
  18         {\r
  19             new Polynomial( //[0]: 1 <= x < 2\r
  20                 4.7019003296642450457550391704e-1,\r
  21                 -5.4129916560282193120142282350e-4,\r
  22                 3.2466045488989207057247716494e-2,\r
  23                 -1.5551574967900995569013129661e-3,\r
  24                 1.2991371652846980287617056420e-3,\r
  25                 -1.2476738559919909068871115854e-4,\r
  26                 4.2164948133638407169423314936e-5,\r
  27                 -5.5944555649744218134732620248e-6,\r
  28                 1.2391258283306184550290288637e-6,\r
  29                 -1.8629332352031266041799875695e-7,\r
  30                 3.3508489330180889238527074083e-8,\r
  31                 -5.1453759994172323811411458769e-9,\r
  32                 8.2778054007509566713316121080e-10,\r
  33                 -1.2419584612081402943997565209e-10,\r
  34                 1.8628225410158766603724913483e-11,\r
  35                 -2.6916642541339090011279256698e-12,\r
  36                 3.8321798025710167662264220379e-13,\r
  37                 -5.2985573947392393836199793226e-14,\r
  38                 7.2214593891835651145391699616e-15,\r
  39                 -1.0552316968012789176241975783e-15,\r
  40                 1.3848365749828243613160210837e-16\r
  41                 ),\r
  42             new Polynomial( //[1]: 2 <= x < 4\r
  43                 0.28246919208993532,//2.8246919208993525100243945507e-1,\r
  44                 1.0295429295263581287408466969e-4,\r
  45                 2.6352174925955979906965931429e-2,\r
  46                 -2.0840262584872572824603082922e-3,\r
  47                 1.6107317263933041407504226538e-3,\r
  48                 -2.6413257091217055241174636048e-4,\r
  49                 9.3720096524725604112989062675e-5,\r
  50                 -2.0394452777501772979594808178e-5,\r
  51                 5.4818843984699154429750547408e-6,\r
  52                 -1.2769242863171188612896799307e-6,\r
  53                 3.0855108190100308040276886416e-7,\r
  54                 -7.1110449954162265412516889674e-8,\r
  55                 1.6277854940656666077194365644e-8,\r
  56                 -3.6388350444865150940102551512e-9,\r
  57                 8.0125451384130244729739875906e-10,\r
  58                 -1.7333139256038089278659882792e-10,\r
  59                 3.6913251414732819339098560977e-11,\r
  60                 -7.7391112861863152629855204763e-12,\r
  61                 1.5985771117554783386097643078e-12,\r
  62                 -3.2558812235239520370430316160e-13,\r
  63                 6.5474397002885039401566856717e-14,\r
  64                 -1.2908539064863451579569117003e-14,\r
  65                 2.4663418649675516649805424542e-15,\r
  66                 -5.0102905606542145531985288605e-16,\r
  67                 1.1637629874200535139854357149e-16,\r
  68                 -1.7702372520803693617903759643e-17\r
  69                 ),\r
  70             new Polynomial( //[2]: 4 <= x < 8\r
  71                 2.6609916698112647641230096605e-1,\r
  72                 1.0052997692961233689565319255e-1,\r
  73                 0.063819586227166714,//6.3819586227166696717161903119e-2,\r
  74                 1.6510111747337704484097091109e-2,\r
  75                 6.7745775475928297665975539362e-3,\r
  76                 1.2613799574673550665257600345e-3,\r
  77                 4.6314050254419785813179921962e-4,\r
  78                 5.4816124629655721757719475954e-5,\r
  79                 2.5069740734847835425277961803e-5,\r
  80                 7.6494811710229433994592919775e-7,\r
  81                 1.3069891449561428965827894398e-6,\r
  82                 -1.1211722505095090109921762128e-7,\r
  83                 7.7627072370605473720667209070e-8,\r
  84                 -1.5067764445724103557487307330e-8,\r
  85                 5.3921138750029687859330065848e-9,\r
  86                 -1.3591247804716958040748207366e-9,\r
  87                 4.0341965879492537779635899830e-10,\r
  88                 -1.0884398300756496354533131958e-10,\r
  89                 3.0363217123340332343882922235e-11,\r
  90                 -8.2439513886349055414081918841e-12,\r
  91                 2.2756305305902268859393373330e-12,\r
  92                 -6.1246789459726400033402401708e-13,\r
  93                 1.4333003386843162197015125496e-13,\r
  94                 -3.7815369812690037058307445308e-14,\r
  95                 1.7401564599128446023009328203e-14,\r
  96                 -4.5875278594765951588791904401e-15\r
  97                 ),\r
  98             new Polynomial( //[3]: 8 <= x < 16\r
  99                 4.6854221014893762619745618301e-2,\r
 100                 -1.5511450952249084104775598329e-2,\r
 101                 5.1178905303441648577722059649e-3,\r
 102                 -1.6829798529769531049332840624e-3,\r
 103                 5.5160777130644620308102757353e-4,\r
 104                 -1.8020184996880086091090661904e-4,\r
 105                 5.8678514133519807272309470851e-5,\r
 106                 -1.9045977453678276537004114523e-5,\r
 107                 6.1623270905278656151589976968e-6,\r
 108                 -1.9875419915464967477768116933e-6,\r
 109                 6.3904356643246431611609250453e-7,\r
 110                 -2.0483278667408340516487849749e-7,\r
 111                 6.5453904816333940501972558207e-8,\r
 112                 -2.0852118199579889358624923805e-8,\r
 113                 6.6229050456699429842283832512e-9,\r
 114                 -2.0972627311983894339808366077e-9,\r
 115                 6.6237821648144864340143708663e-10,\r
 116                 -2.0853485442305791192334927397e-10,\r
 117                 6.5162153190849185531320932146e-11,\r
 118                 -2.0380108335416826794297357698e-11,\r
 119                 6.6463403803777970165860627602e-12,\r
 120                 -2.0777605497301746376639943675e-12,\r
 121                 4.6614566309518191650253075708e-13,\r
 122                 -1.4074713019480904125485292343e-13,\r
 123                 1.0831110710204089408511450325e-13,\r
 124                 -3.3906934201035073659517665140e-14\r
 125                 ),\r
 126             new Polynomial( //[4]: 16 <= x < 28\r
 127                 2.8262089312268317148835811995e-2,\r
 128                 -5.4152623609484073274229170763e-3,\r
 129                 1.6241399412435775154107409724e-3,\r
 130                 -5.0744105591774430405597196162e-3,\r
 131                 -1.7616157737971578130497741968e-2,\r
 132                 2.3083427166129394399962687197e-2,\r
 133                 1.1164742403248709989628167597e-1,\r
 134                 -3.5273994341997794105095275190e-2,\r
 135                 -3.0869428535741202544195312887e-1,\r
 136                 -1.3020536214165453885942308553e-2,\r
 137                 4.3332671745834241123676787595e-1,\r
 138                 9.4385604975718760529994282754e-2,\r
 139                 -3.0238469418692319645741714650e-1,\r
 140                 -9.5995136295717931771565803096e-2,\r
 141                 8.3313274670729987166299100438e-2,\r
 142                 3.1593997744278794785645602547e-2\r
 143                 )\r
 144         };\r
 145         #endregion\r
 146 \r
 147         #region public methods\r
 148         public static double Erf(double x)\r
 149         {\r
 150             double a = Math.Abs(x);\r
 151             if (a >= 1.0)\r
 152             {\r
 153                 return 1.0 - Erfc(x);\r
 154             }\r
 155             double y;\r
 156             double aa = a * a;\r
 157             if (a < 0.0625)\r
 158             {\r
 159                 double c1 = a * aa;\r
 160                 double c2 = c1 * aa;\r
 161                 double c3 = c2 * aa;\r
 162                 double c4 = c3 * aa;\r
 163 \r
 164                 double s1 = c1 * -315.0;// (3*5*7*9 / -3 / 1!);\r
 165                 double s2 = s1 * 0 + c2 * 189.0 / 2.0;// (3*5*7*9 / +5 / 2!)\r
 166                 double s3 = s2 * 0 + c3 * -45.0 / 2.0;// (3*5*7*9 / -7 / 3!)\r
 167                 double s4 = s3 * 0 + c4 * 35.0 / 8.0;// (3*5*7*9 / 9 / 4!)\r
 168                 y = a + (s1 + s2 + s3 + s4) / 945.0; //(3*5*7*9);\r
 169             }\r
 170             else\r
 171             {\r
 172                 double c = a;\r
 173                 double s = 0;\r
 174                 double b = 1.0;\r
 175                 for (double i = -1.0; i > -20.0; i-=1.0)\r
 176                 {\r
 177                     c *= aa / i;\r
 178                     b += 2.0;\r
 179                     s += c / b;\r
 180                 }\r
 181                 y = a + s;\r
 182             }\r
 183             y *= 1.1283791670955126;     // 2 / pi^0.5\r
 184             return Math.Sign(x) > 0 ? y : -y;\r
 185         }\r
 186 \r
 187         public static double Erfc(double x)\r
 188         {\r
 189             double a = Math.Abs(x);\r
 190             if (a < 1.0)\r
 191             {\r
 192                 return 1.0 - Erf(x);\r
 193             }\r
 194             if (a > 27.226017111108501)\r
 195             {\r
 196                 return x > 0 ? 0 : 2.0;\r
 197             }\r
 198             double y = 0.0;\r
 199 \r
 200             if (a < 2.0)\r
 201             {\r
 202                 y = Math.Exp(-1.1688327773406599 * a * a);\r
 203                 y *= erfcTable[0].GetValueAt(a + a - 3.0);\r
 204             }\r
 205             else if (a < 4.0)\r
 206             {\r
 207                 decimal b = -1.05068636145441m * (decimal)a * (decimal)a;\r
 208                 y = Math.Exp((double)b);\r
 209                 y *= erfcTable[1].GetValueAt(a - 3.0);\r
 210             }\r
 211             else if (a < 8.0)\r
 212             {\r
 213                 decimal b = -1.0292687483818601m * (decimal)a * (decimal)a;\r
 214                 y = Math.Exp((double)b);\r
 215                 y *= erfcTable[2].GetValueAt(a * 0.5 - 3);\r
 216             }\r
 217             else if (a < 16.0)\r
 218             {\r
 219                 y = Math.Exp(-a * a);\r
 220                 y *= erfcTable[3].GetValueAt(a * 0.25 - 3);\r
 221             }\r
 222             else\r
 223             {\r
 224                 y = Math.Exp(-a * a);\r
 225                 y *= erfcTable[4].GetValueAt(a * 0.125 - 3);\r
 226             }\r
 227             return x > 0 ? y : 2.0 - y;\r
 228         }\r
 229         #endregion\r
 230     }\r
 231 }\r