在网上发现的Ulam螺旋的描述中，如何将复数绘制为方向矢量？

Question

有问题的代码可以在here中找到

import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

def make_spiral(R):
    """Generates the spiral up to radius R.
    Returns a list of complex numbers (x, y coordinates)."""
    spiral = [(0+0j)]
    for r in range(R+1):
        corner = r - 1j * r 
        side_len = 2 * r 
        current_pos = corner
        for side, direction in zip(range(4), [1j, -1, -1j, 1]):
            for step in range(side_len):
                current_pos += direction
                spiral.append(current_pos)
    return spiral

fig, axes = plt.subplots(figsize=(15, 5), ncols=3)
for ax, r in zip(axes.ravel(), [2, 3, 4]):
    spiral = make_spiral(r)
    ax.plot(np.real(spiral), np.imag(spiral), '-o')

​

​

我倾向于认为1j意味着沿着y轴向上，根据复数的通常笛卡尔表示，这会使每个右下角的“角”点沿着每个螺旋转弯的右侧向上。-1可能会将螺旋从右向左滑动。依此类推，直到完成螺旋的一圈。

但是我不明白为什么在r = 1的时候，从r - 1j * r开始会产生第一个螺旋。或者应该生成[(0, 1j), (1, -1), (2, (-0-1j)), (3, 1)]的zip(range(4), [1j, -1, -1j, 1])如何编码需要向上、向左、向下和向右的步数。

我试图通过查找for side, direction in zip...来理解Python是如何处理这些命令的，但是我只能得到链接的文档作为我的搜索的唯一返回。我看到下一行for step对于理解这一点至关重要。

Answer 1

当然，如果你了解Python，这个问题是愚蠢的，但这里是伪代码/注释代码，适用于像我这样会说破碎的Python的人：

def make_spiral(R):
    """Generates the spiral up to a max number of steps in the last coil of R.
    Returns a list of complex numbers (x, y coordinates)."""
    spiral = [(0+0j)] # Starts at the (0,0) center.
    for r in range(R+1):
    # Because Py starts counting at 0 the first iteration below (range(side_len)))
    # will be lost (side_len = 2*r = 0). So we need to extend to range(R+1)
        corner = r - 1j * r 
    # This corner line repositions the starting point in each coil one diagonal step 
    #... farther away at the top iteration. This corner, however, is not included 
    #... in the output at this point. See (*) 
        side_len = 2 * r 
    # In the actual spiral we go +1 R, +i up, -2 L, -2i down, +3 R, +3i up, -4 L,...
    # odd, odd, even, even, etc...
    # But starting at the diagonal corner (one below ad to the right)
    #... the # of steps by side will be even and equal, doubling each coil. (**)
        current_pos = corner # Starts at the corner... 
        for side, direction in zip(range(4), [1j, -1, -1j, 1]):
        # "side" is a counter or dummy running for the sides of the coil:
        # 0 = R vert; 1 = UP horiz; 2 = L vert: 3 = DOWN horiz
        # This is not actually used in the calculations.
        # "direction" encodes the motion along C in each step.
        # list(zip(range(4), [1j, -1, -1j, 1])):
        # [(0, 1j), (1, -1), (2, (-0-1j)), (3, 1)]
            for step in range(side_len):
            # In, for example, side_len =2, the range is 0 and 1'
            # but 0 is not the number 0, but just one iteration. So...
                current_pos += direction
                spiral.append(current_pos)
            # ...each direction will be added twice (once for 0, and once for 1):
            # 1j added twice will bring the corner first to (1,0), then to (1,1).
            # -1 twice will bring the position first to (0,1), then to (-1,1).
            # -1j twice will bring the position first to (0,-1), then to (-1,-1).
            # 1 twice will bring the position first to (0,-1), then to (1,-1).
            # Notice that (1,-1) is the corner. This time we save it in the output.
    return spiral

(*)

​

​

(**)

​

例如，make_spiral R =2的输出为：

[0j, (1+0j), (1+1j), 1j, (-1+1j), (-1+0j), (-1-1j), -1j, (1-1j), (2-1j), (2+0j), (2+1j), (2+2j), (1+2j), 2j, (-1+2j), (-2+2j), (-2+1j), (-2+0j), (-2-1j), (-2-2j), (-1-2j), -2j, (1-2j), (2-2j)]

现在很容易将其转换为R(没有进行优化时间的尝试)：

spiral <- function(R){
  spiral <- 0 + 0i
    for (r in 1:R){
      corner = r - 1i * r 
      side_len = 2 * r
      current_pos = corner
        for (dir in c(1i, -1, -1i, 1)){
          for (step in seq(side_len)){
            current_pos <- current_pos + dir
            spiral <- c(spiral, current_pos)
          }
        }
    }
  spiral
} 


get_primes <- function(n_min, n_max){
  options(scipen=999)
  result = vector()
  for (x in seq(max(n_min,2), n_max)){
    has_factor <- F
    for (p in seq(2, ceiling(sqrt(x)))){
      if(x %% p == 0) has_factor <- T
      if(has_factor == T) break
    }
    if(has_factor==F) result <- c(result,x)
  }
  result
}



n = 120
sp = spiral(n)
int = 1: length(sp)
prime = get_primes(1,length(sp))
dat = data.frame(int,sp)

ulnam = dat[dat$int %in% prime,]
par(mar=c(.1,.1,.1,.1))
plot(ulnam$sp, asp=1, pch=19, cex=.01, 
     axes=F, xlab="", ylab="NULL", col="slateblue4")

​

​