Imaginary numbers?

[Nexion]

I wanted to test to see how well LOVE could handle a fractal, or if Lua could even support it; but I can't find any way to do imaginary numbers with Lua.
There is no math.i or anything to represent math.sqrt(-1), and that would error the program anyway.
Can anyone think of a way to get imaginary numbers?

[muku]

Well, implement them yourself. You could simply represent them as a table like

Code: Select all

x = { r=2, i=-1}

and the standard operations on complex numbers are easy enough to implement (consult Wikipedia if you are unsure). You could also give them a metatable with the appropriate functions so that you can simply write stuff like x+y or x*y for complex numbers x and y.

Of course, the simplest way is to use a finished module

[Nexion]

Awesome! Thank you!
Does anyone know how the deeper parts of the Mandelbrot Set work? From what I've read it iterates through all the pixels and colors them accordingly, but the equation is Z = Z^2 * C, which escapes me simply due to I have no idea how to actually implement that into the code itself.

[rude]

I just [1] used the Mandelbrot set on a Beowulf cluster ... and I've already forgotten everything, appearently.

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <complex.h>
#include <mpi.h>
#include "plot.h"


#define THRESHOLD 1000.0
#define MAXITER 1500
#define SCALE 30000.0


double complex
    global_lower = -0.512  - 1.024*I,
    global_upper =  1.024  + 1.024*I;

/*
 * Estimate the divergence rate of our function for a given complex point
 */
double
test_point ( double complex c )
{
    double complex z_new, z, z_prev;
    double difference;

    z = z_prev = c;

    // Iterate until divergence or max. number of iterations
    for ( int iter=0; iter<MAXITER; iter++ )
    {
        // Update the value by the formula
        z_new = z*z + 0.56667 - (z_prev * 0.5);

        // Buffer the changes
        z_prev = z;
        z = z_new;

        // Test for "fast enough" divergence - if this is found, we can end
        // the iteration already, assuming that the sequence shortly will
        // disappear into infinity
        difference = cabs(z) - cabs(z_prev);
        if ( difference > THRESHOLD )
            break;
    }

    // Scale the difference value to make a nicer plot
    difference /= SCALE;
    return (difference > 1.0) ? 1.0 : difference;
}


/*
 * This function should evaluate the points in the rectangle bounded by
 * the lower/upper variables, discretized into a grid of w by h points.
 * When this is called, the 'data' pointer should already be appropriately
 * allocated to hold these w by h points.
 *
 * Note that while this array is one-dimensional, there is a convenience
 * macro AT(x,y) in 'plot.h' to address it as a 2D array (writing
 * data[AT(x,y)]).
 */
void
compute_region (
    double complex lower, double complex upper, int w, int h, double *data
)
{
   double complex delta = upper - lower;
   double real = creal(delta)/(double)w;
   double imag = cimag(delta)/(double)h;
   
   for(int i=0;i<w;i++)
   {
   	for(int j=0;j<h;j++)
   	{
   		double complex c = real * i + imag * j * I;
   		data[i*ROWS+j] = test_point( lower + c);
   	}
   }
}


int
main ( int argc, char **argv )
{
    double
        *img, *local_img,         // Global and local image parts
        t_start, t_finish, t_elapsed, max_elapsed;  // Timing variables

    complex double
        local_lower, local_upper, delta;    // Bounds and step size of domain

    int
        rank, size, blocksize;  // Used to compute problem partitioning

    // Initialize MPI, get world rank and size
    MPI_Init ( &argc, &argv );
    MPI_Comm_rank ( MPI_COMM_WORLD, &rank );
    MPI_Comm_size ( MPI_COMM_WORLD, &size );

    // Rank 0 hosts the complete image
    if ( rank == 0 )
        img = calloc ( ROWS * COLS, sizeof(double) );

	// YOUR CODE HERE	
	
	// Assign local variables.
	delta = (creal(global_upper) - creal(global_lower))/(double)size + 0 * I;
	local_lower = global_lower + delta * (double)rank;
	local_upper = creal(local_lower+delta) + cimag(global_upper) * I;

	// Allocate memory.
	int memsize = (ROWS * COLS)/size;
	local_img = calloc(memsize, sizeof(double));
	
	// Compute region
	t_start = MPI_Wtime();
	compute_region(local_lower, local_upper, COLS/size, COLS, local_img);
	t_finish = MPI_Wtime();
	
	t_elapsed = t_finish - t_start;	

	// Gather at 0
	MPI_Gather(local_img, memsize, MPI_DOUBLE, img + memsize * rank, memsize, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	// Reduce
	MPI_Reduce(&t_elapsed, &max_elapsed, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

    // Rank 0 reports greatest time spent in computation, writes final image
    if ( rank == 0 )
    {
        printf ( "%lf [s] spent on computation\n", max_elapsed );
        write_image ( img );
        free ( img );
    }
    
    // Free the local image.
    free(local_img);
    
    MPI_Finalize ();
    exit ( EXIT_SUCCESS );
}

test_point() is what you want to look at.

[1]: Just = ~6 months ago

EDIT: Might be Julia set, but who cares.