在椭圆内寻找整数点的最有效算法

Question

在椭圆内寻找整数点的最有效算法

3

我正在尝试在各种3D椭圆内找到所有整数点。

我希望我的程序可以接受一个整数N，并计算形如ax^2 + by^2 + cz^2 = n的椭圆内的所有格点数量，其中a、b、c是固定整数，n介于1和N之间。该程序应返回N个元组，形式为（n，numlatticePointsWithinEllipse n）。

我目前是通过计算椭圆ax^2 + by^2 + cz^2 = m上的点数，其中m介于0和n之间，并对m求和来实现的。我最初只考虑x、y和z都为正数的情况，然后通过稍后排列它们的符号来添加负数。

理想情况下，我希望在几个小时内达到N = 1,000,000+的数字。

以x ^ 2 + y ^ 2 + 3z ^ 2 = N为具体例子，这是我目前使用的Haskell代码：

import System.Environment

isqrt :: Int -> Int
isqrt 0 = 0
isqrt 1 = 1
isqrt n = head $ dropWhile (\x -> x*x > n) $ iterate (\x -> (x + n `div` x) `div` 2) (n `div` 2)

latticePointsWithoutNegatives :: Int -> [[Int]]
latticePointsWithoutNegatives 0 = [[0,0,0]]
latticePointsWithoutNegatives n = [[x,y,z] | x<-[0.. isqrt n], y<- [0.. isqrt (n - x^2)], z<-[max 0 (isqrt ((n-x^2 -y^2) `div` 3))], x^2 +y^2 + z^2 ==n]

latticePoints :: Int -> [[Int]]
latticePoints n = [ zipWith (*) [x1,x2,x3] y |  [x1,x2,x3] <- (latticePointsWithoutNegatives n), y <- [[a,b,c] | a <- (if x1 == 0 then [0] else [-1,1]), b<-(if x2 == 0 then [0] else [-1,1]), c<-(if x3 == 0 then [0] else [-1,1])]]

latticePointsUpTo :: Int -> Int
latticePointsUpTo n = sum [length (latticePoints x) | x<-[0..n]]

listResults :: Int -> [(Int, Int)]
listResults n = [(x, latticePointsUpTo x) | x<- [1..n]]

main = do
    args <- getArgs
    let cleanArgs = read (head args)
    print (listResults cleanArgs)

我已经使用编译器进行了编译。

ghc -O2 latticePointsTest

但是使用PowerShell的"Measure-Command"命令，我得到了以下结果：

Measure-Command{./latticePointsTest 10}
TotalMilliseconds : 12.0901

Measure-Command{./latticePointsTest 100}
TotalMilliseconds : 12.0901

 Measure-Command{./latticePointsTest 1000}
TotalMilliseconds : 31120.4503

如果我们再向上推一个数量级，就会进入几天的时间尺度，而不是几小时或几分钟。

我的算法有什么根本性问题吗？是什么核心原因导致我的代码不适用于大规模数据？任何指导都将不胜感激。我可能还想在“latticePoints”和“latticePointsUpTo”之间处理数据，因此不能完全依赖巧妙的数论计数技巧-我需要保留底层元组。

- MadMonty

我还没有完全阅读问题，但为什么不直接将isqrt定义为floor . sqrt . fromIntegral呢？ - jub0bs

考虑使用 arithmoi 包中的工具。 - recursion.ninja

肯定有一些通过记忆化进行动态编程的空间。如果存在 ax^2 + by^2 + cz^2 = n/q，那么你知道 qax^2 + qby^2 + qcz^2 = n 存在。这种多重记忆技术将适用于集合 [2..sqrt n] 中的所有数字 q。高效的优化可能会退化为应用数论定理。这可能会在 math.SE 上得到更好的接受。 - recursion.ninja

1

@Jubobs 因为对于一个平方数，sqrt 可能会给出略小于整数的结果，然后 floor 会将其向下舍入到下面的数字。当然，你可以加入一个保护来检查它是否这样做，但我更愿意将所有内容都保持为整数。 - MadMonty

@recursion.ninja 这是一个不错的技巧，它肯定会产生很多答案，但可能存在解a、b、c，使得a^2 + b^2 + c^2 = n（n为合数），其中a、b、c没有公因数。有些方程在这个意义上是可乘的，所以你只需要担心n是质数幂的情况，但我不确定对于一般的椭圆是否显而易见。 - MadMonty

@user1440894，请看我在答案结尾处对 floor 和 sqrt 的评论。 - ErikR

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- ErikR · Accepted Answer

以下是我建议尝试的一些方法：

isqrt 对于您所处理的数值范围来说不够高效。可以直接使用浮点数的 sqrt 函数：

isqrt = floor $ sqrt ((fromIntegral n) :: Double)

或者，您可以在列表综合中使用类似于以下逻辑而不是计算整数平方根的方法：

x <- takeWhile (\x -> x*x <= n) [0..],
y <- takeWhile (\y -> y*y <= n - x*x) [0..]

此外，我建议使用表达式x*x代替x^2。

最后，为什么不用以下方式计算解的数量：

sols a b c n =
  length [ () | x <- takeWhile (\x -> a*x*x <= n) [0..]
              , y <- takeWhile (\y -> a*x*x+b*y*y <= n) [0..]
              , z <- takeWhile (\z -> a*x*x+b*y*y+c*z*z <= n) [0..]
         ]

这并不完全计算出你想要的答案，因为它没有考虑正负解，但你可以很容易地修改它来计算你的答案。思路是使用一个列表推导式，而不是迭代各种值 n 并求和。

最后，我认为在这种情况下使用 floor 和 sqrt 计算整数平方根是完全安全的。此代码通过使用 sqrt 查找 (x*x) == x 对于所有 x <= 3037000499 的整数平方根：

testAll :: Int -> IO ()
testAll n =
  print $ head [ (x,a) | x <- [n,n-1 .. 1], let a = floor $ sqrt (fromIntegral (x*x) :: Double), a /= x ]

main = testAll 3037000499

注意，我正在64位的GHC上运行此代码 - 否则只需使用Int64而不是Int，因为在任何情况下Double都是64位的。只需要一分钟左右即可验证。

这表明，如果y<=3037000499^2，则取sqrt y的floor永远不会得出错误的答案。