Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- from itertools import* #reikalinga funkcijai starmap
- from math import* # to get pi
- from numbers import* #to solve equations with complex numbers
- #SKAICIU TEORIJA
- def primes(n, sieve=[False,False,True]): #by-need calculation of primes; not includes n
- #TEST1
- #sieve: [1,1,1,1,1,1,1,1,1,1]=[1]*10
- #possible primes [0,1,2,3,4,5,6,7,8,9]
- #given sieve [0,0,1,1,0]
- #a=3
- #check from 0 to sqrt(len(sieve))
- #last five entries were checked dividing by primes 0 to a=3
- #sieve is given as list of bools to save memory
- #TEST2
- #primes(20) returns [2, 3, 5, 7, 11, 13, 17, 19], sieve is [False, False, True]
- #after that primes(30) returns [2, 3, 5, 7, 11, 13, 17, 19, 23, 29], taken sieve is list of bools with len=20
- #print sieve - get a primary sieve
- if len(sieve)<n:
- a=len(sieve)
- sieve+=[True]*(n-a)
- i=1
- while i<=int((n-1)**0.5):
- i+=1
- if sieve[i]:
- for j in range(max(2*i, a+(i-a)%i),n,i): sieve[j]=False
- return [i for i in range(n) if sieve[i]]
- from time import*
- def factor(num):
- product=[]
- t=time()
- if num==0: return 0
- if num<0: return None
- for n in primes(int(num**0.5)+1):
- while num%n==0:
- if time()-t>2: break
- product.append(n)
- num/=n
- if num==1: break
- if num<>1: return product+[num]
- return product
- def divisors(num):
- if num<1: return None
- def multiply(array):
- k=1
- for n in array: k*=n
- return k
- M=tuple(factor(num))
- return list(map(lambda x: multiply(x), product(*[[n**i for i in range(M.count(n)+1)] for n in set(M)])))
- #POLINOMU OPERACIJOS, kur a+bx+cx^2+...=[a,b,c...]
- def divide(P1,P2):
- residue=[]
- Pol=[n for n in P1]
- if len(P2)<2: return
- from time import*
- t=time()
- while len(P2)<=len(Pol):
- if time()-t>2:
- break
- print 'Serious error in polynomial division'
- coef=Pol[0]/float(P2[0])
- residue.append(coef)
- del(Pol[0])
- for i in range(len(P2)-1):
- Pol[i]-=coef*P2[i+1]
- if Pol==[0]: return residue
- def roots(pol):
- pol=[float(n) for n in pol]
- deg=max([i for i in range(len(pol)) if pol[i]<>0])
- if deg==0: return None
- elif deg==1: return [-pol[0]/pol[1]]
- elif deg==2:
- a,b,c=pol[2],pol[1],pol[0]
- D=b*b-4*a*c
- if D<0: return [(-b-(0+1j)*(-D)**0.5)/(2*a), (-b+(0+1j)*(-D)**0.5)/(2*a)]
- else: return [(-b-D**0.5)/(2*a), (-b+D**0.5)/(2*a)]
- #Given ax^3+bx^2+cx+d=0. Reduce to x^3+ax^2+bx+c=0. Substitute x=y-a/3.
- #Reduce to y^3+yp+q=0. Substitute y=z-p/3z. Got z^6+qz^3-p^3/27=0
- elif deg==3: #kinda mystics why half of solutions can be omitted
- a,b,c=pol[2]/pol[3],pol[1]/pol[3],pol[0]/pol[3]
- p,q=b-a*a/3, c+2*a*a*a/27-a*b/3
- unity=[cos(2*pi*i/3)+sin(2*pi*i/3)*(0+1j) for i in range(3)]
- solve_z=(roots([-p*p*p/27, q, 1])[0]+0j)**(1./3) #solved for z
- zs=[solve_z*unity[i] for i in range(3)]
- ys=[z-p/(3*z) for z in zs]
- return [y-a/3 for y in ys]
- elif deg==4: #Reduce to y^4=py^2+qy+r
- a,b,c,d=pol[3]/pol[4],pol[2]/pol[4],pol[1]/pol[4], pol[0]/pol[4]
- p,q,r=3*a*a/8-b, a*b/2-a*a*a/8-c,3*a**4/256-a*a*b/16+a*c/4-d
- #print p,q,r, [4*p*r-q*q,8*r,4*p,8]
- D=roots([4*p*r-q*q,8*r,4*p,8]) #solving cubics for a constant :)
- if D[0].imag<D[1].imag: t=D[0].real
- else: t=D[1].real #choosing real root
- pre=((2*t+p+0j)**0.5).real
- return [x-a/4 for x in roots([t+pre*q/(2*(2*t+p)),pre,1])+roots([t-pre*q/(2*(2*t+p)),-pre,1])]
- else:
- print "I can't solve for big degrees cuz I don't know how to approximate complex roots efficiently"
- return None
- #MATRICU OPERACIJOS
- def div(a,b): #norim, kad nebutu stufo, tokio kaip 0.0, o butu vietoj to sveikieji
- if float(a)/b==a/b: return a/b
- return float(a)/b
- def gcd(*args): #rasime skaiciu DBD
- def gcd_of_two(a,b):
- a,b=abs(a),abs(b)
- while a>0:
- if b>a: a,b=b,a
- if a==1: return 1
- a-=(a/b)*b
- return b
- return reduce(lambda x,y: gcd_of_two(x,y), args)
- def T(matrix): return zip(*matrix) #transponuota matrica
- def mul(*args):
- def simple_mul(matrix_a, matrix_b):#matricu daugyba
- 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)
- return reduce(lambda x,y: simple_mul(x,y), args)
- def add(matrix_a, matrix_b, m1=1, m2=1):
- #sudetis m1*matrix_a+m2*matrix_b
- 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))]
- 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
- def unit(length, diagonal=[]):
- if diagonal==[]: diagonal=[1]*length
- return [tuple(int(i==j)*diagonal[i] for i in range(length)) for j in range(length)] #vienetine matrica
- def inverse(matrix): #atvirkstines matricos algoritmas
- #EXAMPLE: [[a1,a2,x1,x2],[b1,b2,y1,y2]] ->[[a3,a4,a5,a6],[x3,x4,x5,x6]], nes tai ekvivalentu zingsniui
- # /a1 a2\|/1 0 x1 x2\ _______\ /1 0\|/a3 a4 a5 a6\
- # \b1 b2/|\0 1 y1 y2/ / \0 1/|\x3 x4 x5 x6/
- #rows can be longer than len of matrix; that case last terms means solutions to several systems of equations
- M=[[m for m in n] for n in matrix] #SIDE EFFECTS ARE ELIMINATED IN THIS ROW
- U=unit(len(M))
- 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)]
- #eiluciu keitimo tvarka, kai reikia suvesti i trikampi, po to i diagolini pavidala
- for i in range(len(M)):
- U[i]=U[i]+tuple(M[i][len(matrix):len(M[0])])
- #eiluciu sudejimo eiliskumas
- try:
- for n in order:
- i,j=n[0],n[1]
- if M[i][i]==0: # jei pasitaike nulinis elementas ant istrizaines, keiciame matricos eilutes. Jei neradome, ka sukeisti - nera atvirkstines matricos
- p=0
- for k in range(i+1,len(M)):
- if M[k][i]<>0:
- M[i],M[k]=M[k],M[i]
- U[i],U[k]=U[k],U[i]
- p=1
- break
- if p==0: return [M[i] for i in range(len(M)) if M[i][i]<>0] #bug we'll use later
- if M[j][i]<>0: #suvedimo i trikampi pavidala vienas zingsnis, kai uznulinama kazkuri eilute
- c=M[i][i]*M[j][i]/gcd(M[i][i], M[j][i])
- c1,c2=c/M[i][i],c/M[j][i]
- M[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(M[j],M[i])))
- U[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(U[j],U[i])))
- return map(lambda i: tuple(map(lambda x: div(x, M[i][i]),U[i])), range(len(U)))
- except IndexError: print 'Nesiskaiciuoja inverse, per daug eiluciu matricoje'
- #pradine matrica pasikeite i matrica su skaiciais ant istrizaines; jos istrizaine normalizavus, graziname U reiksme
- def rank(M): #gudras metodas, kur klaidos atveju inverse funkcijoje grazinamas tiesiskai nepriklausomu eiluciu matrica
- if len(M)>len(M[0]): return len(inverse(T(M))) #transponuojam, nes inverse taikosi tik siuo atveju
- else: return len(inverse(M))
- def det(matrix): #determinanto algoritmas - padaryta pagal inverse funkcija
- D=[]
- M=[[m for m in n] for n in matrix]
- order=[(i,j) for i in range(len(M)-1) for j in range(i+1,len(M))]
- for n in order:
- i,j=n[0],n[1]
- if M[i][i]==0: # jei pasitaike nulinis elememtas ant istrizaines, keiciame matricos eilutes. Jei neradome, ka sukeisti - det=0
- p=0
- for k in range(i+1,len(M)):
- if M[k][i]<>0:
- M[i],M[k]=M[k],M[i]
- p=1
- break
- if p==0: return 0
- if M[j][i]<>0: #suvedimo i trikampi pavidala vienas zingsnis, kai uznulinama kazkuri eilute
- c=M[i][i]*M[j][i]/gcd(M[i][i], M[j][i])
- c1,c2=c/M[i][i],c/M[j][i]
- D.append(-c2)
- M[j]=tuple(map(lambda x: x[1]*c1-x[0]*c2, zip(M[j],M[i])))
- if D<>[]: D=reduce(lambda x,y:x*y,D)
- else: D=1
- return reduce(lambda x,y: x*y, map(lambda i: M[i][i], range(len(M))))/D
- def solve_system(matrix): return [list(n) for n in zip(*inverse(matrix))[len(matrix):]]
- def chr_pol(M): #charakteristinis daugianaris, rastas aprepianciuju minoru pagalba
- def minors(M,k):
- L=len(M)
- #first row is elements of matrix with connection of order=[(M[n],M[n]) for n in second_row]
- #second_row says cyclic permutations with len=k of the sequence with len=L
- 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)]]
- return [sum([det(n) for n in minors(M,k)])*(-1)**k for k in range(len(M),0,-1)]+[1]
- def eigenvalues(M): return roots(chr_pol(M))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement