Gaussian sampler crashes randomly with ArrayIndexOutOfBoundsException

[wasowski]

The following invocation in console, spire 0.18, scala 3.3.1:

val rng = spire.random.rng.SecureJava.apply

spire.random.Gaussian (0.0, 1.0)
  .toLazyList (rng)
  .take(500)
  .toList

Call it a few times, and eventually I am getting:

scala> spire.random.Gaussian(0.0, 1.0).toLazyList(rng).take(500).toList
java.lang.ArrayIndexOutOfBoundsException: Index -1 out of bounds for length 128
  at spire.random.Ziggurat$.loop$1(Ziggurat.scala:83)
  at spire.random.Ziggurat$.nfix(Ziggurat.scala:92)
  at spire.random.Ziggurat$.rnor(Ziggurat.scala:52)
  at spire.random.GaussianInstances.spire$random$GaussianInstances$$anon$2$$_$apply$$anonfun$2(Gaussian.scala:50)
  at spire.random.DistFromGen.apply(Dist.scala:177)
  at spire.random.Dist.toLazyList$$anonfun$1(Dist.scala:127)
  at scala.collection.immutable.LazyList$.$anonfun$continually$1(LazyList.scala:1220)
  at scala.collection.immutable.LazyList.scala$collection$immutable$LazyList$$state$lzycompute(LazyList.scala:259)
  at scala.collection.immutable.LazyList.scala$collection$immutable$LazyList$$state(LazyList.scala:252)
  at scala.collection.immutable.LazyList.isEmpty(LazyList.scala:269)
  at scala.collection.immutable.LazyList.$anonfun$takeImpl$1(LazyList.scala:693)
  at scala.collection.immutable.LazyList.scala$collection$immutable$LazyList$$state$lzycompute(LazyList.scala:259)
  at scala.collection.immutable.LazyList.scala$collection$immutable$LazyList$$state(LazyList.scala:252)
  at scala.collection.immutable.LazyList.isEmpty(LazyList.scala:269)
  at scala.collection.immutable.LazyList$LazyIterator.hasNext(LazyList.scala:1250)
  at scala.collection.immutable.List.prependedAll(List.scala:155)
  at scala.collection.IterableOnceOps.toList(IterableOnce.scala:1288)
  at scala.collection.IterableOnceOps.toList$(IterableOnce.scala:1288)
  at scala.collection.AbstractIterable.toList(Iterable.scala:933)
  ... 42 elided

It seems to happen faster and more often if I increase the number 500.

[wasowski]

I have dug into this deeper. The exception is thrown by the implementation of the Ziggurat algorithm, where some suspicious changes have been made. In the current version, line 83 fails if iz==0:

spire/core/src/main/scala/spire/random/Ziggurat.scala

Lines 74 to 83 in 4a70e3a

	if (iz == 0) {
	while ({
	x = -log(g.nextDouble()) * r1
	y = -log(g.nextDouble())
	y + y < x * x
	})
	return if (hz > 0) r + x else -r - x
	}
	if (fn(iz) + g.nextDouble() * (fn(iz - 1) - fn(iz)) < exp(-.5 * x * x)) return x

The code in lines 74-81 seems to be aimed at preventing that iz==0, but this only works if the while loop is entered. This is in general a rather unusual piece of code, as it acts like an if-condition not a wihle-loop. It is unlikely to have been meant this way.

I dug into the history, and indeed, in this d31279d, more specifically here

d31279d#diff-e1c16b40daf0209587b057492d4d8f5a81ce2d3e1de5e1f1085fdfbd3f4a0474L73-R80

if (iz == 0) {
         while ({
           ...
-        }) () 
-        return if (hz > 0) r + x else -r - x
+        })
+         return if (hz > 0) r + x else -r - x
}

The empty parentheses after the loop condition have been removed and the return line indented. This has changed the semantics from the loop not having a body, and then returning after exit, to an if-condition that returns if satisfied and not otherwise. This will lead to a crash if y + y >= x * x in line 83

This seems to be a large commit with mechanical clean up changes, which increases my suspicion that the change was not intended. The details of this library implementation are alien to me, but I will try to see whether I am able to patch, build and test this hypothesis.