Merge remote-tracking branch 'origin/master' into natalia/qlpt_feature

CMakeLists.txt

@@ -54,7 +54,7 @@ IF(YORICK_SUPPORT)
   ENDIF((EXISTS ${NETCDFCPP_INCLUDE_PATH}/netcdf AND ${NETCDFCPP_LIBRARY} MATCHES "netcdf_c\\+\\+4") OR (INTERNAL_NETCDF))
 ENDIF(YORICK_SUPPORT)
 
-find_program(CYTHON cython3 cython)
+find_program(CYTHON NAMES cython3 cython)
 find_library(ZLIB_LIBRARY z)
 find_library(DL_LIBRARY dl)
 find_library(MATH_LIBRARY m)

builder/build-wheels.sh

@@ -1,6 +1,11 @@
 #!/bin/bash
 set -e -x
 
+CC=cc
+CXX=c++
+
+export CC CXX
+
 # Install a system package required by our library
 #yum install -y atlas-devel
 yum install -y cmake3 gsl-devel zlib-devel
@@ -8,7 +13,7 @@ yum install -y cmake3 gsl-devel zlib-devel
 ln -fs /usr/bin/cmake3 /usr/bin/cmake
 
 
-ALL_PYTHON="cp37-cp37m cp38-cp38" 
+ALL_PYTHON="cp36-cp36m cp37-cp37m cp38-cp38" 
 
 # Compile wheels
 for pkg in $ALL_PYTHON; do

builder/start.sh

@@ -9,4 +9,4 @@ if ! [ -e ${d}/setup.py ] ; then
   exit 1
 fi
 
-podman run -ti  --rm -e PLAT=manylinux2014_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2014_x86_64 /io/builder/build-wheels.sh
+podman run -ti  --rm -e PLAT=manylinux2010_x86_64 -v ${d}:/io:Z quay.io/pypa/manylinux2010_x86_64 /io/builder/build-wheels.sh

conda/build.sh

@@ -0,0 +1 @@
+$PYTHON setup.py install

conda/conda_build_config.yaml

@@ -0,0 +1,7 @@
+python:
+  - 3.7
+  - 3.8
+
+numpy:
+  - 1.11
+  - 1.19

conda/meta.yaml

@@ -0,0 +1,38 @@
+package:
+  name: cosmotool
+  version: "1.0.0a7"
+
+source:
+  git_rev: 7fce73e
+  git_url: https://bitbucket.org/glavaux/cosmotool
+
+requirements:
+  build:
+    - numpy >=1.11
+    - gcc_linux-64
+    - gxx_linux-64
+    - python
+    - setuptools
+    - cython
+    - healpy
+    - numexpr
+    - cffi
+    - pyfftw
+    - gsl
+    - h5py
+
+  run:
+    - numpy
+    - python
+    - healpy
+    - numexpr
+    - cffi
+    - pyfftw
+    - h5py
+
+test:
+  imports:
+    - cosmotool
+
+about:
+  home: https://bitbucket.org/glavaux/cosmotool

external/external_build.cmake

@@ -3,7 +3,7 @@ include(FindOpenMP)
 OPTION(ENABLE_OPENMP "Set to Yes if Healpix and/or you need openMP" OFF)
 
 SET(FFTW_URL "http://www.fftw.org/fftw-3.3.3.tar.gz" CACHE STRING "URL to download FFTW from")
-SET(EIGEN_URL "http://bitbucket.org/eigen/eigen/get/3.2.10.tar.gz" CACHE STRING "URL to download Eigen from")
+SET(EIGEN_URL "https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2" CACHE STRING "URL to download Eigen from")
 SET(GENGETOPT_URL "ftp://ftp.gnu.org/gnu/gengetopt/gengetopt-2.22.5.tar.gz" CACHE STRING "URL to download gengetopt from")
 SET(HDF5_URL "https://support.hdfgroup.org/ftp/HDF5/releases/hdf5-1.8/hdf5-1.8.18/src/hdf5-1.8.18.tar.bz2" CACHE STRING "URL to download HDF5 from")
 SET(NETCDF_URL "ftp://ftp.unidata.ucar.edu/pub/netcdf/netcdf-4.5.0.tar.gz" CACHE STRING "URL to download NetCDF from")
@@ -255,8 +255,8 @@ ELSE (INTERNAL_BOOST)
        set(Boost_VERSION ${TMP_MAJOR}.${TMP_MINOR}.${TMP_PATCHLEVEL})
     ENDIF()
     if (${Boost_VERSION} VERSION_GREATER_EQUAL 1.70.0)
-      set(Boost_DEP Boost::headers)
+      set(Boost_DEP Boost::boost)
-      set(Boost_TARGET Boost::headers)
+      set(Boost_TARGET Boost::boost)
     endif()
   endif()
   if (NOT Boost_FOUND)
@@ -291,13 +291,19 @@ IF(INTERNAL_GSL)
   set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
   SET(cosmotool_DEPS ${cosmotool_DEPS} gsl)
 ELSE(INTERNAL_GSL)
-  find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h)
+  IF (NOT DEFINED GSL_LIBRARY OR NOT GSL_LIBRARY)
-  find_library(GSL_LIBRARY gsl)
+    find_library(GSL_LIBRARY gsl)
-  find_library(GSLCBLAS_LIBRARY gslcblas)
+  ENDIF()
+  IF (NOT DEFINED GSL_INCLUDE_PATH OR NOT GSL_INCLUDE_PATH)
+    find_path(GSL_INCLUDE_PATH NAMES gsl/gsl_blas.h)
+  ENDIF()
+  IF (NOT DEFINED GSL_CBLAS_LIBRARY OR NOT GSL_CBLAS_LIBRARY)
+    find_library(GSLCBLAS_LIBRARY gslcblas)
+  ENDIF()
 
   set(GSL_LIBRARIES ${GSL_LIBRARY} ${GSLCBLAS_LIBRARY})
-
 ENDIF(INTERNAL_GSL)
+
 mark_as_advanced(GSL_LIBRARY GSLCBLAS_LIBRARY GSL_INCLUDE_PATH)
 
 
@@ -350,7 +356,7 @@ ENDIF(INTERNAL_FFTW)
 IF (INTERNAL_EIGEN)
   ExternalProject_Add(eigen
      URL ${EIGEN_URL}
-     URL_HASH SHA256=04f8a4fa4afedaae721c1a1c756afeea20d3cdef0ce3293982cf1c518f178502 
+     URL_HASH MD5=b9e98a200d2455f06db9c661c5610496
      PREFIX ${BUILD_PREFIX}/eigen-prefix
      CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${EXT_INSTALL}
        -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}

python/CMakeLists.txt

@@ -30,6 +30,10 @@ IF(CYTHON)
     COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx
     DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_project.pyx ${CMAKE_CURRENT_SOURCE_DIR}/project_tool.hpp ) 
 
+  add_custom_command(
+    OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp
+    COMMAND ${CYTHON} --cplus -o ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx
+    DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/_cosmomath.pyx ) 
 ENDIF(CYTHON)
 
 
@@ -38,12 +42,13 @@ add_library(_cosmo_power MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_power.cpp)
 add_library(_cosmo_cic MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmo_cic.cpp)
 add_library(_fast_interp MODULE ${CMAKE_CURRENT_BINARY_DIR}/_fast_interp.cpp)
 add_library(_project MODULE ${CMAKE_CURRENT_BINARY_DIR}/_project.cpp)
+add_library(_cosmomath MODULE ${CMAKE_CURRENT_BINARY_DIR}/_cosmomath.cpp)
 target_include_directories(_cosmotool PRIVATE ${PYTHON_INCLUDES}) 
 target_include_directories(_cosmo_power PRIVATE ${PYTHON_INCLUDES}) 
 target_include_directories(_cosmo_cic PRIVATE ${PYTHON_INCLUDES}) 
 target_include_directories(_fast_interp PRIVATE ${PYTHON_INCLUDES}) 
 target_include_directories(_project PRIVATE ${PYTHON_INCLUDES}) 
-
+target_include_directories(_cosmomath PRIVATE ${PYTHON_INCLUDES})
 
 SET(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Bsymbolic-functions")
 if(APPLE)
@@ -53,10 +58,11 @@ endif()
 target_link_libraries(_cosmotool PRIVATE ${CosmoTool_local}  ${GSL_LIBRARIES})
 target_link_libraries(_cosmo_power PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})
 target_link_libraries(_cosmo_cic PRIVATE ${CosmoTool_local}  ${GSL_LIBRARIES})
+target_link_libraries(_cosmomath PRIVATE ${CosmoTool_local} ${GSL_LIBRARIES})
 target_link_libraries(_project )
 target_link_libraries(_fast_interp PRIVATE ${CosmoTool_local} )
 
-SET(ct_TARGETS   _cosmotool _project _cosmo_power _cosmo_cic _fast_interp )
+SET(ct_TARGETS   _cosmotool _project _cosmo_power _cosmo_cic _fast_interp _cosmomath)
 
 if (Boost_FOUND)
   message(STATUS "Building bispectrum support (path = ${Boost_INCLUDE_DIRS})")

python/_cosmomath.pyx

@@ -0,0 +1,38 @@
+import numpy as np
+cimport numpy as np
+
+np.import_array()
+np.import_ufunc()
+
+cdef extern from "sys/types.h":
+   ctypedef np.int64_t int64_t
+
+cdef extern from "special_math.hpp" namespace "CosmoTool":
+   T log_modified_bessel_first_kind[T](T v, T z) nogil except +
+
+
+cdef np.PyUFuncGenericFunction loop_func[1]
+cdef char input_output_types[3]
+cdef void *elementwise_funcs[1]
+
+loop_func[0] = np.PyUFunc_dd_d
+
+input_output_types[0] = np.NPY_DOUBLE
+input_output_types[1] = np.NPY_DOUBLE
+input_output_types[2] = np.NPY_DOUBLE
+
+elementwise_funcs[0] = <void*>log_modified_bessel_first_kind[double]
+
+log_modified_I = np.PyUFunc_FromFuncAndData(
+    loop_func,
+    elementwise_funcs,
+    input_output_types,
+    1, # number of supported input types
+    2, # number of input args
+    1, # number of output args
+    0, # `identity` element, never mind this
+    "log_modified_bessel_first_kind", # function name
+    "log_modified_bessel_first_kind(v: Float, z: Float) -> Float", # docstring
+    0 # unused
+    )
+

python/_cosmotool.pyx

@@ -513,7 +513,8 @@ def writeGadget(str filename, object simulation):
   simdata.TotalNumPart = NumPart
   simdata.NumPart = NumPart
 
-  cxx_writeGadget(filename, &simdata)
+  filename_b = bytes(filename, 'utf-8')
+  cxx_writeGadget(filename_b, &simdata)
 
 def loadRamses(str basepath, int snapshot_id, int cpu_id, bool doublePrecision = False, bool loadPosition = True, bool loadVelocity = False, bool loadId = False, bool loadMass = False):
   """ loadRamses(basepath, snapshot_id, cpu_id, doublePrecision = False, loadPosition = True, loadVelocity = False)

python/cosmotool/__init__.py

@@ -3,6 +3,7 @@ from ._project import *
 from ._cosmo_power import *
 from ._cosmo_cic import *
 from ._fast_interp import *
+from ._cosmomath import *
 from .grafic import writeGrafic, writeWhitePhase, readGrafic, readWhitePhase
 from .borg import read_borg_vol
 from .cic import cicParticles

python/cosmotool/simu.py

@@ -98,12 +98,12 @@ def simpleWriteGadget(filename, positions, boxsize=1.0, Hubble=100, Omega_M=0.30
   
   s.positions = positions
   
-  if velocities:
+  if velocities is not None:
     s.velocities = velocities
   else:
     s.velocities = [np.zeros(positions[0].size,dtype=np.float32)]*3
     
-  if identifiers:
+  if identifiers is not None:
     s.identifiers = identifiers
   else:
     s.identifiers = np.arange(positions[0].size, dtype=np.int64)

sample/CMakeLists.txt

@@ -114,6 +114,9 @@ if (Boost_FOUND)
       add_dependencies(graficToDensity ${all_deps})
     endif()
   endif()
+
+  add_executable(test_special test_special.cpp)
+  target_link_libraries(test_special ${tolink})
 endif (Boost_FOUND)
 
 IF (ENABLE_OPENMP AND YORICK_SUPPORT)

sample/test_special.cpp

@@ -0,0 +1,12 @@
+#include <iostream>
+#include "special_math.hpp"
+
+int main()
+{
+
+    std::cout << CosmoTool::log1p_exp(2.0) << std::endl;
+
+    std::cout << CosmoTool::log_modified_bessel_first_kind(100.0, 0.1) << std::endl;
+    return 0;
+}
+    

setup.py

@@ -156,9 +156,8 @@ class BuildCMakeExt(build_ext):
 
         cython_code = os.path.join(str(build_dir.absolute()),'mycython')
         with open(cython_code, mode="wt") as ff:
-          ff.write(f"#!{sys.executable}\n"
+          ff.write("#!/bin/sh\n"
-                    "from Cython.Compiler.Main import setuptools_main\n"
+                   f"{sys.executable} -c 'from Cython.Compiler.Main import setuptools_main; setuptools_main()' $@")
-                    "setuptools_main()")
         os.chmod(cython_code, stat.S_IXUSR|stat.S_IWUSR|stat.S_IRUSR|stat.S_IRGRP)
 
         # Now that the necessary directories are created, build
@@ -169,6 +168,7 @@ class BuildCMakeExt(build_ext):
         # Below is just an example set of arguments for building Blender as a Python module
 
         self.spawn(['cmake', '-H'+SOURCE_DIR, '-B'+self.build_temp,
+                    '-DCMAKE_C_COMPILER=' + os.environ["CC"], "-DCMAKE_CXX_COMPILER=" + os.environ["CXX"],
                     '-DENABLE_OPENMP=ON','-DINTERNAL_BOOST=ON','-DINTERNAL_EIGEN=ON',
                     '-DINTERNAL_HDF5=ON','-DINTERNAL_NETCDF=ON',
                     '-DBUILD_PYTHON=ON', '-DINSTALL_PYTHON_LOCAL=OFF',
@@ -218,7 +218,7 @@ class BuildCMakeExt(build_ext):
 CosmoTool_extension = CMakeExtension(name="cosmotool")
 
 setup(name='cosmotool',
-      version='1.0.0a5',
+      version='1.1.0',
       packages=["cosmotool"],
       package_dir={'cosmotool': 'python/cosmotool'},
       install_requires=['numpy','cffi','numexpr','pyfftw','h5py'],

src/special_math.hpp

@@ -0,0 +1,231 @@
+// Original code derived from Boost and is distributed here
+// under the Boost license (licenses/boost-license.txt)
+//    Copyright (c) 2006 Xiaogang Zhang
+//    Copyright (c) 2007, 2017 John Maddock
+// Secondary code copyright by its author and is distributed here
+// under the BSD-3 license (LICENSE.md). Derived from
+// stan/math/prim/fun/log_modified_bessel_first_kind.hpp
+
+#ifndef __COSMOTOOL_SPECIAL_MATH_HPP
+#define __COSMOTOOL_SPECIAL_MATH_HPP
+
+#include "algo.hpp"
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/tools/rational.hpp>
+#include <cmath>
+#include <limits>
+
+// Taken and adapted from
+// https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/log_modified_bessel_first_kind.hpp
+
+namespace CosmoTool {
+
+template <typename T> T log1p_exp(T x) {
+  if (x > T(0)) {
+    return x + std::log1p(std::exp(-x));
+  }
+  return std::log1p(std::exp(x));
+}
+
+template <typename T> T multiply_log(T a, T b) {
+  if (a == 0 && b == 0)
+    return 0;
+  return a * std::log(b);
+}
+
+template <typename T> T inf() { return std::numeric_limits<T>::infinity(); }
+
+template <typename T> T log_sum_exp(T const a, T const b) {
+  if (a == -inf<T>()) {
+    return b;
+  }
+  if (a == inf<T>() && b == inf<T>()) {
+    return inf<T>();
+  }
+  if (a > b) {
+    return a + log1p_exp(b - a);
+  }
+  return b + log1p_exp(a - b);
+}
+
+/* Log of the modified Bessel function of the first kind,
+ * which is better known as the Bessel I function. See
+ * modified_bessel_first_kind.hpp for the function definition.
+ * The derivatives are known to be incorrect for v = 0 because a
+ * simple constant 0 is returned.
+ *
+ * @tparam T common type for calculation
+ * @param v Order, can be a non-integer but must be at least -1
+ * @param z Real non-negative number
+ * @throws std::domain_error if either v or z is NaN, z is
+ * negative, or v is less than -1
+ * @return log of Bessel I function
+ */
+template <typename T> T log_modified_bessel_first_kind(T const v, T const z) {
+  using boost::math::tools::evaluate_polynomial;
+  using std::log;
+  using std::pow;
+  using std::sqrt;
+  static const double LOG_TWO = std::log(2.0);
+  static const double EPSILON = std::numeric_limits<double>::epsilon();
+  static const double TWO_PI = 2.0 * boost::math::constants::pi<double>();
+
+  if (z == 0) {
+    if (v == 0) {
+      return 0.0;
+    }
+    if (v > 0) {
+      return -std::numeric_limits<T>::infinity();
+    }
+    return std::numeric_limits<T>::infinity();
+  }
+  if (std::isinf(z)) {
+    return z;
+  }
+  if (v == 0) {
+    // Modified Bessel function of the first kind of order zero
+    // we use the approximating forms derived in:
+    // "Rational Approximations for the Modified Bessel Function of the
+    // First Kind -- I0(x) for Computations with Double Precision"
+    // by Pavel Holoborodko, see
+    // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
+    // The actual coefficients used are [Boost's] own, and extend
+    // Pavel's work to precisions other than double.
+
+    if (z < 7.75) {
+      // Bessel I0 over[10 ^ -16, 7.75]
+      // Max error in interpolated form : 3.042e-18
+      // Max Error found at double precision = Poly : 5.106609e-16
+      //                                       Cheb : 5.239199e-16
+      static const double P[] = {
+          1.00000000000000000e+00, 2.49999999999999909e-01,
+          2.77777777777782257e-02, 1.73611111111023792e-03,
+          6.94444444453352521e-05, 1.92901234513219920e-06,
+          3.93675991102510739e-08, 6.15118672704439289e-10,
+          7.59407002058973446e-12, 7.59389793369836367e-14,
+          6.27767773636292611e-16, 4.34709704153272287e-18,
+          2.63417742690109154e-20, 1.13943037744822825e-22,
+          9.07926920085624812e-25};
+      return log1p_exp(multiply_log(2.0, z) - log(4.0) +
+                       log(evaluate_polynomial(P, 0.25 * square(z))));
+    }
+    if (z < 500) {
+      // Max error in interpolated form : 1.685e-16
+      // Max Error found at double precision = Poly : 2.575063e-16
+      //                                       Cheb : 2.247615e+00
+      static const double P[] = {
+          3.98942280401425088e-01,  4.98677850604961985e-02,
+          2.80506233928312623e-02,  2.92211225166047873e-02,
+          4.44207299493659561e-02,  1.30970574605856719e-01,
+          -3.35052280231727022e+00, 2.33025711583514727e+02,
+          -1.13366350697172355e+04, 4.24057674317867331e+05,
+          -1.23157028595698731e+07, 2.80231938155267516e+08,
+          -5.01883999713777929e+09, 7.08029243015109113e+10,
+          -7.84261082124811106e+11, 6.76825737854096565e+12,
+          -4.49034849696138065e+13, 2.24155239966958995e+14,
+          -8.13426467865659318e+14, 2.02391097391687777e+15,
+          -3.08675715295370878e+15, 2.17587543863819074e+15};
+      return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z);
+    }
+    // Max error in interpolated form : 2.437e-18
+    // Max Error found at double precision = Poly : 1.216719e-16
+    static const double P[] = {3.98942280401432905e-01, 4.98677850491434560e-02,
+                               2.80506308916506102e-02, 2.92179096853915176e-02,
+                               4.53371208762579442e-02};
+    return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z);
+  }
+  if (v == 1) { // WARNING: will not autodiff for v = 1 correctly
+    // modified from Boost's bessel_i1_imp in the double precision case
+    // see credits above in the v == 0 case
+    if (z < 7.75) {
+      // Bessel I0 over[10 ^ -16, 7.75]
+      // Max error in interpolated form: 5.639e-17
+      // Max Error found at double precision = Poly: 1.795559e-16
+
+      static const double P[] = {
+          8.333333333333333803e-02, 6.944444444444341983e-03,
+          3.472222222225921045e-04, 1.157407407354987232e-05,
+          2.755731926254790268e-07, 4.920949692800671435e-09,
+          6.834657311305621830e-11, 7.593969849687574339e-13,
+          6.904822652741917551e-15, 5.220157095351373194e-17,
+          3.410720494727771276e-19, 1.625212890947171108e-21,
+          1.332898928162290861e-23};
+      T a = square(z) * 0.25;
+      T Q[3] = {1, 0.5, evaluate_polynomial(P, a)};
+      return log(z) + log(evaluate_polynomial(Q, a)) - LOG_TWO;
+    }
+    if (z < 500) {
+      // Max error in interpolated form: 1.796e-16
+      // Max Error found at double precision = Poly: 2.898731e-16
+
+      static const double P[] = {
+          3.989422804014406054e-01,  -1.496033551613111533e-01,
+          -4.675104253598537322e-02, -4.090895951581637791e-02,
+          -5.719036414430205390e-02, -1.528189554374492735e-01,
+          3.458284470977172076e+00,  -2.426181371595021021e+02,
+          1.178785865993440669e+04,  -4.404655582443487334e+05,
+          1.277677779341446497e+07,  -2.903390398236656519e+08,
+          5.192386898222206474e+09,  -7.313784438967834057e+10,
+          8.087824484994859552e+11,  -6.967602516005787001e+12,
+          4.614040809616582764e+13,  -2.298849639457172489e+14,
+          8.325554073334618015e+14,  -2.067285045778906105e+15,
+          3.146401654361325073e+15,  -2.213318202179221945e+15};
+      return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z);
+    }
+    // Max error in interpolated form: 1.320e-19
+    // Max Error found at double precision = Poly: 7.065357e-17
+    static const double P[] = {
+        3.989422804014314820e-01, -1.496033551467584157e-01,
+        -4.675105322571775911e-02, -4.090421597376992892e-02,
+        -5.843630344778927582e-02};
+    return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z);
+  }
+  if (z > 100) {
+    // Boost does something like this in asymptotic_bessel_i_large_x
+    T lim = pow((square(v) + 2.5) / (2 * z), 3) / 24;
+    if (lim < (EPSILON * 10)) {
+      T s = 1;
+      T mu = 4 * square(v);
+      T ex = 8 * z;
+      T num = mu - 1;
+      T denom = ex;
+      s -= num / denom;
+      num *= mu - 9;
+      denom *= ex * 2;
+      s += num / denom;
+      num *= mu - 25;
+      denom *= ex * 3;
+      s -= num / denom;
+      s = z - log(sqrt(z * TWO_PI)) + log(s);
+      return s;
+    }
+  }
+
+  T log_half_z = log(0.5 * z);
+  T lgam = (v > -1) ? lgamma(v + 1.0) : 0;
+  T lcons = (2.0 + v) * log_half_z;
+  T out;
+  if (v > -1) {
+    out = log_sum_exp(v * log_half_z - lgam, lcons - lgamma(v + 2));
+    lgam += log1p(v);
+  } else {
+    out = lcons;
+  }
+
+  double m = 2;
+  double lfac = 0;
+  T old_out;
+  do {
+    lfac += log(m);
+    lgam += log(v + m);
+    lcons += 2 * log_half_z;
+    old_out = out;
+    out = log_sum_exp(out, lcons - lfac - lgam); // underflows eventually
+    m++;
+  } while (out > old_out || out < old_out);
+  return out;
+}
+
+} // namespace CosmoTool
+
+#endif