Commit 867ce173 authored by Frederic GOUALARD's avatar Frederic GOUALARD 💬
Browse files

Changed License file to comply with license of the project

parent a8e6746d
This diff is collapsed.
2020-03-24 Frédéric Goualard <Frederic.Goualard@univ-nantes.fr>
* src/float64.c (fpngl_nan64): corrected bug whereby the
function could also return \pm\infty.
2020-03-13 Frédéric Goualard <Frederic.Goualard@univ-nantes.fr>
* src/types.h: removed file. Replaced by frng_t.h
No preview for this file type
......@@ -51,7 +51,7 @@ More information about the FPNGlib library can be found at the project homepage,
@uref{https://gitlab.univ-nantes.fr/goualard-f/fpnglib}.
This manual assumes a basic knowledge of the IEEE~754 standard and refers to the internal representation of floating-point numbers in several places. Readers unfamiliar with it may find beneficial to at least read the @url{https://en.wikipedia.org/wiki/IEEE_754,Wikipedia entry} for the standard first.
This manual assumes a basic knowledge of the IEEE@tie{}754 standard and refers to the internal representation of floating-point numbers in several places. Readers unfamiliar with the standard may find beneficial to at least read its @url{https://en.wikipedia.org/wiki/IEEE_754,Wikipedia entry} first.
The FPNGlib library contains code from third parties, specifically:
@itemize @bullet
......@@ -471,6 +471,8 @@ The number of integer Random Number Generators augments regularly with new versi
@float Table,tab:irngs
@caption{Integer RNGs available in FPNGlib@tie{}@value{VERSION} and their characteristics.}
@tex\hrule
@end tex
@multitable @columnfractions .4 .15 .15 .15 .15
@headitem Name @tab Size @tab Min. @tab Max. @tab Period
@item @code{drand48_lcg64} @tab 64 @tab 0 @tab @math{2^{48}-1} @tab ?
......@@ -484,6 +486,8 @@ The number of integer Random Number Generators augments regularly with new versi
@item @code{xor4096iv32} @tab 32 @tab 0 @tab @math{2^{32}-1} @tab ?
@item @code{xor4096iv64} @tab 64 @tab 0 @tab @math{2^{64}-1} @tab ?
@end multitable
@tex\hrule
@end tex
@end float
@subsection 32 Bits Integer RNGs
......@@ -659,8 +663,233 @@ Return a random 32 bits integer according to the discrete distribution @var{dist
@chapter Floating-Point Random Generators
@c ======================================
The standard libraries of most programming languages offers only one algorithm to compute random IEEE@tie{}754 floating-point numbers. Usually, it computes floats by dividing random integers by some constant. This imparts a specific structure to the fractional parts of the floats, which may adversely impact programs that use them.
FPNGlib gives the user the choice of the method used to produce floats. All float RNGs are either of the type @code{fpngl_frng32_t} or @code{fpngl_frng64_t}; they all must propose the same functionalities:
@deftypefn {Library Function} void fpngl_frng32_delete (fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} void fpngl_frng64_delete (fpngl_frng64_t* @var{rng})
Reclaim the memory allocated to the RNG @var{rng}. If it was constructed with an
integer RNG, the memory for that one is also reclaimed since float RNGs ``own'' the
integer RNGs they are made of.
@end deftypefn
@deftypefn {Library Function} float fpngl_frng32_min (fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} float fpngl_frng32_max (fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} double fpngl_frng64_min (fpngl_frng64_t* @var{rng})
@deftypefnx {Library Function} double fpngl_frng64_max (fpngl_frng64_t* @var{rng})
Return, respectively, the minimum and the maximum possible values drawn by the
floating-point RNG.
In the literature, many FRNGs are described as returning values in the domain
@math{[0,1)}. Technically, their actual domain is not an interval but a set of discrete values. It is then always possible to describe the domain with closed brackets. The functions @code{*_min} and @code{*_max} return the smallest and largest values in this set.
@end deftypefn
@deftypefn {Library Function} float fpngl_frng32_nextf32 (fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} double fpngl_frng64_nextf64 (fpngl_frng64_t* @var{rng})
Return the next random floating-point number uniformly drawn from, respectively,
the domain @code{[fpngl_frng32_min(rng), fpngl_frng32_max(rng)]}, or the domain
@code{[fpngl_frng64_min(rng), fpngl_frng64_max(rng)]}.
@end deftypefn
@deftypefn {Library Function} float fpngl_frng32_next_arrayf32 @
(fpngl_frng32_t* @var{rng}, float *@var{T}, uint32_t @var{n})
@deftypefnx {Library Function} double fpngl_frng64_next_arrayf64 @
(fpngl_frng64_t* @var{rng}, double *@var{T}, uint32_t @var{n})
Fill the array @var{T} with @var{n} random floating-point numbers drawn with the FRNG
@var{rng}.
@strong{Caution}: some FRNGs may allocate a temporary array of size @math{2@times@var{n}} to implement this function.
@end deftypefn
@deftypefn {Library Function} {const char*} fpngl_fprng32_name @
(fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} {const char*} fpngl_fprng64_name @
(fpngl_frng64_t* @var{rng})
Return the string corresponding to the name of the FRNG (@xref{Existing float RNGs}).
@end deftypefn
@deftypefn {Library Function} uint32_t fpngl_fprng32_seed @
(fpngl_frng32_t* @var{rng})
@deftypefnx {Library Function} uint64_t fpngl_fprng64_seed @
(fpngl_frng64_t* @var{rng})
Return the integer used to initialize the FRNG.
@end deftypefn
@node Existing float RNGs
@section Existing float RNGs
@c -------------------------
@node Single Precision FRNGs
@subsection Single Precision FRNGs
@c ...............................
@deftypefn {@file{fpnglib/xorgens.h}} fpngl_frng32_t* fpngl_xor4096rv32 @
(uint32_t @var{seed})
??
@end deftypefn
@node Double Precision FRNGs
@subsection Double Precision FRNGs
@c ...............................
@deftypefn {@file{fpnglib/frng64_division.h}} fpngl_frng64_t* fpngl_matlabp5 (uint64_t @var{seed})
Double precision FRNG used by MATLAB until MATLAB 5.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_division.h}} fpngl_frng64_t* fpngl_drand48bsd @
(uint64_t @var{seed})
Implementation of the POSIX function @code{drand48} in FreeBSD.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_division.h}} fpngl_frng64_t* fpngl_mupad @
(uint64_t @var{seed})
Generator used in Maple 10 and MuPAD Pro 3.1.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_division_k.h}} fpngl_frng64_t* fpngl_div53 @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
Draw a double precision number in the domain @math{@{0@}@cup@{2^{-53},1-2^{-53}@}} by
dividing an integer generated with the IRNG @var{irng} by @math{2^{53}}.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_division_k.h}} fpngl_frng64_t* fpngl_div32 @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
Draw a double precision number in the domain @math{@{0@}@cup@{2^{-32},1-2^{-32}@}} by
dividing an integer generated with the IRNG @var{irng} by @math{2^{32}}.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_drand48gnu.h}} fpngl_frng64_t* fpngl_drand48gnu @
(uint64_t @var{seed})
Implementation of the POSIX @code{drand48} function as implemented in GNU @command{gcc}.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_fog97.h}} fpngl_frng64_t* fpngl_fog97 @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
Computation of a double precision number in @math{[0,1-2^{-52}]}
by generating a 52 bits random
fractional part, concatenating it with an exponent of @math{0} to obtain
a float in @math{[1,2-2^{-52}]}, and then subtracting @math{1.0} from it.
@end deftypefn
@deftypefn {@file{fpnglib/frng64_fog05.h}} fpngl_frng64_t* fpngl_fog05 @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/frng64_golang.h}} fpngl_frng64_t* fpngl_golang @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/frng64_java.h}} fpngl_frng64_t* fpngl_java @
(uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/frng64_KGW.h}} fpngl_frng64_t* fpngl_KGW @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/frng64_lecuyer_simard.h}} fpngl_frng64_t* fpngl_lecuyer_simard @
(fpngl_irng_t* @var{irng}, uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/frng64_rationalLCG10.h}} fpngl_frng64_t* fpngl_rationalLCG10 @
(uint64_t @var{seed})
??
@end deftypefn
@deftypefn {@file{fpnglib/xorgens.h}} fpngl_frng64_t* fpngl_xor4096rv64 @
(uint64_t @var{seed})
??
@end deftypefn
@code Drawing floats in a class
@section Drawing floats in a class
@c -------------------------------
@deftypefn {@file{fpnglib/float64.h}} double fpngl_subnormal64 (fpngl_irng_t* @var{rng})
Return a random double precision subnormal number drawn uniformly from the domain
@math{(-@code{fpngl_lambda64()}, @code{fpngl_lambda64()})}.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} double fpngl_zero64 (fpngl_irng_t* @var{rng})
Return zero in double precision as @math{+0.0} or @math{-0.0} with even probability.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} double fpngl_inf64 (fpngl_irng_t* @var{rng})
Return the double precision value reserved for infinity as @math{+@infty} or
@math{-@infty} with even probability.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} double fpngl_normal64 (fpngl_irng_t* @var{rng})
Return a random double precision number drawn uniformly in the domain
@math{[-@code{fpngl_max64()}, -@code{fpngl_lambda64()}]@cup[@code{fpngl_lambda64()}, @code{fpngl_max64()}]}.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} double fpngl_nan64 (fpngl_irng_t* @var{rng})
Return randomly one of the @math{2^{53}-2} possible @emph{Not-A-Number} double precision
values.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} double fpngl_float64 @
(fpngl_irng_t* @var{rng},@
fpngl_sign_t @var{s},@
uint32_t @var{minexp}, uint32_t @var{maxexp},@
uint64_t @var{minfrac}, uint64_t @var{maxfrac},@
uint64_t @var{andmask},@
uint64_t @var{ormask})
Return a random double precision number whose IEEE@tie{}754 structure (sign, exponent, fractional part) obeys the constraints given as parameters:
@itemize
@item@var{s}: possible sign. May be @code{fpngl_positive}, @code{fpngl_negative}, or @code{fpngl_whatever} if both positive and negative signs are allowed;
@item@var{minexp}: minimum exponent allowed. This is the biased value, which means you must add @math{1023} to the actual value. Possible values are in the domain
@math{[0,2047]};
@item@var{maxexp}: maximum exponent allowed. This is the biased value, which means you must add @math{1023} to the actual value. Possible values are in the domain
@math{[0,2047]};
@item@var{minfrac}: minimum fractional part expressed as a 52 bits unsigned integer
between @code{0x0000000000000000} and @code{0x000fffffffffffff};
@item@var{maxfrac}: maximum fractional part expressed as a 52 bits unsigned integer
between @code{0x0000000000000000} and @code{0x000fffffffffffff};
@item@var{andmask}: bitwise @code{AND} mask applied to the integer obtained by
concatenating the random sign, biased exponent, and fractional part. Should be
@code{fpngl_noand64}@footnote{The @code{fpngl_noand64} constant has value @code{0xffffffffffffffff} and is defined in @file{fpnglib/constants64.h}.} if no mask is to be applied;
@item@var{ormask}: bitwise @code{OR} mask applied to the result of the bitwise AND. Should be @code{fpngl_noor64}@footnote{The code{fpngl_noor64} constant has value @code{0x0000000000000000} and is defined in @file{fpnglib/constants64.h}.} if no mask is needed.
@end itemize
Note that, for better performances, the consistency of the constraints are not tested.
@end deftypefn
@node Drawing floats non-uniformly
@section Drawing floats non-uniformly
@c ----------------------------------
The FPNGlib library proposes a simple way to draw floating-point numbers from one of
the five classes presented in @ref{Drawing floats in a class} according to some discrete probability distribution. The user first defines the probability distribution for each of the classes in the order ``@code{zero, subnormal, normal, infinite, nan}'' with the function @code{fpngl_class_float64_new}; it is then possible to draw floats following that distribution with the function @code{fpngl_float64_distrib}. This functionality is just a facility offered around the more general method presented in @xref{Non-Uniform Distributions}.
@deftypefn {@file{fpnglib/float64.h}} fpngl_distribution_t* fpngl_class_float64_new@
(fpngl_irng_t* @var{irng}, const double @var{P}[static 5])
Create an @code{fpngl_distribution_t} object representing the discrete probability distribution defined by the array @var{P} of at least 5 values:
@itemize
@item @code{P[0]}: probability of returning @math{@pm0.0};
@item @code{P[1]}: probability of returning a subnormal number;
@item @code{P[2]}: probability of returning a normal number;
@item @code{P[3]}: probability of returning @math{@pm@infty};
@item @code{P[4]}: probability of returning an @emph{NaN}.
@end itemize
The sum of the first five values in @var{P} must be equal to 1. Note that, inside the classes 1 and 2 (subnormal and normal numbers), the floats are still drawn uniformly at random.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} fpngl_distribution_t* fpngl_class_float64_new@
(fpngl_irng_t* @var{irng}, const double @var{P}[static 5])
Create an @code{fpngl_distribution_t} object representing the discrete probability distribution defined by the array @var{P} of at least 5 values.
@end deftypefn
@node IEEE 754 and the Floating-Point Unit
......@@ -681,6 +910,7 @@ and @file{fpnglib/constants64.h} ---for the double precision.
@deftypevrx Constant uint32_t fpngl_emax32
@deftypevrx Constant uint32_t fpngl_emax64
Smallest and largest unbiased exponents for single and double precision float.
@noindent @strong{Example}
@cartouche
@verbatiminclude snippets/exponent.c
......@@ -694,13 +924,16 @@ Double precision exponent in [-1022, 1023]
@deftypevr Constant uint32_t fpngl_t32
@deftypevrx Constant uint32_t fpngl_t64
Number of bits in the significand of a single or double precision float.
Number of bits in the significand of a single or double precision float. The @dfn{significand} of a float corresponds to the fractional part together with the hidden bit before the binary point.
@noindent @strong{Example:}
@cartouche
@verbatiminclude snippets/significand.c
@end cartouche
@noindent @strong{Output:}
@verbatim
Number of bits with single precision: 24
Number of bits with double precision: 53
Number of digits with single precision: 7
Number of digits with double precision: 15
@end verbatim
......@@ -856,9 +1089,26 @@ If @var{v} is an NaN, the function returns the same NaN. Otherwise, it returns e
The output is @code{1} since the smallest positive double precision float is separated from the largest negative float by @code{0}.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} uint64_t fpngl_distance_float64 @
(double @var{a}, double @var{b})
Return the number of floats (bounds included) from value @var{a} to value @var{b}. Return
@code{1} if @code{@var{a}==@var{b}} and @code{0} if either @var{a} or var{b} are not
finite floats (infinities or NaNs).
The function is only defined for @code{@var{a} <= @var{b}}.
@end deftypefn
@deftypefn {@file{fpnglib/float64.h}} uint64_t fpngl_signed_distance_float64 @
(double @var{a}, double @var{b})
Return the number of floats (bounds included) from value @var{a} to value @var{b}. Return
@code{1} if @code{@var{a}==@var{b}} and @code{0} if either @var{a} or var{b} are not
finite floats (infinities or NaNs). Return a negative value if @code{@var{a} < @var{b}}.
@end deftypefn
@section Manipulating the FPU
@c ------------------------
@node Extending the Library
@chapter Extending the Library
@c ===========================
......
......@@ -115,7 +115,7 @@ noinst_PROGRAMS = significand$(EXEEXT) exponent$(EXEEXT) \
mu$(EXEEXT) lambda$(EXEEXT) nan$(EXEEXT) inf$(EXEEXT) \
max$(EXEEXT) subnormal$(EXEEXT) float64$(EXEEXT) \
minmaxfrac64$(EXEEXT) noandor$(EXEEXT) distrib$(EXEEXT) \
irng$(EXEEXT) biased_dice_2$(EXEEXT) biased_dice_3$(EXEEXT)
irng$(EXEEXT) biased_dice_2$(EXEEXT)
subdir = doc/snippets
ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
am__aclocal_m4_deps = $(top_srcdir)/m4/ax_check_enable_debug.m4 \
......@@ -148,11 +148,6 @@ biased_dice_2_OBJECTS = biased_dice_2.$(OBJEXT)
biased_dice_2_LDADD = $(LDADD)
biased_dice_2_DEPENDENCIES = $(top_builddir)/src/libfpnglib.la \
$(am__DEPENDENCIES_1)
biased_dice_3_SOURCES = biased_dice_3.c
biased_dice_3_OBJECTS = biased_dice_3.$(OBJEXT)
biased_dice_3_LDADD = $(LDADD)
biased_dice_3_DEPENDENCIES = $(top_builddir)/src/libfpnglib.la \
$(am__DEPENDENCIES_1)
debug_SOURCES = debug.c
debug_OBJECTS = debug.$(OBJEXT)
debug_LDADD = $(LDADD)
......@@ -297,16 +292,16 @@ AM_V_CCLD = $(am__v_CCLD_@AM_V@)
am__v_CCLD_ = $(am__v_CCLD_@AM_DEFAULT_V@)
am__v_CCLD_0 = @echo " CCLD " $@;
am__v_CCLD_1 =
SOURCES = biased_dice.c biased_dice_2.c biased_dice_3.c debug.c dice.c \
distrib.c exponent.c float64.c fsimple.c inf.c irng.c \
isimple.c lambda.c max.c mcpi.c minmaxfrac64.c mu.c nan.c \
nextafter64.c noandor.c previous64.c significand.c subnormal.c \
unitroundoff.c warning.c
DIST_SOURCES = biased_dice.c biased_dice_2.c biased_dice_3.c debug.c \
dice.c distrib.c exponent.c float64.c fsimple.c inf.c irng.c \
isimple.c lambda.c max.c mcpi.c minmaxfrac64.c mu.c nan.c \
nextafter64.c noandor.c previous64.c significand.c subnormal.c \
unitroundoff.c warning.c
SOURCES = biased_dice.c biased_dice_2.c debug.c dice.c distrib.c \
exponent.c float64.c fsimple.c inf.c irng.c isimple.c lambda.c \
max.c mcpi.c minmaxfrac64.c mu.c nan.c nextafter64.c noandor.c \
previous64.c significand.c subnormal.c unitroundoff.c \
warning.c
DIST_SOURCES = biased_dice.c biased_dice_2.c debug.c dice.c distrib.c \
exponent.c float64.c fsimple.c inf.c irng.c isimple.c lambda.c \
max.c mcpi.c minmaxfrac64.c mu.c nan.c nextafter64.c noandor.c \
previous64.c significand.c subnormal.c unitroundoff.c \
warning.c
am__can_run_installinfo = \
case $$AM_UPDATE_INFO_DIR in \
n|no|NO) false;; \
......@@ -508,10 +503,6 @@ biased_dice_2$(EXEEXT): $(biased_dice_2_OBJECTS) $(biased_dice_2_DEPENDENCIES) $
@rm -f biased_dice_2$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(biased_dice_2_OBJECTS) $(biased_dice_2_LDADD) $(LIBS)
biased_dice_3$(EXEEXT): $(biased_dice_3_OBJECTS) $(biased_dice_3_DEPENDENCIES) $(EXTRA_biased_dice_3_DEPENDENCIES)
@rm -f biased_dice_3$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(biased_dice_3_OBJECTS) $(biased_dice_3_LDADD) $(LIBS)
debug$(EXEEXT): $(debug_OBJECTS) $(debug_DEPENDENCIES) $(EXTRA_debug_DEPENDENCIES)
@rm -f debug$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(debug_OBJECTS) $(debug_LDADD) $(LIBS)
......@@ -608,7 +599,6 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/biased_dice.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/biased_dice_2.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/biased_dice_3.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/debug.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/dice.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/distrib.Po@am__quote@
......
......@@ -6,6 +6,10 @@
int main(void)
{
printf("Number of bits with single precision: %u\n",
fpngl_t32);
printf("Number of bits with double precision: %u\n",
fpngl_t64);
printf("Number of digits with single precision: %u\n",
(uint32_t)trunc(fpngl_t32*log10(2)));
printf("Number of digits with double precision: %u\n",
......
@set UPDATED 12 March 2020
@set UPDATED 17 March 2020
@set UPDATED-MONTH March 2020
@set EDITION 0.0.1
@set VERSION 0.0.1
@set UPDATED 12 March 2020
@set UPDATED 17 March 2020
@set UPDATED-MONTH March 2020
@set EDITION 0.0.1
@set VERSION 0.0.1
......@@ -29,7 +29,7 @@
/*
Defining constants according to:
Accuracy and Stability of Numerical Algorithms.
Nicholas J. Higham. SIAM
Nicholas J. Higham. SIAM, 2002. ISBN: 978-0-89871-521-7
*/
// Double precision format **************************************
......
......@@ -46,7 +46,9 @@ static const double fpngl_u64 = 0.5*DBL_EPSILON;
static const double fpngl_lambda64 = DBL_MIN;
// Smallest positive double precision number
static const double fpngl_mu64 = DBL_MIN*DBL_EPSILON;
// Not-a-Number
static const double fpngl_NaN64 = 0.0/0.0;
// Infinity
static const double fpngl_infinity64 = 1.0/0.0;
// Largest positive finite double precision number
static const double fpngl_max64 = DBL_MAX;
......
......@@ -25,7 +25,7 @@
#include <fpnglib/constants64.h>
#include <fpnglib/float64.h>
#include <fpnglib/irange.h>
#include <fpnglib/types.h>
fpngl_irng_t *fpngl_distribution_rng_internal(fpngl_distribution_t *dd);
......@@ -160,7 +160,10 @@ double fpngl_nan64(fpngl_irng_t *rng)
{
fpngl_uintf64_t di;
uint64_t sign = randsign(rng);
uint64_t fract = randfrac(rng);
uint64_t fract;
do {
fract = randfrac(rng);
} while (fract == 0);
di.ui = (sign << 63) | ((uint64_t)0x7ff << 52) | fract;
return di.d;
......@@ -174,6 +177,8 @@ double fpngl_float64(fpngl_irng_t *rng,
uint64_t andmask,
uint64_t ormask)
{
assert(minexp<=2047 && maxexp<=2047);
assert(minfrac<= 0x000fffffffffffff && maxfrac<=0x000fffffffffffff);
fpngl_uintf64_t di;
uint64_t sign;
......@@ -218,6 +223,7 @@ double fpngl_float64_distrib(fpngl_distribution_t *fpd)
uint64_t fpngl_distance_float64(double a, double b)
{
assert(a <= b);
if (isinf(a) || isinf(b) || isnan(a) || isnan(b)) {
return 0; // Error
}
......
......@@ -26,7 +26,7 @@
#include <fpnglib/fpngl_config.h>
#include <stdint.h>
#include <fpnglib/irng_t.h>
#include <fpnglib/types.h>
#include <fpnglib/frng_t.h>
#include <fpnglib/discrete_distribution.h>
#ifdef __cplusplus
......@@ -46,12 +46,14 @@ extern "C" {
double fpngl_previous64(double v, uint64_t n);
// Return a random subnormal double precision number
// in (-fpngl_lambda64, -0.0]\cup[0.0, fpngl_lambda64)
double fpngl_subnormal64(fpngl_irng_t *rng);
// Return +0.0 or -0.0 randomly with even probability
double fpngl_zero64(fpngl_irng_t *rng);
// Return +oo or -oo randomly with even probability
double fpngl_inf64(fpngl_irng_t *rng);
// Return a normal double precision number in [fngl_lambda64,fpngl_max64]
// Return a normal double precision number in
// [-fpngl_max64, -fngl_lambda64]\cup[fngl_lambda64, -fpngl_max64]
double fpngl_normal64(fpngl_irng_t *rng);
// Return a nan
double fpngl_nan64(fpngl_irng_t *rng);
......
......@@ -28,6 +28,8 @@
struct fpngl_frng32_t {
void *state;
const char *name;
float min;
float max;
float(*nextf32)(void*);
void (*next_arrayf32)(void *state, float *T, uint32_t n);
void (*delete)(void*);
......@@ -36,6 +38,7 @@ struct fpngl_frng32_t {
fpngl_frng32_t *fpngl_frng32_new(const char *name,
void *state,
float min, float max,
float (*nextf32)(void*),
void (*next_arrayf32)(void *state,
float *T, uint32_t n),
......@@ -50,6 +53,8 @@ fpngl_frng32_t *fpngl_frng32_new(const char *name,
rng->next_arrayf32 = next_arrayf32;
rng->delete = delete;
rng->seed = seed;
rng->min = min;
rng->max = max;
return rng;
}
......@@ -79,3 +84,14 @@ void fpngl_frng32_next_arrayf32(fpngl_frng32_t *frng, float *T, uint32_t n)
{
frng->next_arrayf32(frng->state,T,n);
}
float fpngl_frng32_min(fpngl_frng32_t *frng)
{
return frng->min;
}
float fpngl_frng32_max(fpngl_frng32_t *frng)
{
return frng->max;
}
......@@ -47,6 +47,7 @@ extern "C" {
*/
fpngl_frng32_t *fpngl_frng32_new(const char *name,
void *state,
float min, float max,
float (*nextf32)(void*),
void (*next_arrayf32)(void *state,
float *T, uint32_t n),
......@@ -66,9 +67,15 @@ extern "C" {
*/
void fpngl_frng32_next_arrayf32(fpngl_frng32_t *frng, float *T, uint32_t n);
// Return the name of the FRNG
const char* fpngl_frng32_name(fpngl_frng32_t *frng);
// Return the minimum value of the domain in which floats are drawn
float fpngl_frng32_min(fpngl_frng32_t *frng);
// Return the maximum value of the domain in which floats are drawn
float fpngl_frng32_max(fpngl_frng32_t *frng);
// Return the seed used by the FRNG. May be `0` if the underlying IRNG
// was initialized with an array and not a unique integer.
uint32_t fpngl_frng32_seed(fpngl_frng32_t *frng);
......
......@@ -93,7 +93,7 @@ fpngl_frng64_t *fpngl_KGW(fpngl_irng_t *irng, uint64_t seed)
state->truncHM = trunc(M/2);
state->eps0 = (M-1+state->truncHM)/(2*sqM);
state->eps1 = 1-(2*M - 1 - state->truncHM)/(2*sqM);
return fpngl_frng64_new("KGW", state,
return fpngl_frng64_new("KGW", state, 0.0, 0x1.fffffffffffffp-1,
(double (*)(void*))KGW_nextf64,
(void (*)(void*, double*, uint32_t))KGW_next_arrayf64,
(void (*)(void*))KGW_delete,
......
......@@ -60,6 +60,7 @@ static void frng64_delete(frng_division_state_t *frngstate)
}
fpngl_frng64_t *fpngl_bydivision_new(const char *name,
double min, double max,
fpngl_irng_t *irng,
uint64_t denominator)
{
......@@ -72,7 +73,7 @@ fpngl_frng64_t *fpngl_bydivision_new(const char *name,
if ((uint64_t)frngstate->denominator != denominator) {
FPNGL_WARNING("Warning: %lu rounded to %.0f\n",denominator,frngstate->denominator);
}
return fpngl_frng64_new(name, frngstate,
return fpngl_frng64_new(name, frngstate, min, max,
(double (*)(void*))nextf64,
(void (*)(void*, double*, uint32_t))next_arrayf64,
(void (*)(void*))frng64_delete,
......@@ -81,19 +82,19 @@ fpngl_frng64_t *fpngl_bydivision_new(const char *name,
fpngl_frng64_t *fpngl_matlabp5(uint64_t seed)
{
return fpngl_bydivision_new("matlabp5",
return fpngl_bydivision_new("matlabp5", 0.0, 0x1.fffffffffffffp-1,
fpngl_irng_new64(fpngl_minstd64(seed)),1UL<<31);
}
fpngl_frng64_t *fpngl_drand48bsd(uint64_t seed)
{
return fpngl_bydivision_new("drand48bsd",
return fpngl_bydivision_new("drand48bsd", 0.0, 0x1.fffffffffffffp-1,
fpngl_irng_new64(fpngl_drand48_lcg64(seed)),1UL<<48);
}
fpngl_frng64_t *fpngl_mupad(uint64_t seed)
{
return fpngl_bydivision_new("mupad",
return fpngl_bydivision_new("mupad", 0.0, 0x1.fffffffffffffp-1,
fpngl_irng_new64(fpngl_mupad_lcg64(seed)),
0xe8d4a50ff5UL);
}
......
......@@ -38,9 +38,10 @@ extern "C" {
The FRNG created owns the IRNG passed as second parameter and is responsible
for its destruction.
*/
fpngl_frng64_t *fpngl_bydivision_new(const char *name,
fpngl_irng_t *irng,
uint64_t denominator);
fpngl_frng64_t *fpngl_bydivision_new(const char *name, double min,
double max,
fpngl_irng_t *irng,
uint64_t denominator);
/*
Double precision floating-point generator used by MATLAB until MATLAB 5.
......
<