1818-Problem 1 (Dual Numbers).sagews

Problem 1 (Dual Numbers)

The square of a small number is very small indeed; e.g., squaring $0.0001$ we get $0.00000001$. Let $\epsilon$ be an idealized small quantity with the property that $\epsilon\neq 0$ and $\epsilon^2=0$. "Numbers" of the form $a+b\epsilon$ with $a$ and $b$ real are called dual numbers. (Compare with the complex numbers, which are of the form $a+bi$ with $i^2=-1$). Dual numbers obey the rules of ordinary algebra, subject to the condition that $\epsilon^2 = 0$; for example, $$(a+b\epsilon) + (c+d\epsilon) = (a+c) + (b+d)\epsilon$$

a. Make a class for calculating with dual numbers. A skeleton is provided below: where it says "pass" you fill in the appropriate definition. Show examples of how each of the overloaded operator works. (Optional: Provide a _latex_ function, if you know LaTeX.)

class DualNum(SageObject):

    
    def __init__(self, a, b=0):
        self.a = RR(a)
        self.b = RR(b)

    def _repr_(self):
        return "DualNum(%s,%s)" % (self.a, self.b)

    def __neg__(self):
        pass

    def __add__(self, other):
        pass

    def __sub__(self, other):
        pass
        
    def __mul__(self, other):
        pass

    def __pow__(self, n):
        pass

 

class DualNum(SageObject):
    
    def __init__(self, a, b=0):
        self.a = a
        self.b = b
    
    def _latex_(self):
        return str(self.a)+"+"+str(self.b)+"\\epsilon"
        
        #"\\frac{"+str(self.num)+"}{"+str(self.den)+"}"
    
    def _repr_(self):
        return "DualNum(%s,%s)" % (self.a, self.b)
        
    def __neg__(self):
        return DualNum(-1*self.a,-1*self.b)
        
    def __add__(self, other):
        if(other in RR):
            other=DualNum(other)
        return DualNum(self.a+other.a,self.b+other.b)
    
    def __radd__(self, other):
        if(other in RR):
            other=DualNum(other)
        return DualNum(self.a+other.a,self.b+other.b)
        
    def __sub__(self, other):
        if(other in RR):
            other=DualNum(other)
        return DualNum(self.a-other.a,self.b-other.b)
        
    def __rsub__(self, other):
        if(other in RR):
            other=DualNum(other)
        return DualNum(self.a-other.a,self.b-other.b)
    
    def __mul__(self, other):
        if(other in RR):
            other=DualNum(other)
        newalpha=self.a*other.a
        newbeta=self.a*other.b+self.b*other.a
        return DualNum(newalpha,newbeta)

    def __rmul__(self, other):
        if(other in RR):
            other=DualNum(other)
        newalpha=self.a*other.a
        newbeta=self.a*other.b+self.b*other.a
        return DualNum(newalpha,newbeta)
        
    def __div__(self, other):
        if(other in RR):
            other=DualNum(other)
        newalpha=self.a/other.a
        newbeta=(self.b*other.a-self.a*other.b)/power(other.a,2)
        return DualNum(newalpha,newbeta)
    
    def __rdiv__(self,other):
        return DualNum(other/self.a, other/self.b)
        
    def __invert__(self):
        return DualNum(1/self.a,-self.b/power(self.a,2))
        
    def __lt__(self,other):
        if(self.a < other.a):
            return True
        if(self.a == other.a):
            if(self.b<other.b):
                return True
        return false
        
    def __gt__(self,other):
        if(self.a > other.a):
            return True
        if(self.a == other.a):
            if(self.b>other.b):
                return True
        return false
        
    def __le__(self,other):
        if(self.a < other.a):
            return True
        if(self.a==other.a):
            if(self.b<=other.b):
                return True
        return false
        
    def __ge__(self,other):
        if(self.a > other.a):
            return True
        if(self.a==other.a):
            if(self.b>=other.b):
                return True
        return false
            
    def __eq__(self,other):
        if(self.a == other.a):
            if(self.b == other.b):
                return True
        return false

    def __pow__(self, n):
        if(n==0):
            return DualNum(1,0)
        raised=1
        for i in range(n):
            raised*=self
        return DualNum(raised.a,raised.b)

pretty_print(3/2)

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{3}{2}

f=DualNum(-2,6)
g=DualNum(1,5)
print f
print g

DualNum(-2,6) DualNum(1,5)

pretty_print(f+g)
print(g+f)

\newcommand{\Bold}[1]{\mathbf{#1}}-1+11\epsilon

DualNum(-1,11)

f-g

DualNum(2.00000000000000,2.00000000000000)

f*g

DualNum(-2,-4)

f/g

DualNum(-2,16)

power(f,2)

DualNum(4,-24)

~f

halló DualNum(-1/2,-3/2)

3/f

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_12.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("My9m"),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpOGaHgC/___code___.py", line 3, in <module> exec compile(ur'_sage_const_3 /f' + '\n', '', 'single') File "", line 1, in <module> File "element.pyx", line 1491, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11901) File "coerce.pyx", line 707, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6079) File "coerce.pyx", line 1178, in sage.structure.coerce.CoercionModel_cache_maps.get_action (sage/structure/coerce.c:10901) File "coerce.pyx", line 1344, in sage.structure.coerce.CoercionModel_cache_maps.discover_action (sage/structure/coerce.c:12577) File "ring.pyx", line 2618, in sage.rings.ring.is_Ring (sage/rings/ring.c:13901) File "/usr/local/sage2/local/lib/python2.6/site-packages/sage/categories/category.py", line 472, in __contains__ c = x.category() TypeError: descriptor 'category' of 'sage.structure.sage_object.SageObject' object needs an argument

b. Make sure that you can add, multiply and substract reals from a dual number correctly. (Simply amend the code above, but provide examples of usage in all cases here below.)

print f+2
print 2+f

DualNum(5.00000000000000,4.00000000000000) DualNum(5.00000000000000,4.00000000000000)

print f-2
print 2-f

DualNum(1.00000000000000,4.00000000000000) DualNum(1.00000000000000,4.00000000000000)

print f*3
print 3*f

DualNum(9.00000000000000,12.0000000000000) DualNum(9.00000000000000,12.0000000000000)

c. Assume that $c\neq 0$. What dual number does $(a+b\epsilon)/(c+d\epsilon)$ evaluate to? Add the definitions of __div__ and __invert__ to the class DualNum that reflects what you found.

Answer:

Division.

$$\frac{a+b \epsilon}{c+d \epsilon}=\frac{(a+b \epsilon)(c-d \epsilon)}{(c+d \epsilon)(c-d \epsilon)}=\frac{ac - ad \epsilon + bc \epsilon - bd \epsilon^2}{c^2-d^2\epsilon^2} = \frac{ac+\epsilon(bc - ad)}{c^2} =  \frac{a}{c}+\frac{bc - ad}{c^2}\epsilon$$

print f/g
print f/3
#print 3/f # works fine on sage.ru.is

DualNum(-2,16) DualNum(-2/3,2)

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_11.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("cHJpbnQgZi9nCnByaW50IGYvMwpwcmludCAzL2Y="),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpmdcRnb/___code___.py", line 5, in <module> exec compile(ur'print _sage_const_3 /f' + '\n', '', 'single') File "", line 1, in <module> File "element.pyx", line 1491, in sage.structure.element.RingElement.__div__ (sage/structure/element.c:11901) File "coerce.pyx", line 707, in sage.structure.coerce.CoercionModel_cache_maps.bin_op (sage/structure/coerce.c:6079) File "coerce.pyx", line 1178, in sage.structure.coerce.CoercionModel_cache_maps.get_action (sage/structure/coerce.c:10901) File "coerce.pyx", line 1344, in sage.structure.coerce.CoercionModel_cache_maps.discover_action (sage/structure/coerce.c:12577) File "ring.pyx", line 2618, in sage.rings.ring.is_Ring (sage/rings/ring.c:13901) File "/usr/local/sage2/local/lib/python2.6/site-packages/sage/categories/category.py", line 472, in __contains__ c = x.category() TypeError: descriptor 'category' of 'sage.structure.sage_object.SageObject' object needs an argument

Invert. We know that the inverse of $a+b\epsilon$ is $\frac{1}{a+b \epsilon}$ because $a+b\epsilon \cdot \frac{1}{a+b \epsilon}=1$. Let´s simplify the inverse.

$$\frac{1}{a+b \epsilon}\frac{a-b \epsilon}{a-b \epsilon} = \frac{a-b\epsilon}{a^2-b^2\epsilon^2}=\frac{1}{a}-\frac{b}{a^2}\epsilon$$

~f

DualNum(0.333333333333333,0.444444444444444)

d. Find a meaningful way of overloading the comparison operators <,<=,==, etc and implement it.

print f<g
print g<f
print g<g

False True False

print f<=g
print g<=f
print g<=g
print g>g
print g>=g

True False True False True

print f == g
print f == f

False True

e. Solve the second degree equation $x^2+x+\epsilon=0$ in dual numbers.

Answer:

First let´s insert an arbitrary dual number in for $x$.

$$(a+b\epsilon)^2+(a+b\epsilon)+\epsilon=0$$

$$a^2+2ab\epsilon+b^2\epsilon^2+a+b\epsilon+\epsilon=0$$

We know that $\epsilon^2=0$...

$$(a^2+a)+(2ab+b+1)\epsilon=0$$

Now we want our dual number to be equal to zero the real number has to be zero and the factor in front of $\epsilon$. Let´s first find out the possible values of $a$.

$$a^2+a=0$$

$$a(a+1)=0$$

So we see that either $a=0$ or $a=-1$. Now let´s find out the possible values of $b$ given these values of $a$. We have the equation:

$$2ab+b+1=0$$

We will first use $a=0$ and obtain.

$$b+1=0$$

$$b=-1$$.

So one of our dual numbers is $0-\epsilon$ or simply $-\epsilon$. Now we will use $a=-1$:

$$-2b+b+1=0$$

$$-b+1=0$$

$$b=1$$

There we have another solution which is $-1+\epsilon$. It´s easy to verify the solutions, the first solution

$$-\epsilon^2-\epsilon+\epsilon=0$$

$$0+(-1+1)\epsilon=0$$

$$0=0$$

And the second solution

$$(-1+\epsilon)^2+(-1+\epsilon)+\epsilon=0$$

$$1-2\epsilon+\epsilon^2-1+\epsilon+\epsilon=0$$

$$(1+0-1)+(-2+1+1)\epsilon=0$$

$$0=0$$

f. (Bonus) Write a function that solves the second degree equation $Ax^2+Bx+C=0$ where $A,B,C$ are dual numbers. The function needs to find or show all dual numbers $x$ that fulfil a given equation.

Problem 2. (Quaternions)

The quaternions are a generalization of the complex numbers. They have the form $$a+bi+cj+dk$$ where $a,b,c,d$ are real numbers and $i,j,k$ fulfill the following equation:

$$i^2=j^2+k^2=ijk=-1.$$

Note that here the order matters in multiplication: $ij=k$ but $ji=-k.

You can read more on Quaterions on Mathworld and Wikipedia.

Write a class for Quaternions and implement all the meaningful operations. Give good examples of how each operator is used.

def prstr(i):
    if(i<0):
        return str(i)
    else:
        return "+"+str(i)

class Quaternions(SageObject):
    
    def __init__(self, a=0, b=0,c=0,d=0):
        self.a = a
        self.b = b
        self.c = c
        self.d = d
                
    def _repr_(self):
        return "Quaternions(%s,%s,%s,%s)" % (self.a, self.b, self.c, self.d)
        
    def _latex_(self):
        return str(self.a)+prstr(self.b)+"i"+prstr(self.c)+"j"+prstr(self.d)+"k"
        
    def __neg__(self):
        return Quaternions(-self.a,-self.b,-self.c,-self.d)
        
    def __add__(self, other):
        if(other in RR):
            other=Quaternions(other)
        return Quaternions(self.a+other.a,self.b+other.b,self.c+other.c,self.d+other.d)
        
    def __radd__(self, other):
        if(other in RR):
            other=Quaternions(other)
        return Quaternions(self.a+other.a,self.b+other.b,self.c+other.c,self.d+other.d)
        
    def __sub__(self, other):
        if(other in RR):
            other=Quaternions(other)
        return Quaternions(self.a-other.a,self.b-other.b,self.c-other.c,self.d-other.d)
        
    def __rsub__(self, other):
        if(other in RR):
            other=Quaternions(other)
        return Quaternions(other.a-self.a,other.b-self.b,other.c-self.c,other.d-self.d)
            
    def __mul__(self, other):
        if(other in RR):
            other=Quaternions(other)
        a=self.a*other.a-self.b*other.b-self.c*other.c-self.d*other.d
        b=self.a*other.b+self.b*other.a+self.c*other.d-self.d*other.c
        c=self.a*other.c-self.b*other.d+self.c*other.a+self.d*other.b
        d=self.a*other.d+self.b*other.c-self.c*other.b+self.d*other.a
        return Quaternions(a,b,c,d)
    
    def __rmul__(self, other):
        if(other in RR):
            other=Quaternions(other)
        a=self.a*other.a-self.b*other.b-self.c*other.c-self.d*other.d
        b=self.a*other.b+self.b*other.a+self.d*other.c-self.c*other.d
        c=self.a*other.c-self.d*other.b+self.c*other.a+self.b*other.d
        d=self.a*other.d+self.c*other.b-self.b*other.c+self.d*other.a
        return Quaternions(a,b,c,d)
        
    def __div__(self,other):
        if(other in RR):
            other=Quaternions(other)
        conj_other=Quaternions(other.a,-other.b,-other.c,-other.d)
        numer=conj_other*self
        denom=other.a*other.a+other.b*other.b+other.c*other.c+other.d*other.d
        return Quaternions(numer.a/denom, numer.b/denom, numer.c/denom, numer.d/denom)
        
    def __rdiv__(self,other):
        if(other in RR):
            other=Quaternions(other)
        conj_self=Quaternions(self.a,-self.b,-self.c,-other.d)
        numer=conj_self*other
        denom=self.a*self.a+self.b*self.b+self.c*self.c+self.d*self.d
        return Quaternions(numer.a/denom, numer.b/denom, numer.c/denom, numer.d/denom)
    
    def __invert__(self):
        conj_self=Quaternions(self.a,-self.b,-self.c,-self.d)
        length=sqrt(self.a*self.a+self.b*self.b+self.c*self.c+self.d*self.d)
        return Quaternions(conj_self.a/(length*length),conj_self.b/(length*length),conj_self.c/(length*length),conj_self.d/(length*length))
        
    def __lt__(self,other):
        return "This operation doesn´t make sense."
        
    def __gt__(self,other):
        return "This operation doesn´t make sense."
        
    def __le__(self,other):
        return "This operation doesn´t make sense."
        
    def __ge__(self,other):
        return "This operation doesn´t make sense."                
    
    def __eq__(self,other):
        if(self.a == other.a):
            if(self.b == other.b):
                if(self.c == other.c):
                    if(self.d == other.d):
                        return True
        return false
    
    def __pow__(self, n):
        if(n==0):
            return Quaternions(1,0,0,0)
        raised=1
        for i in range(n):
            raised*=self
        return Quaternions(raised.a,raised.b,raised.c,raised.d)

f=Quaternions(-2,5,-7,0)
h=Quaternions(-2,-5,7,0) # conjugate of f
g=Quaternions(1,2,3,-4)
pretty_print(f) # it does print 0k, maybe better if it was -2+5i-7j
pretty_print(g)

\newcommand{\Bold}[1]{\mathbf{#1}}-2+5i-7j+0k

\newcommand{\Bold}[1]{\mathbf{#1}}1+2i+3j-4k

print f+g
print 3+g
print g+3
print 3-g
print g-3
print f-g
print f*g

Quaternions(-1,7,-4,-4) Quaternions(4,2,3,-4) Quaternions(4,2,3,-4) Quaternions(2,-2,-3,4) Quaternions(-2,2,3,-4) Quaternions(-3,3,-10,4) Quaternions(9,29,7,37)

print g*f
print 4*f
print f*4
print ~f*f
print f^3
print f/g
#print (f/5)*5  # works fine on sage.ru.is
#print (5/f)*f  # works fine on sage.ru.is
print f+h
print f==g
print f<g #I think it is important to overload all the comparison operators even though most of them do not make sense, otherwise we might be using other comparison operators and get some strange result
print f>g
print f>=g
print f==f

Quaternions(9,-27,-33,-21) Quaternions(-8,20,-28,0) Quaternions(-8,20,-28,0) Quaternions(1,0,0,0) Quaternions(436,-310,434,0) Quaternions(-13/30,37/30,19/30,7/10) Quaternions(-4,0,0,0) False This operation doesn´t make sense. This operation doesn´t make sense. This operation doesn´t make sense. True

Problem 3 (Problem 79 on Project Euler).

A common security method used for online banking is to ask the user for three random characters from a passcode. For example, if the passcode was 531278, they may asked for the 2nd, 3rd, and 5th characters; the expected reply would be: 317.

The text file, keylog.txt, contains fifty successful login attempts.

Given that the three characters are always asked for in order, analyse the file so as to determine the shortest possible secret passcode of unknown length.

v=[[3,1,9],[6,8,0],[1,8,0],[6,9,0],[1,2,9],[6,2,0],[7,6,2],[6,8,9],[7,6,2],[3,1,8],[3,6,8],[7,1,0],[7,2,0],[7,1,0],[6,2,9],[1,6,8],[1,6,0],[6,8,9],[7,1,6],[7,3,1],[7,3,6],[7,2,9],[3,1,6],[7,2,9],[7,2,9],[7,1,0],[7,6,9],[2,9,0],[7,1,9],[6,8,0],[3,1,8],[3,8,9],[1,6,2],[2,8,9],[1,6,2],[7,1,8],[7,2,9],[3,1,9],[7,9,0],[6,8,0],[8,9,0],[3,6,2],[3,1,9],[7,6,0],[3,1,6],[7,2,9],[3,8,0],[3,1,9],[7,2,8],[7,1,6]]

# The algorithm is pretty straight forward.
# It goes through all the samples, finds some number that does not have any number to the
# left of it, then it goes through all the samples and deletes those that have the number
# to the left. Repeat until no items left in v.

def find_code(v):
    code=[]
    while(any_items(v)):
        for i in range(len(v)):
            if(len(v[i])>0):
                if(leftmost(v[i][0],v)):
                    code.append(v[i][0])
                    del_left(v[i][0],v)
                    break
    return code

#finds an item that does not have anything on its left side
def leftmost(k,v):
    for i in range(len(v)):
        if(len(v[i])>0):
            if(v[i][0]!=k):
                for j in range(1,len(v[i])):
                    if(v[i][j]==k):
                        return false
    return True

# deletes the leftmost item in v[i] if it is k
def del_left(k,v):
    for i in range(len(v)):
        if(len(v[i])>0):
            if(v[i][0]==k):
                del(v[i][0])
                
def any_items(v): #checks if there are any items in v
    for i in range(len(v)):
        if(len(v[i])>0):
            return True
    return false

find_code(v)

[7, 3, 1, 6, 2, 8, 9, 0]

Problem 4 (Problem 71 of project Euler)

Consider the fraction, n/d, where n and d are positive integers. If nd and HCF(n,d)=1, it is called a reduced proper fraction.

If we list the set of reduced proper fractions for d 8 in ascending order of size, we get:

1/8, 1/7, 1/6, 1/5, 1/4, 2/7, 1/3, 3/8, 2/5, 3/7, 1/2, 4/7, 3/5, 5/8, 2/3, 5/7, 3/4, 4/5, 5/6, 6/7, 7/8

It can be seen that 2/5 is the fraction immediately to the left of 3/7.

By listing the set of reduced proper fractions for d 1,000,000 in ascending order of size, find the numerator of the fraction immediately to the left of 3/7.

Answer:

The function loops through all denominators from $1$ to $n$ and checks just the numerator that gives us the value that is to the left of $\frac{3}{7}$. The only thing that really needs explanation is how we obtain the numerator. Well we know that our number has to be

$$\frac{n}{d}<\frac{3}{7}$$

$$n < \frac{3}{7}\cdot d$$

We want the numerator to be less than $\frac{3}{7}$ so we take the floor of the right side:

$$n< \lfloor \frac{3}{7} \cdot d\rfloor$$

Then we can be sure that we have the numerator of the fraction that is next to the left from $\frac{3}{7}$ given the denominator $d$.

def proper(n):
    prox=1
    num=0
    den=0
    for d in [1..n]:
        n = floor( (3/7) * d )
        curr_prox=3/7-n/d
        if(curr_prox > 0):
            if(curr_prox<prox):
                prox=curr_prox
                num=n
                den=d
    return num/den

As we can see below the numerator we are looking for is $428570$.

time proper(1000000)

428570/999997 Time: CPU 8.42 s, Wall: 8.42 s