Mathematica - simplifying ulam spiral code

Author

simplifying ulam spiral code

zak

2005-04-23, 3:59 am

hi
in the site:
http://teachers.crescentschool.org/...ral/spiral.html
there is a mathematica code for drawing ULAM'S SPIRAL
the code is:

dat = Table[i, {i, 20000}];
x := 0; y := 0; num := 1; shift := 1;
point := {x, y}; list = {}
While[num <= 20000,
i := 1;
While[i <= shift,
x = x + 1;
If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
i++;]; i := 1; While[i <= shift,
y = y + 1;
If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
i++;]; shift++;
i := 1;
While[i <= shift,
x = x - 1;
If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
i++;]; i := 1;
While[i <= shift,
y = y - 1;
If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
i++;]; shift++;]
ListPlot[Partition[list, 2], PlotStyle -> PointSize[0.007],
AspectRatio -> Automatic, Axes -> True];

or see text in:
http://sr1.mytempdir.com/15902

Please i hope someone transform it to a functional programming code.
thanks
zak

Peter Pein

2005-04-24, 8:58 am

Daniel Lichtblau

2005-04-24, 8:58 am

zak wrote:
> hi
> in the site:
> http://teachers.crescentschool.org/...ral/spiral.html
> there is a mathematica code for drawing ULAM'S SPIRAL
> the code is:
>
> dat = Table[i, {i, 20000}];
> x := 0; y := 0; num := 1; shift := 1;
> point := {x, y}; list = {}
> While[num <= 20000,
> i := 1;
> While[i <= shift,
> x = x + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1; While[i <= shift,
> y = y + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;
> i := 1;
> While[i <= shift,
> x = x - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1;
> While[i <= shift,
> y = y - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;]
> ListPlot[Partition[list, 2], PlotStyle -> PointSize[0.007],
> AspectRatio -> Automatic, Axes -> True];
>
>
> or see text in:
> http://sr1.mytempdir.com/15902
>
>
> Please i hope someone transform it to a functional programming code.
> thanks
> zak

The below possibilities are not in functional programming style. I'm not
convinced that would be a good approach for this problem. These all
agree with one another but I did not check whether they agree with your
code above; there may be some fiddling to do with the initial i and j
values. Also I note that the code above can (and will) run off the end
of the data, giving several messages to that effect.

ulamSpiral1[len_] := Module[
{dat=Range[len], x=0, y=0, num=1, shift=1, j=1,
xyincr={{1,0},{0,1},{-1,0},{0,-1}}, shiftincr = {0,1,0,1}},
Reap[While[num<=len,
Do[If[num>len,Break[]]; {x,y} += xyincr[[Mod[j,4,1]]];
If[PrimeQ[dat[[num]]], Sow[{x,y}]]; ++num;
,{i,shift}];
shift += shiftincr[[Mod[j,4,1]]];
j++;
]][[2,1]]
]

ulamSpiral2[len_] := Module[
{dat=Range[len], x=0, y=0, shift=1, j=0, i=0, big={1000,1000},
xyincr={{1,0},{0,1},{-1,0},{0,-1}}, shiftincr={0,1,0,1}},
Select[Table[
i++;
If [i>=shift, i = 1; shift += shiftincr[[Mod[j,4,1]]]; j++;];
{x,y} += xyincr[[Mod[j,4,1]]];
If[PrimeQ[dat[[num]]], {x,y}, big],
{num,len}], Abs[First[#]]<1000&]
]

ulamSpiral3 = Compile[{{len,_Integer}}, Module[
{dat=Range[len], x=0, y=0, shift=1, i=0, j=0,
xincr = {1,0,-1,0}, yincr = {0,1,0,-1}, shiftincr={0,1,0,1}},
Select[Table[
i++;
If [i>=shift, i = 1; shift += shiftincr[[Mod[j,4,1]]]; j++;];
x += xincr[[Mod[j,4,1]]]; y += yincr[[Mod[j,4,1]]];
If[PrimeQ[dat[[num]]], {x,y}, {1000,1000}],
{num,len}], Abs[First[#]]<1000&]
],
{{PrimeQ[_],True|False}}];

In[31]:= Timing[list1 = ulamSpiral1[20000];]
Out[31]= {0.410938 Second, Null}

In[32]:= Timing[list2 = ulamSpiral2[20000];]
Out[32]= {0.470928 Second, Null}

In[33]:= Timing[list3 = ulamSpiral3[20000];]
Out[33]= {0.100985 Second, Null}

In[34]:= list1===list2===list3
Out[34]= True

Daniel Lichtblau
Wolfram Research

DrBob

2005-04-24, 8:59 am

I wouldn't call this "functional", but it's an improvement:

x = 0;
test = list = First@Last@Reap@
While[(shift = 1 - 2x; num = shift(shift - 1)) ≤ 20000,
Do[PrimeQ[i] && Sow@{x +
i - num, x}, {i, num + 1, num + shift}];
num = shift^2;
Do[PrimeQ[
i] && Sow@{x + shift, x + i - num}, {i, num + 1, num + shift}];
num = shift(shift + 1);
Do[PrimeQ[i] && Sow@{x + shift - i + num, x +
shift}, {i, num + 1, num + shift + 1}];
num = (shift + 1)^2;
Do[PrimeQ[i] && Sow@{x - 1, x + shift - i + num}, {i, num +
1, num + shift + 1}];
x--];
ListPlot[list, PlotStyle -> PointSize[
0.007], AspectRatio -> Automatic, Axes -> True];

Bobby

On Sat, 23 Apr 2005 01:16:49 -0400 (EDT), zak <chocolatez@gmail.com> wrote:

> hi
> in the site:
> http://teachers.crescentschool.org/...ral/spiral.html
> there is a mathematica code for drawing ULAM'S SPIRAL
> the code is:
>
> dat = Table[i, {i, 20000}];
> x := 0; y := 0; num := 1; shift := 1;
> point := {x, y}; list = {}
> While[num <= 20000,
> i := 1;
> While[i <= shift,
> x = x + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1; While[i <= shift,
> y = y + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;
> i := 1;
> While[i <= shift,
> x = x - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1;
> While[i <= shift,
> y = y - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;]
> ListPlot[Partition[list, 2], PlotStyle -> PointSize[0.007],
> AspectRatio -> Automatic, Axes -> True];
>
>
> or see text in:
> http://sr1.mytempdir.com/15902
>
>
> Please i hope someone transform it to a functional programming code.
> thanks
> zak
>
>
>
>

--
DrBob@bigfoot.com

DrBob

2005-04-24, 8:59 am

This is somewhat "functional" and agrees (I think) with the picture at the link (depending on where we start, which isn't actually revealed).

Clear[side, ell]
next = {{0, -1}, {1, 0}};
initial = {start = {0, 0}, {1, 0}, 1};
side[i_]@{corner_, direction_, count_} :=
{Nest[Sow[# + direction] &, corner, count],
next.direction, count + i}
ell@{corner_, direction_, count_} :=
side[1]@side[0]@{corner, direction, count}
limit = Ceiling[(1 + Sqrt[1 + 4*20000])/2];
list = Pick[#, PrimeQ@Range@Length@#] &
@Prepend[First@Last@Reap@Nest[ell, initial, limit], start];
ListPlot[list, PlotStyle -> PointSize[0.007],
AspectRatio -> Automatic, Axes -> True];

The posted code gives {1,1} and {0,1} as the first two colored squares, which means the direction followed on the second side is {-1,0}, since {1,1}+{-1,0}=={0,1}. The picture at the link would have {0,1} as the direction taken on the second side.

Possibly the posted code winds clockwise or plots x and y reversed, or something.

Bobby

On Sat, 23 Apr 2005 01:16:49 -0400 (EDT), zak <chocolatez@gmail.com> wrote:

> hi
> in the site:
> http://teachers.crescentschool.org/...ral/spiral.html
> there is a mathematica code for drawing ULAM'S SPIRAL
> the code is:
>
> dat = Table[i, {i, 20000}];
> x := 0; y := 0; num := 1; shift := 1;
> point := {x, y}; list = {}
> While[num <= 20000,
> i := 1;
> While[i <= shift,
> x = x + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1; While[i <= shift,
> y = y + 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;
> i := 1;
> While[i <= shift,
> x = x - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; i := 1;
> While[i <= shift,
> y = y - 1;
> If[PrimeQ[dat[[num]]] == True, list = Flatten[{list, point}]]; num++;
> i++;]; shift++;]
> ListPlot[Partition[list, 2], PlotStyle -> PointSize[0.007],
> AspectRatio -> Automatic, Axes -> True];
>
>
> or see text in:
> http://sr1.mytempdir.com/15902
>
>
> Please i hope someone transform it to a functional programming code.
> thanks
> zak
>
>
>
>

--
DrBob@bigfoot.com

Carl K. Woll

2005-04-25, 4:02 am

"zak" <chocolatez@gmail.com> wrote in message
news:d4cmlr$39e$1@smc.vnet.net...
> hi
> in the site:
> http://teachers.crescentschool.org/...ral/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

Carl K. Woll

2005-04-26, 4:06 am

"Carl K. Woll" <carlw@u.washington.edu> wrote in message
news:d4hvsv$1jp$1@smc.vnet.net...
> "zak" <chocolatez@gmail.com> wrote in message
> news:d4cmlr$39e$1@smc.vnet.net...
> [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
>
>

I see now that I started with 1 at the origin, and others started with 0 at
the origin. Hence, a small modification of my code will return results that
agree with e.g., Daniel Lichtblau.

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

Basically, I substituted k+1 for k. I also got rid of Power by using n(n+1)
to evaluate n^2+n and I by using Mod to find (-1)^n. These latter changes
increased the speed of coords by about 33%.

Carl Woll

DrBob

2005-04-26, 4:06 am

Brilliant!

Bobby

On Mon, 25 Apr 2005 01:30:51 -0400 (EDT), Carl K. Woll <carlw@u.washington.edu> wrote:

> "zak" <chocolatez@gmail.com> wrote in message
> news:d4cmlr$39e$1@smc.vnet.net...
> [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
>
>
>
>
>

--
DrBob@bigfoot.com

Paul Abbott

2005-04-27, 4:01 am

In article <d4cmlr$39e$1@smc.vnet.net>, zak <chocolatez@gmail.com>
wrote:

> in the site:
> http://teachers.crescentschool.org/...ral/spiral.html
> there is a mathematica code for drawing ULAM'S SPIRAL
> ....
>
> Please i hope someone transform it to a functional programming code.

Ulam's spiral appeared in the second issue of The Mathematica Journal,
1(2): 39 (1990). It uses functional programming and complex numbers.
I've updated that code here:

- Spiral on an Integer Lattice

The following code creates a spiral on an integer lattice:

IntegerSpiral[n_] := {Re[#],Im[#]}& /@
Fold[Join[#1, Last[#1]+I^#2 Range[#2/2]]&, {0}, Range[n]]

For example,

spiral = IntegerSpiral[60];

Indexing the points on the spiral, select those points that have prime
index (Pick is new in 5.1):

primes = Pick[spiral, PrimeQ[Range[Length[spiral]]]];

Draw the spiral and display those with prime index as dots:

spiralplot = ListPlot[spiral, PlotJoined->True, Axes -> False,
AspectRatio->Automatic, Epilog -> {PointSize[0.02], Point /@ primes}];

Some patterns in the distribution of the primes are apparent.

--
Paul Abbott Phone: +61 8 6488 2734
School of Physics, M013 Fax: +61 8 6488 1014
The University of Western Australia (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul@physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul

Paul Abbott

2005-04-27, 4:01 am

Further to my previous message, improved code for Ulam's spiral appeared
in the third issue of The Mathematica Journal, 1(3): 57 (1991). It used

IntegerSpiral[n_] := {Re[#],Im[#]}& /@
Fold[Join[#1, Last[#1]+I^#2 Range[#2/2]]&, {0}, Range[n]]

and explained that

the idea is that the length of the limbs of the spiral follow the
pattern 1,1,2,2,3,3,4,4,5,5,... and that each limb is at a 90 degree
angle to the previous limb (multiplication by I in the complex plane).
This is also a nontrivial use of Fold in the sense that the second
argument to Fold is really used.

Also, instead of using Pick, to select those points that have prime
index from, say

spiral = IntegerSpiral[60];

one can use

primes = spiral[[ Prime[Range[PrimePi[Length[spiral]]]] ]];

Cheers,
Paul

--
Paul Abbott Phone: +61 8 6488 2734
School of Physics, M013 Fax: +61 8 6488 1014
The University of Western Australia (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul@physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul

DrBob

2005-04-27, 4:01 am

Here's an admittedly trivial simplification:

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

"c = 2Mod[n, 2, -1] + 1" seemed a bit tortured, to me.

Bobby

On Tue, 26 Apr 2005 01:33:03 -0400 (EDT), Carl K. Woll <carlw@u.washington.edu> wrote:

> "Carl K. Woll" <carlw@u.washington.edu> wrote in message
> news:d4hvsv$1jp$1@smc.vnet.net...
>
> I see now that I started with 1 at the origin, and others started with 0 at
> the origin. Hence, a small modification of my code will return results that
> agree with e.g., Daniel Lichtblau.
>
> coords[k_List] := Module[{n, a, b, c},
> n = Floor[Sqrt[k + .1]];
> a = k - n(n + 1) ;
> c = 2Mod[n, 2, -1] + 1;
> b = Quotient[2n + 1 - c, 4];
> c Transpose[{Quotient[Abs[a] + a, 2] - b, Quotient[Abs[a] - a, 2] - b}]]
>
> Basically, I substituted k+1 for k. I also got rid of Power by using n(n+1)
> to evaluate n^2+n and I by using Mod to find (-1)^n. These latter changes
> increased the speed of coords by about 33%.
>
> Carl Woll
>
>
>
>
>

--
DrBob@bigfoot.com

Sponsored Links