Mandelbrot and Julia sets C program

I find that Mandelbrot and Julia sets are one of the more artistic results of mathematics. I made a fast C program with an assembler function which computes the points of a set. It runs on ix86 or compatible processor, using the stack of the floating point unit. As the polynomial variables are all stored in the FPU stack, without I/O to RAM, the computation is very fast.

Download

Download the code of mj, the program to calculate the sets, and mj2ppm, the program to colorize the raw data output of mj, and make a ppm image file: mj.tar.gz

How to use it

For Mandelbrot set:

$ ./mj N_of_horizontal_points(unsigned int) N_of_vertical_points(unsigned int)
N_of_maximum_iterations(unsigned int) Squared_radius_of_divergence(double)
X(double) Y(double) Side_length(double) > raw.file

For Julia set:

$ ./mj N_of_horizontal_points(unsigned int) N_of_vertical_points(unsigned int)
N_of_maximum_iterations(unsigned int) Squared_radius_of_divergence(double)
X(double) Y(double) Side_length(double) X_origin(double) Y_origin(double)
> raw.file

To convert raw data into ppm file:

$ ./mj2ppm N_of_min_treshold(int) < raw.file > output.ppm

Example:

$ ./mj  200 100 500 500.0 -1.5 1.0 2.0 | ./mj2ppm 5 > out.ppm

Will give this result:

Here comes the main function code of mj

  mj 3.5.2  Copyright (C) 2007 Marco Bisetto < marco [at] folgorante.net>
  This program comes with ABSOLUTELY NO WARRANTY;
  This is free software, and you are welcome to redistribute it
  under certain conditions;

  Released under the terms of GNU GPLv3.

static inline unsigned int punto(struct mjpt ss)
{
  unsigned int n;
  
  asm ("\tfninit\n");

  asm (
       "fld %%st(1)\n\t"
       "fmul %%st(0)\n\t"
       "fld %%st(1)\n\t"
       "fmul %%st(0)\n\t"
       "movl $0,%%ecx\n"
     "ritor:\t"
       "fsubrp %%st,%%st(1)\n\t"
       "fadd %%st(3)\n\t"
       "fxch %%st(2)\n\t"
       "fmulp %%st,%%st(1)\n\t"
       "fadd %%st(0)\n\t"
       "fadd %%st(3)\n\t"
       "fld %%st(1)\n\t"
       "fmul %%st(0)\n\t"
       "fld %%st(1)\n\t"
       "fmul %%st(0)\n\t"
       "fld %%st(0)\n\t"
       "fadd %%st(2)\n\t"
       "fcomp %%st(7)\n\t"
       "fnstsw %%ax\n\t"
       "and $0x41, %%ah\n\t"
       "jz esci\n\t"
       "incl %%ecx\n\t"
       "cmpl %%ebx,%%ecx\n\t"
       "jb ritor\n"
     "esci:\n"
       : "=c" (n)
       : "b" (ss.ma), "f" (ss.mam2), "f" (ss.yc), "f" (ss.xc),
         "f" (ss.xi), "f" (ss.yi)
       : "%eax"
       );

  asm ("\tfinit\n");
  
  return n;
}