很明顯,你沒有卡在第23號,但它將需要大約120小時來評估它。請參閱http://oeis.org/A005179
這是一個APL NARS2000工作區,可以加快速度。 第一個100只需要1/10秒。
)WSID
C:\Users\Ab\AppData\Roaming\NARS2000\workspaces\combfactor
)FNS
combfactor factors listws save
)VARS
COMBFACTOR DESCRIBE DIV RESULTSfor100
∇ z ← combfactor n;f;p
[1] f ← (⊂,n),factors n ⋄ ⍝ find all factors combinations of n
[2] p ← ⍸0π⍳¯2π32 ⋄ ⍝ The first 32 primes numbers
[3] z ← ((⍴¨f)↑¨⊂⍟p)+.רf-1 ⋄ ⍝ give ratios of all combinations
[4] z ← ∊¯1+((⌊/z)=z)/f ⋄ ⍝ get the combination with minimum ratio
[5] z ← (0x+(⍴z)↑p)×.*z ⋄ ⍝ evaluate p_1^(f_1-1) * p_2^(f_2-1) * ... * p_n^(f_n-1)
∇
∇ z ← {fmax} factors n;f;d
[1] :if 0 = ⎕NC 'fmax' ⋄ fmax ← 2*63 ⋄ :endif
[2] z ← ⍬ ⋄ f ← ⌊fmax⌊n÷2
[3] :while f ≥ 2
[4] :if 0 = f∣n
[5] d ← n÷f ⋄ :if d ∧.≤ f,fmax ⋄ z ,← ⊂f,d ⋄ :endif
[6] z ,← f,¨f factors d
[7] :endif ⋄ f -← 1
[8] :endwhile
∇
COMBFACTOR
RESULTSfor100 ← ⍪(,¨⍳100),¨combfactor¨,¨⍳100 ⋄ ⍝ 0.1 second
DESCRIBE
DESCRIBE
⍝ def factors(number, max_factor=sys.maxint):
⍝ result = []
⍝
⍝ factor = min(number/2, max_factor)
⍝ while factor >= 2:
⍝ if number % factor == 0:
⍝ divisor = number/factor
⍝
⍝ if divisor <= factor and divisor <= max_factor:
⍝ result.append([factor, divisor])
⍝
⍝ result.extend([factor] + item for item in factors(divisor, factor))
⍝
⍝ factor -= 1
⍝
⍝ return result
⍝
⍝ print factors(12) # -> [[6, 2], [4, 3], [3, 2, 2]]
⍝ print factors(24) # -> [[12, 2], [8, 3], [6, 4], [6, 2, 2], [4, 3, 2], [3, 2, 2, 2]]
⍝ print factors(157) # -> []
DIV
divisors←{z[⍋z←∊∘.×/1,¨(∪π⍵)*¨⍳¨∪⍦π⍵]}
RESULTSfor100
1 1
2 2
3 4
4 6
5 16
6 12
7 64
8 24
9 36
10 48
11 1024
12 60
13 4096
14 192
15 144
16 120
17 65536
18 180
19 262144
20 240
21 576
22 3072
23 4194304
24 360
25 1296
26 12288
27 900
28 960
29 268435456
30 720
31 1073741824
32 840
33 9216
34 196608
35 5184
36 1260
37 68719476736
38 786432
39 36864
40 1680
41 1099511627776
42 2880
43 4398046511104
44 15360
45 3600
46 12582912
47 70368744177664
48 2520
49 46656
50 6480
51 589824
52 61440
53 4503599627370496
54 6300
55 82944
56 6720
57 2359296
58 805306368
59 288230376151711744
60 5040
61 1152921504606846976
62 3221225472
63 14400
64 7560
65 331776
66 46080
67 73786976294838206464
68 983040
69 37748736
70 25920
71 1180591620717411303424
72 10080
73 4722366482869645213696
74 206158430208
75 32400
76 3932160
77 746496
78 184320
79 302231454903657293676544
80 15120
81 44100
82 3298534883328
83 4835703278458516698824704
84 20160
85 5308416
86 13194139533312
87 2415919104
88 107520
89 309485009821345068724781056
90 25200
91 2985984
92 62914560
93 9663676416
94 211106232532992
95 21233664
96 27720
97 79228162514264337593543950336
98 233280
99 230400
100 45360
這裏是採取62微秒 C程序!
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define min(a,b) (((a) < (b)) ? (a) : (b))
#define MaxInt 0x7fffffffffffffff
#define ll long long
int primes[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61,
67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131};
ll minimums[16]; // contains minimum for each level
ll ipow(ll base, int exp) {
ll result = 1;
while (exp) {
if (exp & 1) result *= base;
exp >>= 1; base *= base;
} return result;
}
ll factors(ll number, ll max_factor, int level) {
ll result = MaxInt;
ll factor = min(number/2, max_factor);
while (factor >= 2) {
if (number % factor == 0) {
ll divisor = number/factor;
ll tempf = ipow(primes[level], factor-1);
if (divisor <= factor && divisor <= max_factor) {
ll tempd = ipow(primes[level+1], divisor-1);
minimums[level] = min(minimums[level], tempf * tempd);
}
ll fct = factors(divisor, factor, level+1);
if (fct < MaxInt)
minimums[level] = min(minimums[level], tempf * fct);
result = minimums[level];
minimums[level+1] = MaxInt;
}
factor -= 1;
}
return result;
}
ll fact(int number) {
for (int level = 0; level < 16; level++) minimums[level] = MaxInt;
ll res = factors(number, MaxInt, 0);
if (res < MaxInt) return res;
else return 0;
}
int main(int argc, char *argv[]) {
int N = 100;
if (argc > 1) N = atoi(argv[1]);
ll res[N];
clock_t Start = clock();
for(int n = 1; n <= 10000; n++)
for (int i = 1; i <= N; i++) res[i] = fact(i);
printf("\n%0.6f second(s)\n", (clock() - Start)/1000.0/10000.0);
for (int i = 1; i <= N; i++)
if (res[i] > 0) printf("%d %lld\n", i, res[i]);
else if (i > 64) printf("%d 2 power %d\n", i, i-1);
else printf("%d %lld\n", i, ipow(2, i-1));
return 0;
}