''' Created on 06.08.2016 @author: wn ''' import time import random startTime = 0 lastTime = 0 def log(msg, reset=False): global startTime, lastTime currentTime = time.time() if reset: startTime = currentTime if lastTime == 0: lastTime = currentTime print("%f %f: %s" % (currentTime - startTime, currentTime - lastTime, msg)) lastTime = currentTime def s_permutations(seq): items = [[]] for j in seq: new_items = [] for i, item in enumerate(items): if i % 2: # step up new_items += [item[:i] + [j] + item[i:] for i in range(len(item) + 1)] else: # step down new_items += [item[:i] + [j] + item[i:] for i in range(len(item), -1, -1)] items = new_items return [(tuple(item), -1 if i % 2 else 1) for i, item in enumerate(items)] cachedIndexJ = {} def det(m): log("Calculating determinant started") l = len(m) indexListI = range(1, l+1) if l in cachedIndexJ: indexListJ = cachedIndexJ[l] log("J for %i taken from cache" % l) else: indexListJ = s_permutations(indexListI) cachedIndexJ[l] = indexListJ cnt = len(indexListJ) #print(indexListJ) log("J for %i generated, %i entries" % (l, cnt)) summ = 0 steps = 0 for indexJ in indexListJ: sign = indexJ[1] prod = sign for indexI in indexListI: prod *= m[indexI-1][indexJ[0][indexI-1]-1] steps += 1 summ += prod log("Calculating determinant done, %i steps" % steps) return summ def solve(g): coeffMatrix = g[:-1] coeffDet = det(coeffMatrix) order = len(g)-1 unknowns = [] for i in range(order): varMatrix = [] for j in range(order): if (i == j): varMatrix.append(g[order]) else: varMatrix.append(g[j]) varDet = det(varMatrix) var = varDet / coeffDet unknowns.append(var) return unknowns def generate(order, lowerBound, upperBound): log("Generating system started") matrix = [] for i in range(order+1): vector = [] for j in range(order): vector.append(random.uniform(lowerBound, upperBound)) matrix.append(vector) log("Generating system done") return matrix def verify(m, v): coeffs = m[:-1] results = m[-1] epsilons = [] good = True MAX_EPSILON = 1E-10 worstEpsilon = 0.0 for i in range(len(results)): res = 0.0 for j in range(len(results)): res += coeffs[j][i] * v[j] epsilon = res - results[i] if epsilon > MAX_EPSILON: good = False if abs(epsilon) > worstEpsilon: worstEpsilon = abs(epsilon) epsilons.append(res - results[i]) return (good, worstEpsilon, epsilons) #g = [[1.0,2.0,10.0], [5.0,3.0,11.0], [7.0,4.0,3.0], [8.0,9.0,13.0]] #g = [[599.1290806728277, 80.40433203454268, 615.5120772027384], [187.33506822036273, 663.7490700388757, 677.5997857592452], [956.929845461785, 456.8619830134234, 419.21973493934263], [61.824473951393344, 342.8555419523176, 216.32654519677496]] #g = [[8.019558599768107, 903.9449140455304, 919.612242767639, 753.9105365307668, 59.01030128588536, 920.4023338915356, 202.58884034378954, 431.09289053573997, 339.17491130065537, 313.11544472585905], [317.26651051567023, 727.3387587224795, 152.77343043145507, 408.44599685908923, 155.79032672640025, 425.6401391270683, 864.8657207590146, 298.9063131849389, 430.3786466384899, 425.29127250155807], [865.8564877305874, 669.481784333669, 538.9130960586408, 961.8480073483083, 753.9136870076827, 912.3571914611318, 929.0546428258709, 873.5276061317886, 459.59802608296684, 455.8093236762648], [318.8613025219013, 497.149024070492, 343.8577826382859, 780.6069003672162, 319.67906522418554, 628.3446507753014, 37.526503978800974, 490.48198671795006, 527.0266127960343, 180.8362568347207], [156.8237212311948, 510.1641000612821, 229.2078194439, 637.1024703370609, 972.0604419040393, 571.3909858582223, 184.34543232118418, 479.07607501617986, 382.12667148701297, 150.05251231416116], [297.8478733948842, 981.2792974599201, 271.8121041641507, 990.1956266273802, 645.839991521914, 288.9135934264093, 845.7496599336642, 19.88633912497506, 97.46153065555741, 343.1102890796861], [384.5208626880663, 145.38972551894247, 560.1624430301523, 224.07537517241872, 704.7921876231894, 781.2588369275445, 752.0786768851573, 77.91470589306748, 851.5354388193074, 146.64240712820654], [252.4899031289375, 928.4755757806711, 176.3416092863044, 914.652591937962, 379.5647442648579, 418.6432906860023, 397.44937187401274, 943.8352970973906, 495.8158468453109, 746.6463673774376], [330.80903398935936, 560.2262982840332, 881.0811315973745, 181.77380227358952, 315.50808931102347, 158.5072721277919, 947.7031957044982, 929.3796342482069, 131.55810107335486, 982.8660461523164], [310.9354667285096, 770.115136026887, 334.49254019770905, 236.96386990210726, 918.8431629974137, 893.3379619077041, 431.3668259114833, 285.69624238120986, 64.53926288027478, 799.2392946390154], [20.70259163693389, 306.46447894527165, 420.64284211501723, 37.622373302959474, 336.4417266492545, 247.6600784425259, 215.21800533449152, 426.10683135701044, 783.5055424440518, 817.507304763315]] g = [[196.35648286602913, 255.78320214562078, 254.4957402769741, 248.69977961153933, 382.61065017020854, 62.17175785611484, 210.46041763171675, 792.184861527093, 492.86614662413433, 16.418566038295502], [465.59032006096265, 419.1822656204234, 276.93400418307357, 92.75012020060058, 633.9464333410591, 363.96427427067556, 330.61829018827103, 458.1460490942364, 752.3719990444775, 612.9260145499019], [232.7461763483454, 47.02378341457059, 46.89987060240875, 780.5294679474976, 40.868138953817244, 923.4530470705937, 793.3904771973274, 450.0896614228864, 455.645931783578, 212.43143142136833], [56.8274840159122, 66.92825925761848, 977.0220435032436, 688.2546262621239, 618.2761789189016, 9.551896254965552, 863.7431636261247, 602.4057398468475, 787.3808193741036, 422.22770868828144], [374.75867750739735, 634.3158849453325, 432.12940030241697, 787.5745064269738, 573.6257796824623, 902.975285283053, 890.988258484369, 605.4127514364945, 874.1015867803133, 511.5082947334636], [4.822237509823135, 522.9093499283525, 256.68884816246805, 122.51633985600807, 900.3010320832226, 748.4667048177383, 812.8262957013549, 789.2451719039427, 250.48739944027542, 982.5288431098245], [378.8053926479358, 602.2300484993539, 364.26457312412964, 654.3503288596467, 622.8637580845797, 565.5166875112637, 522.8219206732364, 276.1018531710755, 536.7176663236127, 849.518815976257], [476.8990073351791, 895.5485647979485, 398.2982159500944, 45.21333518900506, 535.759609317032, 55.798948045607254, 802.5995294976824, 599.7112194327568, 817.3722187418982, 291.47777811233624], [486.3382733102173, 596.0130673435432, 442.2598417565975, 149.36250900806724, 592.7509203792288, 47.57333396259733, 385.5156216228283, 246.31479580712156, 708.1101889638189, 896.7067437732082], [676.3032176471872, 896.5079308000833, 523.1513042144225, 582.0153418896383, 609.8731927938023, 741.1781501728287, 21.299661773589772, 941.7391463459769, 556.6728614082933, 627.576228023791], [955.1918839549409, 900.0382980231512, 166.1210836385043, 164.25307272366385, 767.6446810122188, 457.81549580183525, 792.5516181980048, 232.54376961359858, 254.9754925437675, 666.654630723243]] # i = 11 log("Starting ...", reset=True) # g = generate(i, 0.0, 1000.0) v = solve(g) log("Solved") # print(i) print(g) print(v) e = verify(g, v) print(e) print