From 6a562a77146f3987c3f71796f9f05d270255ad26 Mon Sep 17 00:00:00 2001 From: Dave Parker Date: Tue, 7 Nov 2006 09:12:53 +0000 Subject: [PATCH] Fixed building of Windows DLL to allow intra-library loading. Moved foxglynn.c/h to prism.c/h. git-svn-id: https://www.prismmodelchecker.org/svn/prism/prism/trunk@123 bbc10eb1-c90d-0410-af57-cb519fbb1720 --- prism/include/dv.h | 39 ++++--- prism/include/{foxglynn.h => prism.h} | 22 ++-- prism/src/dv/dv.cc | 28 ++--- prism/src/hybrid/Makefile | 2 +- prism/src/hybrid/PH_StochBoundedUntil.cc | 2 +- prism/src/hybrid/PH_StochCumulReward.cc | 2 +- prism/src/hybrid/PH_StochTransient.cc | 2 +- prism/src/mtbdd/Makefile | 2 +- prism/src/mtbdd/PM_StochBoundedUntil.cc | 2 +- prism/src/mtbdd/PM_StochCumulReward.cc | 2 +- prism/src/mtbdd/PM_StochTransient.cc | 2 +- prism/src/prism/{foxglynn.cc => prism.cc} | 129 +++++++++++----------- prism/src/sparse/Makefile | 2 +- prism/src/sparse/PS_StochBoundedUntil.cc | 2 +- prism/src/sparse/PS_StochCumulReward.cc | 2 +- prism/src/sparse/PS_StochTransient.cc | 2 +- 16 files changed, 123 insertions(+), 119 deletions(-) rename prism/include/{foxglynn.h => prism.h} (72%) rename prism/src/prism/{foxglynn.cc => prism.cc} (54%) diff --git a/prism/include/dv.h b/prism/include/dv.h index f9f48d28..32ad2e79 100644 --- a/prism/include/dv.h +++ b/prism/include/dv.h @@ -36,6 +36,13 @@ #include #include +// Flags for building Windows DLLs +#ifdef __MINGW32__ + #define EXPORT __declspec(dllexport) +#else + #define EXPORT +#endif + // constants #define DV_GREATER_THAN_EQUALS 1 @@ -57,22 +64,22 @@ struct DistVector // function prototypes -double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd); -double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd, double *res); -DdNode *double_vector_to_mtbdd(DdManager *ddman, double *vec, DdNode **vars, int num_vars, ODDNode *odd); -DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound, DdNode **vars, int num_vars, ODDNode *odd); -DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound1, double bound2, DdNode **vars, int num_vars, ODDNode *odd); - -void filter_double_vector(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); -double get_first_from_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); -double min_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); -double max_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); -double sum_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); -double sum_double_vector_over_mtbdd(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, ODDNode *odd); - -DistVector *double_vector_to_dist(double *v, int n); -void free_dist_vector(DistVector *&dv); -void free_dv_or_dist_vector(double *&v, DistVector *&dv); +EXPORT double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd, double *res); +EXPORT DdNode *double_vector_to_mtbdd(DdManager *ddman, double *vec, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound1, double bound2, DdNode **vars, int num_vars, ODDNode *odd); + +EXPORT void filter_double_vector(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double get_first_from_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double min_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double max_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double sum_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd); +EXPORT double sum_double_vector_over_mtbdd(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, ODDNode *odd); + +EXPORT DistVector *double_vector_to_dist(double *v, int n); +EXPORT void free_dist_vector(DistVector *&dv); +EXPORT void free_dv_or_dist_vector(double *&v, DistVector *&dv); //------------------------------------------------------------------------------ diff --git a/prism/include/foxglynn.h b/prism/include/prism.h similarity index 72% rename from prism/include/foxglynn.h rename to prism/include/prism.h index 78409811..f68599dc 100644 --- a/prism/include/foxglynn.h +++ b/prism/include/prism.h @@ -1,13 +1,6 @@ //============================================================================== // -// File: foxglynn.h -// Date: 14/3/02 -// Author: Dave Parker + Joachim Meyer-Kayser -// Desc: Header for Fox Glynn Poisson probability computation -// -//------------------------------------------------------------------------------ -// -// Copyright (c) 2002, Dave Parker + Joachim Meyer-Kayser +// Copyright (c) 2006, Dave Parker // // This file is part of PRISM. // @@ -27,8 +20,14 @@ // //============================================================================== -// struct +// Flags for building Windows DLLs +#ifdef __MINGW32__ + #define EXPORT __declspec(dllexport) +#else + #define EXPORT +#endif +// Fox-Glynn wieghts struct typedef struct FoxGlynnWeights { int left; @@ -37,8 +36,7 @@ typedef struct FoxGlynnWeights double *weights; } FoxGlynnWeights; -// prototype - -FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, double accuracy); +// Function prototypes +EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, double accuracy); //------------------------------------------------------------------------------ diff --git a/prism/src/dv/dv.cc b/prism/src/dv/dv.cc index 5a54c13c..5ee7c4a8 100644 --- a/prism/src/dv/dv.cc +++ b/prism/src/dv/dv.cc @@ -50,12 +50,12 @@ static double sum_double_vector_over_mtbdd_rec(DdManager *ddman, double *vec, Dd // if not, a new one is created // in either the case, a pointer to the array is returned -double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd) { return mtbdd_to_double_vector(ddman, dd, vars, num_vars, odd, NULL); } -double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd, double *res) +EXPORT double *mtbdd_to_double_vector(DdManager *ddman, DdNode *dd, DdNode **vars, int num_vars, ODDNode *odd, double *res) { int i, n; @@ -99,7 +99,7 @@ void mtbdd_to_double_vector_rec(DdManager *ddman, DdNode *dd, DdNode **vars, int // converts an array of doubles to an mtbdd -DdNode *double_vector_to_mtbdd(DdManager *ddman, double *vec, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT DdNode *double_vector_to_mtbdd(DdManager *ddman, double *vec, DdNode **vars, int num_vars, ODDNode *odd) { return double_vector_to_mtbdd_rec(ddman, vec, vars, num_vars, 0, odd, 0); } @@ -139,12 +139,12 @@ DdNode *double_vector_to_mtbdd_rec(DdManager *ddman, double *vec, DdNode **vars, // converts an array of doubles to a bdd using a relational operator and one or more bounds -DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound, DdNode **vars, int num_vars, ODDNode *odd) { return double_vector_to_bdd(ddman, vec, rel_op, bound, 0, vars, num_vars, odd); } -DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound1, double bound2, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT DdNode *double_vector_to_bdd(DdManager *ddman, double *vec, int rel_op, double bound1, double bound2, DdNode **vars, int num_vars, ODDNode *odd) { return double_vector_to_bdd_rec(ddman, vec, rel_op, bound1, bound2, vars, num_vars, 0, odd, 0); } @@ -190,7 +190,7 @@ DdNode *double_vector_to_bdd_rec(DdManager *ddman, double *vec, int rel_op, doub // filter vector using a bdd (set elements not in filter to 0) -void filter_double_vector(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT void filter_double_vector(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) { filter_double_vector_rec(ddman, vec, filter, vars, num_vars, 0, odd, 0); } @@ -220,7 +220,7 @@ double filter_double_vector_rec(DdManager *ddman, double *vec, DdNode *filter, D // get value of first element in BDD filter -double get_first_from_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double get_first_from_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) { // should never be called with an empty filter - we trap this case and return -1 if (filter == Cudd_ReadZero(ddman)) return -1; @@ -251,7 +251,7 @@ double get_first_from_bdd_rec(DdManager *ddman, double *vec, DdNode *filter, DdN // compute the minimum value of those in the bdd passed in -double min_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double min_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) { return min_double_vector_over_bdd_rec(ddman, vec, filter, vars, num_vars, 0, odd, 0); } @@ -287,7 +287,7 @@ double min_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *fil // compute the maximum value of those in the bdd passed in -double max_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double max_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) { return max_double_vector_over_bdd_rec(ddman, vec, filter, vars, num_vars, 0, odd, 0); } @@ -323,7 +323,7 @@ double max_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *fil // sums up the elements of a double array - but only those in the bdd passed in -double sum_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double sum_double_vector_over_bdd(DdManager *ddman, double *vec, DdNode *filter, DdNode **vars, int num_vars, ODDNode *odd) { return sum_double_vector_over_bdd_rec(ddman, vec, filter, vars, num_vars, 0, odd, 0); } @@ -359,7 +359,7 @@ double sum_double_vector_over_bdd_rec(DdManager *ddman, double *vec, DdNode *fil // sums up the elements of a double array - multiplied by the corresponding values in the mtbdd passed in -double sum_double_vector_over_mtbdd(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, ODDNode *odd) +EXPORT double sum_double_vector_over_mtbdd(DdManager *ddman, double *vec, DdNode *mult, DdNode **vars, int num_vars, ODDNode *odd) { return sum_double_vector_over_mtbdd_rec(ddman, vec, mult, vars, num_vars, 0, odd, 0); } @@ -392,7 +392,7 @@ double sum_double_vector_over_mtbdd_rec(DdManager *ddman, double *vec, DdNode *m // and an array of short pointers to them. Fails and returns NULL if there are more distinct values // than can be indexed by a short or if we run out of memory. -DistVector *double_vector_to_dist(double *v, int n) +EXPORT DistVector *double_vector_to_dist(double *v, int n) { double *buffer, *tmp; int i, j, num_dist, buffer_size, buffer_inc; @@ -455,7 +455,7 @@ DistVector *double_vector_to_dist(double *v, int n) // free distinct vector struct -void free_dist_vector(DistVector *&dv) +EXPORT void free_dist_vector(DistVector *&dv) { free(dv->dist); free(dv->ptrs); @@ -467,7 +467,7 @@ void free_dist_vector(DistVector *&dv) // delete double array, distinct vector struct, or both -void free_dv_or_dist_vector(double *&v, DistVector *&dv) +EXPORT void free_dv_or_dist_vector(double *&v, DistVector *&dv) { if (v) { free(v); v = NULL; } if (dv) { free_dist_vector(dv); dv = NULL; } diff --git a/prism/src/hybrid/Makefile b/prism/src/hybrid/Makefile index 4694fa79..8821c321 100644 --- a/prism/src/hybrid/Makefile +++ b/prism/src/hybrid/Makefile @@ -48,7 +48,7 @@ $(PRISM_DIR_REL)/$(INCLUDE_DIR)/PrismHybrid.h: $(PRISM_DIR_REL)/$(CLASSES_DIR)/$ ($(JAVAH) -classpath $(PRISM_DIR_REL)/$(CLASSES_DIR) -jni -o $@ $(THIS_DIR).PrismHybrid; touch $@) $(PRISM_DIR_REL)/$(LIB_DIR)/$(LIBPREFIX)prismhybrid$(LIBSUFFIX): $(O_FILES) - $(LD) $(SHARED) $(LDFLAGS) -o $@ $(O_FILES) $(PRISM_DIR_REL)/$(OBJ_DIR)/prism/foxglynn.o $(PRISM_DIR_REL)/$(OBJ_DIR)/dv/dv.o $(LIBRARIES) + $(LD) $(SHARED) $(LDFLAGS) -o $@ $(O_FILES) $(LIBRARIES) $(PRISM_DIR_REL)/$(OBJ_DIR)/$(THIS_DIR)/%.o: %.cc $(CPP) $(CPPFLAGS) -c $< -o $@ $(INCLUDES) diff --git a/prism/src/hybrid/PH_StochBoundedUntil.cc b/prism/src/hybrid/PH_StochBoundedUntil.cc index addcc1f0..c3f0a91f 100644 --- a/prism/src/hybrid/PH_StochBoundedUntil.cc +++ b/prism/src/hybrid/PH_StochBoundedUntil.cc @@ -28,7 +28,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "hybrid.h" #include "PrismHybridGlob.h" diff --git a/prism/src/hybrid/PH_StochCumulReward.cc b/prism/src/hybrid/PH_StochCumulReward.cc index d3ca5575..cc3fcb5f 100644 --- a/prism/src/hybrid/PH_StochCumulReward.cc +++ b/prism/src/hybrid/PH_StochCumulReward.cc @@ -28,7 +28,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "hybrid.h" #include "PrismHybridGlob.h" diff --git a/prism/src/hybrid/PH_StochTransient.cc b/prism/src/hybrid/PH_StochTransient.cc index 91cd24be..2338eeae 100644 --- a/prism/src/hybrid/PH_StochTransient.cc +++ b/prism/src/hybrid/PH_StochTransient.cc @@ -28,7 +28,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "hybrid.h" #include "PrismHybridGlob.h" diff --git a/prism/src/mtbdd/Makefile b/prism/src/mtbdd/Makefile index 23670f1e..10c4c34a 100644 --- a/prism/src/mtbdd/Makefile +++ b/prism/src/mtbdd/Makefile @@ -47,7 +47,7 @@ $(PRISM_DIR_REL)/$(INCLUDE_DIR)/PrismMTBDD.h: $(PRISM_DIR_REL)/$(CLASSES_DIR)/$( ($(JAVAH) -classpath $(PRISM_DIR_REL)/$(CLASSES_DIR) -jni -o $@ $(THIS_DIR).PrismMTBDD; touch $@) $(PRISM_DIR_REL)/$(LIB_DIR)/$(LIBPREFIX)prismmtbdd$(LIBSUFFIX): $(O_FILES) - $(LD) $(SHARED) $(LDFLAGS) -o $@ $(O_FILES) $(PRISM_DIR_REL)/$(OBJ_DIR)/prism/foxglynn.o $(LIBRARIES) + $(LD) $(SHARED) $(LDFLAGS) -o $@ $(O_FILES) $(LIBRARIES) $(PRISM_DIR_REL)/$(OBJ_DIR)/$(THIS_DIR)/%.o: %.cc $(CPP) $(CPPFLAGS) -c $< -o $@ $(INCLUDES) diff --git a/prism/src/mtbdd/PM_StochBoundedUntil.cc b/prism/src/mtbdd/PM_StochBoundedUntil.cc index 967bb76a..3acd25d3 100644 --- a/prism/src/mtbdd/PM_StochBoundedUntil.cc +++ b/prism/src/mtbdd/PM_StochBoundedUntil.cc @@ -27,7 +27,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "PrismMTBDDGlob.h" //------------------------------------------------------------------------------ diff --git a/prism/src/mtbdd/PM_StochCumulReward.cc b/prism/src/mtbdd/PM_StochCumulReward.cc index 176d24d0..aa1c2354 100644 --- a/prism/src/mtbdd/PM_StochCumulReward.cc +++ b/prism/src/mtbdd/PM_StochCumulReward.cc @@ -27,7 +27,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "PrismMTBDDGlob.h" //------------------------------------------------------------------------------ diff --git a/prism/src/mtbdd/PM_StochTransient.cc b/prism/src/mtbdd/PM_StochTransient.cc index d339fc75..c18b3f24 100644 --- a/prism/src/mtbdd/PM_StochTransient.cc +++ b/prism/src/mtbdd/PM_StochTransient.cc @@ -27,7 +27,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "PrismMTBDDGlob.h" //------------------------------------------------------------------------------ diff --git a/prism/src/prism/foxglynn.cc b/prism/src/prism/prism.cc similarity index 54% rename from prism/src/prism/foxglynn.cc rename to prism/src/prism/prism.cc index 1aceb23d..ad4711a6 100644 --- a/prism/src/prism/foxglynn.cc +++ b/prism/src/prism/prism.cc @@ -1,6 +1,6 @@ //============================================================================== // -// Copyright (c) 2002-2004, Joachim Meyer-Kayser, Dave Parker +// Copyright (c) 2002-2006, Dave Parker, Joachim Meyer-Kayser // // This file is part of PRISM. // @@ -21,15 +21,16 @@ //============================================================================== // includes -#include "foxglynn.h" +#include "prism.h" #include #include //------------------------------------------------------------------------------ -// this function was written by joachim meyer-kayser (converted from java) +// Compute poisson probabilities for uniformisation (Fox-Glynn method) +// NB: This function was written by Joachim Meyer-Kayser (converted from Java) -FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, double accuracy) +EXPORT FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, double accuracy) { int m; double q; @@ -40,104 +41,103 @@ FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, doub { double m2 = m; double k; - + if (q_tmax == 0.0) { printf("Overflow: TA parameter qtmax = time * maxExitRate = 0."); } if (q_tmax < 25.0) { - fgw.left = 0; + fgw.left = 0; } if (q_tmax < 400.0) { - // Find right using Corollary 1 with q_tmax=400 - double sqrt2 = sqrt(2.0); + // Find right using Corollary 1 with q_tmax=400 + double sqrt2 = sqrt(2.0); double sqrtl = 20; - double a = 1.0025 * exp (0.0625) * sqrt2; - double b = 1.0025 * exp (0.125/400); //exp (0.0003125) - double startk = 1.0/(2.0 * sqrt2 * 400); - double stopk = sqrtl/(2*sqrt2); - - for (k = startk; k <= stopk; k += 3.0) { - double d = 1.0/(1 - exp ((-2.0/9.0)*(k*sqrt2*sqrtl + 1.5))); - double f = a * d * exp (-0.5*k*k) / (k * sqrt (2.0 * 3.1415926)); - - if (f <= accuracy/2.0) - break; - } - + double a = 1.0025 * exp (0.0625) * sqrt2; + double b = 1.0025 * exp (0.125/400); //exp (0.0003125) + double startk = 1.0/(2.0 * sqrt2 * 400); + double stopk = sqrtl/(2*sqrt2); + + for (k = startk; k <= stopk; k += 3.0) { + double d = 1.0/(1 - exp ((-2.0/9.0)*(k*sqrt2*sqrtl + 1.5))); + double f = a * d * exp (-0.5*k*k) / (k * sqrt (2.0 * 3.1415926)); + + if (f <= accuracy/2.0) + break; + } + if (k > stopk) k = stopk; fgw.right = (int) ceil(m2 + k*sqrt2*sqrtl + 1.5); - } - + } if (q_tmax >= 400.0) { - // Find right using Corollary 1 using actual q_tmax - double sqrt2 = sqrt (2.0); - double sqrtl = sqrt (q_tmax); - double a = (1.0 + 1.0/q_tmax) * exp (0.0625) * sqrt2; - double b = (1.0 + 1.0/q_tmax) * exp (0.125/q_tmax); - double startk = 1.0/(2.0 * sqrt2 * q_tmax); - double stopk = sqrtl/(2*sqrt2); - - for (k = startk; k <= stopk; k += 3.0) { - double d = 1.0/(1 - exp ((-2.0/9.0)*(k*sqrt2*sqrtl + 1.5))); - double f = a * d * exp (-0.5*k*k) / (k * sqrt (2.0 * 3.1415926)); - - if (f <= accuracy/2.0) - break; - } + // Find right using Corollary 1 using actual q_tmax + double sqrt2 = sqrt (2.0); + double sqrtl = sqrt (q_tmax); + double a = (1.0 + 1.0/q_tmax) * exp (0.0625) * sqrt2; + double b = (1.0 + 1.0/q_tmax) * exp (0.125/q_tmax); + double startk = 1.0/(2.0 * sqrt2 * q_tmax); + double stopk = sqrtl/(2*sqrt2); + + for (k = startk; k <= stopk; k += 3.0) { + double d = 1.0/(1 - exp ((-2.0/9.0)*(k*sqrt2*sqrtl + 1.5))); + double f = a * d * exp (-0.5*k*k) / (k * sqrt (2.0 * 3.1415926)); + if (f <= accuracy/2.0) + break; + } + if (k > stopk) k = stopk; - + fgw.right = (int) ceil(m2 + k*sqrt2*sqrtl + 1.5); - } + } if (q_tmax >= 25.0) { - // Find left using Corollary 2 using actual q_tmax - double sqrt2 = sqrt (2.0); - double sqrtl = sqrt (q_tmax); - double a = (1.0 + 1.0/q_tmax) * exp (0.0625) * sqrt2; - double b = (1.0 + 1.0/q_tmax) * exp (0.125/q_tmax); - double startk = 1.0/(sqrt2*sqrtl); - double stopk = (m2 - 1.5)/(sqrt2*sqrtl); - - for (k = startk; k <= stopk; k += 3.0) { - if (b * exp(-0.5*k*k)/(k * sqrt (2.0 * 3.1415926)) <= accuracy/2.0) - break; + // Find left using Corollary 2 using actual q_tmax + double sqrt2 = sqrt (2.0); + double sqrtl = sqrt (q_tmax); + double a = (1.0 + 1.0/q_tmax) * exp (0.0625) * sqrt2; + double b = (1.0 + 1.0/q_tmax) * exp (0.125/q_tmax); + double startk = 1.0/(sqrt2*sqrtl); + double stopk = (m2 - 1.5)/(sqrt2*sqrtl); + + for (k = startk; k <= stopk; k += 3.0) { + if (b * exp(-0.5*k*k)/(k * sqrt (2.0 * 3.1415926)) <= accuracy/2.0) + break; } - + if (k > stopk) k = stopk; - + fgw.left = (int) floor(m2 - k*sqrtl - 1.5); - } - - if (fgw.left < 0) { - fgw.left = 0; - //printf("Weighter: negative left truncation point found. Ignored.\n"); - } + } + + if (fgw.left < 0) { + fgw.left = 0; + //printf("Weighter: negative left truncation point found. Ignored.\n"); + } - q = overflow / (pow(10.0, 10.0) * (fgw.right - fgw.left)); + q = overflow / (pow(10.0, 10.0) * (fgw.right - fgw.left)); } fgw.weights = new double[fgw.right-fgw.left+1]; fgw.weights[m-fgw.left] = q; - + // down for (int j=m; j>fgw.left; j--) { fgw.weights[j-1-fgw.left] = (j/q_tmax) * fgw.weights[j-fgw.left]; } - + //up if (q_tmax < 400) { if (fgw.right > 600) { printf("Overflow: right truncation point > 600."); } - + for (int j=m; j underflow/q) { fgw.weights[j+1-fgw.left] = q * fgw.weights[j-fgw.left]; j++; @@ -146,7 +146,6 @@ FoxGlynnWeights fox_glynn(double q_tmax, double underflow, double overflow, doub fgw.right = j; } } - } else { for (int j=m; j #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "PrismSparseGlob.h" diff --git a/prism/src/sparse/PS_StochCumulReward.cc b/prism/src/sparse/PS_StochCumulReward.cc index f932efdc..bd582ea7 100644 --- a/prism/src/sparse/PS_StochCumulReward.cc +++ b/prism/src/sparse/PS_StochCumulReward.cc @@ -28,7 +28,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "PrismSparseGlob.h" diff --git a/prism/src/sparse/PS_StochTransient.cc b/prism/src/sparse/PS_StochTransient.cc index 5f83e838..6461950d 100644 --- a/prism/src/sparse/PS_StochTransient.cc +++ b/prism/src/sparse/PS_StochTransient.cc @@ -28,7 +28,7 @@ #include #include #include -#include "foxglynn.h" +#include #include "sparse.h" #include "PrismSparseGlob.h"