/* k 1 1 1 1 1 1 n 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 1 3 1 0 0 0 0 0 0 0 0 0 0 0 0 4 0 1 7 6 1 0 0 0 0 0 0 0 0 0 0 0 5 0 1 15 25 10 1 0 0 0 0 0 0 0 0 0 0 6 0 1 31 90 65 15 1 0 0 0 0 0 0 0 0 0 7 0 1 63 301 350 140 21 1 0 0 0 0 0 0 0 0 8 0 1 127 966 1701 1050 266 28 1 0 0 0 0 0 0 0 9 0 1 255 3025 7770 6951 2646 462 36 1 0 0 0 0 0 0 10 0 1 511 9330 34105 42525 22827 5880 750 45 1 0 0 0 0 0 11 0 1 1023 28501 145750 246730 179487 63987 11880 1155 55 1 0 0 0 0 12 0 1 2047 86526 611501 1379400 1323652 627396 159027 22275 1705 66 1 0 0 0 13 0 1 4095 261625 2532530 7508501 9321312 5715424 1899612 359502 39325 2431 78 1 0 0 14 0 1 8191 788970 10391745 40075035 63436373 49329280 20912320 5135130 752752 66066 3367 91 1 0 15 0 1 16383 2375101 42355950 210766920 420693273 408741333 216627840 67128490 12662650 1479478 106470 4550 105 1 coincidence? Two formula's: k S2(n,k) = 1/k! * SIGMA[(-1)^(k-j) * C(k,j) * j^n ] j=0 S2(n,k) = k * S2(n-1,k) + S2(n-1,k-1) Initial value's: S2(j,j) = 1 , S2(j,0) = 0 , S2(0,j) = 0 The first formula: * powers, * / factorials, + - (large) numbers . The second one leads to a recursive solution. It's going top-down, a larger number depends on two smaller ones. Bottom-up seems to be an easier, faster, iterative way. Top-down: S2(5,3) = 3 * S2(4,3) + S2(4,2) S2(4,3) = 3 * S2(3,3) + S2(3,2) S2(4,2) = 2 * S2(3,2) + S2(2,2) S2(3,2) = 2 * S2(2,2) + S2(2 ... ...complicated... Bottom-up: Rearrange table, remove empty (zero) items in colums: _ _ _ _ | 1 | | 1 1 1 1 | | 0 1 | =====> | 0 1 3 6 | | 0 1 1 | |_0 1 7 25_| | . 1 3 1 | | . . 7 6 | |_. . . 25_| 1 1 1 + 2 * 0 = 1 1 + 3 * 0 = 1 0 1 1 + 2 * 1 = 3 3 + 3 * 1 = 6 0 1 1 + 2 * 3 = 7 7 + 3 * 6 = 25 = S2(5,3) ...less complicated... |---------------------------------------------------------------| | S2(n,k) | |------|------|-----------------------------------------|-------| | n | k | Time in ms (Athlon X4, XP, 2GB) | bits | |------|------|---------|---------|----------|----------|-------| | 100 | 10 | 0,25 | | | | 311 | | 100 | 25 | | 0,61 | | | 381 | | 100 | 50 | | | 0,81 | | 338 | | 100 | 89 | | | | 0,27 | 108 | | | | | | | | | | 200 | 20 | 1,5 | | | | 804 | | 200 | 43 | | 3,1 | | | 910 | | 200 | 86 | | | 4,4 | | 838 | | 200 | 172 | | | | 1,6 | 295 | | | | | | | | | | 400 | 40 | 9,4 | | | | 1970 | | 400 | 77 | | 18 | | | 2131 | | 400 | 154 | | | 27 | | 1986 | | 400 | 308 | | | | 14 | 982 | | | | | | | | | | 800 | 80 | 75 | | | | 4663 | | 800 | 138 | | 124 | | | 4900 | | 800 | 276 | | | 183 | | 4617 | | 800 | 552 | | | | 119 | 2744 | | | | | | | | | | 1600 | 160 | 672 | | | | 10770 | | 1600 | 250 | | 1010 | | | 11109 | | 1600 | 500 | | | 1620 | | 10546 | | 1600 | 1000 | | | | 1120 | 6974 | | | | | | | | | | 3200 | 320 | 7230 | | | | 24424 | | 3200 | 456 | | 9930 | | | 24889 | | 3200 | 912 | | | 15300 | | 23765 | | 3200 | 1824 | | | | 11900 | 16881 | |------|------|---------|---------|----------|----------|-------| */ using System; using System.Diagnostics; using Xint = System.Numerics.BigInteger; class Stirling_2nd_kind { private static Xint S2(uint n, uint k) { if (k >= n) return k == n ? 1 : 0; if (k <= 2) return k == 0 ? 0 : k == 1 ? 1 : (Xint.One << (int)(n - 1)) - 1; if (k + 1 == n) return (Xint)(n) * k / 2; n -= k - 1; uint i = 0, j = 3; Xint[] A = new Xint[n]; for (; i < n; i++) A[i] = 1; Xint[] B = new Xint[n]; for (i = 0; i < n; i++) B[i] = (Xint.One << (int)(i + 1)) - 1; do { for (i = 1; i < n; i++) A[i] = A[i - 1] * j + B[i]; j++; if (j > k) return A[n - 1]; for (i = 1; i < n; i++) B[i] = B[i - 1] * j + A[i]; j++; } while (j <= k); return B[n - 1]; } // to do: - use uints/ulongs: - for smaller value's / // - as long as possible // - return row (array) // - return previous/next row from row // - k close to n , go from right to left in table ? // - same/similar trick for: - other sequences ? // - sequences using S2(n,k) ? private static Stopwatch sw = new Stopwatch(); static void Main() { Console.WriteLine("S2(50,17) = " + S2(50, 17)); Console.WriteLine(" bits : " + bL(S2(50, 17))); sw.Restart(); for (int i = 0; i < 1000; i++) S2(50, 17); sw.Stop(); Console.WriteLine(" ms : 0," + sw.ElapsedMilliseconds); Console.ReadLine(); // S2(50,17) = 37645241791600906804871080818625037726247519045 // bits : 155 // ms : 0,134 } private static int bL(Xint U) { byte[] bytes = (U.Sign * U).ToByteArray(); int i = bytes.Length - 1; return i * 8 | bitLengthMostSignificantByte(bytes[i]); } private static int bitLengthMostSignificantByte(byte b) { return b < 08 ? b < 02 ? b < 01 ? 0 : 1 : b < 04 ? 2 : 3 : b < 32 ? b < 16 ? 4 : 5 : b < 64 ? 6 : 7; } }
No comments:
Post a Comment