[Change] Now using normal distibution for accuracy for weapons (clamped).

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;
+}
+

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);
+

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);