Σημειωματάριο Πέμπτης 30 Απρ. 2015

class Polynomial

Φτιάχνουμε μια class Polynomial με σκοπό να παραστήσουμε στον υπολογιστή πολυώνυμα, δηλ. συναρτήσεις μιας μεταβλητής της μορφής \[ f(x) = a_n x^n + a_{n+1}x^{n-1}+\cdots+a_2x^2+a_1 x + a_0, \] όπου οι αριθμοί \(a_0, a_1, \ldots, a_n\) είναι σταθεροί. Μια τέτοια συνάρτηση λέμε ότι έχει βαθμό \(n\) αν ο συντελεστής \(a_n\) δεν είναι 0. Βαθμός ενός πολυωνύμου του \(x\) δηλ. είναι η μέγιστη δύναμη του \(x\) που εμφανίζεται σε αυτό.

__init__

Είναι φανερό ότι ένα πολυώνυμο είναι πλήρως καθορισμένο αν ξέρουμε τους συντελεστές του. Για να δημιουργήσουμε ένα πολυώνυμο λοιπόν περνάμε στη μέθοδο __init__ τη λίστα με τους συντελεστές του πολυωνύμου, με το μεγιστοβάθμιο συντελεστή να είναι ο πρώτος (από την αρχή) όρος της λίστας και με το σταθερό όρο του πολυωνύμου να είναι ο τελευταίος όρος της λίστας. Στην __init__ περνάμε αυτή τη λίστα με τους συντελεστές και αυτοί αποθηκεύονται εσωτερικά στο πεδίο coef του self.

Αυτό γίνεται με την ανάθεση

self.coef = coefs[:]

και όχι με την ανάθεση

self.coef = coefs

λόγω του ότι αν κάναμε το δεύτερο το αποτέλεσμα δε θα ήταν να αντιγραφεί η λίστα coefs στη λίστα self.coef (όπως θέλουμε) αλλά να δείχνει το όνομα self.coef στην ίδια λίστα στη μνήμη με το όνομα coefs. Αυτό θα είχε ως αποτέλεσμα οποιαδήποτε αλλαγή, αργότερα, στη λίστα coefs να αντανακλάται και στη λίστα self.coef πράγμα που, σχεδόν σίγουρα, δεν είναι αυτό που θέλουμε. Για να θυμηθείτε τα σχετικά ανατρέξτε στο Σημειωματάριο της 17ης Φεβ. 2015 (<<Πράξεις ανάθεσης σε μεταβλητές τύπου λίστας>>).

Επειδή όμως δε θέλουμε να έχουμε μηδενικό συντελεστή στο μεγιστοβάθμιο όρο του πολυωνύμου (ώστε να ξέρουμε ποιος είναι πραγματικά ο βαθμός του) διατρέχουμε τη λίστα από τα αριστερά προς τα δεξιά μέχρι να βρούμε το πρώτο μη μηδενικό στοιχείο και διαγράφουμε όσους μηδενικούς συντελεστές έχουμε βρει μέχρι τότε. Αυτό κάνουμε στο loop

for c in self.coef:
    if abs(c)<tol:
        self.coef.pop(0)
    else:
        break

Σχόλιο για το πώς ελέγχουμε αν δύο πραγματικοί αριθμοί είανι ίσοι.

Επειδή μιλάμε για πραγματικούς αριθμούς (στους οποίους το να είναι κάτι ακριβώς ίσο με μηδέν δεν έχει νόημα να το ελέγχει κανείς) ελέχγουμε αν η απόλυτη τιμή του συντελεστή είναι το πολύ tol, το οποίο το έχουμε ορίσει να είναι \(10^{-8}\). Αυτό το κάνουμε γενικά όποτε θέλουμε να ελέγξουμε μια ισότητα πραγματικών αριθμών (στους ακεραίους δεν υπάρχει αυτό το πρόβλημα γιατί αυτοί αναπαριστούνται ακριβώς εσωτερικά στον υπολογιστή): για να ελέγξουμε αν δύο αριθμοί είναι ίσοι ελέγχουμε αν η απόλυτη τιμή της διαφοράς τους είναι μικρότερη από κάποιον πολύ μικρό αριθμό. Το πόσο μικρός πρέπει να είναι αυτός ο αριθμός δεν μπορεί να καθοριστεί μονοσήμαντα και εξαρτάται από το ποιο πρόβλημα πάμε να λύσουμε αλλά μια αρχική τιμή \(10^{-8}\) είναι μια καλή αρχή. Αντίθετα, όταν ελέγχουμε ανισότητες δε χρειάζεται να κάνουμε κάτι τέτοιο.

deg

Η επόμενη μέθοδος που ορίζουμε είναι η deg η οποία μας επιστρέφει το βαθμό του πολυωνύμου.

cf

Η μέθοδος cf(i) μας επιστρέφει τον \(i\)-οστό συντελεστή του πολυωνύμου, ακόμη κι αν το \(i\) είναι μεγαλύτερο από το βαθμό, οπότε επιστρέφει 0. Αυτό μας είναι χρήσιμο παρακάτω όταν προσθέτουμε και πολλαπλασιάζουμε πολυώνυμα.

__call__

Η ειδική μέθοδος __call__ μας επιστρέφει την τιμή του πολυωνύμου στη θέση x. Καλείται με το ίδιο το όνομα του πολυωνύμου σαν αυτό να ήταν συνάρτηση:

p = Polynomial([2, 1, -1])
print p(3)

(Τυπώνει 20.0)

__str__

Η ειδική μέθοδος __str__ επιστρέφει ένα string με μια περιγραφή του πολυωνύμου.

p = Polynomial([2, 0, -1, -1])
print p

Τυπώνει: 2x**3-1x**1-1.

Πρέπει να προσέχουμε το πώς τυπώνουμε πρόσημα, μηδενικούς όρους, κλπ, ώστε το αποτέλεσμα να είναι αισθητικά ικανοποιητικό. Δείτε τα σχόλια μέσα στη μέθοδο για το πώς δουλεύει. Μπορεί φυσικά να γίνει και καλύτερη. Π.χ., όταν ο συντελεστής είναι 1 να μη τυπώνεται.

__add__

Η ειδική μέθοδος __add__ προσθέτει δύο πολυώνυμα μεταξύ τους και επιστρέφει το, επίσης πολυώνυμο, άθροισμά τους. Απλά δημιουργούμε μια λίστα με το άθροισμα των συντελεστών των δύο πολυωνύμων. Ο βαθμός του αθροίσματος είναι το πολύ ο μέγιστος των δύο βαθμών. Η __add__ καλείται ως τελεστής: απλά γράφουμε p+q αντί για p.__add__(q):

p = Polynomial([2, 1, 1])
q = Polynomial([3, 1, 1, 1])
print p+q

Τυπώνει: 3.0x**3+3x**2+2x**1+2

__mul__

Η ειδική μέθοδος __mul__ πολλαπλασιάζει μεταξύ τους δύο πολυώνυμα και επιστρέφει το, επίσης πολυώνυμο, γινόμενό τους. Ο βαθμός του γινομένου είναι το άθροισμα των βαθμών. Εύκολα βλέπει κανείς ότι αν \[ a(x) = a_m x^m + a_{m-1}*x^{m-1}+\cdots+a_1x+a_0, \] \[ b(x) = b_n x^n + b_{n-1}x^{n-1}+\cdots+b_1x+b_0, \] και \(c(x) = a(x)b(x)\), και το πολυώνυμο \(c(x)\) έχει συντελεστές \(c_j\) (\(0\le j \le m+n\)) τότε \[ c_0 = a_0 b_0, \] \[ c_1 = a_0 b_1 + a_1 b_0, \] \[ c_2 = a_0 b_2 + a_1 b_1 + a_2 b_0, \] και, γενικά, \[ c_i = a_0 b_i + a_1 b_{i-1}+ \cdots + a_{i-1} b_1 + a_i b_0. \] Αυτόν τον τελευταίο τύπο υλοποιούμε στο πρόγραμμά μας. Η ακολουθία \(c_n\) που ορίζεται με αυτό τον τρόπο ονομάζεται συνέλιξη (convolution) των δύο ακολουθιών \(a_i, b_i\).

roots

Η μέθοδος roots επιστρέφει μια λίστα με τις ρίζες του πολυωνύμου αν αυτό είναι βαθμού μέχρι 2, αλλιώς τυπώνει μήνυμα σφάλματος και επιστρέφει κενή λίστα.

Για να βρούμε τη λίστα ενεργούμε ανάλογα με το βαθμό.

Αν ο βαθμός είναι 2 και η διακρίνουσα είναι αρνητική τότε καλούμε τη συνάρτηση cmath.sqrt αντί για τη math.sqrt η οποία μας επιστρέφει μια (από τις δύο) μιγαδικές τετραγωνικές ρίζες του αριθμού. Η python μπορεί να χειριστεί μιγαδικούς αριθμούς. Ο μιγαδικός αριθμός \(i = \sqrt{-1}\) συμβολίζεται με 1j και ο μιγαδικός αριθμός με πραγματικό μέρος a και φανταστικό μέρος b συμβολίζεται με a+b*1j ή complex(a,b). Δείτε μια σύντομη περιγραφή των μιγαδικών στην python εδώ.

In [9]:
import cmath # Μαθηματικές συναρτήσεις που δουλεύουν και για μιγαδικούς αριθμούς
import math # Μαθηματικές συναρτήσεις
tol=1e-8

class Polynomial:
    
    def __init__(self, coefs):
        """Δημιουργεί ένα πολυώνυμο"""
        self.coef = coefs[:]
        for c in self.coef:
            if abs(c)<tol:
                self.coef.pop(0)
            else:
                break
    
    def deg(self):
        """Επιστρέφει το βαθμό του πολυωνύμου"""
        return len(self.coef)-1
    
    def cf(self, i):
        """Επιστρέφει τον i-οστό συντελεστή, ακόμη και για i>βαθμού"""
        d = self.deg()
        if i<=d:
            return self.coef[d-i]
        else:
            return 0.

    def __call__(self, x):
        """Επιστρέφει την τιμή του πολυωνύμου στη θέση x"""
        d = self.deg()
        v=0.
        for i in range(d+1):
            v += self.coef[i]*x**(d-i)
        return v

    def __str__(self):
        """Επιστρέφει περιγραφή του πολυωνύνου σε μορφή κειμένου"""
        d = self.deg()
        
        s = "{c}x**{ex}".format(c=self.coef[0], ex=d) # Τυπώνει το μεγιστοβάθμιο όρο
        for i in range(1, d): # Τυπώνει τους υπόλοιπους όρους εκτός το σταθερό
            sign = '+' if self.coef[i]>=0 else '' # Χρειάζεται να τυπωθεί + ή - πριν το μονώνυμο
            if abs(self.coef[i])>tol: # αν ο συντελεστής είναι μηδενικός δεν τυπώνουμε το μονώνυμο
                s = s + "{s}{c}x**{ex}".format(s=sign, c=self.coef[i], ex=d-i)
        sign = '+' if self.coef[len(self.coef)-1]>=0 else '' # Τυπώνουμε το σταθερό όρο
        if abs(self.coef[len(self.coef)-1])>tol:
            ss = "{s}{c}".format(s=sign, c=self.coef[len(self.coef)-1])
            s = s+ss
        return s


    def __add__(self, p):
        """Προσθέτουμε δύο πολυώνυμα και επιστρέφουμε το άθροισμά τους"""
        d = max(self.deg(), p.deg())
        c=[]
        for i in range(d+1):
            c.append(self.cf(d-i) + p.cf(d-i))
        return Polynomial(c)
    
    def __mul__(self, p):
        """Πολλαπλασιάζουμε δύο πολυώνυμα και επιστρέφουμε το γινόμενό τους"""
        d = self.deg()+p.deg()
        c=[]
        for i in range(d+1):
            k = d-i
            s = 0
            for j in range(k+1):
                s += self.cf(j)*p.cf(k-j)
            c.append(s)
        return Polynomial(c)
    
    def roots(self):
        """Αν το πολυώνυμο είναι βαθμού μέχρι 2 τότε επιστρέφουμε τις ρίζες του"""
        if self.deg()==0:
            return []
        if self.deg()==1:
            return -self.cf(0)/(0.+self.cf(1))
        if self.deg()>2:
            print "Error"
            return []
        a = self.cf(2)+0.0
        b = self.cf(1)+0.0
        c = self.cf(0)+0.0
        D = b*b-4.*a*c
        sq = cmath.sqrt(D) if D<0 else math.sqrt(D)
        return [(-b+sq)/(2.*a), (-b-sq)/(2.*a)]
        
        


print Polynomial([1,0,1]).roots()
print Polynomial([-2, 3]).roots()

    
[1j, (-0-1j)]
1.5