找回密码
 欢迎注册
查看: 19958|回复: 35

[擂台] projectEuler 454 -- Diophantine reciprocals

[复制链接]
发表于 2014-3-20 22:10:14 | 显示全部楼层 |阅读模式

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

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

×
设对于给定的 \(L\),等式 \(\D\frac1X + \frac1Y = \frac1N, (0 \lt X \lt Y \leqslant L)\) 成立总共有 \(f(L)\) 种情况 .
举几个简单的例子, 有 \(f(6) =1, f(12)=3, f(1000)=1069\), 求\(f(10^{12})\)

参考链接:  http://projecteuler.net/problem=454
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-3-20 23:23:54 | 显示全部楼层
\[\frac1x+\frac1y=\frac1n \iff x+y\,|\, xy \iff \left. (x,y)\,\middle|\,\frac{x+y}{(x,y)} \right.\]
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-3-21 09:09:02 | 显示全部楼层
不错的条件,   能否继续前进,导出计数公式
ps,好像写反了吧:\[\left. \frac{x+y}{(x,y)}\,\middle|\,(x,y) \right.\]
不过,在此基础上,可以导出这种一般性的解表达:   \(x=km(m+n) ,y=kn(m+n), (m,n)=1, m<n\), 即\(L=kn(m+n)\)

但要计算出$f(10^12)$,还需要进一步的挖掘,优化
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-3-21 12:01:14 | 显示全部楼层
直接计算\[ \sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \left[ \frac{L}{n(n+m)} \right] \]只需要`\mathcal{O}(L)`的时间,我认为这是可以承受的。不过我们也可试着优化一下:
\[\begin{split}
\sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \left[ \frac{L}{n(n+m)} \right]
&=\sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} s=n+1 \\ (n,s)=1 \end{subarray}}^{2n-1} \left[ \frac{L}{ns} \right] \\
&=\sum_{n=1}^{\sqrt{L}} \sum_{s=n+1}^{2n-1} \left[ \frac{L}{ns} \right] \left[ \frac{1}{(n,s)} \right] \\
&=\sum_{n=1}^{\sqrt{L}} \sum_{s=n+1}^{2n-1} \left[ \frac{L}{ns} \right] \sum_{d | (n,s)} \mu(d) \\
&=\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \sum_{\begin{subarray}{cc} n < s \leq 2n-1 \\ d | s \end{subarray}} \left[ \frac{L}{ns} \right] \mu(d) \\
&=\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mu(d) \sum_{\frac{n}{d} < s' \leq \frac{2n-1}{d}} \left[ \frac{L}{nds'} \right] \\
&=\sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mu(d) \left( \psi\left( \left[ \frac{L}{nd} \right],\frac{2n-1}{d}\right) - \psi\left(\left[ \frac{L}{nd} \right],\frac{n}{d}\right) \right)
\end{split}
\]
计算`\displaystyle \psi(x,y) = \sum_{n=1}^{y} \left[ \frac{x}{n} \right]`显然是`\mathcal{O}(\sqrt{x})`的。
于是使用该公式可以做到在`\mathcal{O}(L^{3/4})`的复杂度内计算。

点评

$\frac{2n-1}{d}$是不是得加个取整的方括号[]?  发表于 2014-6-26 09:03
强悍!  发表于 2014-3-21 12:31

评分

参与人数 1威望 +12 金币 +12 贡献 +12 经验 +12 鲜花 +12 收起 理由
wayne + 12 + 12 + 12 + 12 + 12 很给力!

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-3-21 13:34:02 | 显示全部楼层
写了个Mathematica函数, 计算10^7花了40秒钟,再算下去不太现实呢,或许C/C++语言实现要快一些

  1. L=10^4;phi=Sum[Floor[#1/n],{n,1,#2}]&;
  2. Sum[MoebiusMu[j](phi[Floor[L/(j i)],(2i-1)/j]-phi[Floor[L/(j i)],i/j]),{i,1,Sqrt[L]},{j,Divisors[i]}]
复制代码
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-3-21 19:18:47 | 显示全部楼层
wayne 发表于 2014-3-21 13:34
写了个Mathematica函数, 计算10^7花了40秒钟,再算下去不太现实呢,或许C/C++语言实现要快一些


经测试,Mathematica计算`\psi(x,y)`时所消耗的时间是平凡的,为`\mathcal{O}(y)`,这会导致复杂度退化为
\[ \sum_{n=1}^{\sqrt{L}} \sum_{d | n} \mathcal{O} \left( \frac{n}{d} \right) = \mathcal{O}(L) \]
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2014-3-21 19:42:04 | 显示全部楼层
Lwins_G 发表于 2014-3-21 19:18
经测试,Mathematica计算`\psi(x,y)`时所消耗的时间是平凡的,为`\mathcal{O}(y)`,这会导致复杂度退化 ...


可以将两个相减的合成一个, 但这也只是常数级的优化,时间砍了一半:
  1. L=10^6;Sum[MoebiusMu[j]Sum[Floor[L/(j i n)],{n,i/j,(2i-1)/j}],{i,1,Sqrt[L]},{j,Divisors[i]}]-L
复制代码

点评

但每一个值出现多少次怎么求出来  发表于 2014-3-21 19:58
`\left[ \frac{N}{n} \right]`这个函数的取值本身就只有`\mathcal{O}(\sqrt{N})`个。  发表于 2014-3-21 19:56
你说的平方根复杂度是怎么做的  发表于 2014-3-21 19:49
所以需要自己实现`\psi(x,y)`。  发表于 2014-3-21 19:44
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-4-17 14:01:16 | 显示全部楼层
Lwins_G 发表于 2014-3-21 12:01
直接计算\[ \sum_{n=1}^{\sqrt{L}} \sum_{\begin{subarray}{cc} m=1 \\ (n,m)=1 \end{subarray}}^{n-1} \le ...


(x,y)中x是10^12的级别的。

点评

所以?  发表于 2014-4-19 10:45

评分

参与人数 1金币 +20 收起 理由
gxqcn + 20 首帖奖励,欢迎常来。

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-4-17 17:46:32 | 显示全部楼层
1/a+1/b=1/c
等价于
(c-a)*(c-b)=c^2
这样问题是不是简单了???????????
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2014-4-17 17:53:25 | 显示全部楼层
(c-a)*(c-b)=c^2写成(a-c)*(b-c)=c^2这种形式更好
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-3-29 06:41 , Processed in 0.073705 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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