PetriDish: Multi-threading for performance in C#

#define SequentialX
#define PFXX

#define DrawImage

#define UseDelay
#define UseBufferThread
#define FlipFlop
#define StopWhenFirstThreadFinishedX
#define UseBalancing
#define UseTwoBitmapsX

#define WaitForPaint
#define WaitForABitX

#define Radius1X
#define Radius2
#define Radius3X
#define Radius4X

#define UnwindLoop
#define UseWeights

#define CalculatePopulation
#define AddDummyWorkX

using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Drawing;
using System.Drawing.Imaging;
using System.IO;
using System.Threading;
using System.Threading.Tasks;
using System.Runtime.InteropServices;

namespace PetriDish
{
	static class Engine
	{

#if Sequential
		public static readonly int ThreadCount = 1;
#else
		public static readonly int ThreadCount = Environment.ProcessorCount;
		//public static readonly int ThreadCount = 2;
		//public static readonly int ThreadCount = Environment.ProcessorCount - ( Environment.ProcessorCount > 2 ? 1 : 0 );

#endif

		public static readonly int X = ( ThreadCount > 4 ? 800 : 500 );
		public static readonly int Y = ( ThreadCount > 2 ? X : X / 2 );
		public static readonly int XY = X * Y;

#if Radius1
		const int R = 1;
#elif Radius2
		const int R = 2;
#elif Radius3
		const int R = 3;
#elif Radius4
		const int R = 4;
#endif
		const int R2 = ( 2 * R ) + 1;
		const int R2Squared = R2 * R2;

		static int Instability = 5;
		static int InstabilityDivision = ( X / 100 ) + ( Y / 100 );
		static int Instabilities = InstabilityDivision * InstabilityDivision * 3 / 4;

		const int OasisSize = 5;
		const int OasisCount = 5;

		const int BalancingChunk = 100;

		public static Bitmap Bitmap { get { return _Bitmap; } }
		public static readonly object BitmapLock = new object();
		public static event EventHandler Update;
		public static long PixelsPerSecond = 0;


		static Thread _Thread = null;
		static Thread _BufferThread = null;
		static AutoResetEvent _UIEvent = new AutoResetEvent( true );
		static AutoResetEvent _ToggleEvent = new AutoResetEvent( true );
		static AutoResetEvent _BufferEvent = new AutoResetEvent( false );

		// Manual
		static Common.TLS<TLS> _ThreadLocalState = new Common.TLS<TLS>();
		static Thread[] _Threads = null;
		static int[] _ThreadStarts = null;
		static int[] _ThreadEnds = null;
		static Semaphore _ThreadSemaphore = null;
		static AutoResetEvent _ThreadEvent = null;
		static int _ThreadCounter = 0;

#pragma warning disable 429
		static TaskManagerPolicy _TaskManagerPolicy = ThreadCount <= 0 ? new TaskManagerPolicy() : new TaskManagerPolicy( ThreadCount, ThreadCount, ThreadCount );
		static TaskManager _TaskManager = new TaskManager( _TaskManagerPolicy );
#pragma warning restore 429

		static Stopwatch _Stopwatch = new Stopwatch();

		static Rectangle sRAll = new Rectangle( 0, 0, X, Y );

		static Bitmap _Bitmap = new Bitmap( X, Y, PixelFormat.Format32bppArgb );
		static Bitmap _Scratch = new Bitmap( X, Y, PixelFormat.Format32bppArgb );
		static Graphics _BitmapGraphics = Graphics.FromImage( _Bitmap );

		static Int32[] _BufferA = new Int32[ XY ];
		static Int32[] _BufferB = new Int32[ XY ];
		static byte[] _Delay = new byte[ XY ];

		static volatile bool _BufferToggle = false;
		static Int32[] ReadBuffer { get { return ( _BufferToggle ? _BufferA : _BufferB ); } }
		static Int32[] WriteBuffer { get { return ( _BufferToggle ? _BufferB : _BufferA ); } }

		static Weights _Weights = null;

		static Engine()
		{
			InitWeights();
			InitDelay();

			_Stopwatch.Start();

			_BufferThread = new Thread( BufferMaster );
			_BufferThread.Name = "BufferThread";
			_BufferThread.IsBackground = true;
			//_BufferThread.Priority = ThreadPriority.AboveNormal;
			_BufferThread.Start();

			if ( ThreadCount > 0 )
			{
				_Threads = new Thread[ ThreadCount ];
				_ThreadStarts = new int[ ThreadCount ];
				_ThreadEnds = new int[ ThreadCount ];
				_ThreadSemaphore = new Semaphore( 0, ThreadCount );
				_ThreadEvent = new AutoResetEvent( false );

				int block = XY / ThreadCount;
				int start = 0, end = block;

				for ( int i = 0 ; i < ThreadCount ; i++ )
				{
					_ThreadStarts[ i ] = 1 * start;
					_ThreadEnds[ i ] = 1 * ( i < ThreadCount - 1 ? end : XY );

					start = end;
					end += block;

					Thread thread = new Thread( ManualWork );
					thread.Name = "Manual thread #" + i;
					thread.IsBackground = true;
					thread.Priority = ThreadPriority.BelowNormal;
					thread.Start( i );

					_Threads[ i ] = thread;
				}
			}
		}

		//-----------------------------------------------------------------------------------------
		// UI public methods

		public static void Start()
		{
			if ( _Thread != null ) return;

			_Thread = new Thread( Master );
			_Thread.Name = "Manual thread master";
			_Thread.IsBackground = true;
			_Thread.Start();
		}

		public static void Stop()
		{
			if ( _Threads != null )
			{
				foreach ( var thread in _Threads )
					if ( thread != null )
					{
						thread.Abort();
						thread.Join();
					}

				_Threads = null;
			}

			if ( _Thread != null )
			{
				_Thread.Abort();
				_Thread.Join();
				_Thread = null;
			}

			if ( _BufferThread != null )
			{
				_BufferThread.Abort();
				_BufferThread.Join();
				_BufferThread = null;
			}

		}

		public static void Painted()
		{
			_UIEvent.Set();
		}

		//-----------------------------------------------------------------------------------------

		static void BufferMaster()
		{
			for ( ; ; )
				try
				{
					_BufferEvent.WaitOne();

					UpdateUI();
				}
				finally
				{
#if FlipFlop
					_ToggleEvent.Set();
#endif
				}
		}

		static void UpdateUI()
		{

#if DrawImage

			var read = ReadBuffer;


			BitmapData b = null;

#if !UseTwoBitmaps
			lock ( BitmapLock )
#endif
			{
				try
				{
#if UseTwoBitmaps
					b = _Scratch.LockBits( sRAll, ImageLockMode.WriteOnly, _Bitmap.PixelFormat );
#else
					b = _Bitmap.LockBits( sRAll, ImageLockMode.WriteOnly, _Bitmap.PixelFormat );
#endif
				}
				catch ( Exception )
				{
					return;
				}

				try
				{
					Marshal.Copy( read, 0, b.Scan0, read.Length );
				}
				finally
				{
#if UseTwoBitmaps
					_Scratch.UnlockBits( b );
#else
					_Bitmap.UnlockBits( b );
#endif
				}

#if UseTwoBitmaps
				lock ( BitmapLock )
					_BitmapGraphics.DrawImageUnscaled( _Scratch, 0, 0 );
#endif
			}

#endif

			var d = Update; if ( d != null ) d( null, EventArgs.Empty );
		}

		static void Master()
		{
			for ( ; ; )
			{
				TimeSpan start = _Stopwatch.Elapsed;

				Work();

				TimeSpan end = _Stopwatch.Elapsed;
				TimeSpan elapsed = end - start;

				PixelsPerSecond = ( long ) ( XY / elapsed.TotalSeconds );

#if WaitForPaint && WaitForABit
				_UIEvent.WaitOne( 3, false );
#elif WaitForPaint
				_UIEvent.WaitOne();
#endif

				//WaitHandle.WaitAll( new WaitHandle[] { _UIEvent, _ToggleEvent } );

#if UseBufferThread

#if FlipFlop
				_ToggleEvent.WaitOne();
#endif

				_BufferToggle = !_BufferToggle;

				_BufferEvent.Set();

#else

				_BufferToggle = !_BufferToggle;

				UpdateUI();
#endif

#if WaitForABit && !WaitForPaint
				Thread.Sleep( 3 );
#endif
			}
		}

		unsafe static void Work()
		{
#if Sequential
			TLS tls = new TLS();
			fixed ( Int32* read = &ReadBuffer[ 0 ] )
			fixed ( Int32* write = &WriteBuffer[ 0 ] )
				for ( int i = 0 ; i < X * Y ; i += 1 ) Next( i, tls, read, write );
#elif PFX
			fixed ( Int32* read = &ReadBuffer[ 0 ] )
			fixed ( Int32* write = &WriteBuffer[ 0 ] )
				Parallel.For( 0, X * Y, 1, () => _ThreadLocalState.Current,
				( i, state ) => Next( i, state.ThreadLocalState, read, write ),
				null, _TaskManager, TaskCreationOptions.None );
#else
#if StopWhenFirstThreadFinished
			finished = false;
#endif
			_BalancingMaster.Reset();
			_ThreadCounter = ThreadCount;
			_ThreadSemaphore.Release( ThreadCount );
			_ThreadEvent.WaitOne();
#endif
		}

#if StopWhenFirstThreadFinished
		static volatile bool finished = false;
#endif

		//-----------------------------------------------------------------------------------------

		unsafe static void ManualWork( object o )
		{
			int thread = ( int ) o;
			TLS tls = _ThreadLocalState.Current;

			for ( ; ; )
			{
				_ThreadSemaphore.WaitOne();

				tls.Enum.Reset();

				try
				{
					fixed ( Int32* read = &ReadBuffer[ 0 ] )
					fixed ( Int32* write = &WriteBuffer[ 0 ] )
					{
#if !UseBalancing
						int start = _ThreadStarts[ thread ];
						int end = _ThreadEnds[ thread ];

						for ( int i = start ; i < end

#if StopWhenFirstThreadFinished
							&& !finished
#endif

							; i += 1 )

#else

						foreach ( int i in tls.Enum )

#endif

							Next( i, tls, read, write );
					}
				}
				finally
				{
#if StopWhenFirstThreadFinished
					finished = true;
#endif
					if ( Interlocked.Decrement( ref _ThreadCounter ) == 0 ) _ThreadEvent.Set();
				}
			}
		}

		//-----------------------------------------------------------------------------------------
		// TLS

		class TLS
		{
			public Random Random = new Random();
			public float[] Results = new float[ 4 ];
			public byte[] Result = new byte[ 4 ];

			public TLSCore Core = new TLSCore();

			public BalancingEnumerator Enum = new BalancingEnumerator( _BalancingMaster );
		}

		unsafe struct TLSCore
		{
			public fixed Int32 Population[ R2Squared ];
		}


		//-----------------------------------------------------------------------------------------

		unsafe static void Next( int i, TLS tls, Int32* read, Int32* write )
		{
#if AddDummyWork
			for ( int x = 0 ; x < 42000 ; x++ ) ;
#endif

#if UseDelay
			byte d = _Delay[ i ];
			if ( d > 0 && tls.Random.Next( 100 ) < d ) return;
#endif

			byte[] population = Population( read, X, Y, i, tls );

			byte b1 = population[ 0 ];
			byte g1 = population[ 1 ];
			byte r1 = population[ 2 ];
			byte a1 = population[ 3 ];

			int a2 = a1;
			int g2 = ( int ) ( ( g1 * 1.1 ) - ( b1 * 0.1 ) - ( r1 * 0.01 ) + 42 - tls.Random.Next( 85 ) );
			int b2 = ( int ) ( ( b1 * 1.1 ) + ( g1 * 0.1 ) - ( r1 * 0.09 ) - tls.Random.Next( 42 ) );
			int r2 = ( int ) ( ( r1 * 1.1 ) + ( b1 * 0.1 ) - ( g1 * 0.00 ) - tls.Random.Next( 42 ) );

			write[ i ] = ( Cap( b2 ) ) | ( Cap( g2 ) << 8 ) | ( Cap( r2 ) << 16 ) | ( Cap( a2 ) << 24 );
		}

		unsafe static byte[] Population( Int32* argb, int width, int height, int index, TLS tls )
		{
#if CalculatePopulation

			int x = index % width;
			int y = index / width;

			PopulationPopulate( argb, width, height, index, x, y, tls );
			PopulationResults( tls );
			PopulationResult( tls );
			PopulationOasis( width, height, x, y, tls );

#else

			tls.Result[ 0 ] = ( byte ) ( argb[ index ] >> 00 );
			tls.Result[ 1 ] = ( byte ) ( argb[ index ] >> 08 );
			tls.Result[ 2 ] = ( byte ) ( argb[ index ] >> 16 );
			tls.Result[ 3 ] = 0xFF;

#endif

			return tls.Result;
		}

		unsafe static void PopulationPopulate( Int32* argb, int width, int height, int index, int x, int y, TLS tls )
		{
			const Int32 outside = 0x0000C000;

			fixed ( Int32* pop = tls.Core.Population )
			{
				for ( int i = 0 ; i < R2Squared ; i++ )
					pop[ i ] = outside;

				for ( int dx = 0 ; dx < R2 ; dx++ )
				{
					int x2 = x + dx - R;

					if ( x2 < 0 || x2 >= width ) continue;

					for ( int dy = 0 ; dy < R2 ; dy++ )
					{
						int y2 = y + dy - R;

						if ( y2 < 0 || y2 >= height ) continue;

						pop[ dx + ( dy * R2 ) ] = argb[ x2 + ( y2 * width ) ];
					}
				}
			}
		}

		unsafe static void PopulationResults( TLS tls )
		{
			fixed ( Int32* pop = tls.Core.Population )
#if UseWeights
			fixed ( float* weights = _Weights.Vector.Value )
#endif
			{
				float[] results = tls.Results;

				results[ 0 ] = 0;
				results[ 1 ] = 0;
				results[ 2 ] = 0;
				results[ 3 ] = 0;

				for ( int dz = 0 ; dz < 4 - 1 ; dz++ )
				{
					int shift = 8 * dz;
					int mask = 0xFF << shift;

#if !UnwindLoop || Radius4
#if UseWeights
					long total = 0;
					for ( int i = 0 ; i < R2Squared ; i++ )
						total += ( long ) ( ( pop[ i ] & mask ) * weights[ i ] );
					total >>= shift;
					results[ dz ] = total / _Weights.Sum;
#else
					long total = 0;
					for ( int i = 0 ; i < R2Squared ; i++ )
						total += pop[ i ] & mask;
					total >>= shift;
					results[ dz ] = total / R2Squared;
#endif
#elif Radius1
#if UseWeights
					// R == 1
					results[ dz ] =
						(
						(
						( int )
						(
						(
						( ( pop[ 00 ] & mask ) * weights[ 00 ] ) +
						( ( pop[ 01 ] & mask ) * weights[ 01 ] ) +
						( ( pop[ 02 ] & mask ) * weights[ 02 ] )
						)
						 + (
						( ( pop[ 03 ] & mask ) * weights[ 03 ] ) +
						( ( pop[ 04 ] & mask ) * weights[ 04 ] ) +
						( ( pop[ 05 ] & mask ) * weights[ 05 ] )
						)
						 + (
						( ( pop[ 06 ] & mask ) * weights[ 06 ] ) +
						( ( pop[ 07 ] & mask ) * weights[ 07 ] ) +
						( ( pop[ 08 ] & mask ) * weights[ 08 ] )
						)
						) >> shift )
						) / _Weights.Sum;
#elif true
					// R == 1
					results[ dz ] =
						(
						(
						(
						( ( pop[ 00 ] & mask ) ) +
						( ( pop[ 01 ] & mask ) ) +
						( ( pop[ 02 ] & mask ) )
						)
						 + (
						( ( pop[ 03 ] & mask ) ) +
						( ( pop[ 04 ] & mask ) ) +
						( ( pop[ 05 ] & mask ) )
						)
						 + (
						( ( pop[ 06 ] & mask ) ) +
						( ( pop[ 07 ] & mask ) ) +
						( ( pop[ 08 ] & mask ) )
						)
						) >> shift )
						/ R2Squared;
#endif
#elif Radius2
#if UseWeights
					// R == 2
					results[ dz ] =
						(
						(
						( int )
						(
						(
						( pop[ 00 ] & mask ) * weights[ 00 ] +
						( pop[ 01 ] & mask ) * weights[ 01 ] +
						( pop[ 02 ] & mask ) * weights[ 02 ] +
						( pop[ 03 ] & mask ) * weights[ 03 ] +
						( pop[ 04 ] & mask ) * weights[ 04 ]
						)
						 +
						 (
						( pop[ 05 ] & mask ) * weights[ 05 ] +
						( pop[ 06 ] & mask ) * weights[ 06 ] +
						( pop[ 07 ] & mask ) * weights[ 07 ] +
						( pop[ 08 ] & mask ) * weights[ 08 ] +
						( pop[ 09 ] & mask ) * weights[ 09 ]
						)
						 +
						 (
						( pop[ 10 ] & mask ) * weights[ 10 ] +
						( pop[ 11 ] & mask ) * weights[ 11 ] +
						( pop[ 12 ] & mask ) * weights[ 12 ] +
						( pop[ 13 ] & mask ) * weights[ 13 ] +
						( pop[ 14 ] & mask ) * weights[ 14 ]
						)
						 +
						 (
						( pop[ 15 ] & mask ) * weights[ 15 ] +
						( pop[ 16 ] & mask ) * weights[ 16 ] +
						( pop[ 17 ] & mask ) * weights[ 17 ] +
						( pop[ 18 ] & mask ) * weights[ 18 ] +
						( pop[ 19 ] & mask ) * weights[ 19 ]
						)
						 +
						 (
						( pop[ 20 ] & mask ) * weights[ 20 ] +
						( pop[ 21 ] & mask ) * weights[ 21 ] +
						( pop[ 22 ] & mask ) * weights[ 22 ] +
						( pop[ 23 ] & mask ) * weights[ 23 ] +
						( pop[ 24 ] & mask ) * weights[ 24 ]
						)
						)
						) >> shift )
						/ _Weights.Sum;
#else
					// R == 2
					results[ dz ] =
						(
						(
						(
						( ( pop[ 00 ] & mask ) ) +
						( ( pop[ 01 ] & mask ) ) +
						( ( pop[ 02 ] & mask ) ) +
						( ( pop[ 03 ] & mask ) ) +
						( ( pop[ 04 ] & mask ) )
						)
						 +
					(
						( ( pop[ 05 ] & mask ) ) +
						( ( pop[ 06 ] & mask ) ) +
						( ( pop[ 07 ] & mask ) ) +
						( ( pop[ 08 ] & mask ) ) +
						( ( pop[ 09 ] & mask ) )
						)
						 +
					(
						( ( pop[ 10 ] & mask ) ) +
						( ( pop[ 11 ] & mask ) ) +
						( ( pop[ 12 ] & mask ) ) +
						( ( pop[ 13 ] & mask ) ) +
						( ( pop[ 14 ] & mask ) )
						)
						 +
					(
						( ( pop[ 15 ] & mask ) ) +
						( ( pop[ 16 ] & mask ) ) +
						( ( pop[ 17 ] & mask ) ) +
						( ( pop[ 18 ] & mask ) ) +
						( ( pop[ 19 ] & mask ) )
						)
						 +
					(
						( ( pop[ 20 ] & mask ) ) +
						( ( pop[ 21 ] & mask ) ) +
						( ( pop[ 22 ] & mask ) ) +
						( ( pop[ 23 ] & mask ) ) +
						( ( pop[ 24 ] & mask ) )
						)
						) >> shift )
						/ R2Squared;
#endif
#elif Radius3
#if UseWeights
					// R == 3
					results[ dz ] =
						(
						(
						( int )
						(
						(
						( ( pop[ 00 ] & mask ) * weights[ 00 ] ) +
						( ( pop[ 01 ] & mask ) * weights[ 01 ] ) +
						( ( pop[ 02 ] & mask ) * weights[ 02 ] ) +
						( ( pop[ 03 ] & mask ) * weights[ 03 ] ) +
						( ( pop[ 04 ] & mask ) * weights[ 04 ] ) +
						( ( pop[ 05 ] & mask ) * weights[ 05 ] ) +
						( ( pop[ 06 ] & mask ) * weights[ 06 ] )
						)
						 + (
						( ( pop[ 07 ] & mask ) * weights[ 07 ] ) +
						( ( pop[ 08 ] & mask ) * weights[ 08 ] ) +
						( ( pop[ 09 ] & mask ) * weights[ 09 ] ) +
						( ( pop[ 10 ] & mask ) * weights[ 10 ] ) +
						( ( pop[ 11 ] & mask ) * weights[ 11 ] ) +
						( ( pop[ 12 ] & mask ) * weights[ 12 ] ) +
						( ( pop[ 13 ] & mask ) * weights[ 13 ] )
						)
						 + (
						( ( pop[ 14 ] & mask ) * weights[ 14 ] ) +
						( ( pop[ 15 ] & mask ) * weights[ 15 ] ) +
						( ( pop[ 16 ] & mask ) * weights[ 16 ] ) +
						( ( pop[ 17 ] & mask ) * weights[ 17 ] ) +
						( ( pop[ 18 ] & mask ) * weights[ 18 ] ) +
						( ( pop[ 19 ] & mask ) * weights[ 19 ] ) +
						( ( pop[ 20 ] & mask ) * weights[ 20 ] )
						)
						 + (
						( ( pop[ 21 ] & mask ) * weights[ 21 ] ) +
						( ( pop[ 22 ] & mask ) * weights[ 22 ] ) +
						( ( pop[ 23 ] & mask ) * weights[ 23 ] ) +
						( ( pop[ 24 ] & mask ) * weights[ 24 ] ) +
						( ( pop[ 25 ] & mask ) * weights[ 25 ] ) +
						( ( pop[ 26 ] & mask ) * weights[ 26 ] ) +
						( ( pop[ 27 ] & mask ) * weights[ 27 ] )
						)
						 + (
						( ( pop[ 28 ] & mask ) * weights[ 28 ] ) +
						( ( pop[ 29 ] & mask ) * weights[ 29 ] ) +
						( ( pop[ 30 ] & mask ) * weights[ 30 ] ) +
						( ( pop[ 31 ] & mask ) * weights[ 31 ] ) +
						( ( pop[ 32 ] & mask ) * weights[ 32 ] ) +
						( ( pop[ 33 ] & mask ) * weights[ 33 ] ) +
						( ( pop[ 34 ] & mask ) * weights[ 34 ] )
						)
						 + (
						( ( pop[ 35 ] & mask ) * weights[ 35 ] ) +
						( ( pop[ 36 ] & mask ) * weights[ 36 ] ) +
						( ( pop[ 37 ] & mask ) * weights[ 37 ] ) +
						( ( pop[ 38 ] & mask ) * weights[ 38 ] ) +
						( ( pop[ 39 ] & mask ) * weights[ 39 ] ) +
						( ( pop[ 40 ] & mask ) * weights[ 40 ] ) +
						( ( pop[ 41 ] & mask ) * weights[ 41 ] )
						)
						 + (
						( ( pop[ 42 ] & mask ) * weights[ 42 ] ) +
						( ( pop[ 43 ] & mask ) * weights[ 43 ] ) +
						( ( pop[ 44 ] & mask ) * weights[ 44 ] ) +
						( ( pop[ 45 ] & mask ) * weights[ 45 ] ) +
						( ( pop[ 46 ] & mask ) * weights[ 46 ] ) +
						( ( pop[ 47 ] & mask ) * weights[ 47 ] ) +
						( ( pop[ 48 ] & mask ) * weights[ 48 ] )
						)
						)
						) >> shift )
						/ _Weights.Sum;
#else
					// R == 3
					results[ dz ] =
						(
						(
						(
						( ( pop[ 00 ] & mask ) ) +
						( ( pop[ 01 ] & mask ) ) +
						( ( pop[ 02 ] & mask ) ) +
						( ( pop[ 03 ] & mask ) ) +
						( ( pop[ 04 ] & mask ) ) +
						( ( pop[ 05 ] & mask ) ) +
						( ( pop[ 06 ] & mask ) )
						)
						 + (
						( ( pop[ 07 ] & mask ) ) +
						( ( pop[ 08 ] & mask ) ) +
						( ( pop[ 09 ] & mask ) ) +
						( ( pop[ 10 ] & mask ) ) +
						( ( pop[ 11 ] & mask ) ) +
						( ( pop[ 12 ] & mask ) ) +
						( ( pop[ 13 ] & mask ) )
						)
						 + (
						( ( pop[ 14 ] & mask ) ) +
						( ( pop[ 15 ] & mask ) ) +
						( ( pop[ 16 ] & mask ) ) +
						( ( pop[ 17 ] & mask ) ) +
						( ( pop[ 18 ] & mask ) ) +
						( ( pop[ 19 ] & mask ) ) +
						( ( pop[ 20 ] & mask ) )
						)
						 + (
						( ( pop[ 21 ] & mask ) ) +
						( ( pop[ 22 ] & mask ) ) +
						( ( pop[ 23 ] & mask ) ) +
						( ( pop[ 24 ] & mask ) ) +
						( ( pop[ 25 ] & mask ) ) +
						( ( pop[ 26 ] & mask ) ) +
						( ( pop[ 27 ] & mask ) )
						)
						 + (
						( ( pop[ 28 ] & mask ) ) +
						( ( pop[ 29 ] & mask ) ) +
						( ( pop[ 30 ] & mask ) ) +
						( ( pop[ 31 ] & mask ) ) +
						( ( pop[ 32 ] & mask ) ) +
						( ( pop[ 33 ] & mask ) ) +
						( ( pop[ 34 ] & mask ) )
						)
						 + (
						( ( pop[ 35 ] & mask ) ) +
						( ( pop[ 36 ] & mask ) ) +
						( ( pop[ 37 ] & mask ) ) +
						( ( pop[ 38 ] & mask ) ) +
						( ( pop[ 39 ] & mask ) ) +
						( ( pop[ 40 ] & mask ) ) +
						( ( pop[ 41 ] & mask ) )
						)
						 + (
						( ( pop[ 42 ] & mask ) ) +
						( ( pop[ 43 ] & mask ) ) +
						( ( pop[ 44 ] & mask ) ) +
						( ( pop[ 45 ] & mask ) ) +
						( ( pop[ 46 ] & mask ) ) +
						( ( pop[ 47 ] & mask ) ) +
						( ( pop[ 48 ] & mask ) )
						)
						) >> shift )
						/ R2Squared;
#endif
#else

#error Select a radius

#endif
				}
			}
		}

		static void PopulationResult( TLS tls )
		{
			tls.Result[ 0 ] = Cap( tls.Results[ 0 ] );
			tls.Result[ 1 ] = Cap( tls.Results[ 1 ] );
			tls.Result[ 2 ] = Cap( tls.Results[ 2 ] );
			//tls.Result[ 3 ] = Cap( tls.Results[ 3 ] );
			tls.Result[ 3 ] = 0xFF;
		}

		static void PopulationOasis( int width, int height, int x, int y, TLS tls )
		{
			int oasisRegion = width / OasisCount;

			int ox = ( int ) Math.Abs( ( x % oasisRegion ) - ( 0.5 * width / OasisCount ) );
			int oy = ( int ) Math.Abs( ( y % oasisRegion ) - ( 0.5 * height / OasisCount ) );

			byte[] result = tls.Result;
			if ( ox < OasisSize && oy < OasisSize ) result[ 1 ] = Math.Max( result[ 1 ], ( byte ) 42 );
		}

		//-----------------------------------------------------------------------------------------
		// Init

		class Weights
		{
			public float Sum = 0;
			public float[ , ] Weight = new float[ R2, R2 ];

			public Vectors Vector = new Vectors();

			public unsafe struct Vectors
			{
				public fixed float Value[ R2Squared ];
			}
		}

		static void InitWeights()
		{
			_Weights = CalcWeights( 0, 0, 0, 0 );
		}

		unsafe static Weights CalcWeights( int l, int t, int r, int b )
		{
			Weights weights = new Weights();

			fixed ( float* vector = weights.Vector.Value )
			{
				float max = ( float ) Math.Sqrt( 2 * R * R );

				for ( int dx = -R ; dx <= R ; dx++ )
					for ( int dy = -R ; dy <= R ; dy++ )
						if ( l - dx <= R && t - dy <= R && r + dx <= R && b + dy <= R )
						{
							float w1 = ( float ) Math.Sqrt( ( dx * dx ) + ( dy * dy ) );
							float w2 = max;
							if ( w1 > 0 ) w2 = max / w1;

							weights.Weight[ dx + R, dy + R ] = w2;
							vector[ dx + R + ( ( dy + R ) * R2 ) ] = w2;
							weights.Sum += w2;
						}
			}

			return weights;
		}

		static void InitDelay()
		{
			Random random = new Random();

			List<int> list = new List<int>();
			for ( int i = 0 ; i < Instabilities ; i++ ) list.Add( random.Next( InstabilityDivision * InstabilityDivision ) );

			for ( int x = 0 ; x < X ; x++ )
				for ( int y = 0 ; y < Y ; y++ )
				{
					int d = ( x * InstabilityDivision / X ) + ( InstabilityDivision * ( y * InstabilityDivision / Y ) );

					if ( list.Contains( d ) ) _Delay[ x + ( X * y ) ] = ( byte ) ( Instability + random.Next( 2 * Instability ) );
				}
		}

		//-----------------------------------------------------------------------------------------
		// Serialization

		[Serializable]
		class State
		{
			public byte[] Delay;
			public Int32[] Read;
		}

		public static void Save( string filepath )
		{
			var s = new State();
			s.Delay = _Delay;
			s.Read = ReadBuffer;

			using ( var fs = File.OpenWrite( filepath ) )
			{
				var f = new System.Runtime.Serialization.Formatters.Binary.BinaryFormatter();
				f.Serialize( fs, s );
			}
		}

		public static void Load( string filepath )
		{
			using ( var fs = File.OpenRead( filepath ) )
			{
				var f = new System.Runtime.Serialization.Formatters.Binary.BinaryFormatter();
				var s = f.Deserialize( fs ) as State;

				Buffer.BlockCopy( s.Delay, 0, _Delay, 0, Math.Min( s.Delay.Length, _Delay.Length ) );
				Buffer.BlockCopy( s.Read, 0, _BufferA, 0, sizeof( Int32 ) * Math.Min( s.Read.Length, _BufferA.Length ) );
				Buffer.BlockCopy( s.Read, 0, _BufferB, 0, sizeof( Int32 ) * Math.Min( s.Read.Length, _BufferA.Length ) );
			}
		}

		//-----------------------------------------------------------------------------------------
		// Balancing

		static BalancingMaster _BalancingMaster = new BalancingMaster();

		class BalancingMaster
		{
			public static readonly int MAX = XY;

			volatile int _Current = 0;
			object _Key = new object();

			public void Reset()
			{
				_Current = 0;
			}

			public int Next
			{
				get
				{
					lock ( _Key )
					{
						if ( _Current >= MAX ) return -1;

						int current = _Current;
						_Current += BalancingChunk;
						return current;
					}
				}
			}
		}

		class BalancingEnumerator : IEnumerable<int>, IEnumerator<int>
		{
			BalancingMaster _Master = null;
			int _Base = 0;
			int _Index = BalancingChunk;
			bool _Finished = false;

			public BalancingEnumerator( BalancingMaster master )
			{
				_Master = master;
			}

			#region IEnumerable<int> Members

			public IEnumerator<int> GetEnumerator()
			{
				return this;
			}

			#endregion

			#region IEnumerable Members

			System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
			{
				return this;
			}

			#endregion

			#region IEnumerator<int> Members

			public int Current
			{
				get { return _Base + _Index; }
			}

			#endregion

			#region IDisposable Members

			public void Dispose()
			{
			}

			#endregion

			#region IEnumerator Members

			object System.Collections.IEnumerator.Current
			{
				get { return _Base + _Index; }
			}

			public bool MoveNext()
			{
				if ( _Finished ) return false;

				_Index += 1;

				if ( _Base + _Index < BalancingMaster.MAX )
				{
					if ( _Index < BalancingChunk ) return true;

					_Index = 0;
					_Base = _Master.Next;
					if ( _Base >= 0 ) return true;
				}

				_Finished = true;
				return false;
			}

			public void Reset()
			{
				_Base = 0;
				_Index = BalancingChunk;
				_Finished = false;
			}

			#endregion
		}

		//-----------------------------------------------------------------------------------------

		static byte Cap( int i ) { return ( byte ) ( i < 0 ? 0 : i > 255 ? 255 : i ); }
		static byte Cap( float i ) { return ( byte ) ( i < 0 ? 0 : i > 255 ? 255 : i ); }
		static byte Cap( double i ) { return ( byte ) ( i < 0 ? 0 : i > 255 ? 255 : i ); }

		//-----------------------------------------------------------------------------------------
	}
}