using Hybridizer.Runtime.CUDAImports; using System; using System.Diagnostics; using System.Threading.Tasks; namespace Mandelbrottyzeug { public class MandelbrotDouble { private readonly int XPixelsDisplayed; // The width of the display in pixels. private readonly int YPixelsDisplayed; // The height of the display in pixels. private readonly int XPixelsSampled; // The number of unique real values sampled at each zoom level. private readonly int YPixelsSampled; // The number of unique imaginary values sampled at each zoom level. private readonly int PixelsSampled; // The total number of pixels sampled. private readonly int XCenterIndex; // The zero-based index for the centermost* column of sampled pixels. private readonly int YCenterIndex; // The zero-based index for the centermost* row of sampled pixels. // * If the number of sampled pixels along a dimension is even, // the central index is rounded up from the true midpoint. public double ZoomLevel { get; private set; } // The ratio of the width/height of the top zoom level to that of the current zoom level. public double PixelWidth { get; private set; } // The amount by which the real or imaginary value changes across a single pixel. public double XMinValue { get; private set; } // The real value of the center of the leftmost column of sampled pixels. public double XMaxValue { get; private set; } // The real value of the center of the rightmost column of sampled pixels. public double YMinValue { get; private set; } // The imaginary value of the center of the topmost row of sampled pixels. public double YMaxValue { get; private set; } // The imaginary value of the center of the bottommost row of sampled pixels. public double XCenterValue { get; private set; } // The real value of the center of the centermost* column of sampled pixels. public double YCenterValue { get; private set; } // The imaginary value of the center of the centermost* row of sampled pixels. // * If the number of sampled pixels along a dimension is even, // the central index is rounded up from the true midpoint. private const int MaxIter = 10000000; // The maximum number of iterations used in generating the mandelbrot set. public bool Iterated { get; private set; } // True if the Mandelbrot is iterated to obtain EscTime. Otherwise EscTime is null. public bool BoundaryIterated { get; private set; } // True if the Mandelbrot boundary is iterated to obtain EscTime. private readonly double BoundaryFactor; // A variable which determines the thickness of the boundary region (larger being thicker). private readonly bool DrawBoundary; // True if the boundary region is colored solid (can save iterations). private const double ColorCycleFactor = 1; // A variable which determines how quickly the exterior gradient changes (smaller being faster). public double[,] EscTime { get; private set; } // Contains the escape time* for the point at the center of the [i,j]th sample. // * Escape time for an interior point can be interpreted as the // time required to converge to an attracting periodic point. public bool?[,] InInterior { get; private set; } // Element [i,j] is true if the [i,j]th sample point is an interior point. public bool[,] InBoundary { get; private set; } // Element [i,j] is true if the [i,j]th sample point is a boundary point. public MandelbrotDouble(int XPixelsDisplayed, int YPixelsDisplayed, int XPixelsSampled, int YPixelsSampled, int XCenterIndex, int YCenterIndex, double ZoomLevel, double PixelWidth, double XCenterValue, double YCenterValue, bool DrawBoundary, double BoundaryFactor, bool Iterated) { this.XPixelsDisplayed = XPixelsDisplayed; this.YPixelsDisplayed = YPixelsDisplayed; this.XPixelsSampled = XPixelsSampled; this.YPixelsSampled = YPixelsSampled; PixelsSampled = XPixelsSampled * YPixelsSampled; this.XCenterIndex = XCenterIndex; this.YCenterIndex = YCenterIndex; this.ZoomLevel = ZoomLevel; this.PixelWidth = PixelWidth; XMinValue = XCenterValue - XCenterIndex * PixelWidth; XMaxValue = XCenterValue + (XPixelsSampled - 1 - XCenterIndex) * PixelWidth; YMinValue = YCenterValue - YCenterIndex * PixelWidth; YMaxValue = YCenterValue + (YPixelsSampled - 1 - YCenterIndex) * PixelWidth; this.XCenterValue = XCenterValue; this.YCenterValue = YCenterValue; this.BoundaryFactor = BoundaryFactor; this.DrawBoundary = DrawBoundary; this.Iterated = Iterated; BoundaryIterated = DrawBoundary && Iterated; if (Iterated) { EscTime = new double[XPixelsSampled, YPixelsSampled]; InInterior = new bool?[XPixelsSampled, YPixelsSampled]; InBoundary = new bool[XPixelsSampled, YPixelsSampled]; CallGenerate(); } } public void CallGenerate() { if (!Iterated) { Iterated = true; EscTime = new double[XPixelsSampled, YPixelsSampled]; InInterior = new bool?[XPixelsSampled, YPixelsSampled]; InBoundary = new bool[XPixelsSampled, YPixelsSampled]; } BoundaryIterated = DrawBoundary || BoundaryIterated; // cuda.GetDeviceProperties(out cudaDeviceProp prop, 0); // int N = prop.multiProcessorCount * 16; // HybRunner runner = HybRunner.Cuda("Mandelbrottyzeug_CUDA.dll").SetDistrib(N, 128); // dynamic wrapper = runner.Wrap(new Program()); MandelbrotDouble NewMandelbrotDouble = GenerateMandelbrot(this); EscTime = NewMandelbrotDouble.EscTime; InInterior = NewMandelbrotDouble.InInterior; InBoundary = NewMandelbrotDouble.InBoundary; } [EntryPoint("DoubleDirectEntry")] private static MandelbrotDouble GenerateMandelbrot(MandelbrotDouble MyMandelbrotDouble) { int XPixelsSampled = MyMandelbrotDouble.XPixelsSampled; int PixelsSampled = MyMandelbrotDouble.PixelsSampled; int XCenterIndex = MyMandelbrotDouble.XCenterIndex; int YCenterIndex = MyMandelbrotDouble.YCenterIndex; double XCenterValue = MyMandelbrotDouble.XCenterValue; double YCenterValue = MyMandelbrotDouble.YCenterValue; double PixelWidth = MyMandelbrotDouble.PixelWidth; double[,] EscTime = MyMandelbrotDouble.EscTime; bool?[,] InInterior = MyMandelbrotDouble.InInterior; bool[,] InBoundary = MyMandelbrotDouble.InBoundary; Parallel.For(0, PixelsSampled, pixelIndex => { int xIndex = pixelIndex % XPixelsSampled; int yIndex = pixelIndex / XPixelsSampled; double zx = XCenterValue + (xIndex - XCenterIndex) * PixelWidth; double zy = YCenterValue + (yIndex - YCenterIndex) * PixelWidth; (double escTime, bool? inInterior, bool inBoundary) = IterateMandelbrot(zx, zy, MyMandelbrotDouble); EscTime[xIndex, yIndex] = escTime; InInterior[xIndex, yIndex] = inInterior; InBoundary[xIndex, yIndex] = inBoundary; }); MyMandelbrotDouble.EscTime = EscTime; MyMandelbrotDouble.InInterior = InInterior; MyMandelbrotDouble.InBoundary = InBoundary; return MyMandelbrotDouble; } [Kernel] private static (double, bool?, bool) IterateMandelbrot(double z0x, double z0y, MandelbrotDouble MyMandelbrotDouble) { double escRad = 10; double escRadSqrd = escRad * escRad; double znx = z0x; double zny = z0y; // znDeriv1 == The derivative of zn with respect to z0 double znDeriv1x = 1; double znDeriv1y = 0; // innerMax == maximum threshold of znDeriv1 for zn to be considered hyperbolic double innerMax = 1e-6; double innerMaxSqrd = innerMax * innerMax; // znDeriv2 == The derivative of zn with respect to c double BoundaryFactor = MyMandelbrotDouble.BoundaryFactor; double PixelWidth = MyMandelbrotDouble.PixelWidth; double cDeriv2x = BoundaryFactor * PixelWidth; double znDeriv2x = cDeriv2x; double znDeriv2y = 0; bool boundaryPt = false; for (int i = 0; i < MaxIter; i++) { double znxSqrd = znx * znx; double znySqrd = zny * zny; double znModulusSqrd = znxSqrd + znySqrd; double znSqrdRealPt; double znSqrdImagPt; if (!boundaryPt) { double znDeriv1xSqrd = znDeriv1x * znDeriv1x; double znDeriv1ySqrd = znDeriv1y * znDeriv1y; double znDeriv1ModulusSqrd = znDeriv1xSqrd + znDeriv1ySqrd; double znDeriv2xSqrd = znDeriv2x * znDeriv2x; double znDeriv2ySqrd = znDeriv2y * znDeriv2y; double znDeriv2ModulusSqrd = znDeriv2xSqrd + znDeriv2ySqrd; // calculate the modulus of zn + c and znDeriv2 + cDeriv2 // to get rid of the boundaries detected at the bulb centers double znxModSqrd = znx + z0x; znxModSqrd *= znxModSqrd; double znyModSqrd = zny + z0y; znyModSqrd *= znyModSqrd; double znModModulusSqrd = znxModSqrd + znyModSqrd; znModModulusSqrd = znModModulusSqrd * (1 + 0.25 * znModModulusSqrd); double znDeriv2xModSqrd = znDeriv2x + cDeriv2x; znDeriv2xModSqrd *= znDeriv2xModSqrd; double znDeriv2ModModulusSqrd = znDeriv2xModSqrd + znDeriv2ySqrd; if (znDeriv1ModulusSqrd < innerMaxSqrd) { // interior point return (1 - Math.Log(i + 1) / Math.Log(MaxIter + 1), true, false); } else if (znModModulusSqrd < znDeriv2ModModulusSqrd) { // boundary point if (MyMandelbrotDouble.DrawBoundary) { return (0, null, true); } else { boundaryPt = true; } } else if (znModulusSqrd > escRadSqrd) { // exterior point double iMod = i + 1 - Math.Log(0.5 * Math.Log(znModulusSqrd), 2); return (128 * Math.Log(iMod) / ColorCycleFactor, false, false); } znSqrdRealPt = znxSqrd - znySqrd; znSqrdImagPt = znx + zny; znSqrdImagPt *= znSqrdImagPt; znSqrdImagPt -= znModulusSqrd; double znDeriv1SqrdRealPt = znDeriv1xSqrd - znDeriv1ySqrd; double znxPlusDeriv1x = znx + znDeriv1x; double znyPlusDeriv1y = zny + znDeriv1y; double znxPlusDeriv1y = znx + znDeriv1y; double znyPlusDeriv1x = zny + znDeriv1x; znDeriv1x = znxPlusDeriv1x * znxPlusDeriv1x - znDeriv1SqrdRealPt - znyPlusDeriv1y * znyPlusDeriv1y - znSqrdRealPt; znDeriv1y = znxPlusDeriv1y * znxPlusDeriv1y - znDeriv1ModulusSqrd + znyPlusDeriv1x * znyPlusDeriv1x - znModulusSqrd; double znDeriv2SqrdRealPt = znDeriv2xSqrd - znDeriv2ySqrd; double znxPlusDeriv2x = znx + znDeriv2x; double znyPlusDeriv2y = zny + znDeriv2y; double znxPlusDeriv2y = znx + znDeriv2y; double znyPlusDeriv2x = zny + znDeriv2x; znDeriv2x = znxPlusDeriv2x * znxPlusDeriv2x - znDeriv2SqrdRealPt - znyPlusDeriv2y * znyPlusDeriv2y - znSqrdRealPt + cDeriv2x; znDeriv2y = znxPlusDeriv2y * znxPlusDeriv2y - znDeriv2ModModulusSqrd + znyPlusDeriv2x * znyPlusDeriv2x - znModulusSqrd; } else { // point will not satisfy interior point test if (znModulusSqrd > escRadSqrd) { // exterior point double iMod = i + 1 - Math.Log(0.5 * Math.Log(znModulusSqrd), 2); return (128 * Math.Log(iMod) / ColorCycleFactor, false, true); } znSqrdRealPt = znxSqrd - znySqrd; znSqrdImagPt = (znx + zny) * (znx + zny) - znModulusSqrd; } znx = znSqrdRealPt + z0x; zny = znSqrdImagPt + z0y; } // not enough iterates (assume interior) return (0, true, boundaryPt); } public void GenerateMandelbrotBoundary() { if (!Iterated) throw new ArgumentNullException(); BoundaryIterated = true; Parallel.For(0, XPixelsSampled * YPixelsSampled, pixelIndex => { int xIndex = pixelIndex % XPixelsSampled, yIndex = pixelIndex / XPixelsSampled; if (InBoundary[xIndex, yIndex]) { double zx = XCenterValue + (xIndex - XCenterIndex) * PixelWidth; double zy = YCenterValue + (yIndex - YCenterIndex) * PixelWidth; (double escTime, bool? inInterior, bool inBoundary) = IterateMandelbrotBoundary(zx, zy); EscTime[xIndex, yIndex] = escTime; InInterior[xIndex, yIndex] = inInterior; InBoundary[xIndex, yIndex] = inBoundary; } }); } private (double, bool?, bool) IterateMandelbrotBoundary(double z0x, double z0y) { double escRad = 10; double escRadSqrd = escRad * escRad; double znx = z0x; double zny = z0y; for (int i = 0; i < MaxIter; i++) { double znxSqrd = znx * znx; double znySqrd = zny * zny; double znModulusSqrd = znxSqrd + znySqrd; // point will not satisfy interior point test if (znModulusSqrd > escRadSqrd) { // exterior point double iMod = i + 1 - Math.Log(0.5 * Math.Log(znModulusSqrd), 2); return (128 * Math.Log(iMod) / ColorCycleFactor, false, true); } double znSqrdRealPt = znxSqrd - znySqrd; double znSqrdImagPt = znx + zny; znSqrdImagPt *= znSqrdImagPt; znSqrdImagPt -= znModulusSqrd; znx = znSqrdRealPt + z0x; zny = znSqrdImagPt + z0y; } // not enough iterates (assume interior) return (0, true, true); } } }