根据matovitch最新的代码草案(我只浏览了一眼)，我可能会想得太多了，但无论如何。。。在

设A=(A.x，A.y)，B=(B.x，B.y)，其中(A.x，A.y，B.x，B.y)是整数。

然后直线p，AB的垂直平分线穿过

M=(M.x，M.y)=((A.x+B.x)/2，(A.y+B.y)/2)

AB和p的斜率乘积为-1，因此p的斜率为

-(B.x-A.x)/(B.y-A.y)

因此，在点斜率形式下，p的方程是

(y-M.y)/(x-M.x)=(A.x-B.x)/(B.y-A.y)

重新排列，

y*(B.y-A.y)+x*(B.x-A.x)=M.y*(B.y-A.y)+M.x*(B.x-A.x)

=((B.y+A.y)*(B.y-A.y)+(B.x+A.x)*(B.x-A.x))/2

=(B.y^2-A.y^2+B.x ^2-A.x^2)/2

显然，对于任何晶格点(x，y)，y*(B.y-A.y)+x*(B.x-A.x)必须是整数。所以只有当(B.y^2-A.y^2+B.x^2-A.x^2)是偶数时，直线p才会通过晶格点。在

现在(B.y^2-A.y^2+B.x^2-A.x^2)是偶数的&仅当(A.x+B.x+A.y+B.y)是偶数，如果偶数(A.x，A.y，B.x，B.y)是奇数。在下面，我假设(A.x+B.x+A.y+B.y)是偶数。在

让

dx=(B.x-A.x)

dy=(B.y-A.y)

s=(B.y^2-A.y^2+B.x^2-A.x^2)/2

所以p的方程是

y*dy+x*dx=s

因为y，dy，x，dx&s都是整数，所以这个方程是一个线性丢番图方程，而求这类方程解的标准方法是使用extended Euclidean algorithm。我们的方程只有在dx&dy的最大公约数除以s时才有解。幸运的是，在这种情况下是正确的，但我不在这里给出证明。在

设Y，X是Y*dy+X*dx=g的解，其中g是gcd(dx，dy)，即

Y*dy+X*dx=g

Y*dy/g+X*dx/g=1

设dy'=dy/g，dx'=dx/g，s'=s/g，所以

Y*dy'+X*dx'=1

将最后一个p的方程除以g，我们得到

y*dy'+x*dx'=s'

我们现在可以为它构造一个解决方案。

(Y*s')*dy'+(X*s')*dx'=s'

i、 e.，(X*s'，Y*s')是直线上的点阵点。在

我们可以得到这样的所有解决方案：

(Y*s'+k*dx')*dy'+(X*s'-k*dy')*dx'=s'，对于所有整数k

为了将解限制在(0，0)到(W，H)之间，我们需要解决k的这些不等式：

0&lt；=X*s'-k*dy'&lt；=W和0&lt；=Y*s'+k*dx'&lt；=H

我不会在这里展示这些不等式的解决方案；有关详细信息，请参阅下面的代码。在#! /usr/bin/env python

''' Lattice Line

Find lattice points, i.e, points with integer co-ordinates,

on the line that is the perpendicular bisector of the line segment AB,

where A & B are lattice points.

See http://stackoverflow.com/q/31265139/4014959

Written by PM 2Ring 2015.07.08

Code for Euclid's algorithm & the Diophantine solver written 2010.11.27

'''

from __future__ import division

import sys

from math import floor, ceil

class Point(object):

''' A simple 2D point '''

def __init__(self, x, y):

self.x, self.y = x, y

def __repr__(self):

return "Point(%s, %s)" % (self.x, self.y)

def __str__(self):

return "(%s, %s)" % (self.x, self.y)

def euclid(a, b):

''' Euclid's extended algorithm for the GCD.

Returns a list of tuples of (dividend, quotient, divisor, remainder)

'''

if a < b:

a, b = b, a

k = []

while True:

q, r = a // b, a % b

k.append((a, q, b, r))

if r == 0:

break

a, b = b, r

return k

def dio(aa, bb):

''' Linear Diophantine solver

Returns [x, aa, y, bb, d]: x*aa + y*bb = d

'''

a, b = abs(aa), abs(bb)

swap = a < b

if swap:

a, b = b, a

#Handle trivial cases

if a == b:

eqn = [2, a, -1, a]

elif a % b == 0:

q = a // b

eqn = [1, a, 1-q, b]

else:

#Generate quotients & remainders list

z = euclid(a, b)[::-1]

#Build equation from quotients & remainders

eqn = [0, 0, 1, 0]

for v in z[1:]:

eqn = [eqn[2], v[0], eqn[0] - eqn[2]*v[1], v[2]]

#Rearrange & fix signs, if required

if swap:

eqn = eqn[2:] + eqn[:2]

if aa < 0:

eqn[:2] = [-eqn[0], -eqn[1]]

if bb < 0:

eqn[2:] = [-eqn[2], -eqn[3]]

d = eqn[0]*eqn[1] + eqn[2]*eqn[3]

if d < 0:

eqn[0], eqn[2], d = -eqn[0], -eqn[2], -d

return eqn + [d]

def lattice_line(pA, pB, pC):

''' Find lattice points, i.e, points with integer co-ordinates, on

the line that is the perpendicular bisector of the line segment AB,

Only look for points in the rectangle from (0,0) to C

Let M be the midpoint of AB. Then M = ((A.x + B.x)/2, (A.y + B.y)/2),

and the equation of the perpendicular bisector of AB is

(y - M.y) / (x - M.x) = (A.x - B.x) / (B.y - A.y)

'''

nosolutions = 'No solutions found'

dx = pB.x - pA.x

dy = pB.y - pA.y

#Test parity of co-ords to see if there are solutions

if (dx + dy) % 2 == 1:

print nosolutions

return

#Handle horizontal & vertical lines

if dx == 0:

#AB is vertical, so bisector is horizontal

y = pB.y + pA.y

if dy == 0 or y % 2 == 1:

print nosolutions

return

y //= 2

for x in xrange(pC.x + 1):

print Point(x, y)

return

if dy == 0:

#AB is horizontal, so bisector is vertical

x = pB.x + pA.x

if x % 2 == 1:

print nosolutions

return

x //= 2

for y in xrange(pC.y + 1):

print Point(x, y)

return

#Compute s = ((pB.x + pA.x)*dx + (pB.y + pA.y)*dy) / 2

#s will always be an integer since (dx + dy) is even

#The desired line is y*dy + x*dx = s

s = (pB.x**2 - pA.x**2 + pB.y**2 - pA.y**2) // 2

#Find ex, ey, g: ex * dx + ey * dy = g, where g is the gcd of (dx, dy)

#Note that g also divides s

eqn = dio(dx, dy)

ex, ey, g = eqn[::2]

#Divide the parameters of the equation by the gcd

dx //= g

dy //= g

s //= g

#Find lattice limits

xlo = (ex * s - pC.x) / dy

xhi = ex * s / dy

if dy < 0:

xlo, xhi = xhi, xlo

ylo = -ey * s / dx

yhi = (pC.y - ey * s) / dx

if dx < 0:

ylo, yhi = yhi, ylo

klo = int(ceil(max(xlo, ylo)))

khi = int(floor(min(xhi, yhi)))

print 'Points'

for k in xrange(klo, khi + 1):

x = ex * s - dy * k

y = ey * s + dx * k

assert x*dx + y*dy == s

print Point(x, y)

def main():

if len(sys.argv) != 7:

print ''' Find lattice points, i.e, points with integer co-ordinates,

on the line that is the perpendicular bisector of the line segment AB,

where A & B are lattice points with co-ords (xA, yA) & (xB, yB).

Only print lattice points in the rectangle from (0, 0) to (W, H)

Usage:

%s xA yA xB yB W H''' % sys.argv[0]

exit(1)

coords = [int(s) for s in sys.argv[1:]]

pA = Point(*coords[0:2])

pB = Point(*coords[2:4])

pC = Point(*coords[4:6])

lattice_line(pA, pB, pC)

if __name__ == '__main__':

main()

我没有对这段代码进行过广泛的测试，但它似乎能正常工作。：)