htswanalysis/MACS/lib/gsl/gsl-1.11/specfunc/TODO

   1 * Could probably return immediately for exact zeros in 3j,6j,9j
   2 functions.  Easiest to implement for 3j.
   3
   4   Note from Serge Winitzki <serge@cosmos.phy.tufts.edu>:
   5
   6   The package "matpack" (www.matpack.de) includes many special functions,
   7   also the 3j symbols. They refer to some quite complicated numerical
   8   methods using recursion relations to get the right answers for large
   9   momenta, and to 1975-1976 papers by Schulten and Gordon for the
  10   description of the algorithms. The papers can be downloaded for free at
  11   http://www.ks.uiuc.edu/Publications/Papers/
  12
  13   http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU76B
  14   http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75A
  15   http://www.ks.uiuc.edu/Publications/Papers/abstract.cgi?tbcode=SCHU75
  16
  17 * add Fresnel Integrals to specfunc.  See TOMS 723 + 2 subsequent
  18 errata.
  19
  20 * make mode variables consistent in specfunc -- some seem to be
  21 unnecessary from performance point of view since the speed difference
  22 is negligible.
  23
  24 * From: "Alexander Babansky" <babansky@mail.ru>
  25 To: "Brian Gough" <bjg@network-theory.co.uk>
  26 Subject: Re: gsl-1.2
  27 Date: Sun, 3 Nov 2002 14:15:15 -0500
  28
  29 Hi Brian,
  30 May I suggest you to add another function to gsl-1.2 ?
  31 It's a modified Ei(x) function:
  32
  33 Em(x)=exp(-x)*Ei(x);
  34
  35 As u  might know, Ei(x) raises as e^x on the negative interval.
  36 Therefore, Ei(100) is very very large.
  37 But Ei(100)*exp(-100) = 0.010;
  38
  39 Unfortunately, if u try x=800 u'll get overflow in Ei(800).
  40 but Ei(800)*exp(-800) should be around 0.0001;
  41
  42 Modified function Em(x) is used in cos, sin integrals such as:
  43 int_0^\infinity dx sin(bx)/(x^2+z^2)=(1/2z)*(Em(bz)-Em(-bz));
  44
  45 int_0^\infinity dx x cos(bx)/(x^2+z^2)=(1/2)*(Em(bz)+Em(-bz));
  46
  47 One of possible ways to add it to the library is:
  48 Em(x) = - PV int_0^\infinity e^(-t)/(t+x) dt
  49
  50 Sincerely,
  51 Alex
  52
  53 DONE: Wed Nov  6 13:06:42 MST 2002 [GJ]
  54
  55
  56 ----------------------------------------------------------------------
  57
  58 The following should be finished before a 1.0 level release.
  59
  60 * Implement the conicalP_sph_reg() functions.
  61   DONE: Fri Nov  6 23:33:53 MST 1998 [GJ]
  62
  63 * Irregular (Q) Legendre functions, at least
  64   the integer order ones. More general cases
  65   can probably wait.
  66   DONE: Sat Nov  7 15:47:35 MST 1998 [GJ]
  67
  68 * Make hyperg_1F1() work right.
  69   This is the last remaining source of test failures.
  70   The problem is with an unstable recursion in certain cases.
  71   Look for the recursion with the variable named "start_pair";
  72   this is stupid hack to keep track of when the recursion
  73   result is going the wrong way for awhile by remembering the
  74   minimum value. An error estimate is amde from that. But it
  75   is just a hack. Somethign must be done abou that case.
  76
  77 * Clean-up Coulomb wave functions. This does not
  78   mean completing a fully controlled low-energy
  79   evaluation, which is a larger project.
  80   DONE: Sun May 16 13:49:47 MDT 1999 [GJ]
  81
  82 * Clean-up the Fermi-Dirac code. The full Fermi-Dirac
  83   functions can probably wait until a later release,
  84   but we should have at least the common j = integer and
  85   j = 1/2-integer cases for the 1.0 release. These
  86   are not too hard.
  87   DONE: Sat Nov  7 19:46:27 MST 1998 [GJ]
  88
  89 * Go over the tests and make sure nothing is left out.
  90
  91 * Sanitize all the error-checking, error-estimation,
  92   algorithm tuning, etc.
  93
  94 * Fill out our scorecard, working from Lozier's
  95   "Software Needs in Special Functions" paper.
  96
  97 * Final Seal of Approval
  98   This section has itself gone through several
  99   revisions (sigh), proving that the notion of
 100   done-ness is ill-defined. So it is worth
 101   stating the criteria for done-ness explicitly:
 102   o interfaces stabilized
 103   o error-estimation in place
 104   o all deprecated constructs removed
 105   o passes tests
 106
 107
 108   - airy.c
 109     INTERFACES:
 110     ERRORESTIM:
 111     DEPRECATED:
 112     PASSTESTS:
 113
 114   - airy_der.c
 115     INTERFACES:
 116     ERRORESTIM:
 117     DEPRECATED:
 118     PASSTESTS:
 119
 120   - airy_zero.c
 121     INTERFACES:
 122     ERRORESTIM:
 123     DEPRECATED:
 124     PASSTESTS:
 125
 126   - atanint.c
 127     INTERFACES:
 128     ERRORESTIM:
 129     DEPRECATED:
 130     PASSTESTS:
 131
 132   - bessel.c
 133     INTERFACES:
 134     ERRORESTIM:
 135     DEPRECATED:
 136     PASSTESTS:
 137
 138   - bessel_I0.c
 139     INTERFACES:
 140     ERRORESTIM:
 141     DEPRECATED:
 142     PASSTESTS:
 143
 144   - bessel_I1.c
 145     INTERFACES:
 146     ERRORESTIM:
 147     DEPRECATED:
 148     PASSTESTS:
 149
 150   - bessel_In.c
 151     INTERFACES:
 152     ERRORESTIM:
 153     DEPRECATED:
 154     PASSTESTS:
 155
 156   - bessel_Inu.c
 157     INTERFACES:
 158     ERRORESTIM:
 159     DEPRECATED:
 160     PASSTESTS:
 161
 162   - bessel_J0.c
 163     INTERFACES:
 164     ERRORESTIM:
 165     DEPRECATED:
 166     PASSTESTS:
 167
 168   - bessel_J1.c
 169     INTERFACES:
 170     ERRORESTIM:
 171     DEPRECATED:
 172     PASSTESTS:
 173
 174   - bessel_Jn.c
 175     INTERFACES:
 176     ERRORESTIM:
 177     DEPRECATED:
 178     PASSTESTS:
 179
 180   - bessel_Jnu.c
 181     INTERFACES:
 182     ERRORESTIM:
 183     DEPRECATED:
 184     PASSTESTS:
 185
 186   - bessel_K0.c
 187     INTERFACES:
 188     ERRORESTIM:
 189     DEPRECATED:
 190     PASSTESTS:
 191
 192   - bessel_K1.c
 193     INTERFACES:
 194     ERRORESTIM:
 195     DEPRECATED:
 196     PASSTESTS:
 197
 198   - bessel_Kn.c
 199     INTERFACES:
 200     ERRORESTIM:
 201     DEPRECATED:
 202     PASSTESTS:
 203
 204   - bessel_Knu.c
 205     INTERFACES:
 206     ERRORESTIM:
 207     DEPRECATED:
 208     PASSTESTS:
 209
 210   - bessel_Y0.c
 211     INTERFACES:
 212     ERRORESTIM:
 213     DEPRECATED:
 214     PASSTESTS:
 215
 216   - bessel_Y1.c
 217     INTERFACES:
 218     ERRORESTIM:
 219     DEPRECATED:
 220     PASSTESTS:
 221
 222   - bessel_Yn.c
 223     INTERFACES:
 224     ERRORESTIM:
 225     DEPRECATED:
 226     PASSTESTS:
 227
 228   - bessel_Ynu.c
 229     INTERFACES:
 230     ERRORESTIM:
 231     DEPRECATED:
 232     PASSTESTS:
 233
 234   - bessel_amp_phase.c
 235     INTERFACES:
 236     ERRORESTIM:
 237     DEPRECATED:
 238     PASSTESTS:
 239
 240   - bessel_i.c
 241     INTERFACES:
 242     ERRORESTIM:
 243     DEPRECATED:
 244     PASSTESTS:
 245
 246   - bessel_j.c
 247     INTERFACES:
 248     ERRORESTIM:
 249     DEPRECATED:
 250     PASSTESTS:
 251
 252   - bessel_k.c
 253     INTERFACES:
 254     ERRORESTIM:
 255     DEPRECATED:
 256     PASSTESTS:
 257
 258   - bessel_olver.c
 259     INTERFACES:
 260     ERRORESTIM:
 261     DEPRECATED:
 262     PASSTESTS:
 263
 264   - bessel_sequence.c
 265     INTERFACES:
 266     ERRORESTIM:
 267     DEPRECATED:
 268     PASSTESTS:
 269
 270   - bessel_temme.c
 271     INTERFACES:
 272     ERRORESTIM:
 273     DEPRECATED:
 274     PASSTESTS:
 275
 276   - bessel_y.c
 277     INTERFACES:
 278     ERRORESTIM:
 279     DEPRECATED:
 280     PASSTESTS:
 281
 282   - bessel_zero.c
 283     INTERFACES:
 284     ERRORESTIM:
 285     DEPRECATED:
 286     PASSTESTS:
 287
 288   - beta.c
 289     INTERFACES:
 290     ERRORESTIM:
 291     DEPRECATED:
 292     PASSTESTS:
 293
 294   - chebyshev.c
 295     INTERFACES:
 296     ERRORESTIM:
 297     DEPRECATED:
 298     PASSTESTS:
 299
 300   - clausen.c
 301     INTERFACES:
 302     ERRORESTIM:
 303     DEPRECATED:
 304     PASSTESTS:
 305
 306   - coulomb.c
 307     INTERFACES:
 308     ERRORESTIM:
 309     DEPRECATED:
 310     PASSTESTS:
 311
 312   - coulomb_bound.c
 313     INTERFACES:
 314     ERRORESTIM:
 315     DEPRECATED:
 316     PASSTESTS:
 317
 318   - coupling.c
 319     INTERFACES:
 320     ERRORESTIM:
 321     DEPRECATED:
 322     PASSTESTS:
 323
 324   - dawson.c
 325     INTERFACES:
 326     ERRORESTIM:
 327     DEPRECATED:
 328     PASSTESTS:
 329
 330   - debye.c
 331     INTERFACES:
 332     ERRORESTIM:
 333     DEPRECATED:
 334     PASSTESTS:
 335
 336   - dilog.c
 337     INTERFACES:
 338     ERRORESTIM:
 339     DEPRECATED:
 340     PASSTESTS:
 341
 342   - elementary.c
 343     INTERFACES:
 344     ERRORESTIM:
 345     DEPRECATED:
 346     PASSTESTS:
 347
 348   - ellint.c
 349     INTERFACES:
 350     ERRORESTIM:
 351     DEPRECATED:
 352     PASSTESTS:
 353
 354   - elljac.c
 355     INTERFACES:
 356     ERRORESTIM:
 357     DEPRECATED:
 358     PASSTESTS:
 359
 360   - erfc.c
 361     INTERFACES:
 362     ERRORESTIM:
 363     DEPRECATED:
 364     PASSTESTS:
 365
 366   - exp.c
 367     INTERFACES:
 368     ERRORESTIM:
 369     DEPRECATED:
 370     PASSTESTS:
 371
 372   - expint.c
 373     INTERFACES:
 374     ERRORESTIM:
 375     DEPRECATED:
 376     PASSTESTS:
 377
 378   - expint3.c
 379     INTERFACES:
 380     ERRORESTIM:
 381     DEPRECATED:
 382     PASSTESTS:
 383
 384   - fermi_dirac.c
 385     INTERFACES:
 386     ERRORESTIM:
 387     DEPRECATED:
 388     PASSTESTS:
 389
 390   - gamma.c
 391     INTERFACES:
 392     ERRORESTIM:
 393     DEPRECATED:
 394     PASSTESTS:
 395
 396   - gamma_inc.c
 397     INTERFACES:
 398     ERRORESTIM:
 399     DEPRECATED:
 400     PASSTESTS:
 401
 402   - gegenbauer.c
 403     INTERFACES:
 404     ERRORESTIM:
 405     DEPRECATED:
 406     PASSTESTS:
 407
 408   - hyperg.c
 409     INTERFACES:
 410     ERRORESTIM:
 411     DEPRECATED:
 412     PASSTESTS:
 413
 414   - hyperg_0F1.c
 415     INTERFACES:
 416     ERRORESTIM:
 417     DEPRECATED:
 418     PASSTESTS:
 419
 420   - hyperg_1F1.c
 421
 422   - hyperg_2F0.c
 423     INTERFACES:
 424     ERRORESTIM:
 425     DEPRECATED:
 426     PASSTESTS:
 427
 428   - hyperg_2F1.c
 429     INTERFACES:
 430     ERRORESTIM:
 431     DEPRECATED:
 432     PASSTESTS:
 433
 434   - hyperg_U.c
 435     INTERFACES:
 436     ERRORESTIM:
 437     DEPRECATED:
 438     PASSTESTS:
 439
 440   - laguerre.c
 441     INTERFACES:
 442     ERRORESTIM:
 443     DEPRECATED:
 444     PASSTESTS:
 445
 446   - legendre_H3d.c
 447     INTERFACES:
 448     ERRORESTIM:
 449     DEPRECATED:
 450     PASSTESTS:
 451
 452   - legendre_Qn.c
 453     INTERFACES:
 454     ERRORESTIM:
 455     DEPRECATED:
 456     PASSTESTS:
 457
 458   - legendre_con.c
 459     INTERFACES:
 460     ERRORESTIM:
 461     DEPRECATED:
 462     PASSTESTS:
 463
 464   - legendre_poly.c
 465     INTERFACES:
 466     ERRORESTIM:
 467     DEPRECATED:
 468     PASSTESTS:
 469
 470   - log.c
 471     INTERFACES:
 472     ERRORESTIM:
 473     DEPRECATED:
 474     PASSTESTS:
 475
 476   - poch.c
 477     INTERFACES:
 478     ERRORESTIM:
 479     DEPRECATED:
 480     PASSTESTS:
 481
 482   - poly.c
 483     INTERFACES:
 484     ERRORESTIM:
 485     DEPRECATED:
 486     PASSTESTS:
 487
 488   - pow_int.c
 489     INTERFACES:
 490     ERRORESTIM:
 491     DEPRECATED:
 492     PASSTESTS:
 493
 494   - psi.c
 495     INTERFACES:
 496     ERRORESTIM:
 497     DEPRECATED:
 498     PASSTESTS:
 499
 500   - result.c
 501     INTERFACES:
 502     ERRORESTIM:
 503     DEPRECATED:
 504     PASSTESTS:
 505
 506   - shint.c
 507     INTERFACES:
 508     ERRORESTIM:
 509     DEPRECATED:
 510     PASSTESTS:
 511
 512   - sinint.c
 513     INTERFACES:
 514     ERRORESTIM:
 515     DEPRECATED:
 516     PASSTESTS:
 517
 518   - synchrotron.c
 519     INTERFACES:
 520     ERRORESTIM:
 521     DEPRECATED:
 522     PASSTESTS:
 523
 524   - transport.c
 525     INTERFACES:
 526     ERRORESTIM:
 527     DEPRECATED:
 528     PASSTESTS:
 529
 530   - trig.c
 531     INTERFACES:
 532     ERRORESTIM:
 533     DEPRECATED:
 534     PASSTESTS:
 535
 536   - zeta.c
 537     INTERFACES:
 538     ERRORESTIM:
 539     DEPRECATED:
 540     PASSTESTS:
 541
 542
 543
 544 ----------------------------------------------------------------------
 545
 546 The following are important but probably will
 547 not see completion before a 1.0 level release.
 548
 549 * Incomplete Fermi-Dirac functions.
 550   Other Fermi-Dirac functions, including the
 551   generic 1/2-integer case, which was not done.
 552
 553 * Implement the low-energy regime for the Coulomb
 554   wave functions. This is fairly well understood in
 555   the recent literature but will require some
 556   detailed work. Specifically this means creating
 557   a drop-in replacement for coulomb_jwkb() which
 558   is controlled and extensible.
 559
 560 * General Legendre functions (at least on the cut).
 561   This subsumes the toroidal functions, so we need not
 562   consider those separately. SLATEC code exists (originally
 563   due to Olver+Smith).
 564
 565 * Characterize the algorithms. A significant fraction of
 566   the code is home-grown and it should be reviewed by
 567   other parties.
 568
 569
 570 ----------------------------------------------------------------------
 571
 572 The following are extra features which need not
 573 be implemented for a version 1.0 release.
 574
 575 * Spheroidal wave functions.
 576
 577 * Mathieu functions.
 578
 579 * Weierstrass elliptic functions.
 580
 581
 582 ----------------------------------------------------------------------
 583
 584 Improve accuracy of ERF
 585
 586 NNTP-Posting-Date: Thu, 11 Sep 2003 07:41:42 -0500
 587 From: "George Marsaglia" <geo@stat.fsu.edu>
 588 Newsgroups: comp.lang.c
 589 References: <t4J7b.18514$98.4310@nwrddc03.gnilink.net>
 590 Subject: Re: When (32-bit) double precision isn't precise enough
 591 Date: Thu, 11 Sep 2003 08:41:40 -0400
 592 X-Priority: 3
 593 X-MSMail-Priority: Normal
 594 X-Newsreader: Microsoft Outlook Express 6.00.2800.1158
 595 X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2800.1165
 596 Message-ID: <wq2dnQBikNwb8P2iU-KYvg@comcast.com>
 597 Lines: 265
 598 NNTP-Posting-Host: 68.35.247.101
 599 X-Trace: sv3-4YY+jkhhdeQvGKAREa99vDBFHJoKVqVBdUTSuRxA71OwlgxX0uUFnKYs54FlnUs0Xb6BRngKigkd75d!tKin8l8rAQKylaP+4vzTI3AO33bivOw1lKDZUUtXe4lUMW1qn+goUp/Pfksstg==
 600 X-Complaints-To: abuse@comcast.net
 601 X-DMCA-Complaints-To: dmca@comcast.net
 602 X-Abuse-and-DMCA-Info: Please be sure to forward a copy of ALL headers
 603 X-Abuse-and-DMCA-Info: Otherwise we will be unable to process your complaint properly
 604 X-Postfilter: 1.1
 605
 606
 607 Why most of those who deal with the normal integral in probability
 608 theory are still stuck with the historical baggage of the error function
 609 is a puzzle to me, as is the poor quality of the results one gets from
 610 standard library implementations of erf().  (One of the most common
 611  is based on ALGORITHM AS66, APPL. STATIST.(1973) Vol.22, .424 by HILL,
 612  which gives only 6-8 digit accuracy).
 613
 614 Here is a listing of my method:
 615
 616 /*
 617 Marsaglia Complementary Normal Distribution Function
 618    cPhi(x) = integral from x to infinity of exp(-.5*t^2)/sqrt(2*pi), x<15
 619  15-digit accuracy for x<15, returns 0 for x>15.
 620 #include <math.h>
 621 */
 622
 623 double cPhi(double x){
 624 long double v[]={0.,.65567954241879847154L,
 625 .42136922928805447322L,.30459029871010329573L,
 626 .23665238291356067062L,.19280810471531576488L,
 627 .16237766089686746182L,.14010418345305024160L,
 628 .12313196325793229628L,.10978728257830829123L,
 629 .99028596471731921395e-1L,.90175675501064682280e-1L,
 630 .82766286501369177252e-1L,.76475761016248502993e-1L,
 631 .71069580538852107091e-1L,.66374235823250173591e-1L};
 632 long double h,a,b,z,t,sum,pwr;
 633 int i,j;
 634       if(x>15.) return (0.);
 635       if(x<-15.) return (1.);
 636         j=fabs(x)+1.;
 637         z=j;
 638         h=fabs(x)-z;
 639         a=v[j];
 640         b=z*a-1.;
 641         pwr=1.;
 642         sum=a+h*b;
 643            for(i=2;i<60;i+=2){
 644            a=(a+z*b)/i;
 645            b=(b+z*a)/(i+1);
 646            pwr=pwr*h*h;
 647            t=sum;
 648            sum=sum+pwr*(a+h*b);
 649            if(sum==t) break; }
 650       sum=sum*exp(-.5*x*x-.91893853320467274178L);
 651       if(x<0.) sum=1.-sum;
 652       return ((double) sum);
 653                       }
 654 */
 655  end of listing
 656 */
 657
 658 The method is based on defining phi(x)=exp(-x^2)/sqrt(2pi) and
 659
 660        R(x)=cPhi(x)/phi(x).
 661
 662 The function R(x) is well-behaved  and terms of its Taylor
 663 series are readily obtained by a two-term recursion.   With an accurate
 664 representation of R(x) at ,say, x=0,1,2,...,15, a simple evaluation
 665 of the Taylor series at intermediate points provides up to
 666 15 digits of accuracy.
 667 An article describing the method will be in the new version of
 668 my Diehard CDROM.   A new version of the Diehard tests
 669 of randomness (but not yet the new DVDROM) is at
 670    http://www.csis.hku.hk/~diehard/
 671
 672
 673  George Marsaglia