From a60c84465f55f8bbc3a8c3e684254e3df1cc1946 Mon Sep 17 00:00:00 2001
From: Allanis <allanis@saracraft.net>
Date: Fri, 16 Aug 2013 13:45:42 +0100
Subject: [PATCH] [Change] Now using normal distibution for accuracy for
weapons (clamped).
---
src/rng.c | 125 +++++++++++++++++++++++++++++++++++++++++++++++++++
src/rng.h | 7 +++
src/weapon.c | 8 ++--
3 files changed, 136 insertions(+), 4 deletions(-)
diff --git a/src/rng.c b/src/rng.c
index 114d42d..3e57db6 100644
--- a/src/rng.c
+++ b/src/rng.c
@@ -2,12 +2,14 @@
#include <math.h>
#include <unistd.h>
#include <time.h>
+#include <errno.h>
#ifdef LINUX
#include <sys/time.h>
#include <fcntl.h>
#endif
#include <SDL.h>
+#include "lephisto.h"
#include "rng.h"
#include "log.h"
@@ -114,3 +116,126 @@ double randfp(void) {
return m / m_div;
}
+/*
+ * Calculate N(x) where N is the normal distribution.
+ *
+ * Approximates to a power series:
+ *
+ * N(x) = 1 - n(x)*(b1*t + b2*t^2 + b3*t^3 + b4*t^4 + b5*t^5) + Err
+ * where t = 1 / (1 + 0.2316419*x)
+ *
+ * Maximum absolute error is 7.5e^-8.
+ */
+double Normal(double x) {
+ double t;
+ double series;
+ const double b1 = 0.319381530;
+ const double b2 = -0.356563782;
+ const double b3 = 1.781477937;
+ const double b4 = -1.821255978;
+ const double b5 = 1.330274429;
+ const double p = 0.2316419;
+ const double c = 0.39894228;
+
+ t = 1. / (1. + p * FABS(x));
+ series = (1. - c * exp( -x * x / 2.) * t *
+ (t*(t*(t*(t*b5 + b4) + b3) + b2) + b1));
+
+ return (x > 0.) ? 1. - series : series;
+}
+
+/*
+ * Lower tail quantile for standard disribution function.
+ *
+ * This function returns an approximation of the inverse cumulative
+ * standard normal distribution system. I.e., given P, it returns
+ * an approximation to the X satisfying P = Pr{Z <= X} Where Z is a
+ * random variable from the standard normal distribution.
+ *
+ * The algorithm uses a minimax approximation by rational functions and
+ * the result has a relative error whose absolute value is less than
+ * 1.15e-9.
+ */
+/* Coefficients in rational approximations. */
+static const double a[] = {
+ -3.969683028665376e+01,
+ 2.209460984245205e+02,
+ -2.759285104469687e+02,
+ 1.383577518672690e+02,
+ -3.066479806614716e+01,
+ 2.506628277459239e+00
+};
+
+static const double b[] = {
+ -5.447609879822406e+01,
+ 1.615858368580409e+02,
+ -1.556989798598866e+02,
+ 6.680131188771972e+01,
+ -1.328068155288572e+01
+};
+
+static const double c[] = {
+ -7.784894002430293e-03,
+ -3.223964580411365e-01,
+ -2.400758277161838e+00,
+ -2.549732539343734e+00,
+ 4.374664141464968e+00,
+ 2.938163982698783e+00
+};
+
+static const double d[] = {
+ 7.784695709041462e-03,
+ 3.224671290700398e-01,
+ 2.445134137142996e+00,
+ 3.754408661907416e+00
+};
+
+#define LOW 0.02425
+#define HIGH 0.97575
+
+double NormalInverse( double p ) {
+ double x, e, u, q, r;
+
+ /* Check for errors. */
+ errno = 0;
+ if((p < 0) || (p > 1)) {
+ errno = EDOM;
+ return 0.;
+ }
+ else if(p == 0.) {
+ errno = ERANGE;
+ return -HUGE_VAL /* minus "infinity". */;
+ }
+ else if(p == 1.) {
+ errno = ERANGE;
+ return HUGE_VAL /* "infinity". */;
+ }
+ /* Use different aproximations for different parts. */
+ else if(p < LOW) {
+ /* Rational approximation for lower region. */
+ q = sqrt(-2*log(p));
+ x = (((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
+ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
+ }
+ else if(p > HIGH) {
+ /* Rational approximation for upper region. */
+ q = sqrt(-2*log(1-p));
+ x = -(((((c[0]*q+c[1])*q+c[2])*q+c[3])*q+c[4])*q+c[5]) /
+ ((((d[0]*q+d[1])*q+d[2])*q+d[3])*q+1);
+ }
+ else {
+ /* Rational approximation for central region. */
+ q = p - 0.5;
+ r = q*q;
+ x = (((((a[0]*r+a[1])*r+a[2])*r+a[3])*r+a[4])*r+a[5])*q /
+ (((((b[0]*r+b[1])*r+b[2])*r+b[3])*r+b[4])*r+1);
+ }
+
+ /* Full machine precision. */
+ e = 0.5 * erfc(-x / M_SQRT2) - p;
+ u = e * 2.5066282746310002 /* sqrt(2*pi) */ * exp((x*x)/2);
+ x = x - u/(1 + x*u/2);
+
+ return x;
+}
+
diff --git a/src/rng.h b/src/rng.h
index ba9fc12..355ab2c 100644
--- a/src/rng.h
+++ b/src/rng.h
@@ -4,7 +4,14 @@
#define RNG_SANE(L,H) ((int)L + (int)((double)(H-L+1) * randfp())) /* L <= RNG <= H */
#define RNGF() (randfp()) /* 0. <= RNGF <= 1. */
+/* Init. */
void rng_init(void);
+
+/* Random functions. */
unsigned int randint(void);
double randfp(void);
+/* Probability functions. */
+double Normal(double x);
+double NormalInverse(double p);
+
diff --git a/src/weapon.c b/src/weapon.c
index 89df1e8..d9f32fe 100644
--- a/src/weapon.c
+++ b/src/weapon.c
@@ -426,11 +426,11 @@ static Weapon* weapon_create(const Outfit* outfit, const double dir, const Vec2*
vect_cset(&v, x, y);
rdir = VANGLE(v);
- rdir += RNG(-outfit->u.blt.accuracy/2.,
- outfit->u.blt.accuracy/2.)/180.*M_PI;
} else /* Fire straight. */
- rdir = dir + RNG(-outfit->u.blt.accuracy/2.,
- outfit->u.blt.accuracy/2.)/180.*M_PI;
+ rdir = dir;
+
+ rdir += NormalInverse(RNGF()*0.9 + 0.05) /* Get rid of extreme values. */
+ * outfit->u.blt.accuracy/2. * 1./180.0*M_PI;
if((rdir > 2.*M_PI) || (rdir < 0.)) rdir = fmod(rdir, 2.*M_PI);