Commits

Ben Bass committed dfc0337

make filenames consistent with project name

Comments (0)

Files changed (9)

 out.ppm: make_ppm mandel.dat
 	./make_ppm
   
-mandel.dat: mandelmap
-	time ./mandelmap
+mandel.dat: brotmap
+	time ./brotmap
 
-make_ppm: make_ppm.cc mandelmap.h
+make_ppm: make_ppm.cc brotmap.h
 	g++ make_ppm.cc -o make_ppm
   
-mandelmap: mandelmap.cc mandelmap.cc mandelmap.h
-	g++ mandelmap.cc -Werror -O3 -o mandelmap -lc -lpthread
+brotmap: brotmap.cc brotmap.h
+	g++ brotmap.cc -Werror -O3 -o brotmap -lc -lpthread
 	
 clean:
-	rm -f mandelmap make_ppm out.ppm out_image
+	rm -f brotmap make_ppm out.ppm out_image
 
 superclean: clean
 	rm -f mandel_*.dat
+// brotmap - Ben Bass 2010-2012
+
+#include <cstdio>
+#include <math.h>  // for NAN. Not in the C++ library...
+#include <cstddef>
+#include "brotmap.h"
+
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <sys/mman.h>
+
+#include <pthread.h>
+
+#ifndef _POSIX_MAPPED_FILES
+#error "no mapped file support"
+#endif
+
+#ifndef MAP_NOCACHE
+#define MAP_NOCACHE 0
+#endif
+
+long mpoint(pinfo* p);
+
+long inside_points = 0;
+
+// these are the dimensions of the fractal we are plotting.
+const FLOAT MIN_X = -2.0;
+const FLOAT MAX_X = 1.0;
+const FLOAT MIN_Y = -1.5;
+const FLOAT MAX_Y = 1.5;
+const FLOAT BAILOUT = 4.0;
+
+// mpoint(real, image, ...) -> 1 if mset else 0
+inline long mpoint(FLOAT r,
+                  FLOAT i,
+                  pinfo* p,
+                  int old_max_iter=0)
+{
+    if (isnan(p->x)) {
+        // we've already processed this entry -
+        // and it's not part of the Mset
+        return 0;
+    }
+    FLOAT x2, y2;
+    FLOAT x = p->x;
+    FLOAT y = p->y;
+    int k = MAX_ITER;
+
+#ifdef CARDIOID_OPT
+    // based on 'Optimizations' section in wikipedia
+    // http://en.wikipedia.org/wiki/Mandelbrot_set#Optimizations
+    double roff = (r-0.25);
+    double isq = i*i;
+    double dist = roff*roff + isq;
+    if (dist * (dist + roff) < 0.25*isq)
+    {
+        return 1;
+    }
+
+    roff = r+1.0;
+    if (roff*roff + isq < 1.0/16.0)
+    {
+        return 1;
+    }
+#endif
+    // this is the core 'inner loop'.
+    // It should be fast, which is why it doesn't
+    // have lots of whitespace in it ;)
+    //
+    // of course this is massively parallelizable
+    do {
+        x2 = x*x;
+        y2 = y*y;
+        if ((x2+y2) >= BAILOUT) break;
+        y = 2.0*x*y+i;
+        x = x2-y2+r;
+    } while (--k);
+
+    if (!k)
+    {
+        // we ran out of iterations. Up to MAX_ITER,
+        // this point seems to be part of the mset.
+        // Store current values of x and y so we can
+        // carry on later.
+        p->x = x;
+        p->y = y;
+        return 1;
+    }
+    else
+    {
+        // this is definitely not part of the mset
+        // (rounding errors excepted). We want to store
+        // that fact, together with an iteration count
+        // and the magnitude of the calculation at that
+        // point (so we can do smooth renderings)
+        p->x = NAN;
+        p->itercount = (old_max_iter + (MAX_ITER - k));
+        return 0;
+    }
+}
+
+pinfo* FPTR_START;
+FLOAT STEP_SIZE;
+pthread_mutex_t acc_lock;
+
+void* worker_start(void* arg)
+{
+    int t_num = *reinterpret_cast<int*>(arg);
+    FLOAT t_start = t_num*STEP_SIZE;
+
+    long local_inside = 0;
+
+    // get the number of points to skip per row
+    ptrdiff_t row_skip = (MAX_X-MIN_X)/STEP_SIZE;
+    // get starting point in the mmapped file - t_num rows in
+    pinfo* local_fptr = FPTR_START+(row_skip*t_num);
+
+    for (FLOAT y=MIN_Y + t_start; y<MAX_Y; y+=(NUM_THREADS*STEP_SIZE))
+    {
+        for (FLOAT x=MIN_X; x<MAX_X; x+=STEP_SIZE)
+        {
+            pinfo* p=local_fptr++;
+            local_inside += mpoint(x, y, p, MAX_ITER);
+        }
+        // skip ahead the relevant number of lines
+        local_fptr += (row_skip * (NUM_THREADS-1));
+    }
+    pthread_mutex_lock(&acc_lock);
+    inside_points += local_inside;
+    pthread_mutex_unlock(&acc_lock);
+}
+
+int main()
+{
+    // for multiples of 2^-n (where n is < the number of mantissa digits)
+    // floating point can exactly represent the values. Work out the
+    // such a value by repeated dividing by two.
+    //
+    // TODO: this is all constant, so use a bit of TMP on it. Just for fun
+    FLOAT step = 1.0;
+    unsigned long long invstep = 1;
+    for (int i=0; i<BINARY_DIGITS; i++)
+    {
+        step /= 2;
+        invstep *= 2;
+    }
+    const FLOAT plot_area = (MAX_Y-MIN_Y)*(MAX_X-MIN_X);  // hopefully an integer.
+    // following assumes we're drawing a square.
+    const unsigned long long total_points = plot_area * invstep * invstep;
+
+    // we will write the data to a file mandel_{BD}.dat where BD
+    // is the 'n' in the step size 2^-n
+    // a 'single precision' full map (23 bits mantissa) would take
+    // up nearly 10 petabytes (2^50) of storage. Much better use
+    // than most of the data out there, me thinks.
+    const unsigned long long mapsize = total_points*sizeof(pinfo) + HEADER_LEN;
+    char fname[21];
+    snprintf(fname, 20, "mandel_%d.dat", BINARY_DIGITS);
+    int fd = open(fname, O_RDWR|O_CREAT, 0644);
+    if (ftruncate(fd, mapsize) != 0) {
+        printf("ftruncate failed\n");
+        return 1;
+    }
+
+    // we will use mmap to write to the data file. It makes things wonderful.
+    void* mapping;
+    if (NULL == (mapping = mmap(0, mapsize,
+            PROT_READ|PROT_WRITE,
+            MAP_FILE|MAP_SHARED|MAP_NOCACHE,
+            fd, 0)))
+    {
+        printf("mmap failed\n");
+        return 1;
+    }
+    close(fd);  // we no longer need the FD when we have the mapping
+
+    void* sourceBuffer = mapping;  // for unmmaping later
+    brotfile_header* fheader = (brotfile_header*)(mapping);
+    pinfo* fptr = (pinfo*)((char*)mapping+HEADER_LEN);
+
+    printf("mapped file %d at %p %llu bytes\n", fd, fptr, mapsize);
+    printf("Current MAX_ITER: %d\n", fheader->max_iter);
+
+    pthread_t t_handle[NUM_THREADS];
+    int t_num[NUM_THREADS];
+
+    // setup global read-only vars //BADBADBAD (but could be worse)
+    FPTR_START = fptr;
+    STEP_SIZE = step;
+
+    printf("Using %d threads\n", NUM_THREADS);
+
+    pthread_mutex_init(&acc_lock, NULL);
+    for (int nt=0; nt < NUM_THREADS; nt++)
+    {
+        t_num[nt] = nt;
+        pthread_create(&t_handle[nt],
+                NULL, worker_start, reinterpret_cast<void*>(&t_num[nt]));
+    }
+    for (int nt=0; nt < NUM_THREADS; nt++)
+    {
+        pthread_join(t_handle[nt], NULL);
+    }
+
+    fheader->max_iter += MAX_ITER;
+    fheader->binary_size = BINARY_DIGITS;
+
+    if (sourceBuffer) {
+        munmap(sourceBuffer, mapsize);
+    }
+    // print total pixels, 'inside' pixels, and area of mset based on
+    // these figures.
+    // Should be approx 1.50648
+    // https://www.fractalus.com/kerry/articles/area/mandelbrot-area.html
+    printf("%llu %lu (%.6f)\n", total_points, inside_points, ((plot_area*inside_points)/float(total_points)));
+}
+// MandelMap - Ben Bass 2010-2012
+
+#ifndef MANDEL_MAP_H
+#define MANDEL_MAP_H
+
+#define CARDIOID_OPT
+
+typedef double FLOAT;
+
+const int BINARY_DIGITS = 11;
+const int NUM_THREADS = 2;
+const int MAX_ITER=1024;
+
+struct pinfo
+{
+    FLOAT x;
+    union {
+      FLOAT y;
+      long itercount;
+    };
+};
+
+
+const int HEADER_LEN=4096;  // probable page size for happiness
+struct brotfile_header
+{
+    FLOAT startx;
+    FLOAT starty;
+    int binary_size;
+    int max_iter;
+};
+
+
+#endif
 #include <stdio.h>
 #include <math.h>
-#include "mandelmap.h"
+#include "brotmap.h"
 
 const int SIZE= (sizeof(pinfo)*3*(1<<BINARY_DIGITS-4));
 

mandelmap.cc

-// MandelMap - Ben Bass 2010-2012
-
-#include <cstdio>
-#include <math.h>  // for NAN. Not in the C++ library...
-#include <cstddef>
-#include "mandelmap.h"
-
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <unistd.h>
-#include <sys/mman.h>
-
-#include <pthread.h>
-
-#ifndef _POSIX_MAPPED_FILES
-#error "no mapped file support"
-#endif
-
-#ifndef MAP_NOCACHE
-#define MAP_NOCACHE 0
-#endif
-
-long mpoint(pinfo* p);
-
-long inside_points = 0;
-
-// these are the dimensions of the fractal we are plotting.
-const FLOAT MIN_X = -2.0;
-const FLOAT MAX_X = 1.0;
-const FLOAT MIN_Y = -1.5;
-const FLOAT MAX_Y = 1.5;
-const FLOAT BAILOUT = 4.0;
-
-// mpoint(real, image, ...) -> 1 if mset else 0
-inline long mpoint(FLOAT r,
-                  FLOAT i,
-                  pinfo* p,
-                  int old_max_iter=0)
-{
-    if (isnan(p->x)) {
-        // we've already processed this entry -
-        // and it's not part of the Mset
-        return 0;
-    }
-    FLOAT x2, y2;
-    FLOAT x = p->x;
-    FLOAT y = p->y;
-    int k = MAX_ITER;
-
-#ifdef CARDIOID_OPT
-    // based on 'Optimizations' section in wikipedia
-    // http://en.wikipedia.org/wiki/Mandelbrot_set#Optimizations
-    double roff = (r-0.25);
-    double isq = i*i;
-    double dist = roff*roff + isq;
-    if (dist * (dist + roff) < 0.25*isq)
-    {
-        return 1;
-    }
-
-    roff = r+1.0;
-    if (roff*roff + isq < 1.0/16.0)
-    {
-        return 1;
-    }
-#endif
-    // this is the core 'inner loop'.
-    // It should be fast, which is why it doesn't
-    // have lots of whitespace in it ;)
-    //
-    // of course this is massively parallelizable
-    do {
-        x2 = x*x;
-        y2 = y*y;
-        if ((x2+y2) >= BAILOUT) break;
-        y = 2.0*x*y+i;
-        x = x2-y2+r;
-    } while (--k);
-
-    if (!k)
-    {
-        // we ran out of iterations. Up to MAX_ITER,
-        // this point seems to be part of the mset.
-        // Store current values of x and y so we can
-        // carry on later.
-        p->x = x;
-        p->y = y;
-        return 1;
-    }
-    else
-    {
-        // this is definitely not part of the mset
-        // (rounding errors excepted). We want to store
-        // that fact, together with an iteration count
-        // and the magnitude of the calculation at that
-        // point (so we can do smooth renderings)
-        p->x = NAN;
-        p->itercount = (old_max_iter + (MAX_ITER - k));
-        return 0;
-    }
-}
-
-pinfo* FPTR_START;
-FLOAT STEP_SIZE;
-pthread_mutex_t acc_lock;
-
-void* worker_start(void* arg)
-{
-    int t_num = *reinterpret_cast<int*>(arg);
-    FLOAT t_start = t_num*STEP_SIZE;
-
-    long local_inside = 0;
-
-    // get the number of points to skip per row
-    ptrdiff_t row_skip = (MAX_X-MIN_X)/STEP_SIZE;
-    // get starting point in the mmapped file - t_num rows in
-    pinfo* local_fptr = FPTR_START+(row_skip*t_num);
-
-    for (FLOAT y=MIN_Y + t_start; y<MAX_Y; y+=(NUM_THREADS*STEP_SIZE))
-    {
-        for (FLOAT x=MIN_X; x<MAX_X; x+=STEP_SIZE)
-        {
-            pinfo* p=local_fptr++;
-            local_inside += mpoint(x, y, p, MAX_ITER);
-        }
-        // skip ahead the relevant number of lines
-        local_fptr += (row_skip * (NUM_THREADS-1));
-    }
-    pthread_mutex_lock(&acc_lock);
-    inside_points += local_inside;
-    pthread_mutex_unlock(&acc_lock);
-}
-
-int main()
-{
-    // for multiples of 2^-n (where n is < the number of mantissa digits)
-    // floating point can exactly represent the values. Work out the
-    // such a value by repeated dividing by two.
-    //
-    // TODO: this is all constant, so use a bit of TMP on it. Just for fun
-    FLOAT step = 1.0;
-    unsigned long long invstep = 1;
-    for (int i=0; i<BINARY_DIGITS; i++)
-    {
-        step /= 2;
-        invstep *= 2;
-    }
-    const FLOAT plot_area = (MAX_Y-MIN_Y)*(MAX_X-MIN_X);  // hopefully an integer.
-    // following assumes we're drawing a square.
-    const unsigned long long total_points = plot_area * invstep * invstep;
-
-    // we will write the data to a file mandel_{BD}.dat where BD
-    // is the 'n' in the step size 2^-n
-    // a 'single precision' full map (23 bits mantissa) would take
-    // up nearly 10 petabytes (2^50) of storage. Much better use
-    // than most of the data out there, me thinks.
-    const unsigned long long mapsize = total_points*sizeof(pinfo) + HEADER_LEN;
-    char fname[21];
-    snprintf(fname, 20, "mandel_%d.dat", BINARY_DIGITS);
-    int fd = open(fname, O_RDWR|O_CREAT, 0644);
-    if (ftruncate(fd, mapsize) != 0) {
-        printf("ftruncate failed\n");
-        return 1;
-    }
-
-    // we will use mmap to write to the data file. It makes things wonderful.
-    void* mapping;
-    if (NULL == (mapping = mmap(0, mapsize,
-            PROT_READ|PROT_WRITE,
-            MAP_FILE|MAP_SHARED|MAP_NOCACHE,
-            fd, 0)))
-    {
-        printf("mmap failed\n");
-        return 1;
-    }
-    close(fd);  // we no longer need the FD when we have the mapping
-
-    void* sourceBuffer = mapping;  // for unmmaping later
-    brotfile_header* fheader = (brotfile_header*)(mapping);
-    pinfo* fptr = (pinfo*)((char*)mapping+HEADER_LEN);
-
-    printf("mapped file %d at %p %llu bytes\n", fd, fptr, mapsize);
-    printf("Current MAX_ITER: %d\n", fheader->max_iter);
-
-    pthread_t t_handle[NUM_THREADS];
-    int t_num[NUM_THREADS];
-
-    // setup global read-only vars //BADBADBAD (but could be worse)
-    FPTR_START = fptr;
-    STEP_SIZE = step;
-
-    printf("Using %d threads\n", NUM_THREADS);
-
-    pthread_mutex_init(&acc_lock, NULL);
-    for (int nt=0; nt < NUM_THREADS; nt++)
-    {
-        t_num[nt] = nt;
-        pthread_create(&t_handle[nt],
-                NULL, worker_start, reinterpret_cast<void*>(&t_num[nt]));
-    }
-    for (int nt=0; nt < NUM_THREADS; nt++)
-    {
-        pthread_join(t_handle[nt], NULL);
-    }
-
-    fheader->max_iter += MAX_ITER;
-    fheader->binary_size = BINARY_DIGITS;
-
-    if (sourceBuffer) {
-        munmap(sourceBuffer, mapsize);
-    }
-    // print total pixels, 'inside' pixels, and area of mset based on
-    // these figures.
-    // Should be approx 1.50648
-    // https://www.fractalus.com/kerry/articles/area/mandelbrot-area.html
-    printf("%llu %lu (%.6f)\n", total_points, inside_points, ((plot_area*inside_points)/float(total_points)));
-}

mandelmap.h

-// MandelMap - Ben Bass 2010-2012
-
-#ifndef MANDEL_MAP_H
-#define MANDEL_MAP_H
-
-#define CARDIOID_OPT
-
-typedef double FLOAT;
-
-const int BINARY_DIGITS = 11;
-const int NUM_THREADS = 2;
-const int MAX_ITER=1024;
-
-struct pinfo
-{
-    FLOAT x;
-    union {
-      FLOAT y;
-      long itercount;
-    };
-};
-
-
-const int HEADER_LEN=4096;  // probable page size for happiness
-struct brotfile_header
-{
-    FLOAT startx;
-    FLOAT starty;
-    int binary_size;
-    int max_iter;
-};
-
-
-#endif
+
+brotmap
+-------
+
+*Precalculating the Mandelbrot set - because it's there*
+
+This is mostly a learning exercise for me.
+
+Key things:
+ - use of anonymous unions and NAN for efficient storage of the
+   mandel_X.dat files (depite the fact that they are rather large still...)
+   This was my original 'cool idea' / excuse for starting this.
+ - use of mmap in a 64 bit address space for basically *doing anything*.
+   Nice.
+ - use of pthreads to utilise multi-core processors
+
+And of course, the key thing:
+ - using a computer for it's intended purpose: computing, taking all the
+   resources it possibly can in the process :-)
+
+For the future:
+ - distributed parallelism
+ - not only incremental in terms of iterations, but also incremental in
+   spatial resolution (which will probably require a change in the file
+   format, or something like a group of 4 mandel_X.dats  which are
+   postprocessed to create a mandel_(X+1).dat, by storing an offset
+   of zero or 2^-(X+1) in each of the real and imaginary axes...
+ - an actual command line interface, with options and everything

readme.txt

-
-brotmap / mandelmap
--------------------
-
-This has two names, because ... well, I can't quite remember :-)
-
-It's rather on the messy side, I'll tidy it up at some point.
-The existence of this readme is a step in the right direction :-)
-
-This is mostly a learning exercise for me.
-
-Key things:
- - use of anonymous unions and NAN for efficient storage of the
-   mandel_X.dat files (depite the fact that they are rather large still...)
-   This was my original 'cool idea' / excuse for starting this.
- - use of mmap in a 64 bit address space for basically *doing anything*.
-   Nice.
- - use of pthreads to utilise multi-core processors
-
-And of course, the key thing:
- - using a computer for it's intended purpose: computing, taking all the
-   resources it possibly can in the process :-)
-
-For the future:
- - distributed parallelism
- - not only incremental in terms of iterations, but also incremental in
-   spatial resolution (which will probably require a change in the file
-   format, or something like a group of 4 mandel_X.dats  which are
-   postprocessed to create a mandel_(X+1).dat, by storing an offset
-   of zero or 2^-(X+1) in each of the real and imaginary axes...
- - an actual command line interface, with options and everything
 #include <stdio.h>
-#include "mandelmap.h"
+#include "brotmap.h"
 
 int mpoint(FLOAT r, FLOAT i, int maxit)
 {