对角线遍历，或按最小和进行遍历

uk702 发表于 2020-11-16 10:39:32

本帖最后由 uk702 于 2020-11-16 15:16 编辑

多维循环时，不知有没有什么好办法按照类似对角进行遍历，比如说三维时，按 i,j,k=(1,1,1) 、(1,1,2)、(1, 2, 1)、(2, 1, 1)、(1, 1, 3)、(1, 3, 1)、(3, 1, 1)、(1, 2, 2)、(2, 1, 2）、(2, 2, 1）、(1, 1, 4）... 这种方式进行遍历，也就是说，按照 i+j+k 的值从小到大进行遍历。

这样的好处是，比如说求解一些不定方程时，可以方便地将 “最小解” 先算出来。

===================================================
好像搞复杂了，以四维循环为例，这样子似就满足要求：
for a=1:n
for x=1:a
for y=1:a-x
for z=1:a-x-y
k = a-x-y-z
...

uk702 发表于 2020-11-16 12:41:43

本帖最后由 uk702 于 2020-11-16 14:12 编辑

经实践，如果维数比较大的话，比如 4 维，预先将对角线（最小和）组成的数组储存起来的方法不可行，太容易 out of memory 了。
不太讲究的话，先将两维的最小和数组排列储存起来，再按 2X2 的方式进行遍历，最后再按 4 个和最小的形式重新整理。

参考代码如下：
a = [(i,j) for i=1:1000 for j=1:1000];
sort(a, by=sum)

dx = 10000.0
for x in a, y in a
global dx
(u, v), (w, k) = x, y
dt = abs((u+v*sqrt(w))/k - pi)
if dt < dx
dx = dt
println((u, v, w, k), " ", dx)
end
end

uk702 发表于 2020-11-16 13:45:40

本帖最后由 uk702 于 2020-11-16 13:53 编辑

julia 语言中，应该可以用 partitions() 函数，写个两重循环。
其它语言可能类似。

using Combinatorics
for n=1:1000
for x in collect(partitions(n, 4))
...

注意：这种输出的 x 是从大到小，并且去重的。

.·.·. 发表于 2020-11-22 21:41:07

uk702 发表于 2020-11-16 13:45
julia 语言中，应该可以用 partitions() 函数，写个两重循环。
其它语言可能类似。

collect(partitions(70,8))大概也会爆炸吧……
一共binomial(77,7)=2404808340种可能
试了试，Rust直接遍历并计算8个数取异或的结果的和，哪怕不开优化也至多需要105.728389385s（每个数字的取值范围是区间的整数，和为70，不假设数字之间的大小关系）
共计算了2404808340个数字，和是5220216392

打开全部优化我敢算到100
const N:i32=100;
fn main(){
let now=std::time::Instant::now();
let mut ans=0u64;
let remain=N;
for a1 in 0..=remain{
   let remain=remain-a1;
   for a2 in 0..=remain{
         let remain=remain-a2;
         for a3 in 0..=remain{
            let remain=remain-a3;
            for a4 in 0..=remain{
               let remain=remain-a4;
               for a5 in 0..=remain{
                     let remain=remain-a5;
                     for a6 in 0..=remain{
                        let remain=remain-a6;
                        for a7 in 0..=remain{
                           //a8=remain-a7;
                           ans+=(((a1^a2)^(a3^a4))^((a5^a6)^(a7^(remain-a7)))).count_ones() as u64;
                        }
                     }
               }
            }
         }
   }
}
println!("{}, cost={:?}",ans,now.elapsed());
}
输出：64995065512, cost=54.381065001s
总感觉调用collect(partitions(100,8))算这个会死

uk702 发表于 2020-11-23 07:22:43

本帖最后由 uk702 于 2020-11-23 07:30 编辑

.·.·. 发表于 2020-11-22 21:41
collect(partitions(70,8))大概也会爆炸吧……
一共binomial(77,7)=2404808340种可能
试了试，Rust直接 ...

虽然我后来也发现使用 partitions() 通常都是无谓的多余，但为了发扬无聊的较真精神，还是写了段程序验证一下。
结论：总的来说不具备可行性。

看了 partitions() 的代码，也无它，其实就是执行几重循环一个个的减，再加上储存操作，当 n 变得稍大，储存开销就变得太恐怖了。

using Combinatorics

function f(n)
a = collect(partitions(n, 8))
ans = 0
for x in a
a1, a2, a3, a4, a5, a6, a7, a8 = x
ans += xor(xor(xor(a1, a2), xor(a3, a4)), xor(xor(a5, a6), xor(a7, a8)))
end
println(ans)
end

@time f(100)
29892426
0.124802 seconds (930.93 k allocations: 134.944 MiB, 53.93% gc time)

@time f(160)
901994896
5.254283 seconds (39.89 M allocations: 5.498 GiB, 59.36% gc time)

@time f(300)
估计至少要10分钟

@time f(1000)
Out of memory

.·.·. 发表于 2020-11-23 15:31:39

本帖最后由 .·.·. 于 2020-11-23 18:26 编辑

uk702 发表于 2020-11-23 07:22
虽然我后来也发现使用 partitions() 通常都是无谓的多余，但为了发扬无聊的较真精神，还是写了段程序验 ...

你的partitions(100,8)假定了a1~a8的大小关系（并且保证了a1~a8大于0）
如果这样的话，Rust其实可以算得更快一些
fn calcN(N:i32){
print!("N={},",N);
let now=std::time::Instant::now();
let mut ans=0u64;
let remain=N;
for a1 in 1..=remain>>3{
   let remain=remain-a1;
   for a2 in a1..=remain/7{
         let remain=remain-a2;
         for a3 in a2..=remain/6{
            let remain=remain-a3;
            for a4 in a3..=remain/5{
               let remain=remain-a4;
               for a5 in a4..=remain>>2{
                     let remain=remain-a5;
                     for a6 in a5..=remain/3{
                        let remain=remain-a6;
                        for a7 in a6..=remain>>1{
                           let a8=remain-a7;
                           ans+=(((a1^a2)^(a3^a4))^((a5^a6)^(a7^a8))) as u64;
                        }
                     }
               }
            }
         }
   }
}
println!("{}, cost={:?}",ans,now.elapsed());
}
fn main(){
calcN(100);
calcN(160);
calcN(300);
calcN(400);
calcN(500);
calcN(600);
calcN(700);
calcN(800);
calcN(900);
calcN(1000);
}
测速：
N=100,29892426, cost=1.616831ms
N=160,901994896, cost=30.246651ms
N=300,109638857854, cost=1.737976853s
N=400,1260273347474, cost=11.636545741s
N=500,6722928203618, cost=51.702947067s
N=600,25125831176186, cost=177.460155807s
N=700,92934978535684, cost=510.614317959s
N=800,297372913223336, cost=1269.206469242s
N=900,725774746541996, cost=2830.599574914s
N=1000,1614149691866148, cost=5760.437703116s

uk702 发表于 2020-11-23 19:13:32

本帖最后由 uk702 于 2020-11-23 19:14 编辑

先佩服一下你的刨根问底的精神。这个帖是我犯的一个错误，搞了这么无聊的一出算法来。

function calcN(n::Int64)
t1 = time()
ans::Int64 = 0

for a1=1:div(n, 8)
   r1 = n - a1
   for a2=a1:div(r1, 7)
         r2 = r1 - a2
         for a3=a2:div(r2, 6)
            r3 = r2 - a3
            for a4=a3:div(r3, 5)
               r4 = r3 - a4
               for a5=a4:div(r4, 4)
                     r5 = r4 - a5
                     for a6=a5:div(r5, 3)
                        r6 = r5 - a6
                        for a7=a6:div(r6, 2)
                           a8=r6 - a7
                           ans += xor(xor(xor(a1, a2), xor(a3, a4)), xor(xor(a5, a6), xor(a7, a8)))
                        end
                     end
               end
            end
         end
   end
end
println("n=$n, ans = $ans, cost = $(time() -t1)")
end

calcN(100);
calcN(160);
calcN(300);
calcN(400);
calcN(500);
calcN(600);
calcN(700);
calcN(800);
calcN(900);
calcN(1000);

$ julia wl.jl
n=100, ans = 29892426, cost = 0.0009999275207519531s
n=160, ans = 901994896, cost = 0.013000011444091797s
n=300, ans = 109638857854, cost = 0.5759999752044678s
n=400, ans = 1260273347474, cost = 3.486999988555908s
n=500, ans = 6722928203618, cost = 14.36300015449524s
n=600, ans = 25125831176186, cost = 46.539000034332275s

页: [1]

数学研发论坛's Archiver

对角线遍历，或按最小和进行遍历