Skip to content
Commits on Source (4)
language: c
compiler:
- gcc
- clang
matrix:
include:
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-4.6']
env:
- COMPILER=gcc-4.6
script: ./autogen.sh && ./configure && make && make check
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-4.7']
env:
- COMPILER=gcc-4.7
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-4.8']
env:
- COMPILER=gcc-4.8
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-4.9']
env:
- COMPILER=gcc-4.9
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-5']
env:
- COMPILER=gcc-5
- os: linux
compiler: gcc
addons:
apt:
sources: ['ubuntu-toolchain-r-test']
packages: ['gcc-6']
env:
- COMPILER=gcc-6
- os: linux
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-precise-3.5']
packages: ['clang-3.5']
env:
- COMPILER=clang-3.5
- os: linux
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-precise-3.6']
packages: ['clang-3.6']
env:
- COMPILER=clang-3.6
- os: linux
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-precise-3.7']
packages: ['clang-3.7']
env:
- COMPILER=clang-3.7
- os: linux
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-precise-3.8']
packages: ['clang-3.8']
env:
- COMPILER=clang-3.8
- os: linux
dist: trusty
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-trusty-3.9']
packages: ['clang-3.9']
env:
- COMPILER=clang-3.9
- os: linux
dist: trusty
compiler: clang
addons:
apt:
sources: ['ubuntu-toolchain-r-test', 'llvm-toolchain-trusty-4.0']
packages: ['clang-4.0']
env:
- COMPILER=clang-4.0
script: ./autogen.sh && CC=$COMPILER ./configure && make && make check
......@@ -2,6 +2,12 @@
All notable changes to `mptp` will be documented in this file.
This project adheres to [Semantic Versioning](http://semver.org/).
## [0.2.4] - 2018-05-14
### Fixed
- If we do not manage to generate a random starting delimitation with the
wanted number of species (randomly chosen), we use the currently generated
delimitation instead.
## [0.2.3] - 2017-07-25
### Fixed
- Replaced hsearch which was causing problems on APPLE with custom hashtable
......
......@@ -78,9 +78,9 @@ where `DIR` is the directory where bash autocompletion is stored. You can use
and the documentation, use the following commands:
```bash
wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.2/mptp-src-0.2.2.tar.gz
tar zxvf mptp-src-0.2.2.tar.gz
cd mptp-src-0.2.2
wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.4/mptp-src-0.2.4.tar.gz
tar zxvf mptp-src-0.2.4.tar.gz
cd mptp-src-0.2.4
./configure
make
make install # as root, or run sudo make install
......@@ -110,12 +110,12 @@ To use the pre-compiled binary, download the appropriate executable for your
system using the following commands if you are using a Linux system:
```bash
wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.2/mptp-0.2.2-linux-x86_64.tar.gz
tar zxvf mptp-0.2.2-linux-x86_64.tar.gz
wget https://github.com/Pas-Kapli/mptp/releases/download/v0.2.4/mptp-0.2.4-linux-x86_64.tar.gz
tar zxvf mptp-0.2.4-linux-x86_64.tar.gz
```
You will now have the binary distribution in a folder called
`mptp-0.2.2-linux-x86_64` in which you will find three subfolders `bin`, `man`
`mptp-0.2.4-linux-x86_64` in which you will find three subfolders `bin`, `man`
and `doc`. We recommend making a copy or a symbolic link to the mptp binary
`bin/mptp` in a folder included in your `$PATH`, and a copy or a symbolic link
to the mptp man page `man/mptp.1` in a folder included in your `$MANPATH`. The
......
......@@ -2,13 +2,15 @@
# Process this file with autoconf to produce a configure script.
AC_PREREQ([2.63])
AC_INIT([mptp], [0.2.3], [Tomas.Flouri@h-its.org])
AC_INIT([mptp], [0.2.4], [t.flouris@ucl.ac.uk])
AM_INIT_AUTOMAKE([subdir-objects])
AC_LANG([C])
AC_CONFIG_SRCDIR([src/mptp.c])
AC_CONFIG_HEADERS([config.h])
AC_CANONICAL_HOST
${CFLAGS=""}
# Checks for programs.
AC_PROG_CC
AC_PROG_RANLIB
......@@ -72,6 +74,18 @@ AS_IF([test "x$enable_pdfman" != "xno"], [
fi
])
# Various OS-related dependencies
case "${host_os}" in
linux*)
;;
cygwin*|mingw*)
AC_CHECK_LIB([psapi],[GetProcessMemoryInfo])
;;
darwin*)
;;
esac
AM_CONDITIONAL(HAVE_PS2PDF, test "x${have_ps2pdf}" = "xyes")
AM_PROG_CC_C_O
......
mptp (0.2.4-1) UNRELEASED; urgency=medium
* New upstream version
-- Andreas Tille <tille@debian.org> Tue, 29 May 2018 08:54:44 +0200
mptp (0.2.3-1) unstable; urgency=medium
* New upstream version
......
......@@ -7,12 +7,12 @@ Description: -mtune=native should not be used
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -5,7 +5,7 @@ libparse_rtree_a_SOURCES = parse_rtree.y
@@ -4,7 +4,7 @@ libparse_utree_a_SOURCES = parse_utree.y
libparse_rtree_a_SOURCES = parse_rtree.y lex_rtree.l
noinst_LIBRARIES = libparse_utree.a libparse_rtree.a
-AM_CFLAGS=-I${srcdir} -O3 -mtune=native -Wall -Wsign-compare -g ${LIBS}
+AM_CFLAGS=-I${srcdir} -O3 -Wall -Wsign-compare -g ${LIBS}
-AM_CFLAGS=-I${srcdir} -O3 -mtune=native -Wall -Wsign-compare -g
+AM_CFLAGS=-I${srcdir} -O3 -Wall -Wsign-compare -g
AM_YFLAGS = -d -p `${SED} -n 's/.*_\(.*\)/\1_/p' <<<"$*"`
AM_LFLAGS = -o lex.yy.c
.\" -*- coding: utf-8 -*-
.\" ============================================================================
.TH mptp 1 "January 31, 2017" "mptp 0.2.2" "USER COMMANDS"
.TH mptp 1 "May 14, 2018" "mptp 0.2.4" "USER COMMANDS"
.\" ============================================================================
.SH NAME
mptp \(em single-locus species delimitation
......@@ -354,5 +354,10 @@ incorrectly formatted trees.
.TP
.BR v0.2.3\~ "released July 25th, 2017"
Replaced hsearch() with custom hashtable. Fixed minor output error messages.
.TP
.BR v0.2.4\~ "released May 14th, 2018"
If we do not manage to generate a random starting delimitation with the wanted
number of species (randomly chosen), we use the currently generated
delimitation instead.
.RE
.LP
......@@ -4,8 +4,7 @@ libparse_utree_a_SOURCES = parse_utree.y lex_utree.l
libparse_rtree_a_SOURCES = parse_rtree.y lex_rtree.l
noinst_LIBRARIES = libparse_utree.a libparse_rtree.a
AM_CFLAGS=-I${srcdir} -O3 -mtune=native -Wall -Wsign-compare -g ${LIBS}
AM_CFLAGS=-I${srcdir} -O3 -mtune=native -Wall -Wsign-compare -g
AM_YFLAGS = -d -p `${SED} -n 's/.*_\(.*\)/\1_/p' <<<"$*"`
AM_LFLAGS = -o lex.yy.c
......
......@@ -258,7 +258,7 @@ static void mcmc_finalize(rtree_t * root,
free(densities);
}
static void dp_recurse(rtree_t * node, int method)
static void dp_recurse(rtree_t * node, long method)
{
int k,j;
......@@ -830,7 +830,7 @@ void aic_mcmc(rtree_t * tree,
/* throw a coin to decide whether to convert a coalescent root to a
speciation or the other way round */
rand_double = erand48(rstate);
rand_double = mptp_erand48(rstate);
int speciation = (rand_double >= 0.5) ? 1 : 0;
if ((speciation && crnodes_count) || (snodes_count == 0))
......@@ -844,7 +844,7 @@ void aic_mcmc(rtree_t * tree,
/* select a coalescent root, split it into two coalescent nodes */
rand_long = nrand48(rstate);
rand_long = mptp_nrand48(rstate);
long r = rand_long % crnodes_count;
rtree_t * node = crnodes[r];
......@@ -923,7 +923,7 @@ void aic_mcmc(rtree_t * tree,
}
/* decide whether to accept or reject proposal */
rand_double = erand48(rstate);
rand_double = mptp_erand48(rstate);
if (rand_double <= a)
{
/* accept */
......@@ -987,7 +987,7 @@ void aic_mcmc(rtree_t * tree,
/ \ / \
CR * * CR C * * C */
rand_long = nrand48(rstate);
rand_long = mptp_nrand48(rstate);
long r = rand_long % snodes_count;
rtree_t * node = snodes[r];
......@@ -1063,7 +1063,7 @@ void aic_mcmc(rtree_t * tree,
}
/* decide whether to accept or reject proposal */
rand_double = erand48(rstate);
rand_double = mptp_erand48(rstate);
if (rand_double <= a)
{
/* accept */
......
......@@ -23,47 +23,56 @@
unsigned long arch_get_memused()
{
#ifdef _WIN32
PROCESS_MEMORY_COUNTERS pmc;
GetProcessMemoryInfo(GetCurrentProcess(),
&pmc,
sizeof(PROCESS_MEMORY_COUNTERS));
return pmc.PeakWorkingSetSize;
#else
struct rusage r_usage;
getrusage(RUSAGE_SELF, & r_usage);
#if defined __APPLE__
# ifdef __APPLE__
/* Mac: ru_maxrss gives the size in bytes */
return (unsigned long)(r_usage.ru_maxrss);
return r_usage.ru_maxrss;
# else
/* Linux: ru_maxrss gives the size in kilobytes */
return (unsigned long)r_usage.ru_maxrss * 1024;
return r_usage.ru_maxrss * 1024;
# endif
#endif
}
unsigned long arch_get_memtotal()
{
#if defined(_SC_PHYS_PAGES) && defined(_SC_PAGESIZE)
#ifdef _WIN32
long phys_pages = sysconf(_SC_PHYS_PAGES);
long pagesize = sysconf(_SC_PAGESIZE);
if ((phys_pages == -1) || (pagesize == -1))
fatal("Cannot determine amount of RAM");
// sysconf(3) notes that pagesize * phys_pages can overflow, such as
// when long is 32-bits and there's more than 4GB RAM. Since vsearch
// apparently targets LP64 systems like x86_64 linux, this will not
// arise in practice on the intended platform.
if (pagesize > LONG_MAX / phys_pages)
return LONG_MAX;
else
return (unsigned long)pagesize * (unsigned long)phys_pages;
MEMORYSTATUSEX ms;
ms.dwLength = sizeof(MEMORYSTATUSEX);
GlobalMemoryStatusEx(&ms);
return ms.ullTotalPhys;
#elif defined(__APPLE__)
int mib [] = { CTL_HW, HW_MEMSIZE };
int64_t ram = 0;
long ram = 0;
size_t length = sizeof(ram);
if(-1 == sysctl(mib, 2, &ram, &length, NULL, 0))
if(sysctl(mib, 2, &ram, &length, NULL, 0) == -1)
fatal("Cannot determine amount of RAM");
return ram;
#elif defined(_SC_PHYS_PAGES) && defined(_SC_PAGESIZE)
long phys_pages = sysconf(_SC_PHYS_PAGES);
long pagesize = sysconf(_SC_PAGESIZE);
if ((phys_pages == -1) || (pagesize == -1))
fatal("Cannot determine amount of RAM");
return pagesize * phys_pages;
#else
struct sysinfo si;
......@@ -73,3 +82,76 @@ unsigned long arch_get_memtotal()
#endif
}
long arch_get_cores()
{
#ifdef _WIN32
SYSTEM_INFO si;
GetSystemInfo(&si);
return si.dwNumberOfProcessors;
#else
return sysconf(_SC_NPROCESSORS_ONLN);
#endif
}
void arch_get_user_system_time(double * user_time, double * system_time)
{
*user_time = 0;
*system_time = 0;
#ifdef _WIN32
HANDLE hProcess = GetCurrentProcess();
FILETIME ftCreation, ftExit, ftKernel, ftUser;
ULARGE_INTEGER ul;
GetProcessTimes(hProcess, &ftCreation, &ftExit, &ftKernel, &ftUser);
ul.u.HighPart = ftUser.dwHighDateTime;
ul.u.LowPart = ftUser.dwLowDateTime;
*user_time = ul.QuadPart * 100.0e-9;
ul.u.HighPart = ftKernel.dwHighDateTime;
ul.u.LowPart = ftKernel.dwLowDateTime;
*system_time = ul.QuadPart * 100.0e-9;
#else
struct rusage r_usage;
getrusage(RUSAGE_SELF, & r_usage);
* user_time = r_usage.ru_utime.tv_sec * 1.0
+ r_usage.ru_utime.tv_usec * 1.0e-6;
* system_time = r_usage.ru_stime.tv_sec * 1.0
+ r_usage.ru_stime.tv_usec * 1.0e-6;
#endif
}
void arch_srandom()
{
/* initialize pseudo-random number generator */
unsigned int seed = opt_seed;
if (seed == 0)
{
#ifdef _WIN32
srand(GetTickCount());
#else
int fd = open("/dev/urandom", O_RDONLY);
if (fd < 0)
fatal("Unable to open /dev/urandom");
if (read(fd, & seed, sizeof(seed)) < 0)
fatal("Unable to read from /dev/urandom");
close(fd);
srandom(seed);
#endif
}
else
{
#ifdef _WIN32
srand(seed);
#else
srandom(seed);
#endif
}
}
long arch_random()
{
#ifdef _WIN32
return rand();
#else
return random();
#endif
}
......@@ -193,7 +193,7 @@ static void link_sequences(rtree_t * root, char ** headers, char ** sequence, lo
rtree_query_tipnodes(root, tipnodes);
/* create a libc hash table of size tip_count */
hashtable_t * ht = hashtable_create(root->leaves);
hashtable_t * ht = hashtable_create((unsigned long)(root->leaves));
/* populate a libc hash table with tree tip labels */
for (i = 0; i < root->leaves; ++i)
......
......@@ -186,7 +186,6 @@ void dp_ptp(rtree_t * tree, long method)
dp_vector_t * vec = tree->vector;
if (method == PTP_METHOD_MULTI)
{
max = vec[0].score_multi;
double min_aic_score = aic(vec[0].score_multi, vec[0].species_count, tree->leaves+2);
for (i = 1; i < tree->edge_count; i++)
{
......@@ -201,6 +200,7 @@ void dp_ptp(rtree_t * tree, long method)
}
}
}
max = vec[best_index].score_multi;
}
else
{
......@@ -223,8 +223,10 @@ void dp_ptp(rtree_t * tree, long method)
tree->edge_count,
2 * tree->leaves - 2);
printf("Score Null Model: %.6f\n", tree->coal_logl);
if (method == PTP_METHOD_SINGLE)
fprintf(stdout, "Best score for single coalescent rate: %.6f\n",
vec[best_index].score_single);
else
fprintf(stdout, "Best score for multi coalescent rate: %.6f\n",
vec[best_index].score_multi);
}
......
......@@ -582,7 +582,7 @@ void fillheader()
"%s %s_%s, %1.fGB RAM, %ld cores",
PROG_NAME, PROG_VERSION, PROG_ARCH,
arch_get_memtotal() / 1024.0 / 1024.0 / 1024.0,
sysconf(_SC_NPROCESSORS_ONLN));
arch_get_cores());
}
void show_header()
......
......@@ -22,7 +22,8 @@
#define _GNU_SOURCE
#include <assert.h>
#include <search.h>
#include <ctype.h>
#include <fcntl.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
......@@ -33,10 +34,11 @@
#include <limits.h>
#include <locale.h>
#include <math.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <stdbool.h>
#include <time.h>
#ifdef HAVE_CONFIG_H
#include "config.h"
......@@ -51,12 +53,40 @@
#define PROG_NAME PACKAGE
#define PROG_VERSION PACKAGE_VERSION
#ifdef __APPLE__
#define PROG_ARCH "macosx_x86_64"
#ifdef __PPC__
#ifdef __LITTLE_ENDIAN__
#define PROG_CPU "ppc64le"
#else
#error "Big endian ppc64 CPUs not supported"
#endif
#else
#define PROG_ARCH "linux_x86_64"
#define PROG_CPU "x86_64"
#endif
#ifdef __APPLE__
#define PROG_OS "osx"
#include <sys/resource.h>
#include <sys/sysctl.h>
#endif
#ifdef __linux__
#define PROG_OS "linux"
#include <sys/resource.h>
#include <sys/sysinfo.h>
#endif
#ifdef _WIN32
#define PROG_OS "win"
#include <windows.h>
#include <psapi.h>
#endif
#define PROG_ARCH PROG_OS "_" PROG_CPU
#define PLL_FAILURE 0
#define PLL_SUCCESS 1
#define PLL_LINEALLOC 2048
......@@ -206,7 +236,7 @@ typedef struct hashtable_s
typedef struct pair_s
{
char * label;
int index;
size_t index;
} pair_t;
/* macros */
......@@ -289,9 +319,10 @@ char * xstrchrnul(char *s, int c);
char * xstrdup(const char * s);
char * xstrndup(const char * s, size_t len);
long getusec(void);
void show_rusage(void);
FILE * xopen(const char * filename, const char * mode);
void random_init(unsigned short * rstate, long seedval);
double mptp_erand48(unsigned short * rstate);
long mptp_nrand48(unsigned short * rstate);
/* functions in mptp.c */
......@@ -366,6 +397,7 @@ void lca_destroy(void);
unsigned long arch_get_memused(void);
unsigned long arch_get_memtotal(void);
long arch_get_cores(void);
/* functions in dp.c */
......
......@@ -155,7 +155,7 @@ void multirun(rtree_t * root, long method)
/* generate one seed for each run */
seeds = (long *)xmalloc((size_t)opt_mcmc_runs * sizeof(long));
for (i = 0; i < opt_mcmc_runs; ++i)
seeds[i] = nrand48(global_xsubi);
seeds[i] = mptp_nrand48(global_xsubi);
if (opt_mcmc_runs == 1)
seeds[0] = opt_seed;
......@@ -174,7 +174,7 @@ void multirun(rtree_t * root, long method)
across all MCMC runs */
double * combined_val;
combined_val = (double *)xmalloc((size_t)(root->leaves-1) * sizeof(double));
memset(combined_val,0,(root->leaves-1)*sizeof(double));
memset(combined_val,0,(unsigned long)(root->leaves-1)*sizeof(double));
rtree_t ** inner_node_list = (rtree_t **)xmalloc((size_t)(root->leaves-1) *
sizeof(rtree_t *));
......@@ -291,7 +291,6 @@ void multirun(rtree_t * root, long method)
double mean, var, stdev, avg_stdev = 0;
for (i = 0; i < support_count; ++i)
{
int j;
mean = var = stdev = 0;
for (j = 0; j < opt_mcmc_runs; ++j)
mean += support[j][i];
......
......@@ -51,7 +51,7 @@ static int cb_node_select(rtree_t * node)
}
/* otherwise, we just throw a coin and select one of the two cases */
rand_double = erand48(g_rstate);
rand_double = mptp_erand48(g_rstate);
if (rand_double >= 0.5)
{
/* don't select */
......@@ -86,7 +86,7 @@ double random_delimitation(rtree_t * root,
max_species = root->max_species_count;
g_rstate = rstate;
rand_long = nrand48(rstate);
rand_long = mptp_nrand48(rstate);
if (!root->max_species_count)
species_count = (rand_long % root->leaves) + 1;
else
......@@ -116,6 +116,13 @@ double random_delimitation(rtree_t * root,
free(inner_node_list);
assert(count <= species_count);
if (count < species_count)
{
/* TODO: This fixes issue #82, but we should implement a better, non-biased
way of generatng random starting delimitations */
species_count = count;
}
assert(count == species_count);
*delimited_species = species_count;
......
......@@ -205,7 +205,7 @@ static void rtree_traverse_recursive(rtree_t * node,
return;
}
rand_double = erand48(rstate);
rand_double = mptp_erand48(rstate);
if (rand_double >= 0.5)
{
rtree_traverse_recursive(node->left, cbtrav, index, rstate, outbuffer);
......@@ -450,7 +450,7 @@ static rtree_t ** rtree_tipstring_nodes(rtree_t * root,
sizeof(rtree_t *));
/* create a hashtable of tip labels */
hashtable_t * ht = hashtable_create(root->leaves);
hashtable_t * ht = hashtable_create((unsigned long)(root->leaves));
for (i = 0; i < (unsigned int)(root->leaves); ++i)
{
......
......@@ -21,6 +21,11 @@
#include "mptp.h"
#define MPTP_RAND48_MULT_0 (0xe66d)
#define MPTP_RAND48_MULT_1 (0xdeec)
#define MPTP_RAND48_MULT_2 (0x0005)
#define MPTP_RAND48_ADD (0x000b)
static const char * progress_prompt;
static unsigned long progress_next;
static unsigned long progress_size;
......@@ -144,23 +149,6 @@ long getusec(void)
return tv.tv_sec * 1000000 + tv.tv_usec;
}
void show_rusage()
{
struct rusage r_usage;
getrusage(RUSAGE_SELF, & r_usage);
fprintf(stderr, "Time: %.3fs (user)", r_usage.ru_utime.tv_sec * 1.0 + (double) r_usage.ru_utime.tv_usec * 1.0e-6);
fprintf(stderr, " %.3fs (sys)", r_usage.ru_stime.tv_sec * 1.0 + r_usage.ru_stime.tv_usec * 1.0e-6);
#if defined __APPLE__
/* Mac: ru_maxrss gives the size in bytes */
fprintf(stderr, " Memory: %.0fMB\n", r_usage.ru_maxrss * 1.0e-6);
#else
/* Linux: ru_maxrss gives the size in kilobytes */
fprintf(stderr, " Memory: %.0fMB\n", r_usage.ru_maxrss * 1.0e-3);
#endif
}
FILE * xopen(const char * filename, const char * mode)
{
FILE * out = fopen(filename, mode);
......@@ -177,3 +165,42 @@ void random_init(unsigned short * rstate, long seedval)
rstate[1] = seedval & 0xffffl;
rstate[2] = seedval >> 16;
}
static void mptp_randomize(unsigned short * rstate)
{
unsigned long accu;
unsigned short temp[2];
accu = MPTP_RAND48_MULT_0 * (unsigned long)rstate[0] +
MPTP_RAND48_ADD;
temp[0] = (unsigned short)accu; /* lower 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += (unsigned long)MPTP_RAND48_MULT_0 * (unsigned long)rstate[1] +
(unsigned long)MPTP_RAND48_MULT_1 * (unsigned long)rstate[0];
temp[1] = (unsigned short)accu; /* middle 16 bits */
accu >>= sizeof(unsigned short) * 8;
accu += MPTP_RAND48_MULT_0 * rstate[2] +
MPTP_RAND48_MULT_1 * rstate[1] +
MPTP_RAND48_MULT_2 * rstate[0];
rstate[0] = temp[0];
rstate[1] = temp[1];
rstate[2] = (unsigned short)accu;
}
double mptp_erand48(unsigned short * rstate)
{
mptp_randomize(rstate);
return ldexp((double)rstate[0], -48) +
ldexp((double)rstate[1], -32) +
ldexp((double)rstate[2], -16);
}
long mptp_nrand48(unsigned short * rstate)
{
mptp_randomize(rstate);
return ((long)rstate[2] << 15) + ((long)rstate[1] >> 1);
}