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

uk702

多维循环时，不知有没有什么好办法按照类似对角进行遍历，比如说三维时，按 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

经实践，如果维数比较大的话，比如 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

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

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


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

.·.·.

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

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

打开全部优化我敢算到100

1. const N:i32=100;
   
2. fn main(){
   
3.     let now=std::time::Instant::now();
   
4.     let mut ans=0u64;
   
5.     let remain=N;
   
6.     for a1 in 0..=remain{
   
7.         let remain=remain-a1;
   
8.         for a2 in 0..=remain{
   
9.             let remain=remain-a2;
   
10.             for a3 in 0..=remain{
    
11.                 let remain=remain-a3;
    
12.                 for a4 in 0..=remain{
    
13.                     let remain=remain-a4;
    
14.                     for a5 in 0..=remain{
    
15.                         let remain=remain-a5;
    
16.                         for a6 in 0..=remain{
    
17.                             let remain=remain-a6;
    
18.                             for a7 in 0..=remain{
    
19.                                 //a8=remain-a7;
    
20.                                 ans+=(((a1^a2)^(a3^a4))^((a5^a6)^(a7^(remain-a7)))).count_ones() as u64;
    
21.                             }
    
22.                         }
    
23.                     }
    
24.                 }
    
25.             }
    
26.         }
    
27.     }
    
28.     println!("{}, cost={:?}",ans,now.elapsed());
    
29. }
    

复制代码

输出：

1. 64995065512, cost=54.381065001s

复制代码

总感觉调用collect(partitions(100,8))算这个会死

uk702

.·.·. 发表于 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

.·.·.

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

你的partitions(100,8)假定了a1~a8的大小关系（并且保证了a1~a8大于0）
如果这样的话，Rust其实可以算得更快一些

1. fn calcN(N:i32){
   
2.     print!("N={},",N);
   
3.     let now=std::time::Instant::now();
   
4.     let mut ans=0u64;
   
5.     let remain=N;
   
6.     for a1 in 1..=remain>>3{
   
7.         let remain=remain-a1;
   
8.         for a2 in a1..=remain/7{
   
9.             let remain=remain-a2;
   
10.             for a3 in a2..=remain/6{
    
11.                 let remain=remain-a3;
    
12.                 for a4 in a3..=remain/5{
    
13.                     let remain=remain-a4;
    
14.                     for a5 in a4..=remain>>2{
    
15.                         let remain=remain-a5;
    
16.                         for a6 in a5..=remain/3{
    
17.                             let remain=remain-a6;
    
18.                             for a7 in a6..=remain>>1{
    
19.                                 let a8=remain-a7;
    
20.                                 ans+=(((a1^a2)^(a3^a4))^((a5^a6)^(a7^a8))) as u64;
    
21.                             }
    
22.                         }
    
23.                     }
    
24.                 }
    
25.             }
    
26.         }
    
27.     }
    
28.     println!("{}, cost={:?}",ans,now.elapsed());
    
29. }
    
30. fn main(){
    
31.     calcN(100);
    
32.     calcN(160);
    
33.     calcN(300);
    
34.     calcN(400);
    
35.     calcN(500);
    
36.     calcN(600);
    
37.     calcN(700);
    
38.     calcN(800);
    
39.     calcN(900);
    
40.     calcN(1000);
    
41. }

复制代码

测速：

1. N=100,29892426, cost=1.616831ms
   
2. N=160,901994896, cost=30.246651ms
   
3. N=300,109638857854, cost=1.737976853s
   
4. N=400,1260273347474, cost=11.636545741s
   
5. N=500,6722928203618, cost=51.702947067s
   
6. N=600,25125831176186, cost=177.460155807s
   
7. N=700,92934978535684, cost=510.614317959s
   
8. N=800,297372913223336, cost=1269.206469242s
   
9. N=900,725774746541996, cost=2830.599574914s
   
10. N=1000,1614149691866148, cost=5760.437703116s

复制代码

uk702

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

1. function calcN(n::Int64)
   
2.     t1 = time()
   
3.     ans::Int64 = 0
   
4.    
   
5.     for a1=1:div(n, 8)
   
6.         r1 = n - a1
   
7.         for a2=a1:div(r1, 7)
   
8.             r2 = r1 - a2
   
9.             for a3=a2:div(r2, 6)
   
10.                 r3 = r2 - a3
    
11.                 for a4=a3:div(r3, 5)
    
12.                     r4 = r3 - a4
    
13.                     for a5=a4:div(r4, 4)
    
14.                         r5 = r4 - a5
    
15.                         for a6=a5:div(r5, 3)
    
16.                             r6 = r5 - a6
    
17.                             for a7=a6:div(r6, 2)
    
18.                                 a8=r6 - a7
    
19.                                 ans += xor(xor(xor(a1, a2), xor(a3, a4)), xor(xor(a5, a6), xor(a7, a8)))
    
20.                             end
    
21.                         end
    
22.                     end
    
23.                 end
    
24.             end
    
25.         end
    
26.     end
    
27.     println("n=$n, ans = $ans, cost = $(time() -t1)")
    
28. end
    
29. 
    
30. calcN(100);
    
31. calcN(160);
    
32. calcN(300);
    
33. calcN(400);
    
34. calcN(500);
    
35. calcN(600);
    
36. calcN(700);
    
37. calcN(800);
    
38. calcN(900);
    
39. calcN(1000);   
    
40. 
    
41. $ julia wl.jl
    
42. n=100, ans = 29892426, cost = 0.0009999275207519531s
    
43. n=160, ans = 901994896, cost = 0.013000011444091797s
    
44. n=300, ans = 109638857854, cost = 0.5759999752044678s
    
45. n=400, ans = 1260273347474, cost = 3.486999988555908s
    
46. n=500, ans = 6722928203618, cost = 14.36300015449524s
    
47. n=600, ans = 25125831176186, cost = 46.539000034332275s

复制代码