block({
  include(`largenum-prg.txt'):
  let(gammapi,4*atan(1)):

  func(gammafactorial,{gamma([x]+1)}):
  func(gammarecurrenceerror,{
    abs(gamma([x]+1)-([x]*gamma([x])))
  }):
  func(gammareflectionerror,{
    abs((gamma([x])*gamma(1-[x]))-([gammapi]/sin([gammapi]*[x])))
  }):
  func(betagamma,{
    (gamma([x])*gamma([y]))/gamma([x]+[y])
  }):
  func(gammachoose,{
    gamma([x]+1)/(gamma([y]+1)*gamma(([x]-[y])+1))
  }):
  func(gammarising,{
    gamma([x]+[y])/gamma([x])
  }):
  func(gammapdf,{
    (([x]^([y]-1))*eul(-([x]/[z])))/(gamma([y])*([z]^[y]))
  }):

  print_`Gamma function: real arguments, binary64 arithmetic\n\n':

  print_`1. Direct evaluation and vector application\n':
  print_`gamma(1)     = ' gamma(1) `\n':
  print_`gamma(6)     = ' gamma(6) `  = 5!\n':
  let(gammainputs,<0.25,0.5,1,1.5,2,2.5,5,10>):
  print_`x            = ' [gammainputs] `\n':
  print_`gamma(x)     = ' apply(gamma,[gammainputs]) `\n\n':

  print_`2. Factorials: n! = gamma(n+1)\n':
  print_`0!           = ' gammafactorial(0) `\n':
  print_`1!           = ' gammafactorial(1) `\n':
  print_`5!           = ' gammafactorial(5) `\n':
  print_`10!          = ' gammafactorial(10) `\n':
  print_`0.5!         = ' gammafactorial(0.5) `\n':
  print_`2.5!         = ' gammafactorial(2.5) `\n':
  print_`30! binary64 = ' gammafactorial(30) `\n':
  print_`30! exact    = ' largefactorial(30) `\n\n':

  print_`3. Half-integer values\n':
  print_`gamma(0.5)   = ' gamma(0.5) `\n':
  print_`sqrt(pi)     = ' sqrt([gammapi]) `\n':
  print_`gamma(1.5)   = ' gamma(1.5) `\n':
  print_`sqrt(pi)/2   = ' (sqrt([gammapi])/2) `\n':
  print_`gamma(2.5)   = ' gamma(2.5) `\n':
  print_`3*sqrt(pi)/4 = ' ((3*sqrt([gammapi]))/4) `\n\n':

  print_`4. Negative non-integers and values near zero\n':
  print_`quote(...) preserves signed numbers in labelled output.\n':
  print_`gamma(-0.5)  = ' quote(gamma(-0.5)) `\n':
  print_`-2*sqrt(pi)  = ' quote(-2*sqrt([gammapi])) `\n':
  print_`gamma(-1.5)  = ' quote(gamma(-1.5)) `\n':
  print_`gamma(-2.5)  = ' quote(gamma(-2.5)) `\n':
  print_`gamma(-0.001)= ' quote(gamma(-0.001)) `\n':
  print_`gamma(0.001) = ' gamma(0.001) `\n\n':

  print_`5. Recurrence: gamma(x+1) = x*gamma(x)\n':
  print_`Absolute residuals at <-1.5,-0.5,0.25,1.5,5.25>:\n':
  print_apply(gammarecurrenceerror,<-1.5,-0.5,0.25,1.5,5.25>) `\n\n':

  print_`6. Reflection: gamma(x)*gamma(1-x) = pi/sin(pi*x)\n':
  print_`Absolute residuals at <-1.5,-0.5,0.25,0.75,1.25,3.5>:\n':
  print_apply(gammareflectionerror,<-1.5,-0.5,0.25,0.75,1.25,3.5>) `\n\n':

  print_`7. Beta function: B(a,b) = gamma(a)*gamma(b)/gamma(a+b)\n':
  print_`For positive a and b.\n':
  print_`B(2,3)       = ' betagamma(2,3) `  = 1/12\n':
  print_`B(0.5,0.5)   = ' betagamma(0.5,0.5) `  = pi\n\n':

  print_`8. Binomial coefficients and rising factorials\n':
  print_`C(n,k) = gamma(n+1)/(gamma(k+1)*gamma(n-k+1))\n':
  print_`C(10,3)      = ' gammachoose(10,3) `\n':
  print_`C(2.5,2)     = ' gammachoose(2.5,2) `\n':
  print_`Rising factorial (a)_m = gamma(a+m)/gamma(a), a>0\n':
  print_`(0.5)_4      = ' gammarising(0.5,4) `  = 0.5*1.5*2.5*3.5\n\n':

  print_`9. Gamma probability density\n':
  print_`f(x;k,theta) = x^(k-1)*exp(-x/theta)/(gamma(k)*theta^k)\n':
  print_`Here x>0, shape k>0, scale theta>0; eul is exp in sxript.\n':
  print_`f(1;0.5,2)   = ' gammapdf(1,0.5,2) `\n':
  print_`f(2;2.5,1.5) = ' gammapdf(2,2.5,1.5) `\n\n':

  print_`Domain and numerical limits\n':
  print_`gamma(0), gamma(-1), gamma(-2), ... raise a domain error.\n':
  print_`Gamma results and the derived functions use ordinary floating point.\n':
  print_`Small identity residuals are floating-point rounding errors.\n':
  print_`Use largefactorial for exact nonnegative integer factorials.\n':
  print_`Direct gamma ratios can overflow at large arguments.'
})