Revision 601 trunk/src/df1b2-separable/df3gammp.cpp

df3gammp.cpp (revision 601)
1
/*
1
/**
2 2
 * $Id$
3 3
 *
4 4
 * Author: David Fournier
5 5
 * Copyright (c) 2009-2012 ADMB Foundation
6 6
 */
7
/**
8
 * \file
9
 * Description not yet available.
10
 */
11 7
#include <df1b2fun.h>
12
#include <df13fun.h>
13
#include "mconf.h"
8
#define ITMAX 100
9
#define EPS 1.0e-9
10
//#define EPS 3.0e-7
11
#define FPMIN 1.0e-30
12
double get_values(double x,double y,int print_switch);
14 13

  
15
#ifdef NANS
16
#undef NANS
17
#endif
18
#ifdef INFINITIES
19
#undef INFINITIES
20
#endif
21 14

  
22
double get_values(double x, double y, int print_switch);
23
df3_two_variable igam(const df3_two_variable & _a,
24
		      const df3_two_variable & _x);
25
df3_two_variable igamc(const df3_two_variable & _a,
26
		       const df3_two_variable & _x);
27

  
28
namespace Cephes
15
df1b2variable log_negbinomial_density(double x,const df1b2variable& _xmu, 
16
  const df1b2variable& _xtau)
29 17
{
18
  ADUNCONST(df1b2variable,xmu)
19
  ADUNCONST(df1b2variable,xtau)
20
  init_df3_two_variable mu(xmu);
21
  init_df3_two_variable tau(xtau);
22
  *mu.get_u_x()=1.0;
23
  *tau.get_u_y()=1.0;
24
  if (value(tau)-1.0<0.0)
25
  {
26
    cerr << "tau <=1 in log_negbinomial_density " << endl;
27
    ad_exit(1);
28
  }
29
  df3_two_variable r=mu/(1.e-120+(tau-1.0));
30
  df3_two_variable tmp;
31
  tmp=gammln(x+r)-gammln(r) -gammln(x+1)
32
    +r*log(r)+x*log(mu)-(r+x)*log(r+mu);
33
  df1b2variable tmp1;
34
  tmp1=tmp;
35
  return tmp1;
36
}
30 37

  
31
   int sgngam = 0;
32
   extern const double MAXLOG = 200;
33
   extern const double MAXNUM = 1.7976931348623158E+308;
34
   extern const double LOGPI = 1.14472988584940017414;
35
   extern const double big = 4.503599627370496e15;
36
   extern const double biginv = 2.22044604925031308085e-16;
37
   extern const double MACHEP = 2.22045e-16;
38
   extern const double MYINF = 1.7976931348623158E+308;
38
/** Log gamma function.
39
    \param xx \f$x\f$
40
    \return \f$\ln\bigr(\Gamma(x)\bigl)\f$
39 41

  
40
   df3_two_variable polevl(const df3_two_variable & x, const double *_coef, int N);
41
   df3_two_variable p1evl(const df3_two_variable & x, const double *_coef, int N);
42

  
43
/* Coefficents used for the rational function in gamma(x) */
44
#ifdef UNK
45
   extern const double P[] = {
46
      1.60119522476751861407E-4,
47
      1.19135147006586384913E-3,
48
      1.04213797561761569935E-2,
49
      4.76367800457137231464E-2,
50
      2.07448227648435975150E-1,
51
      4.94214826801497100753E-1,
52
      9.99999999999999996796E-1
53
   };
54

  
55
   extern const double Q[] = {
56
      -2.31581873324120129819E-5,
57
      5.39605580493303397842E-4,
58
      -4.45641913851797240494E-3,
59
      1.18139785222060435552E-2,
60
      3.58236398605498653373E-2,
61
      -2.34591795718243348568E-1,
62
      7.14304917030273074085E-2,
63
      1.00000000000000000320E0
64
   };
65

  
66
#define MAXGAM 171.624376956302725
... This diff was truncated because it exceeds the maximum size that can be displayed.

Also available in: Unified diff