from numpy import * # because arrays are defined in numpy

def Mand(z0, max_steps):
    z = 0j  # no need to specify type.  
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

def Mandelbrot(ext, Nxy, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    data = zeros((Nxy,Nxy)) # initialize a 2D dynamic array
    for i in range(Nxy):
        for j in range(Nxy):
            x = ext[0] + (ext[1]-ext[0])*i/(Nxy-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Nxy-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
    return data
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

data = Mandelbrot([-2,1,-1,1], 500, 1000)

from pylab import *    # plotting library
%matplotlib inline

ext=[-2,1,-1,1]
# pylab's function for displaying 2D image

imshow(data.T,origin='lower',extent=ext);

help(imshow)

Help on function imshow in module matplotlib.pyplot:

imshow(X, cmap=None, norm=None, *, aspect=None, interpolation=None, alpha=None, vmin=None, vmax=None, origin=None, extent=None, interpolation_stage=None, filternorm=True, filterrad=4.0, resample=None, url=None, data=None, **kwargs)
    Display data as an image, i.e., on a 2D regular raster.
    
    The input may either be actual RGB(A) data, or 2D scalar data, which
    will be rendered as a pseudocolor image. For displaying a grayscale
    image set up the colormapping using the parameters
    ``cmap='gray', vmin=0, vmax=255``.
    
    The number of pixels used to render an image is set by the Axes size
    and the *dpi* of the figure. This can lead to aliasing artifacts when
    the image is resampled because the displayed image size will usually
    not match the size of *X* (see
    :doc:`/gallery/images_contours_and_fields/image_antialiasing`).
    The resampling can be controlled via the *interpolation* parameter
    and/or :rc:`image.interpolation`.
    
    Parameters
    ----------
    X : array-like or PIL image
        The image data. Supported array shapes are:
    
        - (M, N): an image with scalar data. The values are mapped to
          colors using normalization and a colormap. See parameters *norm*,
          *cmap*, *vmin*, *vmax*.
        - (M, N, 3): an image with RGB values (0-1 float or 0-255 int).
        - (M, N, 4): an image with RGBA values (0-1 float or 0-255 int),
          i.e. including transparency.
    
        The first two dimensions (M, N) define the rows and columns of
        the image.
    
        Out-of-range RGB(A) values are clipped.
    
    cmap : str or `~matplotlib.colors.Colormap`, default: :rc:`image.cmap`
        The Colormap instance or registered colormap name used to map scalar data
        to colors.
    
        This parameter is ignored if *X* is RGB(A).
    
    norm : str or `~matplotlib.colors.Normalize`, optional
        The normalization method used to scale scalar data to the [0, 1] range
        before mapping to colors using *cmap*. By default, a linear scaling is
        used, mapping the lowest value to 0 and the highest to 1.
    
        If given, this can be one of the following:
    
        - An instance of `.Normalize` or one of its subclasses
          (see :doc:`/tutorials/colors/colormapnorms`).
        - A scale name, i.e. one of "linear", "log", "symlog", "logit", etc.  For a
          list of available scales, call `matplotlib.scale.get_scale_names()`.
          In that case, a suitable `.Normalize` subclass is dynamically generated
          and instantiated.
    
        This parameter is ignored if *X* is RGB(A).
    
    vmin, vmax : float, optional
        When using scalar data and no explicit *norm*, *vmin* and *vmax* define
        the data range that the colormap covers. By default, the colormap covers
        the complete value range of the supplied data. It is an error to use
        *vmin*/*vmax* when a *norm* instance is given (but using a `str` *norm*
        name together with *vmin*/*vmax* is acceptable).
    
        This parameter is ignored if *X* is RGB(A).
    
    aspect : {'equal', 'auto'} or float, default: :rc:`image.aspect`
        The aspect ratio of the Axes.  This parameter is particularly
        relevant for images since it determines whether data pixels are
        square.
    
        This parameter is a shortcut for explicitly calling
        `.Axes.set_aspect`. See there for further details.
    
        - 'equal': Ensures an aspect ratio of 1. Pixels will be square
          (unless pixel sizes are explicitly made non-square in data
          coordinates using *extent*).
        - 'auto': The Axes is kept fixed and the aspect is adjusted so
          that the data fit in the Axes. In general, this will result in
          non-square pixels.
    
    interpolation : str, default: :rc:`image.interpolation`
        The interpolation method used.
    
        Supported values are 'none', 'antialiased', 'nearest', 'bilinear',
        'bicubic', 'spline16', 'spline36', 'hanning', 'hamming', 'hermite',
        'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell',
        'sinc', 'lanczos', 'blackman'.
    
        If *interpolation* is 'none', then no interpolation is performed
        on the Agg, ps, pdf and svg backends. Other backends will fall back
        to 'nearest'. Note that most SVG renderers perform interpolation at
        rendering and that the default interpolation method they implement
        may differ.
    
        If *interpolation* is the default 'antialiased', then 'nearest'
        interpolation is used if the image is upsampled by more than a
        factor of three (i.e. the number of display pixels is at least
        three times the size of the data array).  If the upsampling rate is
        smaller than 3, or the image is downsampled, then 'hanning'
        interpolation is used to act as an anti-aliasing filter, unless the
        image happens to be upsampled by exactly a factor of two or one.
    
        See
        :doc:`/gallery/images_contours_and_fields/interpolation_methods`
        for an overview of the supported interpolation methods, and
        :doc:`/gallery/images_contours_and_fields/image_antialiasing` for
        a discussion of image antialiasing.
    
        Some interpolation methods require an additional radius parameter,
        which can be set by *filterrad*. Additionally, the antigrain image
        resize filter is controlled by the parameter *filternorm*.
    
    interpolation_stage : {'data', 'rgba'}, default: 'data'
        If 'data', interpolation
        is carried out on the data provided by the user.  If 'rgba', the
        interpolation is carried out after the colormapping has been
        applied (visual interpolation).
    
    alpha : float or array-like, optional
        The alpha blending value, between 0 (transparent) and 1 (opaque).
        If *alpha* is an array, the alpha blending values are applied pixel
        by pixel, and *alpha* must have the same shape as *X*.
    
    origin : {'upper', 'lower'}, default: :rc:`image.origin`
        Place the [0, 0] index of the array in the upper left or lower
        left corner of the Axes. The convention (the default) 'upper' is
        typically used for matrices and images.
    
        Note that the vertical axis points upward for 'lower'
        but downward for 'upper'.
    
        See the :doc:`/tutorials/intermediate/imshow_extent` tutorial for
        examples and a more detailed description.
    
    extent : floats (left, right, bottom, top), optional
        The bounding box in data coordinates that the image will fill.
        These values may be unitful and match the units of the Axes.
        The image is stretched individually along x and y to fill the box.
    
        The default extent is determined by the following conditions.
        Pixels have unit size in data coordinates. Their centers are on
        integer coordinates, and their center coordinates range from 0 to
        columns-1 horizontally and from 0 to rows-1 vertically.
    
        Note that the direction of the vertical axis and thus the default
        values for top and bottom depend on *origin*:
    
        - For ``origin == 'upper'`` the default is
          ``(-0.5, numcols-0.5, numrows-0.5, -0.5)``.
        - For ``origin == 'lower'`` the default is
          ``(-0.5, numcols-0.5, -0.5, numrows-0.5)``.
    
        See the :doc:`/tutorials/intermediate/imshow_extent` tutorial for
        examples and a more detailed description.
    
    filternorm : bool, default: True
        A parameter for the antigrain image resize filter (see the
        antigrain documentation).  If *filternorm* is set, the filter
        normalizes integer values and corrects the rounding errors. It
        doesn't do anything with the source floating point values, it
        corrects only integers according to the rule of 1.0 which means
        that any sum of pixel weights must be equal to 1.0.  So, the
        filter function must produce a graph of the proper shape.
    
    filterrad : float > 0, default: 4.0
        The filter radius for filters that have a radius parameter, i.e.
        when interpolation is one of: 'sinc', 'lanczos' or 'blackman'.
    
    resample : bool, default: :rc:`image.resample`
        When *True*, use a full resampling method.  When *False*, only
        resample when the output image is larger than the input image.
    
    url : str, optional
        Set the url of the created `.AxesImage`. See `.Artist.set_url`.
    
    Returns
    -------
    `~matplotlib.image.AxesImage`
    
    Other Parameters
    ----------------
    data : indexable object, optional
        If given, all parameters also accept a string ``s``, which is
        interpreted as ``data[s]`` (unless this raises an exception).
    
    **kwargs : `~matplotlib.artist.Artist` properties
        These parameters are passed on to the constructor of the
        `.AxesImage` artist.
    
    See Also
    --------
    matshow : Plot a matrix or an array as an image.
    
    Notes
    -----
    Unless *extent* is used, pixel centers will be located at integer
    coordinates. In other words: the origin will coincide with the center
    of pixel (0, 0).
    
    There are two common representations for RGB images with an alpha
    channel:
    
    -   Straight (unassociated) alpha: R, G, and B channels represent the
        color of the pixel, disregarding its opacity.
    -   Premultiplied (associated) alpha: R, G, and B channels represent
        the color of the pixel, adjusted for its opacity by multiplication.
    
    `~matplotlib.pyplot.imshow` expects RGB images adopting the straight
    (unassociated) alpha representation.

import time            # timeing
t0 = time.time()
data = Mandelbrot([-2,1,-1,1], 1000, 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  12.632213115692139 s

from numpy import * # because arrays are defined in numpy
from numba import njit  # This is the new line with numba

@njit   # this is an alias for @jit(nopython=True)
def Mand(z0, max_steps):
    z = 0j  # no need to specify type. 
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

@njit
def Mandelbrot2(ext, Nxy, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    data = zeros((Nxy,Nxy)) # initialize a 2D dynamic array
    for i in range(Nxy):
        for j in range(Nxy):
            x = ext[0] + (ext[1]-ext[0])*i/(Nxy-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Nxy-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
    return data
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

import time            # timeing
t0 = time.time()
data = Mandelbrot2(array([-2,1,-1,1]), 1000, 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  1.0402450561523438 s

from numpy import * # because arrays are defined in numpy
from numba import njit  # This is the new line with numba
from numba import prange

@njit   # this is an alias for @jit(nopython=True)
def Mand(z0, max_steps):
    z = 0j  # no need to specify type. 
    # To initialize to complex number, just assign 0j==i*0
    for itr in range(max_steps):
        if abs(z)>2:
            return itr
        z = z*z + z0
    return max_steps

@njit(parallel=True)
def Mandelbrot3(data, ext, max_steps):
    """
    ext[4]    -- array of 4 values [min_x,max_x,min_y,max_y]
    Nxy       -- int number of points in x and y direction
    max_steps -- how many steps we will try at most before we conclude the point is in the set
    """
    Nx,Ny = shape(data) # 2D array should be already allocated we get its size
    for i in range(Nx):
        for j in prange(Ny):    # note that we used prange instead of range.
                                # this switches off parallelization of this loop, so that
                                # only the outside loop over i is parallelized.
            x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.)
            y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.)
            # creating complex number of the fly
            data[i,j] = Mand(x + y*1j, max_steps)  
# data now contains integers. 
# MandelbrotSet has value 1000, and points not in the set have value <1000.

import time            # timeing
data = zeros((1000,1000))
t0 = time.time()
Mandelbrot3(data, array([-2,1,-1,1]), 1000)
t1 = time.time()
print ('clock time: ',t1-t0,'s')

clock time:  0.377485990524292 s

import matplotlib.cm as cm

imshow(-log(data.T), extent=[-2,1,-1,1], cmap=cm.hot, origin='lower');

%%file mandc.cc

#include <iostream> // for printing to stdout
#include <complex>  // for complex numbers
#include <ctime>    // for measuring time
#include <vector>   // for using array/vector
#include <omp.h>    // openMP and measuing time with openMP
using namespace std;// everything from standard template library is available withouth prefix std
// std::cout and std::endl can be called cout, endl, etc...

int Mandelb(const complex<double>& z0, int max_steps)
{
  complex<double> z=0; // specify type where we start using it. (still need to specify type not like python with auto-typing)
  for (int i=0; i<max_steps; i++){ // specify int type inside the loop, so it is local to the loop. useful for reducing bugs.
    if (abs(z)>2.) return i;
    z = z*z + z0;
  }
  return max_steps;
}

int main()
{
  const int Nx = 1000; // constants are named cost rather than parameter
  const int Ny = 1000;
  int max_steps = 1000; // allowed to change later.
  double ext[]={-2,1,-1,1}; // this specifies fixed array of four numbers and initializes it.

  //double mand[Nx][Ny];
  
  vector<int> mand(Nx*Ny);  // allocating one dimensional array (vector) of size Nx*Ny
  // multidimensional dynamic arrays are not standard in C++. One can use pointers to allocate/deallocate memory,
  // and write one's own class interface for that. Or one has to use extension of C++ (blitz++ is excellent).
  // Unfortunately the standard C++ still does not support standard class for that.
  clock_t startTimec = clock();  // cpu time at the start 
  double start = omp_get_wtime(); // wall time at the start

  #pragma omp parallel for
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.); // x in the interval ext[0]...ext[1]
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.); // y in the interval ext[2]...ext[3]
      mand[i*Ny+j] = Mandelb(complex<double>(x,y), max_steps); // storing values in 2D array using 1D array
      //mand[i*Ny+j] = Mandelb(x+y*complex<double>(0.,1.), max_steps); // storing values in 2D array using 1D array
    }
  }

  clock_t endTimec = clock();
  double diffc = double(endTimec-startTimec)/CLOCKS_PER_SEC; // how to get seconds from cpu time
  double diff = omp_get_wtime()-start;                       // openMP time is already in seconds
  
  clog<<"clock time : "<<diffc<<"s"<<" with wall time="<<diff<<"s "<<endl; // printout of time
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
      cout<<x<<" "<<y<<" "<< 1./mand[i*Ny+j] << endl; // prinout of mandelbrot set
    }
  }
}

Overwriting mandc.cc

!{ls -l mandc.cc}

-rw-r--r--  1 haule  staff  2660 Sep  3 17:38 mandc.cc

cmd="g++-12 -fopenmp -Ofast -o mandc mandc.cc"

!{cmd}

!{time mandc > mand.dat}

clock time : 0.720274s with wall time=0.16768s 
mandc > mand.dat  1.24s user 0.92s system 120% cpu 1.802 total

def PlotFile(filename):
    import numpy as np
    import pylab as pl
    from matplotlib.colors import LogNorm

    X,Y,Z = np.loadtxt(filename,unpack=True)
    ext=[X[0],X[-1],Y[0],Y[-1]]
    N = int(np.sqrt(len(Z)))
    dat = Z.reshape(N,N)
    print(ext)
    pl.imshow(dat.T,extent=ext,cmap=cm.hot,norm=LogNorm())   # pylab's function for displaying 2D image

PlotFile('mand.dat')

[-2.0, 1.0, -1.0, 1.0]

%%file mandf.f90

INTEGER Function Mandelb(z0, max_steps)
  IMPLICIT NONE ! Every variable needs to be declared. It is very prudent to use that.
  COMPLEX*16, intent(in) :: z0
  INTEGER, intent(in)    :: max_steps
  ! locals
  COMPLEX*16 :: z
  INTEGER    :: i
  z=0.
  do i=1,max_steps
     if (abs(z)>2.) then ! |f(z)|>2 we know it will explode for large n, hence not part of the set.
        Mandelb = i-1      ! return how many iterations it took to know this is not part of the set.
        return 
     end if
     z = z*z + z0
  end do
  Mandelb = max_steps
  return
END Function Mandelb

program mand
  use omp_lib
  IMPLICIT NONE
  ! external function
  INTEGER :: Mandelb ! Need to declare the external function
  ! locals
  INTEGER :: i, j
  REAL*8  :: x, y
  COMPLEX*16  :: z0
  INTEGER, parameter :: Nx = 1000
  INTEGER, parameter :: Ny = 1000
  INTEGER, parameter :: max_steps = 1000
  REAL*8  :: ext(4) = (/-2., 1., -1., 1./) ! The limits of plotting
  REAL*8  :: mande(Nx,Ny)
  REAL*8  :: start, finish, startw, finishw
  
  call cpu_time(start)
  startw  = OMP_get_wtime()
  
  !$OMP PARALLEL DO  PRIVATE(j,x,y,z0)
  do i=1,Nx
     do j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.) ! x \in [ext(1)...ext(2)]
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.) ! y \in [ext(3)...ext(4)]
        z0 = dcmplx(x,y)
        mande(i,j) = Mandelb(z0, max_steps)
     enddo
  enddo
  !$OMP END PARALLEL DO

  finishw = OMP_get_wtime()
  call cpu_time(finish)
  WRITE(0, '("clock time : ",f6.3,"s  wall time=",f6.3,"s")') finish-start, finishw-startw

  do i=1,Nx
     do j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.)
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.)
        print *, x, y, 1./mande(i,j)
     enddo
  enddo
end program mand

Overwriting mandf.f90

cmd="gfortran -fopenmp -Ofast -o mandf mandf.f90"

!{cmd}

!{time mandf > mandf.dat}

clock time :  0.723s  wall time= 0.169s
mandf > mandf.dat  1.57s user 0.05s system 132% cpu 1.215 total

PlotFile('mandf.dat')

[-2.0, 1.0, -1.0, 1.0]

%%file imanc.cc   

#include "pybind11/pybind11.h"
#include "pybind11/numpy.h"
#include "pybind11/stl.h"
#include <cstdint>

namespace py = pybind11;
using namespace std;

void mand(py::array_t<double>& data, int Nx, int Ny, int max_steps, const vector<int>& ext)
{
  auto dat = data.mutable_unchecked<2>();
  #pragma omp parallel for
  for (int i=0; i<Nx; i++){
    for (int j=0; j<Ny; j++){
      dat(j,i) = max_steps;
      double x = ext[0] + (ext[1]-ext[0])*i/(Nx-1.);
      double y = ext[2] + (ext[3]-ext[2])*j/(Ny-1.);
      complex<double> z0(x,y);
      complex<double> z=0;
      for (int itr=0; itr<max_steps; itr++){
        if (norm(z)>4.){
          dat(j,i) = itr;
          break;
        }
        z = z*z + z0;
      }
    }
  }
}

PYBIND11_MODULE(imanc,m){
  m.doc() = "pybind11 wrap for mandelbrot";
  m.def("mand", &mand);
}

Overwriting imanc.cc

!{ls -l imanc.cc}

-rw-r--r--  1 haule  staff  828 Sep  3 17:39 imanc.cc

#cmd="g++ `python -m pybind11 --includes` -undefined dynamic_lookup -O3 -shared -std=c++11 -fPIC imanc.cc -o imanc.so"
cmd="g++-12 -fopenmp `python -m pybind11 --includes` -undefined dynamic_lookup -O3  -shared -std=c++11 -fPIC imanc.cc -o imanc.so"

!{cmd}

import imanc  # We now import pybind11 module, produced with C++ code

data3 = ones((1000,1000))

t0 = time.time()
imanc.mand(data3, 1000, 1000, 1000, [-2,1,-1,1])
t1 = time.time()

print('pybind11: walltime: ', t1-t0)

pybind11: walltime:  0.20494413375854492

import matplotlib.cm as cm

imshow(-log(data3), extent=ext, cmap=cm.hot)

<matplotlib.image.AxesImage at 0x127bbc510>

%%file mandel.f90

!-----------------------------------------------------
! Produces Mandelbrot plot in the range [-ext[0]:ext[1]]x[ext[2]:ext[3]]
! It uses formula z = z*z + z0 iteratively until
! asb(z) gets bigger than 2 
! (deciding that z0 is not in mandelbrot)
! The value returned is 1/(#-iterations to escape)
!-----------------------------------------------------
SUBROUTINE Mandelb(data, ext, Nx, Ny, max_steps)
  IMPLICIT NONE ! Don't use any implicit names of variables!
  ! Function arguments
  REAL*8, intent(out) :: data(Nx,Ny)
  REAL*8, intent(in)  :: ext(4)              ! [xa,xb,ya,yb]
  INTEGER, intent(in) :: max_steps
  INTEGER, intent(in) :: Nx, Ny
  !f2py integer optional, intent(in)          :: max_steps=1000
  !  !f2py integer intent(hide), depend(data) :: Nx=shape(data,0)  ! it will be hidden automatically
  !  !f2py integer intent(hide), depend(data) :: Ny=shape(data,1)  ! it will be hidden automatically
  ! Local variables
  INTEGER    :: i, j, itt
  COMPLEX*16 :: z0, z
  REAL*8     :: x, y
  data(:,:) = max_steps
  !$OMP PARALLEL DO  PRIVATE(j,x,y,z0,z,itt)
  DO i=1,Nx
     DO j=1,Ny
        x = ext(1) + (ext(2)-ext(1))*(i-1.)/(Nx-1.)
        y = ext(3) + (ext(4)-ext(3))*(j-1.)/(Ny-1.)
        z0 = dcmplx(x,y)
        z=0
        DO itt=1,max_steps
           IF (abs(z)>2.) THEN
              data(i,j) = itt-1 !1./itt         ! result is number of iterations
              EXIT
           ENDIF
           z = z**2 + z0             ! f(z) = z**2+z0 -> z
        ENDDO
        !if (abs(z)<2) data(i,j) = max_steps
     ENDDO
  ENDDO
  !$OMP END PARALLEL DO
  RETURN
END SUBROUTINE Mandelb

Overwriting mandel.f90

!{f2py -c mandel.f90 --f90flags='-fopenmp' -m mandel}

/Users/haule/anaconda3/lib/python3.11/site-packages/numpy/f2py/f2py2e.py:719: VisibleDeprecationWarning: distutils has been deprecated since NumPy 1.26.xUse the Meson backend instead, or generate wrapperswithout -c and use a custom build script
  builder = build_backend(
running build
running config_cc
INFO: unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
INFO: unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
INFO: build_src
INFO: building extension "mandel" sources
INFO: f2py options: []
INFO: f2py:> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandelmodule.c
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11
Reading fortran codes...
	Reading file 'mandel.f90' (format:free)
Post-processing...
	Block: mandel
			Block: mandelb
Applying post-processing hooks...
  character_backward_compatibility_hook
Post-processing (stage 2)...
Building modules...
    Building module "mandel"...
    Generating possibly empty wrappers"
    Maybe empty "mandel-f2pywrappers.f"
        Constructing wrapper function "mandelb"...
          data = mandelb(ext,nx,ny,[max_steps])
    Wrote C/API module "mandel" to file "/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandelmodule.c"
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/fortranobject.c' to sources.
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11' to include_dirs.
copying /Users/haule/anaconda3/lib/python3.11/site-packages/numpy/f2py/src/fortranobject.c -> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11
copying /Users/haule/anaconda3/lib/python3.11/site-packages/numpy/f2py/src/fortranobject.h -> /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11
INFO:   adding '/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.f' to sources.
INFO: build_src: building npy-pkg config files
/Users/haule/anaconda3/lib/python3.11/site-packages/setuptools/_distutils/cmd.py:66: SetuptoolsDeprecationWarning: setup.py install is deprecated.
!!

        ********************************************************************************
        Please avoid running ``setup.py`` directly.
        Instead, use pypa/build, pypa/installer or other
        standards-based tools.

        See https://blog.ganssle.io/articles/2021/10/setup-py-deprecated.html for details.
        ********************************************************************************

!!
  self.initialize_options()
running build_ext
INFO: customize UnixCCompiler
INFO: customize UnixCCompiler using build_ext
INFO: get_default_fcompiler: matching types: '['gnu95', 'nag', 'nagfor', 'absoft', 'ibm', 'intel', 'gnu', 'g95', 'pg']'
INFO: customize Gnu95FCompiler
INFO: Found executable /opt/homebrew/bin/gfortran
INFO: customize Gnu95FCompiler
INFO: customize Gnu95FCompiler using build_ext
INFO: building 'mandel' extension
INFO: compiling C sources
INFO: C compiler: clang -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /Users/haule/anaconda3/include -arch arm64 -fPIC -O2 -isystem /Users/haule/anaconda3/include -arch arm64 -I/usr/local/opt/openjdk@8/include

creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw
creating /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11
INFO: compile options: '-DNPY_DISABLE_OPTIMIZATION=1 -I/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11 -I/Users/haule/anaconda3/lib/python3.11/site-packages/numpy/core/include -I/Users/haule/anaconda3/include/python3.11 -c'
INFO: clang: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandelmodule.c
INFO: clang: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/fortranobject.c
INFO: compiling Fortran sources
INFO: Fortran f77 compiler: /opt/homebrew/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops
Fortran f90 compiler: /opt/homebrew/bin/gfortran -fopenmp -fPIC -O3 -funroll-loops
Fortran fix compiler: /opt/homebrew/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fopenmp -fPIC -O3 -funroll-loops
INFO: compile options: '-I/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11 -I/Users/haule/anaconda3/lib/python3.11/site-packages/numpy/core/include -I/Users/haule/anaconda3/include/python3.11 -c'
INFO: gfortran:f90: mandel.f90
INFO: gfortran:f77: /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.f
INFO: /opt/homebrew/bin/gfortran -Wall -g -Wall -g -undefined dynamic_lookup -bundle /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandelmodule.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/fortranobject.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/mandel.o /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw/src.macosx-11.0-arm64-3.11/mandel-f2pywrappers.o -L/opt/homebrew/Cellar/gcc/14.1.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin23/14 -L/opt/homebrew/Cellar/gcc/14.1.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin23/14/../../.. -L/opt/homebrew/Cellar/gcc/14.1.0_1/bin/../lib/gcc/current/gcc/aarch64-apple-darwin23/14/../../.. -lgfortran -o ./mandel.cpython-311-darwin.so
Removing build directory /var/folders/j8/d9m3r0zx7j37l3ktfl_n1xw00000gn/T/tmpolfvlyhw

import mandel  # importing module created by f2py
import time

ext=[-2,1,-1,1]
Nx = Ny = 1000
max_steps = 1000
tc = time.process_time() # cpu time
tw = time.time()  # wall time
data = mandel.mandelb(ext,Nx,Ny)
print('# wall time : ', time.time()-tw, 's  clock time : ', time.process_time() - tc, 's')

# wall time :  0.21518993377685547 s  clock time :  0.9273019999999974 s

imshow(-log(data.T), extent=ext, cmap=cm.hot)

<matplotlib.image.AxesImage at 0x127b693d0>

Installation¶

Mandelbrot set¶

Optional: Using C++ with `pybind11` to speed up the code¶

Optional installation of C++. These instructions are MAC specific.¶

Optional installation of openMP C++ and gfortran. These instructions are MAC specific.¶

Install Homebrew¶

Install gnu-gcc and gfortran¶

(Optional) Install gnuplot and gsl¶

C++ and fortran with openMP implementation of the same algorithm¶

Return to `pybind11` to speed up the code¶

`f2py` as the alternative to `pybind11`¶

Popularity of programing languages¶

Installation¶

Mandelbrot set¶

Optional: Using C++ with pybind11 to speed up the code¶

Optional installation of C++. These instructions are MAC specific.¶

Optional installation of openMP C++ and gfortran. These instructions are MAC specific.¶

Install Homebrew¶

Install gnu-gcc and gfortran¶

(Optional) Install gnuplot and gsl¶

C++ and fortran with openMP implementation of the same algorithm¶

Return to pybind11 to speed up the code¶

f2py as the alternative to pybind11¶

Popularity of programing languages¶

Optional: Using C++ with `pybind11` to speed up the code¶

Return to `pybind11` to speed up the code¶

`f2py` as the alternative to `pybind11`¶