Programming Forums > Script Programming > Python

UnKnown X · Oct 10th, 2007, 11:18 AM

I have a project where I am to write a programme that illustrates how computers can be used for solving iterative tasks, using Euclid's algorithm as an example, and presenting and explaining it to a class in which most students have little to no knowledge about programming.

Thus, I need to keep it as simple as possible, and I chose to use a high-level language like Python.

Euclid's algorithm went smoothly, but I was additionally requested to add a function for solving linear diophantine equations with two variables and one equation, and I'm not even sure how to begin tackling this. I'm not asking for a solution at all, since this is my task; I only need a way to get started with this.

For example, given the equation 75x + 42y = 3, the greatest common divisor is 3.
33 = 75 - 1 * 42
9 = 42 - 1 * 33
6 = 33 - 3 * 9
3 = 9 - 1 * 6
Having stored all the numbers in the process (in the form of [[a,b,c,d],[a,b,c,d],[a,b,c,d],...]), how can I work my way 'back up'? Do I have to treat the equation as a string and keep replacing on my way up, or is there a simpler way?

Here's the function I wrote to find gdc(a,b) (in Norwegian):

(Toggle Plain Text)

# Finn st�rste felles divisor av tallene a og b.
def Finn_sfd(a, b):
    global stop, eualg

    if a % b == 0:  # Sjekk om 'b' g�r opp i 'a'
        print "%d modulo %d = 0\n" % (a, b)
        if stop == True:
            raw_input("Trykk Enter...\n")
        return b    # Hvis ja, ferdig; send tilbake resultat

    k = a / b       # Finn k (antall ganger 'b' g�r opp i 'a')
    r = a - b * k   # Finn rest
                    # a = k * b + r
    print "\t%d = %d * %d + %d" % (a, k, b, r)
                    #a - k * b = r
    print "\t%d - %d * %d = %d" % (a, k, b, r)

    # sfd(a, b) = sfd(b, r)
    print "sfd(%d, %d) = sfd(%d, %d)\n" % (a, b, b, r)
    if stop == True:
        raw_input("Trykk Enter...\n")
    eualg.append([r, a, k, b]) # Lagre tallene for senere bruk
    return Finn_sfd(b, r)

Sample run of the gcd function (also Norwegian):

(Toggle Plain Text)

Finner sfd(108, 63)...


Trykk Enter...

        108 = 1 * 63 + 45
        108 - 1 * 63 = 45
sfd(108, 63) = sfd(63, 45)

        63 = 1 * 45 + 18
        63 - 1 * 45 = 18
sfd(63, 45) = sfd(45, 18)

        45 = 2 * 18 + 9
        45 - 2 * 18 = 9
sfd(45, 18) = sfd(18, 9)

18 modulo 9 = 0

Arevos · Oct 11th, 2007, 11:53 AM

Quote:

Originally Posted by UnKnown X

Having stored all the numbers in the process (in the form of [[a,b,c,d],[a,b,c,d],[a,b,c,d],...]), how can I work my way 'back up'? Do I have to treat the equation as a string and keep replacing on my way up, or is there a simpler way?

You probably want the Extended Euclidean Algorithm. There's some useful pseudo-code down at the bottom of the linked Wikipedia article.

UnKnown X · Oct 11th, 2007, 11:59 AM

Thanks.

I had the extended Euclidean algorithm in mind for solving this, though I still have to find the two variables before that.

I'll check the article and links from there, though, thanks.

Arevos · Oct 11th, 2007, 12:06 PM

Can't you just ask the user for three numbers, e.g.

 python Syntax (Toggle Plain Text) 
print "ax + bx = c"
a = int(raw_input("a = "))
b = int(raw_input("b = "))
c = int(raw_input("c = "))
 
...
print "ax + bx = c"
a = int(raw_input("a = "))
b = int(raw_input("b = "))
c = int(raw_input("c = "))

...

UnKnown X · Oct 11th, 2007, 12:15 PM

That's what the program does already, but it's ax + by = c, and I have to find (i.e. calculate) a possible solution for x and y, then generalise it and find all possible integral solutions.

ZenMasterJG · Oct 11th, 2007, 12:38 PM

I took a quick stab at it. I'm not done, but heres a strategy I was working on.
First, find the maximum value for x (by minimizing y) and then a max value for y (by minimizing x)

to find the max value of x, write a loop that breaks when a solution is found, set y and x to 1, and increase x until the result is greater then c, then increment y and repeat (hope thats understandable)

To find all solutions, then, you could just bruteforce it and try every combination of numbers less then max x and max y. Theres certainly cleverer ways to do it, thats just as far as i've gotten so far and I have to go to class now

good luck, hope that hint helps.

Arevos · Oct 11th, 2007, 1:00 PM

Quote:

Originally Posted by UnKnown X

That's what the program does already, but it's ax + by = c, and I have to find (i.e. calculate) a possible solution for x and y, then generalise it and find all possible integral solutions.

Well, you're practically there, already!

Let's go over what we have so far. We have the Euclidean function gcd, and the Extended Euclidean function ex_gcd.

Given the equation 75x + 42y = 6, we can deduce the following:

 python Syntax (Toggle Plain Text) 
a = 75
b = 42
c = 6
 
d = gcd(a, b) = gcd(75, 42) = 3
 
(x0, y0) = ex_gcd(a, b) = ex_gcd(75, 42) = (-5, 9)
a = 75
b = 42
c = 6

d = gcd(a, b) = gcd(75, 42) = 3

(x0, y0) = ex_gcd(a, b) = ex_gcd(75, 42) = (-5, 9)

And the output we want is a complete range of solutions like:

(Toggle Plain Text)

x = -10 + 14n
y = 18 - 25n

You basically have all the information you need. You just need to use it to generate the right set of output ranges.

UnKnown X · Oct 11th, 2007, 7:51 PM

I've mostly figured it out. Thanks for the help!

The only things bothering me now are printing every step described in painstaking detail, and finding solutions for positive x and y--I'm on the right track, but it keeps messing up for some numbers, and working for others! I'll figure it out, though...

UnKnown X · Oct 12th, 2007, 4:41 PM

Okay, it's pretty much done now. I've posted a sample run at the bottom of this post (translated to English).

I might post the code later, if anyone's interested, but the code needs a lot more documentation and spiffing up (or down, in some cases) before being presentable.

Also, I realise that my usage of the implication, equivalence and inequality signs near the end is a bit clumsy and hard to read. If you have a better suggestion, do post it.

(Toggle Plain Text)

Input three integers a, b, c in the equation ax + bx = c

a = 108

b = 63

c = 765

------------------------------------------------------------

Finding gcd(108, 63)...

        108 = 1 * 63 + 45
        108 - 1 * 63 = 45
gcd(108, 63) = gcd(63, 45)

        63 = 1 * 45 + 18
        63 - 1 * 45 = 18
gcd(63, 45) = gcd(45, 18)

        45 = 2 * 18 + 9
        45 - 2 * 18 = 9
gcd(45, 18) = gcd(18, 9)

18 modulo 9 = 0

gcd(108, 63) = 9
        108 / 9 = 12 + 0
        63 / 9 = 7 + 0

------------------------------------------------------------

Finding lcm(108, 63)...

a * b = gcd(a, b) * lcm(a, b)
lcm(a, b) = a * b / gcd(a, b)
lcm = 108 * 63 / 9
lcm = 756

lcm(108, 63) = 756

------------------------------------------------------------

Checking if 108x + 63y = 765 has a solution...
        765 modulo 9 = 0

Has solution! Solving 108x + 63y = 765...

108 = 1 * 63 + 45
63 = 1 * 45 + 18
45 = 2 * 18 + 9
18 = 2 * 9 + 0
Returning 1 and -2...
Returning -2 and 3...
Returning 3 and -5...

9 = (3)*108 + (-5)*63

Finding a general solution for 108x + 63y = 765...

x[0] = 3 * (765 / 9) = 3 * 85 = 255
y[0] = -5 * (765 / 9) = -5 * 85 = -425
x = x[0] + b / d * z
y = y[0] - a / d * z
x = 255 + 63 / 9 * z = 255 + 7z
y = -425 - 108 / 9 * z = -425 - 12z
108(255 + 7z) + 63(-425 - 12z) = 765

Finding positive integral solutions...

... for x...
        x >= 0 =>
        255 + 7z >= 0 <=>
        7z >= -255 =>
        z >= -36.429
        z >= -36
... for y...
        y >= 0 =>
        -425 - 12z >= 0 <=>
        -12z <= 425 <=>
        12z <= -425 =>
        z <= -35.417
        z <= -35

There are 2 positive integral solutions for z between -35 and -36.

z = -36
        108(255 + 7 * (-36)) + 63(-425 - 12 * (-36))
        = 108 * (3) + 63 * (7)
        = 765
z = -35
        108(255 + 7 * (-35)) + 63(-425 - 12 * (-35))
        = 108 * (10) + 63 * (-5)
        = 765
------------------------------------------------------------

Done!

Oct 10th, 2007, 11:18 AM	#1
UnKnown X Hobbyist Programmer Join Date: Dec 2005 Location: Sandvika, Norway Posts: 114 Rep Power: 0	Solving Diophantine Equations I have a project where I am to write a programme that illustrates how computers can be used for solving iterative tasks, using Euclid's algorithm as an example, and presenting and explaining it to a class in which most students have little to no knowledge about programming. Thus, I need to keep it as simple as possible, and I chose to use a high-level language like Python. Euclid's algorithm went smoothly, but I was additionally requested to add a function for solving linear diophantine equations with two variables and one equation, and I'm not even sure how to begin tackling this. I'm not asking for a solution at all, since this is my task; I only need a way to get started with this. For example, given the equation 75x + 42y = 3, the greatest common divisor is 3. 33 = 75 - 1 * 42 9 = 42 - 1 * 33 6 = 33 - 3 * 9 3 = 9 - 1 * 6 Having stored all the numbers in the process (in the form of [[a,b,c,d],[a,b,c,d],[a,b,c,d],...]), how can I work my way 'back up'? Do I have to treat the equation as a string and keep replacing on my way up, or is there a simpler way? Here's the function I wrote to find gdc(a,b) (in Norwegian): (Toggle Plain Text) # Finn st�rste felles divisor av tallene a og b. def Finn_sfd(a, b): global stop, eualg if a % b == 0: # Sjekk om 'b' g�r opp i 'a' print "%d modulo %d = 0\n" % (a, b) if stop == True: raw_input("Trykk Enter...\n") return b # Hvis ja, ferdig; send tilbake resultat k = a / b # Finn k (antall ganger 'b' g�r opp i 'a') r = a - b * k # Finn rest # a = k * b + r print "\t%d = %d * %d + %d" % (a, k, b, r) #a - k * b = r print "\t%d - %d * %d = %d" % (a, k, b, r) # sfd(a, b) = sfd(b, r) print "sfd(%d, %d) = sfd(%d, %d)\n" % (a, b, b, r) if stop == True: raw_input("Trykk Enter...\n") eualg.append([r, a, k, b]) # Lagre tallene for senere bruk return Finn_sfd(b, r) # Finn st�rste felles divisor av tallene a og b. def Finn_sfd(a, b): global stop, eualg if a % b == 0: # Sjekk om 'b' g�r opp i 'a' print "%d modulo %d = 0\n" % (a, b) if stop == True: raw_input("Trykk Enter...\n") return b # Hvis ja, ferdig; send tilbake resultat k = a / b # Finn k (antall ganger 'b' g�r opp i 'a') r = a - b * k # Finn rest # a = k * b + r print "\t%d = %d * %d + %d" % (a, k, b, r) #a - k * b = r print "\t%d - %d * %d = %d" % (a, k, b, r) # sfd(a, b) = sfd(b, r) print "sfd(%d, %d) = sfd(%d, %d)\n" % (a, b, b, r) if stop == True: raw_input("Trykk Enter...\n") eualg.append([r, a, k, b]) # Lagre tallene for senere bruk return Finn_sfd(b, r) Sample run of the gcd function (also Norwegian): (Toggle Plain Text) Finner sfd(108, 63)... Trykk Enter... 108 = 1 * 63 + 45 108 - 1 * 63 = 45 sfd(108, 63) = sfd(63, 45) 63 = 1 * 45 + 18 63 - 1 * 45 = 18 sfd(63, 45) = sfd(45, 18) 45 = 2 * 18 + 9 45 - 2 * 18 = 9 sfd(45, 18) = sfd(18, 9) 18 modulo 9 = 0 Finner sfd(108, 63)... Trykk Enter... 108 = 1 * 63 + 45 108 - 1 * 63 = 45 sfd(108, 63) = sfd(63, 45) 63 = 1 * 45 + 18 63 - 1 * 45 = 18 sfd(63, 45) = sfd(45, 18) 45 = 2 * 18 + 9 45 - 2 * 18 = 9 sfd(45, 18) = sfd(18, 9) 18 modulo 9 = 0

Oct 11th, 2007, 12:06 PM	#4
Arevos Programming Guru Join Date: Aug 2005 Location: England Posts: 1,499 Rep Power: 5	Can't you just ask the user for three numbers, e.g. python Syntax (Toggle Plain Text) print "ax + bx = c" a = int(raw_input("a = ")) b = int(raw_input("b = ")) c = int(raw_input("c = ")) ... print "ax + bx = c" a = int(raw_input("a = ")) b = int(raw_input("b = ")) c = int(raw_input("c = ")) ...

Oct 12th, 2007, 4:41 PM	#9
UnKnown X Hobbyist Programmer Join Date: Dec 2005 Location: Sandvika, Norway Posts: 114 Rep Power: 0	Okay, it's pretty much done now. I've posted a sample run at the bottom of this post (translated to English). I might post the code later, if anyone's interested, but the code needs a lot more documentation and spiffing up (or down, in some cases) before being presentable. Also, I realise that my usage of the implication, equivalence and inequality signs near the end is a bit clumsy and hard to read. If you have a better suggestion, do post it. (Toggle Plain Text) Input three integers a, b, c in the equation ax + bx = c a = 108 b = 63 c = 765 ------------------------------------------------------------ Finding gcd(108, 63)... 108 = 1 * 63 + 45 108 - 1 * 63 = 45 gcd(108, 63) = gcd(63, 45) 63 = 1 * 45 + 18 63 - 1 * 45 = 18 gcd(63, 45) = gcd(45, 18) 45 = 2 * 18 + 9 45 - 2 * 18 = 9 gcd(45, 18) = gcd(18, 9) 18 modulo 9 = 0 gcd(108, 63) = 9 108 / 9 = 12 + 0 63 / 9 = 7 + 0 ------------------------------------------------------------ Finding lcm(108, 63)... a * b = gcd(a, b) * lcm(a, b) lcm(a, b) = a * b / gcd(a, b) lcm = 108 * 63 / 9 lcm = 756 lcm(108, 63) = 756 ------------------------------------------------------------ Checking if 108x + 63y = 765 has a solution... 765 modulo 9 = 0 Has solution! Solving 108x + 63y = 765... 108 = 1 * 63 + 45 63 = 1 * 45 + 18 45 = 2 * 18 + 9 18 = 2 * 9 + 0 Returning 1 and -2... Returning -2 and 3... Returning 3 and -5... 9 = (3)108 + (-5)63 Finding a general solution for 108x + 63y = 765... x[0] = 3 * (765 / 9) = 3 * 85 = 255 y[0] = -5 * (765 / 9) = -5 * 85 = -425 x = x[0] + b / d * z y = y[0] - a / d * z x = 255 + 63 / 9 * z = 255 + 7z y = -425 - 108 / 9 * z = -425 - 12z 108(255 + 7z) + 63(-425 - 12z) = 765 Finding positive integral solutions... ... for x... x >= 0 => 255 + 7z >= 0 <=> 7z >= -255 => z >= -36.429 z >= -36 ... for y... y >= 0 => -425 - 12z >= 0 <=> -12z <= 425 <=> 12z <= -425 => z <= -35.417 z <= -35 There are 2 positive integral solutions for z between -35 and -36. z = -36 108(255 + 7 * (-36)) + 63(-425 - 12 * (-36)) = 108 * (3) + 63 * (7) = 765 z = -35 108(255 + 7 * (-35)) + 63(-425 - 12 * (-35)) = 108 * (10) + 63 * (-5) = 765 ------------------------------------------------------------ Done! Input three integers a, b, c in the equation ax + bx = c a = 108 b = 63 c = 765 ------------------------------------------------------------ Finding gcd(108, 63)... 108 = 1 * 63 + 45 108 - 1 * 63 = 45 gcd(108, 63) = gcd(63, 45) 63 = 1 * 45 + 18 63 - 1 * 45 = 18 gcd(63, 45) = gcd(45, 18) 45 = 2 * 18 + 9 45 - 2 * 18 = 9 gcd(45, 18) = gcd(18, 9) 18 modulo 9 = 0 gcd(108, 63) = 9 108 / 9 = 12 + 0 63 / 9 = 7 + 0 ------------------------------------------------------------ Finding lcm(108, 63)... a * b = gcd(a, b) * lcm(a, b) lcm(a, b) = a * b / gcd(a, b) lcm = 108 * 63 / 9 lcm = 756 lcm(108, 63) = 756 ------------------------------------------------------------ Checking if 108x + 63y = 765 has a solution... 765 modulo 9 = 0 Has solution! Solving 108x + 63y = 765... 108 = 1 * 63 + 45 63 = 1 * 45 + 18 45 = 2 * 18 + 9 18 = 2 * 9 + 0 Returning 1 and -2... Returning -2 and 3... Returning 3 and -5... 9 = (3)108 + (-5)63 Finding a general solution for 108x + 63y = 765... x[0] = 3 * (765 / 9) = 3 * 85 = 255 y[0] = -5 * (765 / 9) = -5 * 85 = -425 x = x[0] + b / d * z y = y[0] - a / d * z x = 255 + 63 / 9 * z = 255 + 7z y = -425 - 108 / 9 * z = -425 - 12z 108(255 + 7z) + 63(-425 - 12z) = 765 Finding positive integral solutions... ... for x... x >= 0 => 255 + 7z >= 0 <=> 7z >= -255 => z >= -36.429 z >= -36 ... for y... y >= 0 => -425 - 12z >= 0 <=> -12z <= 425 <=> 12z <= -425 => z <= -35.417 z <= -35 There are 2 positive integral solutions for z between -35 and -36. z = -36 108(255 + 7 * (-36)) + 63(-425 - 12 * (-36)) = 108 * (3) + 63 * (7) = 765 z = -35 108(255 + 7 * (-35)) + 63(-425 - 12 * (-35)) = 108 * (10) + 63 * (-5) = 765 ------------------------------------------------------------ Done!

Thread Tools
Show Printable Version Email this Page
Display Modes
Linear Mode Switch to Hybrid Mode Switch to Threaded Mode

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Win PCs, xboxes, zunes and more by solving puzzles!	bulio	Coder's Corner Lounge	0	Jan 9th, 2007 8:10 PM
solving warning messages, in a program	kalulu	C	4	Nov 25th, 2005 8:41 PM

Oct 11th, 2007, 11:59 AM	#3
UnKnown X Hobbyist Programmer Join Date: Dec 2005 Location: Sandvika, Norway Posts: 114 Rep Power: 0	Thanks. I had the extended Euclidean algorithm in mind for solving this, though I still have to find the two variables before that. I'll check the article and links from there, though, thanks.

Oct 11th, 2007, 12:15 PM	#5
UnKnown X Hobbyist Programmer Join Date: Dec 2005 Location: Sandvika, Norway Posts: 114 Rep Power: 0	That's what the program does already, but it's ax + by = c, and I have to find (i.e. calculate) a possible solution for x and y, then generalise it and find all possible integral solutions.

Oct 11th, 2007, 12:38 PM	#6
ZenMasterJG Hobbyist Programmer Join Date: Nov 2004 Location: Boston, MA Posts: 148 Rep Power: 4	I took a quick stab at it. I'm not done, but heres a strategy I was working on. First, find the maximum value for x (by minimizing y) and then a max value for y (by minimizing x) to find the max value of x, write a loop that breaks when a solution is found, set y and x to 1, and increase x until the result is greater then c, then increment y and repeat (hope thats understandable) To find all solutions, then, you could just bruteforce it and try every combination of numbers less then max x and max y. Theres certainly cleverer ways to do it, thats just as far as i've gotten so far and I have to go to class now good luck, hope that hint helps.

Oct 11th, 2007, 7:51 PM	#8
UnKnown X Hobbyist Programmer Join Date: Dec 2005 Location: Sandvika, Norway Posts: 114 Rep Power: 0	I've mostly figured it out. Thanks for the help! The only things bothering me now are printing every step described in painstaking detail, and finding solutions for positive x and y--I'm on the right track, but it keeps messing up for some numbers, and working for others! I'll figure it out, though...

Currently Active Users Viewing This Thread: 1 (0 members and 1 guests)

Programming Forums > Script Programming > Python

Solving Diophantine Equations