View difference between Paste ID: HRs8YtLi and NXAYTLXV
SHOW: | | - or go back to the newest paste.
1
from itertools import* #reikalinga funkcijai starmap
2
#SKAICIU TEORIJA
3
def primes(n, sieve=[False,False,True]): #by-need calculation of primes; not includes n
4
    #TEST1
5
    #sieve:          [1,1,1,1,1,1,1,1,1,1]=[1]*10
6
    #possible primes [0,1,2,3,4,5,6,7,8,9]
7
    #given sieve     [0,0,1,1,0]
8
    #a=3
9
    #check from 0 to sqrt(len(sieve))
10
    #last five entries were checked dividing by primes 0 to a=3
11
    #sieve is given as list of bools to save memory
12
    #TEST2
13
    #primes(20) returns [2, 3, 5, 7, 11, 13, 17, 19], sieve is [False, False, True]
14
    #after that primes(30) returns [2, 3, 5, 7, 11, 13, 17, 19, 23, 29], taken sieve is list of bools with len=20
15
    #print sieve - get a primary sieve
16
    if len(sieve)<n:
17
        a=len(sieve)
18
        sieve+=[True]*(n-a)
19
        i=1
20
        while i<=int((n-1)**0.5):
21
            i+=1
22
            if sieve[i]:
23
                for j in range(max(2*i, a+(i-a)%i),n,i): sieve[j]=False
24
    return [i for i in range(n) if sieve[i]]
25
from time import*
26
def factor(num):
27
    product=[]
28
    t=time()
29
    if num==0: return 0
30
    if num<0: return None
31
    for n in primes(int(num**0.5)+1):
32
        while num%n==0:
33
            if time()-t>2: break
34
            product.append(n)
35
            num/=n
36
        if num==1: break
37
    if num<>1: return product+[num]
38
    return product
39
def divisors(num):
40
    if num<1: return None
41
    def multiply(array):
42
        k=1
43
        for n in array: k*=n
44
        return k
45
    M=tuple(factor(num))
46
    return list(map(lambda x: multiply(x), product(*[[n**i for i in range(M.count(n)+1)] for n in set(M)])))
47
#POLINOMU OPERACIJOS, kur a+bx+cx^2+...=[a,b,c...]
48-
def interpolation(ordered_pairs):
48+
def divide(P1,P2):
49-
    n=len(ordered_pairs)
49+
    residue=[]
50-
    return solve_system([[m[0]**i for i in range(n)]+[m[1]] for m in ordered_pairs])
50+
    Pol=[n for n in P1]
51
    if len(P2)<2: return
52
    from time import*
53
    t=time()
54
    while len(P2)<=len(Pol):
55
        if time()-t>2:
56
            break
57
            print 'Serious error in polynomial division'
58
        coef=Pol[0]/float(P2[0])
59
        residue.append(coef)
60
        del(Pol[0])
61
        for i in range(len(P2)-1):
62
            Pol[i]-=coef*P2[i+1]
63
    if Pol==[0]: return residue
64
def roots(P):
65
    #mains bugs: finds just integral roots - need quotients (use substitution of polynomial variable to count this)
66
    #need general algorhytm for quadratic and cubic polynomials
67
    #possible liar result if roots not integral - abs(P[0]/P[-1]) is not always exact number
68
    if len(P)<2: return []
69
    guess=divisors(abs(P[0]/P[-1]))
70
    guess+=[-n for n in guess]
71
    p=True
72
    for test in guess:
73
        M=divide(P,[-test,1])
74
        if M<>None:
75
            p=False
76
            return [test]+roots(M)
77
    if p==True: return []
78
#MATRICU OPERACIJOS
79
def div(a,b): #norim, kad nebutu stufo, tokio kaip 0.0, o butu vietoj to sveikieji
80
        if float(a)/b==a/b: return a/b
81
        return float(a)/b
82
def gcd(*args): #rasime skaiciu DBD
83
        def gcd_of_two(a,b):
84
            a,b=abs(a),abs(b)
85
            while a>0:
86
                if b>a: a,b=b,a
87
                if a==1: return 1
88
                a-=(a/b)*b
89
            return b
90
        return reduce(lambda x,y: gcd_of_two(x,y), args)
91
def T(matrix): return zip(*matrix) #transponuota matrica
92
def mul(*args):
93
    def simple_mul(matrix_a, matrix_b):#matricu daugyba
94
        if len(matrix_a[0])==len(matrix_b): return map(lambda x: tuple(sum(starmap(lambda j,k: k*j, zip(x,n))) for n in zip(*matrix_b)),matrix_a)
95
    return reduce(lambda x,y: simple_mul(x,y), args)
96
def add(matrix_a, matrix_b, m1=1, m2=1):
97
    #sudetis m1*matrix_a+m2*matrix_b
98
    return [tuple(m1*matrix_a[i][j]+m2*matrix_b[i][j] for j in range(len(matrix_a[0]))) for i in range(len(matrix_a))]
99
def scalmul(matrix,const): return [tuple(matrix[i][j]*const for j in range(len(matrix[0]))) for i in range(len(matrix))] #matrica dauginta is skaliaro
100
def unit(length, diagonal=[]):
101
    if diagonal==[]: diagonal=[1]*length
102
    return [tuple(int(i==j)*diagonal[i] for i in range(length)) for j in range(length)] #vienetine matrica
103
def inverse(matrix): #atvirkstines matricos algoritmas
104
    #EXAMPLE: [[a1,a2,x1,x2],[b1,b2,y1,y2]] ->[[a3,a4,a5,a6],[x3,x4,x5,x6]], nes tai ekvivalentu zingsniui
105
    # /a1 a2\|/1 0 x1 x2\  _______\     /1 0\|/a3 a4 a5 a6\
106
    # \b1 b2/|\0 1 y1 y2/         /     \0 1/|\x3 x4 x5 x6/
107
    #rows can be longer than len of matrix; that case last terms means solutions to several systems of equations
108
    M=[[m for m in n] for n in matrix] #SIDE EFFECTS ARE ELIMINATED IN THIS ROW
109
    U=unit(len(M))
110
    order=[(i,j) for i in range(len(M)-1) for j in range(i+1,len(M))]+[(i,j) for i in range(len(M)-1,0,-1) for j in range(i-1,-1,-1)]
111
    #eiluciu keitimo tvarka, kai reikia suvesti i trikampi, po to i diagolini pavidala
112
    for i in range(len(M)):
113
        U[i]=U[i]+tuple(M[i][len(matrix):len(M[0])])
114
    #eiluciu sudejimo eiliskumas
115
    try:
116
        for n in order:
117
            i,j=n[0],n[1]
118
            if M[i][i]==0: # jei pasitaike nulinis elementas ant istrizaines, keiciame matricos eilutes. Jei neradome, ka sukeisti - nera atvirkstines matricos
119
                p=0
120
                for k in range(i+1,len(M)):
121
                    if M[k][i]<>0:
122
                        M[i],M[k]=M[k],M[i]
123
                        U[i],U[k]=U[k],U[i]
124
                        p=1
125
                        break
126
                if p==0: return [M[i] for i in range(len(M)) if M[i][i]<>0] #bug we'll use later
127
            if M[j][i]<>0: #suvedimo i trikampi pavidala vienas zingsnis, kai uznulinama kazkuri eilute
128
                c=M[i][i]*M[j][i]/gcd(M[i][i], M[j][i])
129
                c1,c2=c/M[i][i],c/M[j][i]
130
                M[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(M[j],M[i])))
131
                U[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(U[j],U[i])))
132
        return map(lambda i: tuple(map(lambda x: div(x, M[i][i]),U[i])), range(len(U)))
133
    except IndexError: print 'Nesiskaiciuoja inverse, per daug eiluciu matricoje'
134
    #pradine matrica pasikeite i matrica su skaiciais ant istrizaines; jos istrizaine normalizavus, graziname U reiksme
135
def rank(M): #gudras metodas, kur klaidos atveju inverse funkcijoje grazinamas tiesiskai nepriklausomu eiluciu matrica
136
    if len(M)>len(M[0]): return len(inverse(T(M))) #transponuojam, nes inverse taikosi tik siuo atveju
137
    else: return len(inverse(M))
138
def det(matrix): #determinanto algoritmas - padaryta pagal inverse funkcija
139
    D=[]
140
    M=[[m for m in n] for n in matrix]
141
    order=[(i,j) for i in range(len(M)-1) for j in range(i+1,len(M))]
142
    for n in order:
143
        i,j=n[0],n[1]
144
        if M[i][i]==0: # jei pasitaike nulinis elememtas ant istrizaines, keiciame matricos eilutes. Jei neradome, ka sukeisti - det=0
145
            p=0
146
            for k in range(i+1,len(M)):
147
                if M[k][i]<>0:
148
                    M[i],M[k]=M[k],M[i]
149
                    p=1
150
                    break
151
            if p==0: return 0
152
        if M[j][i]<>0: #suvedimo i trikampi pavidala vienas zingsnis, kai uznulinama kazkuri eilute
153
            c=M[i][i]*M[j][i]/gcd(M[i][i], M[j][i])
154
            c1,c2=c/M[i][i],c/M[j][i]
155
            D.append(-c2)
156
            M[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(M[j],M[i])))
157
    if D<>[]: D=reduce(lambda x,y:x*y,D)
158
    else: D=1
159
    return reduce(lambda x,y: x*y, map(lambda i: M[i][i], range(len(M))))/D
160
def solve_system(matrix): return [list(n) for n in zip(*inverse(matrix))[len(matrix):]]
161
def orthoappend(*args):
162
    #kas zingsni sprendziama sistema, kur like nezinomieji(kuriu kiekis yra defect) prilyginti 1
163
    M=[list(n) for n in list(args)]
164
    defect=len(args[0])-len(M)
165
    while defect>0:
166
        N=[M[i][:-defect]+[-sum(M[i][-defect:])] for i in range(len(M))]
167
        M+=[list(solve_system(N)[0])+[1]*defect]
168
        defect-=1
169
    return M
170
def chr_pol(M): #charakteristinis daugianaris, rastas aprepianciuju minoru pagalba
171
    def minors(M,k):
172
        L=len(M)
173
        #first row is elements of matrix with connection of order=[(M[n],M[n]) for n in second_row]
174
        #second_row says cyclic permutations with len=k of the sequence with len=L
175
        return [[[M[i][j] for j in n] for i in n] for n in [[x%L for x in range(i,i+k)] for i in range(L)]]
176
    return [sum([det(n) for n in minors(M,k)])*(-1)**k for k in range(len(M),0,-1)]+[1]
177
def eigenvalues(M): return roots(chr_pol(M))