home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / c / cephes.arc / NDTRI.C < prev    next >
Encoding:
C/C++ Source or Header  |  1987-06-20  |  10.0 KB  |  375 lines
  1. /*                            ndtri.c
  2.  *
  3.  *    Inverse of Normal distribution function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, ndtri();
  10.  *
  11.  * x = ndtri( y );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the argument, x, for which the area under the
  18.  * Gaussian probability density function (integrated from
  19.  * minus infinity to x) is equal to y.
  20.  *
  21.  *
  22.  * For small arguments 0 < y < exp(-2), the program computes
  23.  * z = sqrt( -2.0 * log(y) );  then the approximation is
  24.  * x = z - (log(z) - log(sqrt(2pi)))/z  -  P(-z)/Q(-z).
  25.  * There are two rational functions P/Q, one for 0 < y < 2e-40
  26.  * and the other for y up to exp(-2).  For larger arguments,
  27.  * z = y - 0.5, and  x = sqrt(2pi) (z + z**3 P4(z**2)/Q8(z**2)).
  28.  *
  29.  *
  30.  * ACCURACY:
  31.  *
  32.  *                      Relative error:
  33.  * arithmetic   range         # trials      peak         rms
  34.  *    DEC      0.125, 1        11000       7.6e-17     1.6e-17
  35.  *    DEC      6e-39, 0.135    10000       6.0e-17     1.3e-17
  36.  *    IEEE     0.125, 1         6500       7.1e-16     1.3e-16
  37.  *    IEEE     3e-308, 0.135    5000       4.1e-16     9.5e-17
  38.  *
  39.  *
  40.  * ERROR MESSAGES:
  41.  *
  42.  *   message         condition    value returned
  43.  * ndtri domain       x <= 0        -MAXNUM
  44.  * ndtri domain       x >= 1         MAXNUM
  45.  *
  46.  */
  47.  
  48.  
  49. /*
  50. Cephes Math Library Release 2.0:  April, 1987
  51. Copyright 1984, 1987 by Stephen L. Moshier
  52. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  53. */
  54.  
  55. #include "mconf.h"
  56. extern double MAXNUM;
  57.  
  58. #ifdef UNK
  59. static double P0[] = {
  60. -5.99633501014107891266E1,
  61.  9.80010754185999654266E1,
  62. -5.66762857469070288892E1,
  63.  1.39312609387279678315E1,
  64. -1.23916583867381256906E0
  65. };
  66. static double Q0[] = {
  67. /*  1.00000000000000000000E0, */
  68.  1.95448858338141759078E0,
  69.  4.67627912898881536003E0,
  70.  8.63602421390890584914E1,
  71. -2.25462687854119368911E2,
  72.  2.00260212380060658814E2,
  73. -8.20372256168333333195E1,
  74.  1.59056225126211694148E1,
  75. -1.18331621121330002082E0
  76. };
  77. static double P1[] = {
  78. -8.06295164062928296467E-7,
  79. -1.24053576627339949494E-4,
  80.  6.49932444421283598941E-3,
  81. -1.75578085015445987182E1,
  82.  1.13340546911064577348E3,
  83. -1.55308233703367942177E4,
  84.  5.96478957762231923216E4,
  85. -1.34118762557568719035E5,
  86.  1.33922736098486867793E5,
  87. -7.68504745178723676740E4
  88. };
  89. static double Q1[] = {
  90. /*  1.00000000000000000000E0,*/
  91. -1.47602102962036448989E2,
  92.  4.36545092375268950884E3,
  93. -4.06167070356628991869E4,
  94.  1.54083842563020896864E5,
  95. -3.52535896087944211765E5,
  96.  3.74088548110300571177E5,
  97. -2.40053383557002665939E5,
  98.  3.38930498523655775059E4,
  99.  9.17564445564876422447E2
  100. };
  101. static double P2[] = {
  102.  8.94084651318768121212E-8,
  103.  5.74305670024326318968E-5,
  104.  2.89444514333676761949E-2,
  105. -1.42646079401864981264E1,
  106.  5.65275920555989502947E2,
  107.  1.20932159059563133049E3,
  108. -7.30662837143628743787E4
  109. };
  110. static double Q2[] = {
  111. /* 1.00000000000000000000E0,*/
  112. -9.42696522534933841013E1,
  113.  1.32905713172348184767E3,
  114.  1.30035792204994811143E4,
  115. -1.45436997023716148384E5,
  116.  8.82097696844164863759E4,
  117.  5.20585037820295592936E4
  118. };
  119. /* -log(sqrt(2pi)) */
  120. static double lstpi = -9.18938533204672741780E-1;
  121. /* sqrt(2pi) */
  122. static double s2pi = 2.50662827463100050242E0;
  123. #endif
  124.  
  125. #ifdef DEC
  126. static short P0[] = {
  127. 0141557,0155170,0071360,0120550,
  128. 0041704,0000214,0172417,0067307,
  129. 0141542,0132204,0040066,0156722,
  130. 0041136,0163161,0157276,0007746,
  131. 0140236,0116374,0073666,0051764
  132. };
  133. static short Q0[] = {
  134. /*0040200,0000000,0000000,0000000,*/
  135. 0040372,0026256,0110403,0123707,
  136. 0040625,0122024,0020277,0026661,
  137. 0041654,0134161,0124134,0007243,
  138. 0142141,0073162,0133021,0131371,
  139. 0042110,0041235,0043516,0057766,
  140. 0141644,0011417,0036155,0137305,
  141. 0041176,0076556,0004043,0125427,
  142. 0140227,0073347,0152776,0067251
  143. };
  144. static short P1[] = {
  145. 0133130,0070056,0104154,0123744,
  146. 0135002,0012140,0157600,0156617,
  147. 0036324,0174110,0173614,0163423,
  148. 0141214,0073144,0046674,0102277,
  149. 0042615,0126371,0115133,0034141,
  150. 0143562,0125513,0020627,0172034,
  151. 0044150,0177745,0050627,0027445,
  152. 0144402,0174660,0146676,0041276,
  153. 0044402,0144257,0016074,0151757,
  154. 0144226,0014474,0136400,0065630
  155. };
  156. static short Q1[] = {
  157. /*0040200,0000000,0000000,0000000,*/
  158. 0142023,0115043,0065562,0142562,
  159. 0043210,0065633,0076751,0113124,
  160. 0144036,0124265,0000112,0004502,
  161. 0044426,0074365,0166215,0071345,
  162. 0144654,0021374,0126300,0117727,
  163. 0044666,0124421,0105036,0116362,
  164. 0144552,0066530,0106062,0125647,
  165. 0044004,0062414,0141437,0163711,
  166. 0042545,0062037,0160112,0060416
  167. };
  168. static short P2[] = {
  169. 0032300,0000322,0151003,0072465,
  170. 0034560,0160632,0006537,0051737,
  171. 0036755,0016352,0004722,0136054,
  172. 0141144,0035725,0104425,0147223,
  173. 0042415,0050650,0127261,0010470,
  174. 0042627,0025112,0074134,0057002,
  175. 0144216,0132444,0050300,0111372
  176. };
  177. static short Q2[] = {
  178. /*0040200,0000000,0000000,0000000,*/
  179. 0141674,0105017,0156064,0013764,
  180. 0042646,0020724,0002750,0076546,
  181. 0043513,0027121,0017455,0134707,
  182. 0144416,0003477,0147474,0107614,
  183. 0044254,0044342,0102404,0155206,
  184. 0044113,0055200,0173733,0166505
  185. };
  186. static short lstp[] = {0140153,0037616,0041445,0172645};
  187. #define lstpi *(double *)lstp
  188. static short s2p[] = {0040440,0066230,0177661,0034055};
  189. #define s2pi *(double *)s2p
  190. #endif
  191.  
  192. #ifdef IBMPC
  193. static short P0[] = {
  194. 0x142d,0x0e5e,0xfb4f,0xc04d,
  195. 0xedd9,0x9ea1,0x8011,0x4058,
  196. 0xdbba,0x8806,0x5690,0xc04c,
  197. 0xc1fd,0x3bd7,0xdcce,0x402b,
  198. 0xca7e,0x8ef6,0xd39f,0xbff3
  199. };
  200. static short Q0[] = {
  201. /*0x0000,0x0000,0x0000,0x3ff0,*/
  202. 0x74f9,0xd220,0x4595,0x3fff,
  203. 0xe5b6,0x8417,0xb482,0x4012,
  204. 0x81d4,0x350b,0x970e,0x4055,
  205. 0x365f,0x56c2,0x2ece,0xc06c,
  206. 0xcbff,0xa8e9,0x0853,0x4069,
  207. 0xb7d9,0xe78d,0x8261,0xc054,
  208. 0x7563,0xc104,0xcfad,0x402f,
  209. 0xcdd5,0xfabf,0xeedc,0xbff2
  210. };
  211. static short P1[] = {
  212. 0x94fd,0xd10d,0x0e05,0xbeab,
  213. 0x1bb2,0x1bf0,0x428c,0xbf20,
  214. 0x9ce2,0x1ef1,0x9f09,0x3f7a,
  215. 0x9098,0x89b7,0x8ecc,0xc031,
  216. 0x670c,0x334b,0xb59f,0x4091,
  217. 0xfe84,0x6432,0x5569,0xc0ce,
  218. 0xe5e5,0xaa32,0x1ffc,0x40ed,
  219. 0xc858,0x19b7,0x5f36,0xc100,
  220. 0x9a7e,0xe387,0x5915,0x4100,
  221. 0x0d73,0x97a0,0xc327,0xc0f2
  222. };
  223. static short Q1[] = {
  224. /*0x0000,0x0000,0x0000,0x3ff0,*/
  225. 0x58ae,0x6d6e,0x7344,0xc062,
  226. 0x32cb,0x6fbd,0x0d73,0x40b1,
  227. 0x4128,0xa009,0xd516,0xc0e3,
  228. 0xae5d,0xbd91,0xcf1e,0x4102,
  229. 0x13fb,0x9598,0x845f,0xc115,
  230. 0xd39e,0x3143,0xd522,0x4116,
  231. 0x5575,0x1186,0x4dab,0xc10d,
  232. 0xfcf9,0x9863,0x8ca1,0x40e0,
  233. 0x4c22,0xfc09,0xac83,0x408c
  234. };
  235. static short P2[] = {
  236. 0x6ea7,0x5a40,0x001a,0x3e78,
  237. 0xea7c,0x41ab,0x1c33,0x3f0e,
  238. 0x5785,0x413a,0xa39d,0x3f9d,
  239. 0xb9d2,0xb122,0x877a,0xc02c,
  240. 0x2227,0x15d6,0xaa35,0x4081,
  241. 0x8bc0,0x4f0b,0xe549,0x4092,
  242. 0x125f,0x8a18,0xd6a4,0xc0f1
  243. };
  244. static short Q2[] = {
  245. /*0x0000,0x0000,0x0000,0x3ff0,*/
  246. 0x82ff,0xfb86,0x9141,0xc057,
  247. 0x0fad,0x80bd,0xc43a,0x4094,
  248. 0xb739,0x23e5,0x65ca,0x40c9,
  249. 0x91f2,0xf9e7,0xc0e7,0xc101,
  250. 0x9b51,0x50a0,0x891c,0x40f5,
  251. 0x7da9,0x1efb,0x6b50,0x40e9
  252. };
  253. static short lstp[] = {0xbeb5,0xc864,0x67f1,0xbfed};
  254. #define lstpi *(double *)lstp
  255. static short s2p[] = {0x2706,0x1ff6,0x0d93,0x4004};
  256. #define s2pi *(double *)s2p
  257. #endif
  258.  
  259. #ifdef MIEEE
  260. static short P0[] = {
  261. 0xc04d,0xfb4f,0x0e5e,0x142d,
  262. 0x4058,0x8011,0x9ea1,0xedd9,
  263. 0xc04c,0x5690,0x8806,0xdbba,
  264. 0x402b,0xdcce,0x3bd7,0xc1fd,
  265. 0xbff3,0xd39f,0x8ef6,0xca7e
  266. };
  267. static short Q0[] = {
  268. 0x3fff,0x4595,0xd220,0x74f9,
  269. 0x4012,0xb482,0x8417,0xe5b6,
  270. 0x4055,0x970e,0x350b,0x81d4,
  271. 0xc06c,0x2ece,0x56c2,0x365f,
  272. 0x4069,0x0853,0xa8e9,0xcbff,
  273. 0xc054,0x8261,0xe78d,0xb7d9,
  274. 0x402f,0xcfad,0xc104,0x7563,
  275. 0xbff2,0xeedc,0xfabf,0xcdd5
  276. };
  277. static short P1[] = {
  278. 0xbeab,0x0e05,0xd10d,0x94fd,
  279. 0xbf20,0x428c,0x1bf0,0x1bb2,
  280. 0x3f7a,0x9f09,0x1ef1,0x9ce2,
  281. 0xc031,0x8ecc,0x89b7,0x9098,
  282. 0x4091,0xb59f,0x334b,0x670c,
  283. 0xc0ce,0x5569,0x6432,0xfe84,
  284. 0x40ed,0x1ffc,0xaa32,0xe5e5,
  285. 0xc100,0x5f36,0x19b7,0xc858,
  286. 0x4100,0x5915,0xe387,0x9a7e,
  287. 0xc0f2,0xc327,0x97a0,0x0d73
  288. };
  289. static short Q1[] = {
  290. 0xc062,0x7344,0x6d6e,0x58ae,
  291. 0x40b1,0x0d73,0x6fbd,0x32cb,
  292. 0xc0e3,0xd516,0xa009,0x4128,
  293. 0x4102,0xcf1e,0xbd91,0xae5d,
  294. 0xc115,0x845f,0x9598,0x13fb,
  295. 0x4116,0xd522,0x3143,0xd39e,
  296. 0xc10d,0x4dab,0x1186,0x5575,
  297. 0x40e0,0x8ca1,0x9863,0xfcf9,
  298. 0x408c,0xac83,0xfc09,0x4c22
  299. };
  300. static short P2[] = {
  301. 0x3e78,0x001a,0x5a40,0x6ea7,
  302. 0x3f0e,0x1c33,0x41ab,0xea7c,
  303. 0x3f9d,0xa39d,0x413a,0x5785,
  304. 0xc02c,0x877a,0xb122,0xb9d2,
  305. 0x4081,0xaa35,0x15d6,0x2227,
  306. 0x4092,0xe549,0x4f0b,0x8bc0,
  307. 0xc0f1,0xd6a4,0x8a18,0x125f
  308. };
  309. static short Q2[] = {
  310. 0xc057,0x9141,0xfb86,0x82ff,
  311. 0x4094,0xc43a,0x80bd,0x0fad,
  312. 0x40c9,0x65ca,0x23e5,0xb739,
  313. 0xc101,0xc0e7,0xf9e7,0x91f2,
  314. 0x40f5,0x891c,0x50a0,0x9b51,
  315. 0x40e9,0x6b50,0x1efb,0x7da9
  316. };
  317. static short lstp[] = {
  318. 0xbfed,0x67f1,0xc864,0xbeb5
  319. };
  320. #define lstpi *(double *)lstp
  321. static short s2p[] = {
  322. 0x4004,0x0d93,0x1ff6,0x2706
  323. };
  324. #define s2pi *(double *)s2p
  325. #endif
  326.  
  327.  
  328. double ndtri(y0)
  329. double y0;
  330. {
  331. double d, x, y, y2, x0, x1;
  332. int code;
  333. double polevl(), p1evl(), log(), sqrt();
  334.  
  335. if( y0 <= 0.0 )
  336.     {
  337.     mtherr( "ndtri", DOMAIN );
  338.     return( -MAXNUM );
  339.     }
  340. if( y0 >= 1.0 )
  341.     {
  342.     mtherr( "ndtri", DOMAIN );
  343.     return( MAXNUM );
  344.     }
  345. code = 1;
  346. y = y0;
  347. if( y > (1.0 - 0.13533528323661269189) ) /* 0.135... = exp(-2) */
  348.     {
  349.     y = 1.0 - y;
  350.     code = 0;
  351.     }
  352.  
  353. if( y > 0.13533528323661269189 )
  354.     {
  355.     y = y - 0.5;
  356.     y2 = y * y;
  357.     x = y + y * (y2 * polevl( y2, P0, 4)/p1evl( y2, Q0, 8 ));
  358.     x = x * s2pi; 
  359.     return(x);
  360.     }
  361.  
  362. x = sqrt( -2.0 * log(y) );
  363. x0 = x - (log(x) - lstpi)/x;
  364.  
  365. x = -x;
  366. if( x > -13.5 ) /* y > 2.6602064e-40 */
  367.     x1 = polevl( x, P1, 9 )/p1evl( x, Q1, 9 );
  368. else
  369.     x1 = polevl( x, P2, 6 )/p1evl( x, Q2, 6 );
  370. x = x0 - x1;
  371. if( code != 0 )
  372.     x = -x;
  373. return( x );
  374. }
  375.