Module: Rubystats::SpecialMath

Includes:
NumericalConstants
Included in:
BetaDistribution, BinomialDistribution, ExponentialDistribution, GammaDistribution, LognormalDistribution, NormalDistribution, ProbabilityDistribution, StudentTDistribution, UniformDistribution
Defined in:
lib/rubystats/modules.rb

Overview

Ruby port of SpecialMath.php from PHPMath, which is a port of JSci methods found in SpecialMath.java.

Ruby port by Bryan Donovan bryandonovan.com

Author

Jaco van Kooten

Author

Paul Meagher

Author

Bryan Donovan

Constant Summary

Constants included from NumericalConstants

NumericalConstants::EPS, NumericalConstants::GAMMA, NumericalConstants::GAMMA_X_MAX_VALUE, NumericalConstants::GOLDEN_RATIO, NumericalConstants::LOG_GAMMA_X_MAX_VALUE, NumericalConstants::MAX_FLOAT, NumericalConstants::MAX_ITERATIONS, NumericalConstants::MAX_VALUE, NumericalConstants::PRECISION, NumericalConstants::SQRT2, NumericalConstants::SQRT2PI, NumericalConstants::TWO_PI, NumericalConstants::XMININ

Instance Attribute Summary collapse

Instance Method Summary collapse

Instance Attribute Details

#log_beta_cache_p ⇒ Object (readonly)

Returns the value of attribute log_beta_cache_p.



45
46
47
 
# File 'lib/rubystats/modules.rb', line 45

def log_beta_cache_p
  @log_beta_cache_p
end

#log_beta_cache_q ⇒ Object (readonly)

Returns the value of attribute log_beta_cache_q.



45
46
47
 
# File 'lib/rubystats/modules.rb', line 45

def log_beta_cache_q
  @log_beta_cache_q
end

#log_beta_cache_res ⇒ Object (readonly)

Returns the value of attribute log_beta_cache_res.



45
46
47
 
# File 'lib/rubystats/modules.rb', line 45

def log_beta_cache_res
  @log_beta_cache_res
end

#log_gamma_cache_res ⇒ Object (readonly)

Returns the value of attribute log_gamma_cache_res.



45
46
47
 
# File 'lib/rubystats/modules.rb', line 45

def log_gamma_cache_res
  @log_gamma_cache_res
end

#log_gamma_cache_x ⇒ Object (readonly)

Returns the value of attribute log_gamma_cache_x.



45
46
47
 
# File 'lib/rubystats/modules.rb', line 45

def log_gamma_cache_x
  @log_gamma_cache_x
end

Instance Method Details

#beta(p, q) ⇒ Object

Beta function.

Author

Jaco van Kooten



436
437
438
439
440
441
442
 
# File 'lib/rubystats/modules.rb', line 436

def beta(p, q)
  if p <= 0.0 || q <= 0.0 || (p + q) > LOG_GAMMA_X_MAX_VALUE
    return 0.0
  else
    return Math.exp(log_beta(p, q))
  end
end

#beta_fraction(x, p, q) ⇒ Object

Evaluates of continued fraction part of incomplete beta function. Based on an idea from Numerical Recipes (W.H. Press et al, 1992).

Author

Jaco van Kooten



474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
 
# File 'lib/rubystats/modules.rb', line 474

def beta_fraction(x, p, q)
  c = 1.0
  sum_pq  = p + q
  p_plus  = p + 1.0
  p_minus = p - 1.0
  h = 1.0 - sum_pq * x / p_plus
  if h.abs < XMININ
    h = XMININ
  end
  h     = 1.0 / h
  frac  = h
  m     = 1
  delta = 0.0

  while (m <= MAX_ITERATIONS) && ((delta - 1.0).abs > PRECISION)
    m2 = 2 * m
    # even index for d
    d = m * (q - m) * x / ( (p_minus + m2) * (p + m2))
    h = 1.0 + d * h
    if h.abs < XMININ
      h = XMININ
    end
    h = 1.0 / h
    c = 1.0 + d / c
    if c.abs < XMININ
      c = XMININ
    end
    frac *= (h * c)
    # odd index for d
    d = -(p + m) * (sum_pq + m) * x / ((p + m2) * (p_plus + m2))
    h = 1.0 + d * h
    if h.abs < XMININ
      h = XMININ
    end
    h = 1.0 / h
    c = 1.0 + d / c
    if c.abs < XMININ
      c = XMININ
    end
    delta = h * c
    frac *= delta
    m += 1
  end
  return frac
end

#complementary_error(x) ⇒ Object

Complementary error function. Based on C-code for the error function developed at Sun Microsystems. author Jaco van Kooten



699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
 
# File 'lib/rubystats/modules.rb', line 699

def complementary_error(x)
 # Coefficients for approximation of erfc in [1.25,1/.35]

 eRa = [-9.86494403484714822705e-03,
   -6.93858572707181764372e-01,
   -1.05586262253232909814e01,
   -6.23753324503260060396e01,
   -1.62396669462573470355e02,
   -1.84605092906711035994e02,
   -8.12874355063065934246e01,
   -9.81432934416914548592e00 ]

 eSa = [ 1.96512716674392571292e01,
   1.37657754143519042600e02,
   4.34565877475229228821e02,
   6.45387271733267880336e02,
   4.29008140027567833386e02,
   1.08635005541779435134e02,
   6.57024977031928170135e00,
   -6.04244152148580987438e-02 ]

 # Coefficients for approximation to erfc in [1/.35,28]

 eRb = [-9.86494292470009928597e-03,
   -7.99283237680523006574e-01,
   -1.77579549177547519889e01,
   -1.60636384855821916062e02,
   -6.37566443368389627722e02,
   -1.02509513161107724954e03,
   -4.83519191608651397019e02 ]

 eSb = [ 3.03380607434824582924e01,
   3.25792512996573918826e02,
   1.53672958608443695994e03,
   3.19985821950859553908e03,
   2.55305040643316442583e03,
   4.74528541206955367215e02,
   -2.24409524465858183362e01 ]

 abs_x = (if x >= 0.0 then x else -x end)
 if abs_x < 1.25
   retval = 1.0 - error(abs_x)
 elsif abs_x > 28.0
   retval = 0.0

   # 1.25 < |x| < 28
 else
   s = 1.0/(abs_x * abs_x)
   if abs_x < 2.8571428
     r = eRa[0] + s * (eRa[1] + s *
                       (eRa[2] + s * (eRa[3] + s * (eRa[4] + s *
                                                    (eRa[5] + s * (eRa[6] + s * eRa[7])
                                                    )))))

                                                    s = 1.0 + s * (eSa[0] + s * (eSa[1] + s *
                                       (eSa[2] + s * (eSa[3] + s * (eSa[4] + s *
                        (eSa[5] + s * (eSa[6] + s * eSa[7])))))))

   else
     r = eRb[0] + s * (eRb[1] + s *
                       (eRb[2] + s * (eRb[3] + s * (eRb[4] + s *
          (eRb[5] + s * eRb[6])))))

     s = 1.0 + s * (eSb[0] + s *
                    (eSb[1] + s * (eSb[2] + s * (eSb[3] + s *
          (eSb[4] + s * (eSb[5] + s * eSb[6]))))))
   end
   retval =  Math.exp(-x * x - 0.5625 + r/s) / abs_x
 end
 return ( if x >= 0.0 then retval else 2.0 - retval end )
end

#error(x) ⇒ Object

Error function. Based on C-code for the error function developed at Sun Microsystems.

Author

Jaco van Kooten



629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
 
# File 'lib/rubystats/modules.rb', line 629

def error(x)
  e_efx = 1.28379167095512586316e-01

  ePp = [ 1.28379167095512558561e-01,
    -3.25042107247001499370e-01,
    -2.84817495755985104766e-02,
    -5.77027029648944159157e-03,
    -2.37630166566501626084e-05 ]

  eQq = [ 3.97917223959155352819e-01,
    6.50222499887672944485e-02,
    5.08130628187576562776e-03,
    1.32494738004321644526e-04,
    -3.96022827877536812320e-06 ]

  # Coefficients for approximation to erf in [0.84375,1.25]
  ePa = [-2.36211856075265944077e-03,
    4.14856118683748331666e-01,
    -3.72207876035701323847e-01,
    3.18346619901161753674e-01,
    -1.10894694282396677476e-01,
    3.54783043256182359371e-02,
    -2.16637559486879084300e-03 ]

  eQa = [ 1.06420880400844228286e-01,
    5.40397917702171048937e-01,
    7.18286544141962662868e-02,
    1.26171219808761642112e-01,
    1.36370839120290507362e-02,
    1.19844998467991074170e-02 ]

  e_erx = 8.45062911510467529297e-01

  abs_x = (if x >= 0.0 then x else -x end)
  # 0 < |x| < 0.84375
  if abs_x < 0.84375
    #|x| < 2**-28
    if abs_x < 3.7252902984619141e-9
      retval = abs_x + abs_x * e_efx
    else
      s = x * x
      p = ePp[0] + s * (ePp[1] + s * (ePp[2] + s * (ePp[3] + s * ePp[4])))

      q = 1.0 + s * (eQq[0] + s * (eQq[1] + s *
                                   ( eQq[2] + s * (eQq[3] + s * eQq[4]))))
      retval = abs_x + abs_x * (p / q)
    end
  elsif abs_x < 1.25
    s = abs_x - 1.0
    p = ePa[0] + s * (ePa[1] + s *
                      (ePa[2] + s * (ePa[3] + s *
       (ePa[4] + s * (ePa[5] + s * ePa[6])))))

    q = 1.0 + s * (eQa[0] + s *
                   (eQa[1] + s * (eQa[2] + s *
       (eQa[3] + s * (eQa[4] + s * eQa[5])))))
    retval = e_erx + p / q

  elsif abs_x >= 6.0
    retval = 1.0
  else
    retval = 1.0 - complementary_error(abs_x)
  end
  return (if x >= 0.0 then retval else -retval end)
end

#gamma(_x) ⇒ Object

TODO test this



215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
 
# File 'lib/rubystats/modules.rb', line 215

def gamma(_x)
  return 0 if _x == 0.0

  p0 = 1.000000000190015
  p = {1 => 76.18009172947146,
  2 => -86.50532032941677,
  3 => 24.01409824083091,
  4 => -1.231739572450155,
  5 => 1.208650973866179e-3,
  6 => -5.395239384953e-6}

  y = x = _x
  tmp = x + 5.5
  tmp -= (x + 0.5) * Math.log(tmp)

  summer = p0
  for j in (1 ... 6)
    y += 1
    summer += (p[j] / y)
  end
  return Math.exp(0 - tmp + Math.log(2.5066282746310005 * summer / x))
end

#gamma_fraction(a, x) ⇒ Object

Author

Jaco van Kooten



405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
 
# File 'lib/rubystats/modules.rb', line 405

def gamma_fraction(a, x)
  b  = x + 1.0 - a
  c  = 1.0 / XMININ
  d  = 1.0 / b
  h  = d
  del= 0.0
  an = 0.0
  for i in (1...MAX_ITERATIONS)
    if (del-1.0).abs > PRECISION
      an = -i * (i - a)
      b += 2.0
      d  = an * d + b
      c  = b + an / c
      if c.abs < XMININ
        c = XMININ
        if d.abs < XMININ
          c = XMININ
          d   = 1.0 / d
          del = d * c
          h  *= del
        end
      end
    end
    return Math.exp(-x + a * Math.log(x) - log_gamma(a)) * h
  end
end

#gamma_series_expansion(a, x) ⇒ Object

Author

Jaco van Kooten



389
390
391
392
393
394
395
396
397
398
399
400
401
402
 
# File 'lib/rubystats/modules.rb', line 389

def gamma_series_expansion(a, x)
  ap  = a
  del = 1.0 / a
  sum = del
  (1...MAX_ITERATIONS).each do
    ap += 1
    del *= x / ap
    sum += del
    if del < sum * PRECISION
      return sum * Math.exp(-x + a * Math.log(x) - log_gamma(a))
    end
  end
  return "Maximum iterations exceeded: please file a bug report."
end

#incomplete_beta(x, p, q) ⇒ Object

Incomplete Beta function.

Author

Jaco van Kooten

Author

Paul Meagher

The computation is based on formulas from Numerical Recipes,
Chapter 6.4 (W.H. Press et al, 1992).



452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
 
# File 'lib/rubystats/modules.rb', line 452

def incomplete_beta(x, p, q)
  if x <= 0.0
    return 0.0
  elsif x >= 1.0
    return 1.0
  elsif (p <= 0.0) || (q <= 0.0) || (p + q) > LOG_GAMMA_X_MAX_VALUE
    return 0.0
  else
    beta_gam = Math.exp( -log_beta(p, q) + p * Math.log(x) + q * Math.log(1.0 - x) )
    if x < (p + 1.0) / (p + q + 2.0)
      return beta_gam * beta_fraction(x, p, q) / p
    else
      return 1.0 - (beta_gam * beta_fraction(1.0 - x, q, p) / q)
    end
  end
end

#incomplete_gamma(a, x) ⇒ Object

Incomplete Gamma function. The computation is based on approximations presented in Numerical Recipes, Chapter 6.2 (W.H. Press et al, 1992).

Parameters:

  • a —

    require a>=0

  • x —

    require x>=0

Returns:

  • 0 if x<0, a<=0 or a>2.55E305 to avoid errors and over/underflow

Author:

  • Jaco van Kooten



378
379
380
381
382
383
384
385
386
 
# File 'lib/rubystats/modules.rb', line 378

def incomplete_gamma(a, x)
  if x <= 0.0 || a <= 0.0 || a > LOG_GAMMA_X_MAX_VALUE
    return 0.0
  elsif x < (a + 1.0)
    return gamma_series_expansion(a, x)
  else
    return 1.0-gamma_fraction(a, x)
  end
end

#log_beta(p, q) ⇒ Object 



54
55
56
57
58
59
60
61
62
63
64
65
 
# File 'lib/rubystats/modules.rb', line 54

def log_beta(p,q)
  if p != log_beta_cache_p || q != log_beta_cache_q
    @log_beta_cache_p = p
    @log_beta_cache_q = q
    if (p <= 0.0) || (q <= 0.0) || (p + q) > LOG_GAMMA_X_MAX_VALUE
      @log_beta_cache_res = 0.0
    else
      @log_beta_cache_res = log_gamma(p) + log_gamma(q) - log_gamma(p + q)
    end
  end
  return log_beta_cache_res
end

#log_gamma(x) ⇒ Object 



238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
 
# File 'lib/rubystats/modules.rb', line 238

def log_gamma(x)
  lg_d1 = -0.5772156649015328605195174
  lg_d2 = 0.4227843350984671393993777
  lg_d4 = 1.791759469228055000094023

  lg_p1 = [ 4.945235359296727046734888,
    201.8112620856775083915565, 2290.838373831346393026739,
    11319.67205903380828685045, 28557.24635671635335736389,
    38484.96228443793359990269, 26377.48787624195437963534,
    7225.813979700288197698961 ]

  lg_p2 = [ 4.974607845568932035012064,
    542.4138599891070494101986, 15506.93864978364947665077,
    184793.2904445632425417223, 1088204.76946882876749847,
    3338152.967987029735917223, 5106661.678927352456275255,
    3074109.054850539556250927 ]

  lg_p4 = [ 14745.02166059939948905062,
    2426813.369486704502836312, 121475557.4045093227939592,
    2663432449.630976949898078, 29403789566.34553899906876,
    170266573776.5398868392998, 492612579337.743088758812,
    560625185622.3951465078242 ]

  lg_q1 = [ 67.48212550303777196073036,
    1113.332393857199323513008, 7738.757056935398733233834,
    27639.87074403340708898585, 54993.10206226157329794414,
    61611.22180066002127833352, 36351.27591501940507276287,
    8785.536302431013170870835 ]

  lg_q2 = [ 183.0328399370592604055942,
    7765.049321445005871323047, 133190.3827966074194402448,
    1136705.821321969608938755, 5267964.117437946917577538,
    13467014.54311101692290052, 17827365.30353274213975932,
    9533095.591844353613395747 ]

  lg_q4 = [ 2690.530175870899333379843,
    639388.5654300092398984238, 41355999.30241388052042842,
    1120872109.61614794137657, 14886137286.78813811542398,
    101680358627.2438228077304, 341747634550.7377132798597,
    446315818741.9713286462081 ]

  lg_c  = [ -0.001910444077728,8.4171387781295e-4,
    -5.952379913043012e-4, 7.93650793500350248e-4,
    -0.002777777777777681622553, 0.08333333333333333331554247,
    0.0057083835261 ]

  # Rough estimate of the fourth root of logGamma_xBig
  lg_frtbig = 2.25e76
  pnt68     = 0.6796875

  if x == log_gamma_cache_x
    return log_gamma_cache_res
  end

  y = x
  if y > 0.0 && y <= LOG_GAMMA_X_MAX_VALUE
    if y <= EPS
      res = -Math.log(y)
    elsif y <= 1.5
      # EPS .LT. X .LE. 1.5
      if y < pnt68
        corr = -Math.log(y)
        # xm1 is x-m-one, not x-m-L
        xm1 = y
      else
        corr = 0.0
        xm1 = y - 1.0
      end
      if y <= 0.5 || y >= pnt68
        xden = 1.0
        xnum = 0.0
        for i in (0...8)
          xnum = xnum * xm1 + lg_p1[i]
          xden = xden * xm1 + lg_q1[i]
        end
        res = corr + xm1 * (lg_d1 + xm1 * (xnum / xden))
      else
        xm2 = y - 1.0
        xden = 1.0
        xnum = 0.0
        for i in (0 ... 8)
          xnum = xnum * xm2 + lg_p2[i]
          xden = xden * xm2 + lg_q2[i]
        end
        res = corr + xm2 * (lg_d2 + xm2 * (xnum / xden))
      end
    elsif y <= 4.0
      # 1.5 .LT. X .LE. 4.0
      xm2 = y - 2.0
      xden = 1.0
      xnum = 0.0
      for i in (0 ... 8)
        xnum = xnum * xm2 + lg_p2[i]
        xden = xden * xm2 + lg_q2[i]
      end
      res = xm2 * (lg_d2 + xm2 * (xnum / xden))
    elsif y <= 12.0
      # 4.0 .LT. X .LE. 12.0
      xm4 = y - 4.0
      xden = -1.0
      xnum = 0.0
      for i in (0 ... 8)
        xnum = xnum * xm4 + lg_p4[i]
        xden = xden * xm4 + lg_q4[i]
      end
      res = lg_d4 + xm4 * (xnum / xden)
    else
      # Evaluate for argument .GE. 12.0
      res = 0.0
      if y <= lg_frtbig
        res = lg_c[6]
        ysq = y * y
        for i in (0...6)
          res = res / ysq + lg_c[i]
        end
      end
      res = res/y
      corr = Math.log(y)
      res = res + Math.log(SQRT2PI) - 0.5 * corr
      res = res + y * (corr - 1.0)
    end
  else
    #return for bad arguments
    res = MAX_VALUE
  end
  # final adjustments and return
  @log_gamma_cache_x = x
  @log_gamma_cache_res = res
  return res
end

#orig_gamma(x) ⇒ Object

Gamma function. Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz<BR> Applied Mathematics Division<BR> Argonne National Laboratory<BR> Argonne, IL 60439<BR> <P> References: <OL> <LI>“An Overview of Software Development for Special Functions”, W. J. Cody, Lecture Notes in Mathematics, 506, Numerical Analysis Dundee, 1975, G. A. Watson (ed.), Springer Verlag, Berlin, 1976. <LI>Computer Approximations, Hart, Et. Al., Wiley and sons, New York, 1968. </OL></P><P> From the original documentation: </P><P> This routine calculates the Gamma function for a real argument X. Computation is based on an algorithm outlined in reference 1. The program uses rational functions that approximate the Gamma function to at least 20 significant decimal digits. Coefficients for the approximation over the interval (1,2) are unpublished. Those for the approximation for X .GE. 12 are from reference 2. The accuracy achieved depends on the arithmetic system, the compiler, the intrinsic functions, and proper selection of the machine-dependent constants. </P><P> Error returns:<BR> The program returns the value XINF for singularities or when overflow would occur. The computation is believed to be free of underflow and overflow. </P>

Author

Jaco van Kooten



96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
 
# File 'lib/rubystats/modules.rb', line 96

def orig_gamma(x)
  # Gamma related constants
  g_p = [ -1.71618513886549492533811, 24.7656508055759199108314,
    -379.804256470945635097577, 629.331155312818442661052,
    866.966202790413211295064, -31451.2729688483675254357,
    -36144.4134186911729807069, 66456.1438202405440627855 ]
  g_q = [-30.8402300119738975254353, 315.350626979604161529144,
    -1015.15636749021914166146, -3107.77167157231109440444,
    22538.1184209801510330112, 4755.84627752788110767815,
    -134659.959864969306392456, -115132.259675553483497211 ]
  g_c = [-0.001910444077728, 8.4171387781295e-4, -5.952379913043012e-4,
    7.93650793500350248e-4, -0.002777777777777681622553,
    0.08333333333333333331554247, 0.0057083835261 ]
  fact = 1.0
  i = 0
  n = 0
  y = x
  parity = false
  if y <= 0.0
    # ----------------------------------------------------------------------
    #  Argument is negative
    # ----------------------------------------------------------------------
    y = -(x)
    y1 = y.to_i
    res = y - y1
    if res != 0.0
      if y1 != (((y1*0.5).to_i) * 2.0)
        parity = true
        fact = -M_pi/sin(M_pi * res)
        y += 1
      end
    else
      return MAX_VALUE
    end
  end

  # ----------------------------------------------------------------------
  #  Argument is positive
  # ----------------------------------------------------------------------
  if y < EPS
    # ----------------------------------------------------------------------
    #  Argument .LT. EPS
    # ----------------------------------------------------------------------
    if y >= XMININ
      res = 1.0 / y
    else
      return MAX_VALUE
    end
  elsif y < 12.0
    y1 = y
    #end
    if y < 1.0
      # ----------------------------------------------------------------------
      #  0.0 .LT. argument .LT. 1.0
      # ----------------------------------------------------------------------
      z = y
      y += 1
    else
      # ----------------------------------------------------------------------
      #  1.0 .LT. argument .LT. 12.0, reduce argument if necessary
      # ----------------------------------------------------------------------
      n = y.to_i - 1
      y -= n.to_f
      z = y - 1.0
    end
    # ----------------------------------------------------------------------
    #  Evaluate approximation for 1.0 .LT. argument .LT. 2.0
    # ----------------------------------------------------------------------
    xnum = 0.0
    xden = 1.0
    for i in (0...8)
      xnum = (xnum + g_p[i]) * z
      xden = xden * z + g_q[i]
    end
    res = xnum / xden + 1.0
    if y1 < y
      # ----------------------------------------------------------------------
      #  Adjust result for case  0.0 .LT. argument .LT. 1.0
      # ----------------------------------------------------------------------
      res /= y1
    elsif y1 > y
      # ----------------------------------------------------------------------
      #  Adjust result for case  2.0 .LT. argument .LT. 12.0
      # ----------------------------------------------------------------------
      for i in (0...n)
        res *= y
        y += 1
      end
    end
  else
    # ----------------------------------------------------------------------
    #  Evaluate for argument .GE. 12.0
    # ----------------------------------------------------------------------
    if y <= GAMMA_X_MAX_VALUE
      ysq = y * y
      sum = g_c[6]
      for i in(0...6)
        sum = sum / ysq + g_c[i]
        sum = sum / y - y + log(SQRT2PI)
        sum += (y - 0.5) * log(y)
        res = Math.exp(sum)
      end
    else
      return MAX_VALUE
    end
    # ----------------------------------------------------------------------
    #  Final adjustments and return
    # ----------------------------------------------------------------------
    if parity
      res = -res
      if fact != 1.0
        res = fact / res
        return res
      end
    end
  end
end