数学研发论坛

 找回密码
 欢迎注册
查看: 131|回复: 12

[讨论] 对角线遍历,或按最小和进行遍历

[复制链接]
发表于 2020-11-16 10:39:32 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
本帖最后由 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
                                ...
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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   
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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(每个数字的取值范围是区间[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))算这个会死
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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

点评

话说julia的collect是什么意思?julia能写成类似python的迭代器的形式吗?  发表于 2020-11-23 18:29
julia直接写循环不用partitions要算多久?总觉得julia不应该这么慢的……慢出境界了  发表于 2020-11-23 15:33
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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其实可以算得更快一些
  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
复制代码

点评

我实在没搞明白,分配5G内存然后GC需要5秒钟这么慢吗?  发表于 2020-11-23 18:31
只是想知道julia如果安装这里的算法需要算多久。  发表于 2020-11-23 18:27
没搞懂你想测试什么,比较 rust 的运行速度吗?实在没这个必要吧。 想搞懂 julia 调 partitions() 为什么“慢出境界”吗?大量地分配内存再大量地 gc。  发表于 2020-11-23 16:45
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2020-11-23 19:13:32 | 显示全部楼层
本帖最后由 uk702 于 2020-11-23 19:14 编辑

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


  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. calcN(100);
  30. calcN(160);
  31. calcN(300);
  32. calcN(400);
  33. calcN(500);
  34. calcN(600);
  35. calcN(700);
  36. calcN(800);
  37. calcN(900);
  38. calcN(1000);   

  39. $ julia wl.jl
  40. n=100, ans = 29892426, cost = 0.0009999275207519531s
  41. n=160, ans = 901994896, cost = 0.013000011444091797s
  42. n=300, ans = 109638857854, cost = 0.5759999752044678s
  43. n=400, ans = 1260273347474, cost = 3.486999988555908s
  44. n=500, ans = 6722928203618, cost = 14.36300015449524s
  45. n=600, ans = 25125831176186, cost = 46.539000034332275s
复制代码

点评

n=600, ans = 25125831176186, cost = 34.18349599838257我真不知道julia这么快……感觉是时候研究一下rust出了什么问题or学一门新语言了  发表于 2020-11-23 19:43
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2020-12-3 21:57 , Processed in 0.103370 second(s), 22 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表