[Date Index] [Thread Index] [Author Index]

Re: simplifying ulam spiral code

• To: mathgroup at smc.vnet.net
• Subject: [mg56435] Re: simplifying ulam spiral code
• From: "Carl K. Woll" <carlw at u.washington.edu>
• Date: Mon, 25 Apr 2005 01:30:51 -0400 (EDT)
• Organization: University of Washington
• References: <d4cmlr$39e$1@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

"zak" <chocolatez at gmail.com> wrote in message 
news:d4cmlr$39e$1 at smc.vnet.net...
> hi
> in the site:
> http://teachers.crescentschool.org/weissworld/m3/spiral/spiral.html
> there is a mathematica code for drawing ULAM'S SPIRAL
> the code is:
>
[snip]

I have to confess that I don't understand how zak's code relates to the 
above link. The text at the link says that 1 is placed at the origin, 2 is 
placed to the right of 1, and succeeding integers are placed in a 
counterclockwise spiral. Hence, 3 ought to have the coordinates {1,1}, 4 
ought to have the coordinates {0,1}, etc.

At any rate, it is not too difficult to program a function to determine the 
coordinates of an integer in the above spiral. If one notices that the 
bottom right and top left corners have the integer values n^2+1, then one 
eventually gets

coords[k_Integer] := Module[{n, a, b},
    n = Floor[Sqrt[k - .9]];
    a = k - n^2 - n - 1;
    b = Quotient[2n + 1 - (-1)^n, 4];
    (-1)^n {Quotient[Abs[a] + a, 2] - b, Quotient[Abs[a] - a, 2] - b}]

coords[k_List] := Module[{n, a, b, c},
    n = Floor[Sqrt[k - .9]];
    a = k - n^2 - n - 1;
    c = (-1)^n;
    b = Quotient[2n + 1 - c, 4];
    c Transpose[{Quotient[Abs[a] + a, 2] - b, Quotient[Abs[a] - a, 2] - b}]]

A few comments may be in order. Concentrating on the second function 
definition, I used a number of ideas to speed up it's execution. I used 
Sqrt[k-.9] so that Mathematica is taking square roots of real numbers 
instead of integers, which perhaps surprisingly is much faster. I used 
Sqrt[k-.9] instead of Sqrt[k-1.] to avoid spurious cancellation errors when 
the Floor of the result is evaluated. I wanted to make sure that all my 
arrays were packed, so I used Quotient instead of dividing two integers. 
Even though 2n+1-(-1)^n is always divisible by 4, (2n+1-(-1)^n)/4 is not 
packed even when 2n+1-(-1)^n is packed. Finally, I used (Abs[a]+a)/2 (with 
Quotient instead of using division) to change all negative values in the 
list a to 0 and (Abs[a]-a)/2 to change all positive values in the list a to 
0.

At any rate, using coords on a list of the first million integers takes a 
bit less than 2 seconds on my machine.

Now, we are ready to used coords to find the coordinates of the primes. For 
example, if we are interested in the first million primes:

data = Prime[Range[10^6]]; // Timing
{9.547 Second, Null}

Now, we use coords to get the coordinates of the primes.

pts = coords[data]; // Timing
{2.375 Second, Null}

Applying Developer`ToPackedArray to the data would speed up the pts 
computation by a bit. Looking at the first few points of pts reveals

Take[pts, 10]
{{1, 0}, {1, 1}, {-1, 1}, {-1, -1}, {2, 0}, {2, 2}, {-2, 2}, {-2, 0}, 
{0, -2}, {3, 1}}

which looks correct to me.

Carl Woll 

• Follow-Ups:
  • Re: Re: simplifying ulam spiral code
    • From: DrBob <drbob@bigfoot.com>