(*Whatlayer,i.e.whatsquare,isthepointin?*)
In[]:=
quadlayer[q_?(Positive[#]&&IntegerQ[#]&)]:=Ceiling
Sqrt[q]-1
2

(*Howmanypointsaroundthegivenlayeristhepoint?*)
In[]:=
quadaround[q_?(Positive[#]&&IntegerQ[#]&)]:=q-
2
(2quadlayer[q]-1)
In[]:=
quadpoint[q_?(Positive[#]&&IntegerQ[#]&)]:=quadlayer[q]Sqrt[2]AngleVector
π
2
#1-
1
2
+(#2)AngleVector
π
2
(#1+1)&@@Quiet[QuotientRemainder[quadaround[q],2quadlayer[q]]]
In[]:=
ListPlot[quadpoint[#]&/@Range[400],JoinedTrue,AspectRatioAutomatic]
Out[]=
-5
5
10
-10
-5
5
10
In[]:=
ListPlot[quadpoint[#]&/@(Prime[#]&/@Range[10000]),AspectRatioAutomatic,AxesFalse]
Out[]=
In[]:=
AbsoluteTiming[ListPlot[quadpoint[#]&/@(Prime[#]&/@Range[10000]),AspectRatioAutomatic,AxesFalse]]
Out[]=
0.975403,

In[]:=
ListPlot[quadpoint[#-40]&/@(Prime[#]&/@Range[13,10000]),AspectRatioAutomatic,AxesFalse]
Out[]=
(*Wecanalsoextendthisnotiontootherspiralsaroundtheplane*)
In[]:=
mgonLayer[q_Integer,m_Integer]:=Ceiling
1
2
1+
-8+m+8q
m
-1
In[]:=
mgonLayer[#,6]&/@Range[20]
Out[]=
{0,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3}
In[]:=
maround[q_,m_]:=q-1+
m
2
((mgonLayer[q,m]-1)(mgonLayer[q,m]));
In[]:=
mpoint[q_,m_]:=mgonLayer[q,m]AngleVector
2π
m
(#1-1)+#22Sin
π
m
AngleVector2
π
m
(#1-1)+
π
2
+
π
m
&@@Quiet[QuotientRemainder[maround[q,m],mgonLayer[q,m]]]
In[]:=
ListPlot[mpoint[#,6]&/@Range[1000],AxesFalse,AspectRatioAutomatic,JoinedTrue]
Out[]=