代码之家 › 专栏 › 技术社区 › delfuego

如何找到系列中的第一百万个数字:2 3 4 6 9 13 19 28 42 63…?

series recursion math bash

delfuego · 技术社区 · 7 年前

我的成绩达到3000分需要几分钟,但我需要知道这个系列的第一百万个数字。这个定义是递归的,所以我看不到任何快捷方式,除了计算第一百万个数字之前的所有内容。你怎么能快速计算出这个数列中的百万分之一?

系列DEF

n_{i+1} = \floor{ 3/2 * n_{i} } 和 n_{0}=2 .

有趣的是,只有一个网站根据谷歌列出了这个系列: this one .

太慢的bash代码

#!/bin/bash

function series 
{
        n=$( echo "3/2*$n" | bc -l | tr '\n' ' ' | sed -e 's@\\@@g' -e 's@ @@g' );
                                        # bc gives \ at very large numbers, sed-tr for it
        n=$( echo $n/1 | bc )           #DUMMY FLOOR func
}

n=2
nth=1

while [ true ]; #$nth -lt 500 ];
do
        series $n                        # n gets new value in the function through global value
        echo $nth $n
        nth=$( echo $nth + 1 | bc )     #n++
done

12 回复 | 直到 13 年前

therealhoff 16 年前

您可以通过考虑二进制的问题来轻松地解决这个问题。地板(3/2*i)基本上是右移、截断和添加。在伪代码中:

0b_n[0]   = 10              // the start is 2
0b_n[1]   = 10 + 1(0) = 11  // shift right, chop one bit off and add 
0b_n[i+1] = 0b_n[i] + Truncate(ShiftRight(0b_n[i]))

这对于任何形式的实现都应该是相当快的。

我刚刚在Mathematica中实现了这一点,看起来bitshiftright操作也会将位截断到单元位置之外,这样就可以自动进行处理。这是一条线:

In[1] := Timing[num = Nest[(BitShiftRight[#] + #) &, 2, 999999];]
Out[2] = {16.6022, Null}

16秒,数字打印得很好,但很长:

In[2] := IntegerLength[num]
Out[2] = 176092

In[3] := num
Out[3] = 1963756...123087

全部结果 here .

Xavier Ho 9 年前

你差点就找到了。下次,看看 Online Encyclopedia of Integer Series .

条目如下: http://oeis.org/A061418

     FORMULA    

a(n) = ceiling[K*(3/2)^n] where K=1.08151366859...

The constant K is 2/3*K(3) (see A083286). - Ralf Stephan, May 29, 2003

上面说:

>>> def f(n):
...     K = 1.08151366859
...     return int(ceil(K*(1.5)**n))

清醒测试:

>>> for x in range(1, 10):
...     print f(x)
...     
2
3
4
6
9
13
19
28
42

令人惊叹的!现在100万怎么样:

>>> f(1000000)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "<input>", line 3, in f
OverflowError: (34, 'Result too large')

嗯,我试过了。 :] 但你明白了。

再次编辑: 找到解决方案!见 Timo 或 Lasse V. Karlsen 的答案。

编辑: 使用蒂莫的一点转变想法:

import gmpy
n=gmpy.mpz(2)
for _ in xrange(10**6-1):
    n+=n>>1
print(n)

产量

1963756763…226123087(176092位)

% time test.py > test.out

real    0m21.735s
user    0m21.709s
sys 0m0.024s

Jurney 16 年前

你的脚本这么慢的原因是它正在生成 bc 三次, tr 一次 sed 一次循环。

重写整个内容 公元前 然后做 塞德 最后。我的测试显示 公元前 -唯一的版本是600倍以上的速度。在一个旧的慢系统上花了不到16分钟的时间 公元前 找到第100000个值的版本(仅打印最后一个)。

另外,请注意,您的“floor”函数实际上是“int”。

#!/usr/bin/bc -lq
define int(n) {
    auto oscale
    oscale = scale
    scale = 0
    n = n/1
    scale = oscale
    return n
}

define series(n) {
    return int(3/2*n)
}

n = 2
end = 1000
for (count = 1; count < end; count++ ) {
    n = series(n)
}
print count, "\t", n, "\n"
quit

注意 print 是的扩展和某些版本 公元前 可能没有。如果是这样,只需引用变量本身及其值就可以输出。

现在你可以做 chmod +x series.bc 这样称呼它:

./series.bc | tr -d '\n' | sed 's.\\..g'

polygenelubricants 16 年前

我使用了下面的Java程序:

import java.math.*;

public class Series {
    public static void main(String[] args) {
        BigInteger num = BigInteger.valueOf(2);
        final int N = 1000000;

        long t = System.currentTimeMillis();
        for (int i = 1; i < N; i++) {
            num = num.shiftLeft(1).add(num).shiftRight(1);
        }
        System.out.println(System.currentTimeMillis() - t);
        System.out.println(num);
    }
}

裁剪后的输出:( full output on pastebin )

516380 (milliseconds)
196375676351034182442....29226123087

所以在我的普通机器上花了8.5分钟。我用过 -Xmx128M 但不确定是否真的有必要。

可能有更好的算法,但这总共需要10分钟,包括编写幼稚的实现和运行程序。

样本运行

N=500 (同意) OEIS 061418 )
N=100K ( 1445232038814....522773508 )

Shawn 9 年前

这是一个在我10岁的笔记本电脑上运行的python版本,运行大约需要220秒:

import math;
import timeit;

def code():
  n = 2
  nth = 1

  while nth < 1000000:
    n = (n * 3) >> 1
    nth = nth + 1

  print(n);

t = timeit.Timer(setup="from __main__ import code", stmt="code()");
print(t.timeit(1));

它产生的结果与 this answer 在Pastebin上(也就是说,我验证了它的开始和结束,而不是全部。)

paxdiablo 16 年前

隐马尔可夫模型, bash 不是我用来高速数字处理的。给自己一份 GMP 然后编一个C程序来完成它。

很可能有一个数学公式可以很快地给出给你,但是,在你计算出来的时候,GMP可能会把结果扔给你:—)

Alex Martelli 16 年前

这被确定为序列 A061418 在 sequences 站点(又称“整数序列在线百科全书”);每 the relevant page ,

公式 a(n) =A061419(n)+1 = ceiling[K*(3/2)^n] 在哪里? K=1.08151366859 …常数k是 2/3*K(3) (见A08328)。

并且有一个合适的高精度的库(如已经建议的那样是gmp,或者是mpir,也许还有一个像我的宝宝一样的包装纸 gmpy 对于python)您可以使用封闭形式的公式来更快地计算“系列中的第百万项”等。

通常可以将递归指定的递归放入闭合公式中。对于一个广泛的初学者对这个主题的介绍, Concrete Mathematics (格雷厄姆,克努斯和帕塔什尼克)真的很难打败。

Jerry Coffin 16 年前

通过使用更合适的语言,例如,Scheme,您可能会更近一点。

(define (series n) (if (= n 0) 2
                       (quotient (* 3 (series (- n 1))) 2)))

这将计算 (series 100000) 在我的机器上大约8秒钟。不幸的是, (series 1000000) 即使是一台更新/更快的机器也需要很长时间才能在一分钟内完成。

用大整数库(NTL,在这种情况下)切换到C++:

#include <NTL/zz.h>
#include <fstream>
#include <time.h>
#include <iostream>

int main(int argc, char **argv) {
    NTL::ZZ sum;

    clock_t start = clock();
    for (int i=0; i<1000000; i++) 
        sum = (sum*3)/2;
    clock_t finish = clock();
    std::cout << "computation took: " << double(finish-start)/CLOCKS_PER_SEC << " seconds.\n";

    std::ofstream out("series_out.txt");
    out << sum;

    return 0;
}

这将在我的机器上4分钟35秒内计算1000000个序列。速度足够快几乎相信一台真正快速的新机器至少能在一分钟内完成(是的,我检查了使用轮班而不是乘法/除法时发生的事情——它比较慢)。

不幸的是,其他人提出的封闭式计算似乎没有什么帮助。要使用它,您需要计算常数k到足够的精度。我看不到k的闭式计算,所以这实际上只是将迭代转换为计算k,并且看起来计算k到足够的精度比计算原始系列快一点(如果有的话)。

Charles 16 年前

进去很容易 Pari :

n=2;for(k=1,10^6,n+=n>>1)

这在我的机器上需要14秒。当然,有更快的方法——gmp出现在脑海中——但是为什么要麻烦呢?您将无法从运行时中节省超过10秒的时间,开发时间将按分钟 . :)

小题:原公式中不明确的是第一百万项n _{九十九万九千九百九十九} 期望或n _一百万 ,索引为一百万的数字;我给出了后者,因为我看到前者已经在上面计算过了。

Joe 16 年前

这几乎是一阶递推关系,除了地板,它把事情搞得一团糟。如果你不想要地板,

http://en.wikipedia.org/wiki/Recurrence_relation

另外,不要使用bash。

Billy ONeal 16 年前

在大多数情况下,递归公式需要相当长的时间 Curcumstances,因为它必须维护机器堆栈。为什么不使用动态改为编程?

即(伪代码)

bignum x = 2
for (int i = 1; i < 1000000; i++) {
    x = floor(3.0/2*x)
}

当然,为了获得有意义的结果,您需要一个高精度的数字库。

delfuego 16 年前

我把蒂莫的想法转化为伊利普。它以100失败,给出负数。失败,请看 no BigNums !

(progn
  (let ((a 2)
        (dummy 0))
    (while (< dummy 100)
      (setq a (+ a (lsh a -1)))
      (setq dummy (1+ dummy)))
    (message "%d" a)))
-211190189 #WRONG evalution