in reply to Fractal structure: PDL, memory, time

Try this. It's correct at levels 1, 2, & 3, and Ithink it is correct for all levels, but it gets harder to discern.

I found this really tough to get right, whether through recursion or iteration. I'd given up and moved on before the pattern dawned on me. I think a mathematician may be able to reduce this to a formula, maybe involving some kind of polynomial, but it's beyond my abilities. 

#! perl -slw
use strict;

$|++;

our $DEPTH ||= 3;

sub f{
    use integer;
    my( $x, $y, $depth ) = @_;
    my $bsize = 4**( $depth-1 );
    {
        if( $x == $y ) {
            return $depth;
        }
        elsif( $x/4 == $y/4 ) {
            return $depth-1;
        }
        elsif( $x/$bsize == $y/$bsize ) {
            $x %= 4;
            $y %= 4;
            --$depth;
            redo;
        }
        else {
            $x %= $bsize;
            $y %= $bsize;
            --$depth;
            redo;
        }
    }
}

for my $depth ( $DEPTH ) {
    for my $y ( 0 .. 4**$depth-1 ) {
        my $c = 0;
        for my $x ( 0 .. 4 ** $depth -1 ) {
            printf +( ( ++$c % 4 ) == 0 ? " %2d " : " %2d" ),
                f( $x, $y, $depth );
        }
        print +( $y % 4 ) == 3 ? "\n" : '';
    }
    print '';
}

__END__

[download]
Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
"Science is about questioning the status quo. Questioning authority".
In the absence of evidence, opinion is indistinguishable from prejudice.

Replies are listed 'Best First'.
Re^2: Fractal structure: PDL, memory, time
by FFRANK (Beadle) on Jul 23, 2007 at 01:45 UTC
    Hi BrowserUK,

    I tried also to get your code working, and could'nt. I verified level 4 after your update and it does not yet match - although close ! (if you compare w/ output from PDL fractal code).

    I'll meet a mathematicien and if I find something better, I'll post in this thread.

    I'm also looking at simpler solutions using nDeep arrays.

    Thanks very much,

    Frank

[reply]

      Try these. fc() is just an Inline C version of f().

      I'm pretty sure that they are now correct, but you should verify them yourself. 

      #! perl -slw
use strict;
use Inline C => << '__C__', NAME => "_627288"; # => CLEAN_AFTER_BUILD 
+=> 0;
#include <stdio.h>

int fc( double x, double y, int d ) {
    double bsize = pow( 4, ( d-1 ) );

    if( x == y )
        return d;
    else if( (int)( x / 4 ) == (int)( y / 4 ) )
        return d-1;
    else if( (int)( x / bsize ) == (int)( y / bsize ) )
        return 1 + fc( fmod( x, bsize ), fmod( y, bsize ), d-1 );
    else
        return     fc( fmod( x, bsize ), fmod( y, bsize ), d-1 );
}

__C__

sub f{
    my( $x, $y, $d ) = @_;
    my $bsize = 4**($d-1);

    if( $x == $y ) {
        return $d
    }
    elsif( int( $x / 4 ) == int( $y / 4 ) ) {
        return $d-1;
    }
    elsif( int( $x / $bsize ) == int( $y / $bsize ) ) {
        return 1 + f( $x % $bsize, $y % $bsize, $d-1 );
    }
    else{
        return     f( $x % $bsize, $y % $bsize, $d-1 );
    }
}

      [download]

      Using this I have produced plots of 1000x1000 chunks upto level 31 successfully. At that point the coordinates get so large that they blow the limits of doubles. I'm also fairly sure that you are loosing precision long before level 31. At level 16 everything appears to check out.

      I don't think any single set of 20 lines of code has ever (in the best part of 30 years), given me as much trouble to wrap my brain around as this!

      Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
      "Science is about questioning the status quo. Questioning authority".
      In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]
[d/l]
[select]
        Hi BrowserUK,

        This code is fanstastic! I am impressed.

        I hope to be able to fit it to my needs (indices), which I've been dealing with so far by nesting structures, and not completely satisfied with the runtime, maybe this code will be faster, or something inspired by it...

        Frank :-)

[reply]
      I verified level 4 after your update and it does not yet match - although close !

      Gah. I See where this is headed. I add another ifelse clause to get level 4 right, and then another for level 5...

      This has to be possible recursively! It's amazing how something so complex evolves out of such a simple process. It's also amazing how something with so many symmetries and essentially only one dissymmetry could be so hard to generate.

      Part of the problem is the shear size the things grow to so rapidly. I wrote a simple png plotter that produces a png of a 1000x1000 chunk and allows you to 'move' around and zoom in--very slowly. But as the level increases, it becomes impossible to keep track of where you are, so it doesn't help much.

      Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
      "Science is about questioning the status quo. Questioning authority".
      In the absence of evidence, opinion is indistinguishable from prejudice.
[reply]

In Section Seekers of Perl Wisdom