/* Example: N D SD SSD A024916 1 1 1 1 2 1,2 3 4 3 1,3 4 8 4 1,2,4 7 15 5 1,5 6 21 6 1,2,3,6 12 33 Sum of all divisors of all positive integers <= 10^n. A072692 as a simple table n a(n) time SSD SSD3 0 1 00:00:00.0000061 1 87 00:00:00.0000058 2 8299 00:00:00.0000092 3 823081 00:00:00.0000215 4 82256014 00:00:00.0000555 5 8224740835 00:00:00.0002824 6 822468118437 00:00:00.0007623 7 82246711794796 00:00:00.0022483 8 8224670422194237 00:00:00.0075755 9 822467034112360628 00:00:00.0247891 10 82246703352400266400 00:00:00.1029980 00:00:00.0536221 11 8224670334323560419029 00:00:00.3310292 12 822467033425357340138978 00:00:01.0792704 00:00:00.5139608 13 82246703342420509396897774 00:00:03.3665267 14 8224670334241228180927002517 00:00:10.7148327 00:00:05.1645348 15 822467033424114009326065894639 00:00:34.2567848 16 82246703342411333689227187822414 00:01:49.7741595 00:00:53.1228885 17 8224670334241132270081671519064067 00:05:55.3252935 18 822467033424113219363487627735401433 00:19:50.2083007 19 82246703342411321831834750379375676781 01:12:36.6691244 20 8224670334241132182459113152714495956789 04:38:51.4847697 21 822467033424113218236966661186847524013409 15:19:44.3197506 AMD Athlon II X4 640, 3 GHz, 2 GB, XP N=10^16 time(s) CPU Usage(%) SSD(N) 110 25 SSD3(N) 53 60 FSSD(N) 92 25 FSSD3(N) 41 60 FSSD4(N) 33 80
*/ using System; using System.Diagnostics; using System.Threading.Tasks; using Xint = System.Numerics.BigInteger; class A024916 { static Xint SSD(Xint N) { Xint D = 1, Q = N, S = 0; while (D < Q) { S += D++ * (Q - D); S += Q++ * Q / 2; Q = N / D; } S += Q++ * Q / 2; S -= D-- * D-- * D / 6; return S; } static void Main() { var sw = new Stopwatch(); int p = 0; Xint N = 1, S; L0: sw.Restart(); S = SSD(N); sw.Stop(); Console.WriteLine(p + " " + S + " " + sw.Elapsed); p++; N *= 10; goto L0; } private static Xint[] SA, DA, QA; static Xint SSD3(Xint N) { SA = new Xint[3]; DA = new Xint[3]; QA = new Xint[3]; Parallel.For(0, 3, (int d) => SSDP3(N, d)); Xint S = SA[0] + SA[1] + SA[2], D = DA[0], Q = QA[0]; if (DA[1] < D) { D = DA[1]; Q = QA[1]; } if (DA[2] < D) { D = DA[2]; Q = QA[2]; } S += Q++ * Q / 2; S -= D-- * D-- * D / 6; return S; } static void SSDP3(Xint N, int d) { Xint D = d + 1, Q = N / D, S = 0; while (D < Q) { S += D * (Q - D - 1); S += Q * (Q + 1) / 2; D += 3; Q = N / D; } SA[d] = S; DA[d] = D; QA[d] = Q; } static Xint FSSD(Xint N) { Xint D = 1, Q = N, S = 0; while (D < Q) { S += Q * (Q + 1 + 2 * D) / 2; D += 1; Q = N / D; } S += Q * (Q + 1) / 2; S -= D * D * (D - 1) / 2; return S; } static Xint FSSD3(Xint N) { SA = new Xint[3]; DA = new Xint[3]; QA = new Xint[3]; Parallel.For(0, 3, (int d) => FSSDP3(N, d)); Xint S = SA[0] + SA[1] + SA[2], D = DA[0], Q = QA[0]; if (DA[1] < D) { D = DA[1]; Q = QA[1]; } if (DA[2] < D) { D = DA[2]; Q = QA[2]; } S += Q * (Q + 1) / 2; S -= D * D * (D - 1) / 2; return S; } static void FSSDP3(Xint N, int d) { Xint D = d + 1, Q = N / D, S = 0; while (D < Q) { S += Q * (Q + 1 + 2 * D) / 2; D += 3; Q = N / D; } SA[d] = S; DA[d] = D; QA[d] = Q; } static Xint FSSD4(Xint N) { SA = new Xint[4]; DA = new Xint[4]; QA = new Xint[4]; Parallel.For(0, 4, (int d) => FSSDP4(N, d)); Xint S = SA[0] + SA[1] + SA[2] + SA[3], D = DA[0], Q = QA[0]; if (DA[1] < D) { D = DA[1]; Q = QA[1]; } if (DA[2] < D) { D = DA[2]; Q = QA[2]; } if (DA[3] < D) { D = DA[3]; Q = QA[3]; } S += Q * (Q + 1) / 2; S -= D * D * (D - 1) / 2; return S; } static void FSSDP4(Xint N, int d) { Xint D = d + 1, Q = N / D, S = 0; while (D < Q) { S += Q * (Q + 1 + 2 * D) / 2; D += 4; Q = N / D; } SA[d] = S; DA[d] = D; QA[d] = Q; } }
No comments:
Post a Comment