Alright, let's tackle this problem step by step. We have a 2D matrix A shaped (100,50), and for each starting point, we need to traverse along the vector (d, 1) (where d = (1-X)/Y from your code) until we hit a value ≥1.5 in A, then record the number of grid steps taken along the y-dimension. If we never reach the threshold, we mark it as nan.

First, let's align on the setup to avoid confusion:

• We start with the function A(x,y) = x*y, discretized into a grid using:

  import numpy as np
  x = np.linspace(0.01, 0.99, 100)
  y = np.linspace(0.5, 2, 50)
  X, Y = np.meshgrid(x, y, indexing='ij')
  A = X * Y
  d = (1 - X) / Y  # Simplified from Y*(1-X)/(Y*Y)
  

• For each grid point (i,j) (x-index i, y-index j), we move along a path where each step increases the y-index by 1, and the x-index by d[i,j] (since dy/dz=1 and dx/dz=d). We need to find the smallest number of steps k such that:
  1. The new x-index i + k*d[i,j] stays within [0,99]
  2. The new y-index j + k stays within [0,49]
  3. A[round(i + k*d[i,j]), j + k] ≥ 1.5
• If no valid k exists, set R[i,j] = nan.

1. 高效优先方案（向量化Numpy运算） #

This approach avoids slow Python loops by leveraging Numpy's vectorization, which is critical for handling the (100,50) grid efficiently. We precompute all possible step counts and use boolean masking to find the minimal valid k for each point.

代码实现 #

import numpy as np

# Initialize grid and parameters
x = np.linspace(0.01, 0.99, 100)
y = np.linspace(0.5, 2, 50)
X, Y = np.meshgrid(x, y, indexing='ij')
A = X * Y
d = (1 - X) / Y

# Get grid indices for x and y dimensions
X_idx, Y_idx = np.meshgrid(np.arange(100), np.arange(50), indexing='ij')

# Calculate maximum possible steps for each point (can't exceed y-grid bounds)
max_k = 50 - Y_idx - 1  # j + k < 50 → k ≤ 49 - j
max_possible_k = max_k.max()

# Create a step array broadcastable to (100,50, max_possible_k)
k_array = np.arange(1, max_possible_k + 1)[None, None, :]

# Compute new x/y indices for every possible step count
x_new = X_idx[..., None] + k_array * d[..., None]
y_new = Y_idx[..., None] + k_array

# Mask valid index positions (x within 0-99, y within 0-49)
valid_mask = (x_new >= 0) & (x_new <= 99) & (y_new <= 49)

# Round x indices to nearest integer and clamp to valid range
x_new_idx = np.clip(np.round(x_new).astype(int), 0, 99)

# Get A values for all candidate steps
A_vals = A[x_new_idx, y_new]

# Mask where A meets the threshold AND indices are valid
threshold_mask = (A_vals >= 1.5) & valid_mask

# Find the first valid step for each point
first_k_idx = np.argmax(threshold_mask, axis=-1)
# Convert to actual k value; set to nan if no valid step exists
R_efficient = np.where(np.any(threshold_mask, axis=-1), k_array[0,0,first_k_idx], np.nan)

优势 #

• Blazing fast: Fully vectorized operations eliminate Python loops, making this orders of magnitude faster for large grids.
• Clean maintainability: Uses standard Numpy functions, so the code is readable and easy to scale.

2. 精度优先方案（精确根求解+插值） #

This approach prioritizes accuracy by calculating the exact continuous point where A(x(t), y(t)) = 1.5 (instead of relying on discrete grid steps), then converts that to the number of y-grid steps taken. Since A(x,y)=x*y is a simple quadratic function along our path, we can solve it analytically.

实现思路 #

For each starting point (x0, y0) with direction d=(1-x0)/y0, our path is:

• x(t) = x0 + d*t
• y(t) = y0 + t
• We solve (x0 + d*t)(y0 + t) = 1.5 for the smallest positive t.

This expands to a quadratic equation: d*t² + (x0 + d*y0)*t + (x0*y0 - 1.5) = 0. We use the quadratic formula to find valid t, then check if the solution stays within grid bounds. Finally, we convert t to y-grid steps using the y-dimension's step size.

代码实现 #

import numpy as np

# Initialize grid and parameters (same as before)
x = np.linspace(0.01, 0.99, 100)
y = np.linspace(0.5, 2, 50)
X, Y = np.meshgrid(x, y, indexing='ij')
A = X * Y
d = (1 - X) / Y

# Calculate the step size of the y grid
y_step = (y[-1] - y[0]) / (len(y)-1)

# Coefficients for quadratic equation a*t² + b*t + c = 0
a = d
b = X + d * Y
c = A - 1.5

# Compute discriminant to check for real solutions
discriminant = b**2 - 4*a*c

# Mask valid solutions: discriminant ≥0, positive t, and final x/y within grid bounds
valid_disc = discriminant >= 0
# Use quadratic formula to get the positive root (since t must be >0)
t = (-b + np.sqrt(discriminant)) / (2*a)
x_final = X + d*t
y_final = Y + t
valid_solution = valid_disc & (t > 0) & (x_final >= x[0]) & (x_final <= x[-1]) & (y_final >= y[0]) & (y_final <= y[-1])

# Convert t to number of y-grid steps; set to nan if invalid
R_precise = np.where(valid_solution, t / y_step, np.nan)

# Optional: Round up to integer steps to match the example's integer R values
# R_precise_int = np.where(valid_solution, np.ceil(t / y_step), np.nan)

优势 #

• Maximum accuracy: Doesn't rely on discrete grid sampling—calculates the exact point where the threshold is crossed.
• Flexibility: Can return fractional steps (for precise distance) or rounded integer steps to align with the problem's example.

Key notes:

• For the efficient approach, rounding x_new to the nearest integer is a speed-for-precision trade-off. If you need better accuracy here, you could use linear interpolation between adjacent x grid points instead of rounding.
• Both solutions handle out-of-bounds cases by setting invalid entries to nan, which matches the problem's requirements.

内容的提问来源于stack exchange，提问作者FooBar