Showing posts with label math. Show all posts
Showing posts with label math. Show all posts

Sunday, April 12, 2015

Mandelbrot Set Rendering

Benoit Mandelbrot was a Polish-born mathematician. His most well known work is in the field of fractals. He came up with a function describing a set of complex numbers called the Mandelbrot Set.
The Mandlebrot set describes a set of numbers "c" such that the sequence (c1, c2 .... cx) where cn+1 = (cn*cn)+cn does not approach infinity.
So you're looking at c, c²+c, (c²+c)²+c, ((c²+c)²+c)²+c...

A complex number has two parts. It's a sum of real an imaginary components. Let's say we're looking at A + Bi.
A is the real component.
Bi is the imaginary component.
'i' is a symbol representing the √-1 (square root of -1)
i.e., (Bi)² would be -B².

Mandelbrot's equation for describing the model is z'=z²+c.
Let z be our A + Bi from above.
The next element of the sequence is 
If z=(A+Bi)² then z=A²-B²+2*A*Bi


If plotting, we can us X, Y in place of A, and B.
so, z=X²-Y²+2*X*Y(i)
We have real and imaginary components here.
So Xnew=X²-Y², and Ynew=2*X*Y
the "c" is the starting point for Mandelbrot's equation.
Z=Z²+c, so we need to add X0 and Y0 to the above:

Xnew=X²-Y²+X0
Ynew=2*X*Y+Y0
And from here, we keep iterating on this till we hit set limits to number of iterations, or if X/Y are sufficiently large.

So, for a point X, Y,
Repeat the following until x0²+y0²>4 or until set max iterations

Until (max iterations, or x0²+y0²>4)
  x0=0 and y0=0
  x1=x0*x0-y0*y0+X
  y1=2*x0*y0+Y
  x0=x1
  y0=y1
Color point X,Y based on iterations.




Point X,Y is then associated with the number of iterations to break out of that loop.

Using this we can construct a basic Mandelbrot generator. I did this using Javascript, since it's easy to use any browser to display a canvas element. This works in Firefox and Chrome (has worked for a while) and I've tested it in IE11 and it works there too.




<!DOCTYPE HTML>
<html>
<script language="javascript">
  function rgbToHex(r, g, b) { //thanks stackoverflow
    return "#" + ((1 << 24) + (r << 16) + (g << 8) + b).toString(16).slice(1);
  }
  function putpixel(x,y,r,g,b,ctx){
    ctx.fillStyle = rgbToHex(r,g,b);
    ctx.fillRect(x,y,1,1);
  }
  function paint(cnvname){
    mcanvas = document.getElementById(cnvname);
    ctx = mcanvas.getContext('2d');
    ctx.clearRect(0, 0, 799, 599);
    for (x = 0; x < 799; x++) {
      for (y = 0; y < 599; y++) {
        tx=(x-400)/(200); // tx is "translated x" for graph centering"
        ty=(300-y)/(200); // ty is "translated y" for graph centering"
        xtmp=0;
        ytmp=0;
        xtmp2=0;
        iter=0
        maxiter=1000;
        for (; (iter<maxiter) && (xtmp*xtmp+ytmp*ytmp < 4);) {
          xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
          ytmp=2*xtmp*ytmp+ty;
          xtmp=xtmp2;
          iter++;
        }
        contrast=0.15
        color=Math.floor( (Math.pow(iter,contrast)*255)/Math.pow(maxiter,contrast));
        putpixel(x,y,0,255-color,255-color,ctx);
      }
    }
  }
</script>
<body onload='paint("mcanvas")'>
  <div>
    <canvas id="mcanvas" width="800" height="600"></canvas>
  </div>
</body>
</html>


This is a simple bit of javascript that iterates over the pixels of the canvas created in the body. It applies the function described earlier to determine the max number of iterations to the set.



Pretty right? Of course there's much more to mandelbrot than just the bulb.
We can zoom in on a portion of this by applying a zoom and x/y offset to the above code.

We can add a zoom and offset feature to explore the set.


  for (x = 0; x < 799; x++) {
    for (y = 0; y < 599; y++) {
      zoom=20; xoff=0.6; yoff=0.6;
      tx=-xoff+(x-400)/(zoom*200); // tx is "translated x" for graph centering"
      ty=yoff+(300-y)/(zoom*200); // ty is "translated y" for graph centering"

The above change will produce the following:

This is zoomed in 20x onto the small bulb in the top left portion of the largest bulb. Note the pattern replicates, and there are several tiny bulbs that look a lot like the main one (with slight variations).

Now the fun part - we can color based on x, y or iterations. Apply a sinusoidal function to iterations for red, blue and green, and you can have a never ending palette of color.

e.g.
From the contrast and color line:

//contrast=0.25
//color=Math.floor( (Math.pow(iter,contrast)*255)/Math.pow(maxiter,contrast));
cycle=(maxiter)/(3.14159265359*100);
cr=Math.floor(128+127*Math.sin(0.00+iter/(cycle*1.5000)));
cg=Math.floor(128+127*Math.sin(1.57+iter/(cycle*1.6000)));
cb=Math.floor(128+127*Math.sin(3.14+iter/(cycle*1.5666)));
putpixel(x,y,cr,cg,cb,ctx);

So, what are we saying? A cycle is 2*pi. We're dividing the iterations by 100*pi, so we're looking at at 50 cycles. Red, green and blue values are sinusoidal cycles with different start points - red starts at 128+127*sin(0) (so 128, and increasing). Green is 128+127*sin(1.57) (so 255, 1.57 being half pi, sin(pi/2) is 1) - i.e. max green and heading down to 0. Blue is starting at 128+sin(pi), which means middle intensity and getting darker.
Each sinusoidal cycle has it's own multiplier so each cycle has it's own rate, giving many more colors ans each cycle creates a new palette.

Here's the result:


We're looking at integers since we're only counting max iterations. We can smoothen this by looking at how far past our boundary condition each iteration went.

My script was based on methods described in wikipedia.


log2=Math.log(2)
if (iter<maxiter){ //smoother transitions
  zn=Math.sqrt(xtmp*xtmp+ytmp*ytmp);
  zl=Math.log(Math.log(zn)/log2)/log2;
  iter+=1-zl;
} 






Now consider the interior of the bulb. It's what happens when the iterations have been maxed. We can color the interior by using the same palette. I haven't quite figured out the how to map distance, but we can use the existing xtmp and ytmp variables. Since iterating more doesn't get us greater than the boundary condition, we can assume xtmp and ytmp are small. 

What I tried was several permutations of setting iter to 1/xtmp and 1/ytmp. I settled on using iter=Math.log( (1/(xtmp*xtmp+ytmp*ytmp)) )/(log2);
The idea is that if my values aren't going to break out of the boundary, then they're converging. So I'm using the inverse of the values to dictate what they will be against the existing palette. Not all sub-bulbs show this though.


Here's the overall code for this basic mandelbrot renderer...



<!DOCTYPE HTML>
<html>
<script language="javascript">
  function rgbToHex(r, g, b) { //thanks stackoverflow
    return "#" + ((1 << 24) + (r << 16) + (g << 8) + b).toString(16).slice(1);
  }
  function putpixel(x,y,r,g,b,ctx){
    ctx.fillStyle = rgbToHex(r,g,b);
    ctx.fillRect(x,y,1,1);
  }
  function paint(cnvname){
    mcanvas = document.getElementById(cnvname);
    ctx = mcanvas.getContext('2d');
    ctx.clearRect(0, 0, 799, 599);
    log2=Math.log(2)
    for (x = 0; x < 799; x++) {
      for (y = 0; y < 599; y++) {
        zoom=20; xoff=0.6; yoff=0.6;
        tx=-xoff+(x-400)/(zoom*200); // tx is "translated x" for graph centering"
        ty=yoff+(300-y)/(zoom*200); // ty is "translated y" for graph centering"
        xtmp=0;
        ytmp=0;
        xtmp2=0;
        iter=0
        maxiter=1000; // increase with larger zoom values
        totes=0;
        for (; (iter<maxiter) && (xtmp*xtmp+ytmp*ytmp < 4);) {
          xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
          ytmp=2*xtmp*ytmp+ty;
          xtmp=xtmp2;
          iter++;
        }
        if (iter<maxiter){ //smoother transitions
          zn=Math.sqrt(xtmp*xtmp+ytmp*ytmp);
          zl=Math.log(Math.log(zn)/log2)/log2;
          iter+=1-zl;
        } else { // for internal color of Mandelbulb
          iter=Math.log( (1/(xtmp*xtmp+ytmp*ytmp)) )/(log2);
        }
        cycle=(maxiter)/(3.14159265359*100); //50 cycles from 0-maxiter
        cr=Math.floor(128+127*Math.sin(0.00+iter/(cycle*1.5000)));
        cg=Math.floor(128+127*Math.sin(1.57+iter/(cycle*1.6000)));
        cb=Math.floor(128+127*Math.sin(3.14+iter/(cycle*1.5666)));
        putpixel(x,y,cr,cg,cb,ctx);
      }
    }
  }
</script>
<body onload='paint("mcanvas")'>
  <div>
    <canvas id="mcanvas" width="800" height="600"></canvas>
  </div>
</body>
</html>


Manipulation of the color palette constants, zoom, offsets, and even the equation can produce some beautiul results.

For example,
Set the initial xtmp to tx and ytmp to ty.

change:
xtmp2=xtmp*xtmp-ytmp*ytmp+tx;
ytmp=2*xtmp*ytmp+ty;


to:
xtmp2=xtmp*xtmp-ytmp*ytmp-0.32;
ytmp=2*xtmp*ytmp+0.6575;


Now you're making a Julia plot.
This is another fractal you can zoom in continuously.

Let's look at zooming in on a portion of the Mandelbrot set.
zoom=1.00; xoff=0; yoff=0;


zoom=5; xoff=0.7348; yoff=0.1798;

zoom=25; xoff=0.8; yoff=0.2;



zoom=125; xoff=0.79; yoff=0.158;


zoom=625; xoff=0.79; yoff=0.162;

zoom=3125; xoff=0.7893; yoff=0.163;


At this point we need to increase the maxiter constant. I raised this to 10000.
zoom=9125; xoff=0.78937; yoff=0.16307;


 zoom=15625; xoff=0.78940; yoff=0.163055;


zoom=78225; xoff=0.78940; yoff=0.163055;
Note the zoom - we're at almost 80000x and as you can see there's a lot more structures here!

Friday, September 6, 2013

Euler 3 continued:

Here's the code for Euler 3 in Python (written in 2.7)

Runs quicker than i was expecting. For a large prime like "600851475143987", it completes it's run in under 10 seconds on my ancient core 2 duo work machine.

Since the number is prime, the list of prime factors is just itself.
[600851475143987L]
9.84399986267 seconds...


An unoptimized version that doesn't employ that bit of array arithmetic (2*index+1=possible prime factor) takes about 15 seconds, so it looks quirky, but these optimizations count!


unoptimized prime factor was:

def lowestprimefactor_bysieve(number):
    sieve=[0]*int(round(1+pow(number,0.5)))
    init=2
    k=2
    sieve[k]=2
    while (k < len(sieve)):
        if (sieve[k]!=1):
            sieve[k]=2 #number is prime
            if ((number%k==0)):
                return int(k)
            for idx in xrange(k,len(sieve),k):
                sieve[idx]=1 #mark as a being a factor
        if k>2:
            k+=2
        else:
            k+=1
    return int(number)



Here's the full code:


import sys
import time

def lowestprimefactor_bysieve(number):
    sieve=[0]*int(round(1+pow(number,0.5)/2))
    init=2
    if (number%2==0):
        return 2
    k=1
    sieve[k]=2
    while (k < len(sieve)):
        if (sieve[k]!=1):
            realvalue=k*2+1
            sieve[k]=2 #number is prime
            if ((number%realvalue==0)):
                return int(realvalue)
            for idx in xrange(k,len(sieve),realvalue):
                sieve[idx]=1 #mark as a being a factor
        k+=1
    return int(number)

def primefactorlist_sieve(number):
    lfac=lowestprimefactor_bysieve(number)
    if lfac==number:
        return [lfac]
    return [lfac]+primefactorlist_sieve(number//lfac)


def main():
    #What is the largest prime factor of the number 600851475143 ?
    number=600851475143
    start = time.time()
    print(primefactorlist_sieve(number))
    print(time.time() - start)

if __name__ == '__main__':
    main()

Euler 3 - prime numbers

http://projecteuler.net/problem=3

Hi all, here's where the projects start requiring some thought, and often some research.

What is the largest prime factor of the number 600851475143?



Here's the gist.
it's easier to find the lowest prime factor, divide the number by it, and repeat.

e.g. if a number X is the product of prime factors A x B x C
Then, if the lowest prime factor is A, you now have X/A = B x C
Repeat finding the lowest prime factor of X/A and you'll get B.
Now you have to find the lowest prime factor is (X/A)/B.

X=AxBxC
X/A=(A/A)xBxC
(X/A)/B=(A/A)x(B/B)xC
(X/A)/B =C <= the largest prime factor.

So, how do you find a prime factor? well, for starters, it's a prime number.
We start by looking through a list of numbers - 2, 3, 5, 7, 11, 13, 17 etc.. to see if X is divisible by any.

Lets start by setting boundaries - how high should we count?

If we want the prime factors of 91 (i.e. 7x13).
We wont get past 7 before dividing.  i.e. we'll count up to the lowest prime factor. 

AxB=C. 

If A gets bigger, B must get smaller. So we'll iterate up to B and then divide according to the algorithm. The maximum iterations will then be when A=B and this is the square root of C.

The largest prime factor will be the number itself - i.e. it is prime. Well, how do we know? Once again, we only need to count up to the square root.
e.g. what's the prime factors of 31?
we'll check, 2, 3, 5, 7- we've went past the square root and stop here. Why? Well, since we've past the square root, the other number must be smaller - and if prime would have been counted already! 31/7 ~ 4.43. We've already checked smaller primes, so we know they aren't going to divide it.

The checks look like this:
AxB=C
2x31/2=31...2x15.50=31...A=2, B=15.5
3x31/3=31...3x10.33=31...A=3, B=10.33
5x31/5=31...5x 6.20=31...A=5, B=6.2
7x31/5=31...7x 4.43=31...A=7, B=4.43

So stop checking. Increasing A just makes B cover the number range we already went through, so the number is prime.

So - iterate over primes to find the first prime factor, range of primes from 2 to the square root of the number.

Great - how do you know what IS prime anyway? You're testing primes only right? why not test 4 or 9 and see if out number is divisible by that?

Well, here's where we create an array to keep in mind what's prime or not.
e.g. 
Lets say we wanted to know the prime factors of 221. The Square root is ~14.9, so we only need to check up to there.

01 02 03 04 05 06 07 08 09 10 11 12 13 14 15

00 00 00 00 00 00 00 00 00 00 00 00 00 00 00

Prime numbers and the number 2, are odd. So do a check for 2, then iterate over the odds. no point spending half your time checking even numbers. if it's not divisible by 2, then it's not divisible by any even number.
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
00 xx 00 00 00 00 00 00 00 00 00 00 00 00 00

Now check 3. It's not divisible by 3, so mark off every multiple of 3 as not a factor.
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
00 xx xx 00 00 xx 00 00 xx 00 00 xx 00 00 xx
See what happened? Now I didn't mark off evens before since there's no point of bothering to check

Last check was 3. Add 2, and check 5 - was 5 marked? No? ok! 5 is not a factor, we lose 5 and 10.
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
00 xx xx 00 xx xx 00 00 xx xx 00 xx 00 00 xx

7 Isn't a factor, so we mark out 7 and 14.
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
00 xx xx 00 xx xx xx 00 xx xx 00 xx 00 xx xx


Next odd number is 9 but that's marked. so skip it. next odd number is 11.
01 02 03 04 05 06 07 08 09 10 11 12 13 14 15
00 xx xx 00 xx xx xx 00 xx xx xx xx 00 xx xx


Not a factor.  Next is 13 - that's a factor!

Factors are 13 + the factors of 221/13 - which means it's 13 and factors of 17. We'd do the same check on 17 as illustrated above and return no prime factor found (meaning the number itself is prime).


Now optimizations. We can check 2 initially, and there's no need to store even number flags if they're not being checked.

So our array can be half the size.
00 01 02 03 04 05 06 07 08 ==>N
03 05 07 09 11 13 15 17 19 ==>2N+3

If we're looking at a possible prime factor in the array - say 11, it corresponds to (11-3)/2 = 4 in the index for the flag.

So we're only counting up to the square root of a number, and we're skipping every even number to look for a factor.