API tools faq
paste
Advertisement
UF6

Problem 3

UF6
Feb 28th, 2016
306
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 11.82 KB | None | 0 0
raw download clone embed print report
  1.  
  2. from decimal import getcontext, Decimal
  3. from fractions import Fraction
  4. import sys
  5.  
  6. class LUDecomposition():
  7.     def __init__(self, A):
  8.         self.LU = A.get_array_copy()
  9.         self.m = A.m
  10.         self.n = A.n
  11.         self.pivot_sign = 1
  12.         self.pivot_vector = [i for i in range(self.m)]
  13.  
  14.         LU_row_i = []
  15.         LU_col_j = [0 for i in range(self.m)]
  16.         for j in range(self.n):
  17.             for i in range(self.m):
  18.                 LU_col_j[i] = self.LU[i][j]
  19.             for i in range(self.m):
  20.                 LU_row_i = self.LU[i]
  21.                 s = sum([LU_row_i[k] * LU_col_j[k] for k in range(min(i, j))])
  22.                 LU_row_i[j] -= s
  23.                 LU_col_j[i] -= s
  24.  
  25.             p = j
  26.             for i in range(j + 1, self.m):
  27.                 if abs(LU_col_j[i]) > abs(LU_col_j[p]):
  28.                     p = i
  29.             if p != j:
  30.                 for k in range(self.n):
  31.                     self.LU[p][k], self.LU[j][k] = self.LU[j][k], self.LU[p][k]
  32.                 self.pivot_vector[p], self.pivot_vector[j] = self.pivot_vector[j], self.pivot_vector[p]
  33.                 self.pivot_sign = -self.pivot_sign
  34.              
  35.             if j < self.m and self.LU[j][j] != 0:
  36.                 for i in range(j + 1, self.m):
  37.                     self.LU[i][j] *= Fraction(1, self.LU[j][j])
  38.  
  39.     def solve(self, B):
  40.         assert(B.m == self.m)
  41.         assert(self.nonsingular())
  42.  
  43.         X = B.get_submatrix(self.pivot_vector, 0, B.n - 1)
  44.         for k in range(self.n):
  45.             for i in range(k + 1, self.n):
  46.                 for j in range(B.n):
  47.                     X.entries[i][j] -= self.LU[i][k] * X.entries[k][j]
  48.  
  49.         for k in range(self.n - 1, -1, -1):
  50.             for j in range(B.n):
  51.                 X.entries[k][j] *= Fraction(1, self.LU[k][k])
  52.             for i in range(k):
  53.                 for j in range(B.n):
  54.                     X.entries[i][j] -= self.LU[i][k] * X.entries[k][j]
  55.         return X
  56.  
  57.     def nonsingular(self):
  58.         for j in range(self.n):
  59.             if self.LU[j][j] == 0:
  60.                 return False
  61.         return True
  62.  
  63.  
  64.  
  65. class Matrix():
  66.     def identity(m, n):
  67.         matrix = Matrix(m=m, n=n)
  68.         for i in range(matrix.m):
  69.             for j in range(matrix.n):
  70.                 if i == j:
  71.                     matrix.entries[i][j] = 1
  72.         return matrix
  73.  
  74.     def __init__(self, *args, **kwargs):
  75.         self.entries = None
  76.         self.m = None
  77.         self.n = None
  78.         if 'm' in kwargs and 'n' in kwargs and 'entries' not in kwargs:
  79.             self.m = kwargs['m']
  80.             self.n = kwargs['n']
  81.             self.entries = [[0 for j in range(self.n)] for i in range(self.m)]
  82.         elif 'm' not in kwargs and 'n' not in kwargs and 'entries' in kwargs:
  83.             self.m = len(kwargs['entries'])
  84.             self.n = len(kwargs['entries'][0])
  85.             self.entries = kwargs['entries']
  86.  
  87.     def get_array_copy(self):
  88.         copy = [[0 for j in range(self.n)] for i in range(self.m)]
  89.         for i in range(self.m):
  90.             for j in range(self.n):
  91.                 copy[i][j] = self.entries[i][j]
  92.         return copy
  93.  
  94.     def get_submatrix(self, row_index_array, begin_column_index, end_column_index):
  95.         submatrix = Matrix(m=len(row_index_array), n=end_column_index - begin_column_index + 1)
  96.         for i in range(len(row_index_array)):
  97.             for j in range(begin_column_index, end_column_index + 1):
  98.                 submatrix.entries[i][j - begin_column_index] = self.entries[row_index_array[i]][j]
  99.         return submatrix
  100.  
  101.     def minus(self, B):
  102.         X = Matrix(m=self.m, n=self.n)
  103.         for i in range(self.m):
  104.             for j in range(self.n):
  105.                 X.entries[i][j] = self.entries[i][j] - B.entries[i][j]
  106.         return X
  107.  
  108.     def times(self, scalar):
  109.         X = Matrix(m=self.m, n=self.n)
  110.         for i in range(self.m):
  111.             for j in range(self.n):
  112.                 X.entries[i][j] = scalar * self.entries[i][j]
  113.         return X
  114.  
  115.     def get_sum_of_entries(self):
  116.         result = 0
  117.         for i in range(self.m):
  118.             for j in range(self.n):
  119.                 result += self.entries[i][j]
  120.         return result
  121.  
  122.     def inverse(self):
  123.         return self.solve(Matrix.identity(self.m, self.m))
  124.  
  125.     def solve(self, B):
  126.         assert(self.m == self.n)
  127.         return LUDecomposition(self).solve(B)
  128.  
  129.  
  130.  
  131. class KempnerSeries():
  132.     def __init__(self, T):
  133.         self.T = T
  134.         self.n = len(T) - 1
  135.         self.epsilon = Fraction(1, 10**20)
  136.         self.psi = []
  137.  
  138.     def get_sum(self):
  139.         total_sum = self.get_one_digit_sum()
  140.         total_sum += self.get_two_digits_sum()
  141.         total_sum += self.get_more_digits_sum()
  142.         total_sum += self.get_remaining_sum()
  143.         return total_sum
  144.  
  145.     def get_one_digit_sum(self):
  146.         result = 0
  147.         for s in range(1, 10):
  148.             if self.T[0][s]:
  149.                 result += Fraction(1, s)
  150.         return result
  151.  
  152.     def get_two_digits_sum(self):
  153.         k = 1
  154.         while True:
  155.             v = Matrix(m=self.n, n=1)
  156.             for m1 in range(1, 10):
  157.                 for m2 in range(10):
  158.                     i = self.T[0][m1]
  159.                     if i and self.T[i][m2]:
  160.                         s = m1 * 10 + m2
  161.                         j = self.T[i][m2]
  162.                         v.entries[j - 1][0] += Fraction(1, s**k)
  163.             if k > 1 and v.get_sum_of_entries() < self.epsilon:
  164.                 break
  165.             self.psi.append(v)
  166.             k += 1
  167.         return self.psi[0].get_sum_of_entries()
  168.  
  169.     def get_more_digits_sum(self):
  170.         total_sum = 0
  171.         i = 3
  172.         while len(self.psi) > 1:
  173.             k = 1
  174.             while k <= len(self.psi):
  175.                 v = Matrix(m=self.n, n=1)
  176.                 for l in range(1, self.n + 1):
  177.                     for m in range(10):
  178.                         j = self.T[l][m]
  179.                         if not j:
  180.                             continue
  181.                         for w in range(len(self.psi) - k + 1):
  182.                             v.entries[j - 1][0] += self.get_coefficient_a(k, w, m) * self.psi[k - 1 + w].entries[l - 1][0]
  183.                 if k > 1 and v.get_sum_of_entries() < self.epsilon:
  184.                     self.psi = self.psi[:k - 1]
  185.                     break
  186.                 else:
  187.                     self.psi[k - 1] = v
  188.                 k += 1
  189.             total_sum += self.psi[0].get_sum_of_entries()
  190.             i += 1
  191.         return total_sum
  192.  
  193.     def get_remaining_sum(self):
  194.         b = self.get_vector_b(self.get_matrix_B(self.get_matrix_A()))
  195.         return self.inner_product(b, self.psi[0])
  196.  
  197.     def get_matrix_A(self):
  198.         result = Matrix(m=self.n, n=self.n)
  199.         for m in range(10):
  200.             for l in range(1, self.n + 1):
  201.                 for j in range(1, self.n + 1):
  202.                     if self.T[l][m] == j:
  203.                         result.entries[j - 1][l - 1] += 1
  204.         return result.times(Fraction(1, 10))
  205.    
  206.     def get_matrix_B(self, A):
  207.         identity = Matrix.identity(m=self.n, n=self.n)
  208.         return identity.minus(A).inverse().minus(identity)
  209.  
  210.     def get_vector_b(self, B):
  211.         result = Matrix(m=self.n, n=1)
  212.         for l in range(1, self.n + 1):
  213.             for j in range(1, self.n + 1):
  214.                 result.entries[l - 1][0] += B.entries[j - 1][l - 1]
  215.         return result
  216.  
  217.     def get_coefficient_a(self, k, w, m):
  218.         return Fraction(1, 10**(k + w)) * (-1)**w * self.get_binomial_coefficient(k + w - 1, w) * m**w
  219.  
  220.     def get_binomial_coefficient(self, n, m):
  221.         if m > n - m:
  222.             return self.get_binomial_coefficient(n, n - m)
  223.         result = 1
  224.         i, j = 1, n
  225.         while i <= m:
  226.             result *= Fraction(j , i)
  227.             i, j = i + 1, j - 1
  228.         return result
  229.    
  230.     def inner_product(self, A, B):
  231.         assert(A.m == B.m)
  232.         assert(A.n == 1)
  233.         assert(B.n == 1)
  234.         result = 0
  235.         for i in range(A.m):
  236.             result += A.entries[i][0] * B.entries[i][0]
  237.         return result
  238.  
  239.  
  240.  
  241. class Problem():
  242.     def __init__(self):
  243.         getcontext().prec = 20
  244.  
  245.     def solve(self):
  246.         #self.benchmark()
  247.         self.get()
  248.  
  249.     def get(self):
  250.         T = [
  251.                 [ 0, 2, 3, 4, 5, 6, 7, 8, 9, 10 ],
  252.                 [ 11, 2, 3, 4, 5, 6, 7, 8, 9, 10 ],
  253.                 [ 1, 12, 3, 4, 5, 6, 7, 8, 9, 10 ],
  254.                 [ 1, 2, 13, 4, 5, 6, 7, 8, 9, 10 ],
  255.                 [ 1, 2, 3, 14, 5, 6, 7, 8, 9, 10 ],
  256.                 [ 1, 2, 3, 4, 15, 6, 7, 8, 9, 10 ],
  257.                 [ 1, 2, 3, 4, 5, 16, 7, 8, 9, 10 ],
  258.                 [ 1, 2, 3, 4, 5, 6, 17, 8, 9, 10 ],
  259.                 [ 1, 2, 3, 4, 5, 6, 7, 18, 9, 10 ],
  260.                 [ 1, 2, 3, 4, 5, 6, 7, 8, 19, 10 ],
  261.                 [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 20 ],
  262.                 [ 0, 2, 3, 4, 5, 6, 7, 8, 9, 10 ],
  263.                 [ 1, 0, 3, 4, 5, 6, 7, 8, 9, 10 ],
  264.                 [ 1, 2, 0, 4, 5, 6, 7, 8, 9, 10 ],
  265.                 [ 1, 2, 3, 0, 5, 6, 7, 8, 9, 10 ],
  266.                 [ 1, 2, 3, 4, 0, 6, 7, 8, 9, 10 ],
  267.                 [ 1, 2, 3, 4, 5, 0, 7, 8, 9, 10 ],
  268.                 [ 1, 2, 3, 4, 5, 6, 0, 8, 9, 10 ],
  269.                 [ 1, 2, 3, 4, 5, 6, 7, 0, 9, 10 ],
  270.                 [ 1, 2, 3, 4, 5, 6, 7, 8, 0, 10 ],
  271.                 [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 0 ],
  272.         ]
  273.         print(self.rep(KempnerSeries(T).get_sum()))
  274.  
  275.     def benchmark(self):
  276.         print('get_sum_with_no_digit_nine =>', self.get_sum_with_no_digit_nine())
  277.         print('get_sum_with_no_even_digits =>', self.get_sum_with_no_even_digits())
  278.         print('get_sum_with_no_odd_digits =>', self.get_sum_with_no_odd_digits())
  279.         print('get_sum_with_no_42 =>', self.get_sum_with_no_42())
  280.         print('get_sum_with_no_314 =>', self.get_sum_with_no_314())
  281.         print('get_sum_with_no_even_digits_no_55_no_13579 =>', self.get_sum_with_no_even_digits_no_55_no_13579())
  282.  
  283.     def get_sum_with_no_digit_nine(self):
  284.         T = [
  285.                 [ 0, 1, 1, 1, 1, 1, 1, 1, 1, 0 ],
  286.                 [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 0 ],
  287.         ]
  288.         return self.rep(KempnerSeries(T).get_sum())
  289.  
  290.     def get_sum_with_no_even_digits(self):
  291.         T = [
  292.                 [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 ],
  293.                 [ 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 ],
  294.         ]
  295.         return self.rep(KempnerSeries(T).get_sum())
  296.  
  297.     def get_sum_with_no_odd_digits(self):
  298.         T = [
  299.                 [ 0, 0, 1, 0, 1, 0, 1, 0, 1, 0 ],
  300.                 [ 1, 0, 1, 0, 1, 0, 1, 0, 1, 0 ],
  301.         ]
  302.         return self.rep(KempnerSeries(T).get_sum())
  303.  
  304.     def get_sum_with_no_42(self):
  305.         T = [
  306.                 [ 0, 1, 1, 1, 2, 1, 1, 1, 1, 1 ],
  307.                 [ 1, 1, 1, 1, 2, 1, 1, 1, 1, 1 ],
  308.                 [ 1, 1, 0, 1, 2, 1, 1, 1, 1, 1 ],
  309.         ]
  310.         return self.rep(KempnerSeries(T).get_sum())
  311.  
  312.     def get_sum_with_no_314(self):
  313.         T = [
  314.                 [ 0, 1, 1, 2, 1, 1, 1, 1, 1, 1 ],
  315.                 [ 1, 1, 1, 2, 1, 1, 1, 1, 1, 1 ],
  316.                 [ 1, 3, 1, 2, 1, 1, 1, 1, 1, 1 ],
  317.                 [ 1, 1, 1, 2, 0, 1, 1, 1, 1, 1 ],
  318.         ]
  319.         return self.rep(KempnerSeries(T).get_sum())
  320.  
  321.     def get_sum_with_no_even_digits_no_55_no_13579(self):
  322.         T = [
  323.                 [ 0, 3, 0, 1, 0, 2, 0, 1, 0, 1 ],
  324.                 [ 0, 3, 0, 1, 0, 2, 0, 1, 0, 1 ],
  325.                 [ 0, 3, 0, 1, 0, 0, 0, 1, 0, 1 ],
  326.                 [ 0, 3, 0, 4, 0, 2, 0, 1, 0, 1 ],
  327.                 [ 0, 3, 0, 1, 0, 5, 0, 1, 0, 1 ],
  328.                 [ 0, 3, 0, 1, 0, 0, 0, 6, 0, 1 ],
  329.                 [ 0, 3, 0, 1, 0, 2, 0, 1, 0, 0 ],
  330.         ]
  331.         return self.rep(KempnerSeries(T).get_sum())
  332.  
  333.     def rep(self, fraction):
  334.         return Decimal(fraction.numerator) / Decimal(fraction.denominator)
  335.  
  336. def main():
  337.     problem = Problem()
  338.     problem.solve()
  339.  
  340. if __name__ == '__main__':
  341.     sys.exit(main())
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!