• Jason Rhinelander's avatar
    Use Cholesky decomposition to invert · f89bf6f1
    Jason Rhinelander authored
    The previous version was using a generic LU inverse, which ran into
    numerical stability problems for large q values.  For example,
    
    q=7, c=0, b=(1.40, 1.41, ..., 1.50) looked like:
    
    0.1767019
    0.1741197
    0.1640931
    0.1330334
    0.2003187
    0.1792207
    0.1792029
    0.175216
    0.1836156
    0.2131243
    0.2235138
    
    because the source data files have only slight variation in that range.
    Using a cholesky decomposition inverse fixes this, and gives numbers
    agreeing exactly with JGM's fortran version (which also uses Cholesky
    decomposition):
    
    0.1744341
    0.1780252
    0.181262
    0.1841298
    0.1866257
    0.1887341
    0.1904548
    0.1917758
    0.1926928
    0.1932153
    0.1933418
    
    which agrees precisely with JGM's (trimmed) fracdist output of:
    
     P value = 0.174
     P value = 0.178
     P value = 0.181
     P value = 0.184
     P value = 0.187
     P value = 0.189
     P value = 0.190
     P value = 0.192
     P value = 0.193
     P value = 0.193
     P value = 0.193
    f89bf6f1
Name
Last commit
Last update
boost Loading commit data...
cmake Loading commit data...
doc Loading commit data...
eigen3 Loading commit data...
fracdist Loading commit data...
.gitignore Loading commit data...
.ycm_extra_conf.py Loading commit data...
BUILDING.md Loading commit data...
CHANGELOG.md Loading commit data...
CMakeLists.txt Loading commit data...
DataParser.pm Loading commit data...
LICENSE Loading commit data...
README.md Loading commit data...
build-data.pl Loading commit data...
fdcrit.cpp Loading commit data...
fdpval.cpp Loading commit data...
parse-vals.hpp Loading commit data...