Ο κώδικας είναι γραμμένος στα αγγλικά και αναρτήθηκε πρώτη φορά στο RosettaCode.
Για τη λύση έχει επιλεχθεί η απαλοιφή Gauss όπως ήταν και το απαιτούμενο για το έργο στο RosettaCode.
Μας ενδιαφέρει η ακρίβεια. Για το σκοπό αυτό θα χρησιμοποιήσουμε αριθμούς τύπου Decimal. Πληροφορίες για τον τύπο αριθμών Decimal υπάρχουν εδώ. Ο αριθμός 10@ είναι τύπου decimal λόγω του χαρακτήρα @ που ακολουθεί τα ψηφία. Έτσι το α=1212.121@ θα κάνει την α τύπου Decimal.
Από την γραμμή εντολών μπορούμε να ορίσουμε γενικές μεταβλητές τύπου Decimal.
>α=1212.121@
>? Τύπος$(α)
Decimal
Ακολουθούν δυο παραδείγματα και στα δύο αυτά υπάρχει ο πίνακας a() με δυο σετ τιμών. Διαλέγουμε πριν εκτελέσουμε κάποιο από τα δύο (σκιάζουμε την γραμμή με τα στοιχεία με έναν χαρακτήρα"\" - συνήθως βάζουμε δύο- \\.
Σε καθένα θα γίνει η εύρεση των λύσεων (για απλοποίηση δεχόμαστε ότι τα δεδομένα που βάζουμε έχουν λύση).
Στο τέλος θα γίνει επαλήθευση ως προς ένα βαθμό στρογγυλοποίησης. Η Round() κάνει στρογγυλοποίηση στο μισό, δηλαδή αν ζητάμε στο 3 δεκαδικό τότε αν το τέταρτο είναι από 5 και πάνω, κερδίζουμε ένα χιλιοστό.
Οι Decimal αριθμοί έχουν μέγιστο 28 ψηφία, και μπορούν να μετακινούν το σημείο δεκαδικών, έτσι αν έχουμε έναν αριθμό με ακέραια ψηφία 24 δεν μπορούμε να έχουμε 5 δεκαδικά, αλλά το πολύ τέσσερα. Ομοίως αν έχουμε έναν αριθμό με ακέραια ψηφία 4 δεν μπορούμε να έχουμε 25 δεκαδικά, αλλά το πολύ 24.
Αν παρελπίδα δώσουμε σε μεταβλητή που δέχεται decimal αριθμό που δεν χωράει τότε βγαίνει λάθος. Αν ο Decimal είναι σε πίνακα, τότε δεν βγαίνει λάθος αλλά αλλάζει μορφή γίνεται Double.
Μας ενδιαφέρει το πρόγραμμα να εξαγει αποτελέσματα στην οθόνη και στο πρόχειρο (για να τα επικολλήσουμε κάπου, όπως στο Rosettacode). Για το σκοπό αυτό υπάρχει ένα αντικείμενο Document, το οποίο είναι αλφαριθμητικό με παραγράφους. Σε αυτό τροφοδοτούμε στοιχεία με το "=" τα οποία προστίθενται στο τέλος ( λειτουργία append). Επειδή θέλουμε να βάλουμε και τους πίνακες μέσα έχουμε μια βοηθητική συνάρτηση η οποία έχει εκτός από τον πίνακα που μας ενδιαφέρει (θέλουμε να υποστηρίζει πίνακες μιας και δυο διαστάσεων) και άλλες παραμέτρους. Δέχεται προαιρετικά το μέγεθος περιθωρίου αριστερά (left margin), επίσης προαιρετικά το μέγεθος πεδίου αριθμού, και τέλος επίσης προαιρετικά ένα αλφαριθμητικό που μπορεί να είναι κενό ή να έχει έναν αριθμό. Αν είναι κενό τότε αφήνει την υποδιαστολή ως έχει, αν δεν είναι κενό τότε αν είναι μηδέν δεν δείχνει την υποδιαστολή και στρογγυλεύει το νούμερο στο ακέραιο (με το μισό και πάνω υπέρ της μονάδας). Αν δώσουμε το "2" τότε θα στρογγυλέψει στα δυο δεκαδικά και ακόμα αν δεν έχουμε δεκαδικά θα βάλει δυο μηδενικά!
Στο πρώτο πρόγραμμα χρησιμοποιούμε αποθήκευση σκέτου αριθμού, ενώ στο δεύτερο αποθηκεύουμε κλάσματα.
Μια εντολή form 80,50 μπορεί να μεγαλώσει την ανάλυση της κονσόλας, για να βλέπετε όλα τα αποτελέσματα. (ή Φόρμα 80,80 με ή χωρίς τόνο στην εντολή).
Εδώ στο δεύτερο πρόγραμμα αποθηκεύουμε σε πίνακα δυο στοιχείων κάθε στοιχείο του πίνακα Α, όπου το στοιχείο στη θέση 0 είναι ο αριθμητής και στη θέση 1 ο παρανομαστής. Με αυτό τον τρόπο έχουμε πιο αργό αλγόριθμο αλλά έχει μεγαλύτερη ακρίβεια! (ένα δεκαδικό)
Ο αλγόριθμος για τα κλάσματα παίρνει βελτιώσεις, αλλά για την ώρα μένουμε εδώ (για να μην δυσκολέψει το πράγμα στη κατανόηση).
Παρακάτω ακολουθούν αποτελέσματα με αλλαγές στο acc, και στην επιλογή τιμών στον πίνακα a().
Για τη λύση έχει επιλεχθεί η απαλοιφή Gauss όπως ήταν και το απαιτούμενο για το έργο στο RosettaCode.
Μας ενδιαφέρει η ακρίβεια. Για το σκοπό αυτό θα χρησιμοποιήσουμε αριθμούς τύπου Decimal. Πληροφορίες για τον τύπο αριθμών Decimal υπάρχουν εδώ. Ο αριθμός 10@ είναι τύπου decimal λόγω του χαρακτήρα @ που ακολουθεί τα ψηφία. Έτσι το α=1212.121@ θα κάνει την α τύπου Decimal.
Από την γραμμή εντολών μπορούμε να ορίσουμε γενικές μεταβλητές τύπου Decimal.
>α=1212.121@
>? Τύπος$(α)
Decimal
Ακολουθούν δυο παραδείγματα και στα δύο αυτά υπάρχει ο πίνακας a() με δυο σετ τιμών. Διαλέγουμε πριν εκτελέσουμε κάποιο από τα δύο (σκιάζουμε την γραμμή με τα στοιχεία με έναν χαρακτήρα"\" - συνήθως βάζουμε δύο- \\.
Σε καθένα θα γίνει η εύρεση των λύσεων (για απλοποίηση δεχόμαστε ότι τα δεδομένα που βάζουμε έχουν λύση).
Στο τέλος θα γίνει επαλήθευση ως προς ένα βαθμό στρογγυλοποίησης. Η Round() κάνει στρογγυλοποίηση στο μισό, δηλαδή αν ζητάμε στο 3 δεκαδικό τότε αν το τέταρτο είναι από 5 και πάνω, κερδίζουμε ένα χιλιοστό.
Οι Decimal αριθμοί έχουν μέγιστο 28 ψηφία, και μπορούν να μετακινούν το σημείο δεκαδικών, έτσι αν έχουμε έναν αριθμό με ακέραια ψηφία 24 δεν μπορούμε να έχουμε 5 δεκαδικά, αλλά το πολύ τέσσερα. Ομοίως αν έχουμε έναν αριθμό με ακέραια ψηφία 4 δεν μπορούμε να έχουμε 25 δεκαδικά, αλλά το πολύ 24.
Αν παρελπίδα δώσουμε σε μεταβλητή που δέχεται decimal αριθμό που δεν χωράει τότε βγαίνει λάθος. Αν ο Decimal είναι σε πίνακα, τότε δεν βγαίνει λάθος αλλά αλλάζει μορφή γίνεται Double.
Μας ενδιαφέρει το πρόγραμμα να εξαγει αποτελέσματα στην οθόνη και στο πρόχειρο (για να τα επικολλήσουμε κάπου, όπως στο Rosettacode). Για το σκοπό αυτό υπάρχει ένα αντικείμενο Document, το οποίο είναι αλφαριθμητικό με παραγράφους. Σε αυτό τροφοδοτούμε στοιχεία με το "=" τα οποία προστίθενται στο τέλος ( λειτουργία append). Επειδή θέλουμε να βάλουμε και τους πίνακες μέσα έχουμε μια βοηθητική συνάρτηση η οποία έχει εκτός από τον πίνακα που μας ενδιαφέρει (θέλουμε να υποστηρίζει πίνακες μιας και δυο διαστάσεων) και άλλες παραμέτρους. Δέχεται προαιρετικά το μέγεθος περιθωρίου αριστερά (left margin), επίσης προαιρετικά το μέγεθος πεδίου αριθμού, και τέλος επίσης προαιρετικά ένα αλφαριθμητικό που μπορεί να είναι κενό ή να έχει έναν αριθμό. Αν είναι κενό τότε αφήνει την υποδιαστολή ως έχει, αν δεν είναι κενό τότε αν είναι μηδέν δεν δείχνει την υποδιαστολή και στρογγυλεύει το νούμερο στο ακέραιο (με το μισό και πάνω υπέρ της μονάδας). Αν δώσουμε το "2" τότε θα στρογγυλέψει στα δυο δεκαδικά και ακόμα αν δεν έχουμε δεκαδικά θα βάλει δυο μηδενικά!
Στο πρώτο πρόγραμμα χρησιμοποιούμε αποθήκευση σκέτου αριθμού, ενώ στο δεύτερο αποθηκεύουμε κλάσματα.
Μια εντολή form 80,50 μπορεί να μεγαλώσει την ανάλυση της κονσόλας, για να βλέπετε όλα τα αποτελέσματα. (ή Φόρμα 80,80 με ή χωρίς τόνο στην εντολή).
module
checkit {
Dim Base 1, a(6, 6), b(6)
a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
\\ remove \\ to feed next array
\\ a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
function GaussJordan(a(), b()) {
cols=dimension(a(),1)
rows=dimension(a(),2)
\\ make augmented matrix
Dim Base 1, a(cols, rows)
\\ feed array with rationals
Dim Base 1, b(Len(b()))
for diag=1 to rows {
max_row=diag
max_val=abs(a(diag, diag))
if diag<rows Then {
for ro=diag+1 to rows {
d=abs(a(ro, diag))
if d>max_val then max_row=ro : max_val=d
}
}
\\ SwapRows diag, max_row
if diag<>max_row then {
for i=1 to cols {
swap a(diag, i), a(max_row, i)
}
swap b(diag), b(max_row)
}
invd= a(diag, diag)
if diag<=cols then {
for col=diag to cols {
a(diag, col)/=invd
}
}
b(diag)/=invd
for ro=1 to rows {
d1=a(ro,diag)
d2=d1*b(diag)
if ro<>diag Then {
for col=diag to cols {a(ro, col)-=d1*a(diag, col)}
b(ro)-=d2
}
}
}
=b()
}
Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
\\ defualt no set decimals, can show any number
ex$={
}
const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
if dimension(a())=1 then {
m=each(a())
while m {ex$+=format$(way$,array(m))}
Insert 3, 2 ex$=string$(" ", leftmargin)
=ex$ : Break
}
for i=1 to dimension(a(),1) {
ex1$=""
for j=1 to dimension(a(),2 ) {
ex1$+=format$(way$,a(i,j))
}
Insert 1,2 ex1$=string$(" ", leftmargin)
ex$+=ex1$+{
}
}
=ex$
}
mm=GaussJordan(a(), b())
c=each(mm)
while c {
print array(c)
}
\\ check accuracy
link mm to r()
\\ prepare output document
Document out$={Algorithm using decimals
}+"Matrix A:"+ArrayLines$(a(),,,"2")+{
}+"Vector B:"+ArrayLines$(b(),,,"2")+{
}+"Solution: "+{
}
acc=25
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
p$=format$("Coef. {0::-2}, rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
Print p$
Out$=p$+{
}
next i
Report out$
clipboard out$
}
checkit
Dim Base 1, a(6, 6), b(6)
a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
\\ remove \\ to feed next array
\\ a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
function GaussJordan(a(), b()) {
cols=dimension(a(),1)
rows=dimension(a(),2)
\\ make augmented matrix
Dim Base 1, a(cols, rows)
\\ feed array with rationals
Dim Base 1, b(Len(b()))
for diag=1 to rows {
max_row=diag
max_val=abs(a(diag, diag))
if diag<rows Then {
for ro=diag+1 to rows {
d=abs(a(ro, diag))
if d>max_val then max_row=ro : max_val=d
}
}
\\ SwapRows diag, max_row
if diag<>max_row then {
for i=1 to cols {
swap a(diag, i), a(max_row, i)
}
swap b(diag), b(max_row)
}
invd= a(diag, diag)
if diag<=cols then {
for col=diag to cols {
a(diag, col)/=invd
}
}
b(diag)/=invd
for ro=1 to rows {
d1=a(ro,diag)
d2=d1*b(diag)
if ro<>diag Then {
for col=diag to cols {a(ro, col)-=d1*a(diag, col)}
b(ro)-=d2
}
}
}
=b()
}
Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
\\ defualt no set decimals, can show any number
ex$={
}
const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
if dimension(a())=1 then {
m=each(a())
while m {ex$+=format$(way$,array(m))}
Insert 3, 2 ex$=string$(" ", leftmargin)
=ex$ : Break
}
for i=1 to dimension(a(),1) {
ex1$=""
for j=1 to dimension(a(),2 ) {
ex1$+=format$(way$,a(i,j))
}
Insert 1,2 ex1$=string$(" ", leftmargin)
ex$+=ex1$+{
}
}
=ex$
}
mm=GaussJordan(a(), b())
c=each(mm)
while c {
print array(c)
}
\\ check accuracy
link mm to r()
\\ prepare output document
Document out$={Algorithm using decimals
}+"Matrix A:"+ArrayLines$(a(),,,"2")+{
}+"Vector B:"+ArrayLines$(b(),,,"2")+{
}+"Solution: "+{
}
acc=25
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
p$=format$("Coef. {0::-2}, rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
Print p$
Out$=p$+{
}
next i
Report out$
clipboard out$
}
checkit
Εδώ στο δεύτερο πρόγραμμα αποθηκεύουμε σε πίνακα δυο στοιχείων κάθε στοιχείο του πίνακα Α, όπου το στοιχείο στη θέση 0 είναι ο αριθμητής και στη θέση 1 ο παρανομαστής. Με αυτό τον τρόπο έχουμε πιο αργό αλγόριθμο αλλά έχει μεγαλύτερη ακρίβεια! (ένα δεκαδικό)
Ο αλγόριθμος για τα κλάσματα παίρνει βελτιώσεις, αλλά για την ώρα μένουμε εδώ (για να μην δυσκολέψει το πράγμα στη κατανόηση).
Module
Checkit2 {
Dim Base 1, a(6, 6), b(6)
\\ a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
\\ modules/function to use rational nymbers
Module Global subd(m as array, n as array) { ' change m
link m to m()
link n to n()
if m(0)=0 then return m, 0:=-n(0), 1:=n(1) : exit
if n(0)=0 then exit
return m, 0:=m(0)*(n(1)/m(1))-n(0), 1:=n(1)
}
Function Global Inv(m as array){
link m to m()
if m(0)=0@ then =m : exit
=(m(1), m(0))
}
Function Global mul(m as array, n as array){' nothing change
link m to m()
link n to n()
if n(0)=0 or n(1)=0 then =(0@,0@) : exit
=((m(0)/n(1))*n(0),m(1))
}
Module Global mul(m as array, n as array) { ' change m
link m to m()
link n to n()
if n(0)=0 or n(1)=0 then m=(0@,0@) : exit
return m, 0:=(m(0)/n(1))*n(0)
}
Function Global Res(m as array) {
link m to m()
if m(0)=0@ then =0@: exit
=m(0)/m(1)
}
\\ GaussJordan get arrays byvalue
function GaussJordan(a(), b()) {
Function copypointer(m) { Dim a() : a()=m:=a()}
\\ we can use : def copypointer(a())=a(0),a(1)
cols=dimension(a(),1)
rows=dimension(a(),2)
Dim Base 1, a(cols, rows)
for i=1 to cols : for j=1 to rows : a(i, j)=(a(i, j), 1@) : next j : next i
def d as decimal
for j=1 to rows : b(j)=(b(j), 1@) : next j
for diag=1 to rows {
max_row=diag
max_val=abs(Res(a(diag, diag)))
if diag<rows Then {
for ro=diag+1 to rows {
d=abs(Res(a(ro, diag)))
if d>max_val then max_row=ro : max_val=d
}
}
\\ SwapRows diag, max_row
if diag<>max_row then {
for i=1 to cols {
swap a(diag, i), a(max_row, i)
}
swap b(diag), b(max_row)
}
invd= Inv(a(diag, diag))
if diag<=cols then {
for col=diag to cols {
mul a(diag, col), invd
}
}
mul b(diag), invd
for ro=1 to rows {
\\ work also d1=(a(ro,diag)(0), a(ro,diag)(1))
d1=copypointer(a(ro, diag))
if ro<>diag Then {
for col=diag to cols {subd a(ro, col), mul(d1, a(diag, col))}
subd b(ro), mul(d1, b(diag))
}
}
}
dim base 1, ans(len(b()))
for i=1 to cols {
ans(i)=res(b(i)) \\ : Print b(i) ' print pairs
}
=ans()
}
Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
\\ defualt no set decimals, can show any number
ex$={
}
const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
if dimension(a())=1 then {
m=each(a())
while m {ex$+=format$(way$,array(m))}
Insert 3, 2 ex$=string$(" ", leftmargin)
=ex$ : Break
}
for i=1 to dimension(a(),1) {
ex1$=""
for j=1 to dimension(a(),2 ) {
ex1$+=format$(way$,a(i,j))
}
Insert 1,2 ex1$=string$(" ", leftmargin)
ex$+=ex1$+{
}
}
=ex$
}
mm=GaussJordan(a(), b())
c=each(mm)
while c {
print array(c)
}
\\ check accuracy
link mm to r()
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
Print round(sum-b(i),26), b(i)
next i
\\ check accuracy
Document out$={Algorithm using pair of decimals as rational numbers
}+"Matrix A:"+ArrayLines$(a(),,,"2")+{
}+"Vector B:"+ArrayLines$(b(),,,"2")+{
}+"Solution: "+{
}
acc=26
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
p$=format$("Coef. {0::-2}, rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
Print p$
Out$=p$+{
}
next i
Report out$
clipboard out$
}
Checkit2
Dim Base 1, a(6, 6), b(6)
\\ a(1,1)= 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.63, 0.39, 0.25, 0.16, 0.10, 1.00, 1.26, 1.58, 1.98, 2.49, 3.13, 1.00, 1.88, 3.55, 6.70, 12.62, 23.80, 1.00, 2.51, 6.32, 15.88, 39.90, 100.28, 1.00, 3.14, 9.87, 31.01, 97.41, 306.02
a(1,1)=1.1,0.12,0.13,0.12,0.14,-0.12,1.21,0.63,0.39,0.25,0.16,0.1,1.03,1.26,1.58,1.98,2.49,3.13, 1.06,1.88,3.55,6.7,12.62,23.8, 1.12,2.51,6.32,15.88,39.9,100.28,1.16,3.14,9.87,31.01,97.41,306.02
for i=1 to 6 : for j=1 to 6 : a(i,j)=val(a(i,j)->Decimal) :Next j:Next i
b(1)=-0.01, 0.61, 0.91, 0.99, 0.60, 0.02
for i=1 to 6 : b(i)=val(b(i)->Decimal) :Next i
\\ modules/function to use rational nymbers
Module Global subd(m as array, n as array) { ' change m
link m to m()
link n to n()
if m(0)=0 then return m, 0:=-n(0), 1:=n(1) : exit
if n(0)=0 then exit
return m, 0:=m(0)*(n(1)/m(1))-n(0), 1:=n(1)
}
Function Global Inv(m as array){
link m to m()
if m(0)=0@ then =m : exit
=(m(1), m(0))
}
Function Global mul(m as array, n as array){' nothing change
link m to m()
link n to n()
if n(0)=0 or n(1)=0 then =(0@,0@) : exit
=((m(0)/n(1))*n(0),m(1))
}
Module Global mul(m as array, n as array) { ' change m
link m to m()
link n to n()
if n(0)=0 or n(1)=0 then m=(0@,0@) : exit
return m, 0:=(m(0)/n(1))*n(0)
}
Function Global Res(m as array) {
link m to m()
if m(0)=0@ then =0@: exit
=m(0)/m(1)
}
\\ GaussJordan get arrays byvalue
function GaussJordan(a(), b()) {
Function copypointer(m) { Dim a() : a()=m:=a()}
\\ we can use : def copypointer(a())=a(0),a(1)
cols=dimension(a(),1)
rows=dimension(a(),2)
Dim Base 1, a(cols, rows)
for i=1 to cols : for j=1 to rows : a(i, j)=(a(i, j), 1@) : next j : next i
def d as decimal
for j=1 to rows : b(j)=(b(j), 1@) : next j
for diag=1 to rows {
max_row=diag
max_val=abs(Res(a(diag, diag)))
if diag<rows Then {
for ro=diag+1 to rows {
d=abs(Res(a(ro, diag)))
if d>max_val then max_row=ro : max_val=d
}
}
\\ SwapRows diag, max_row
if diag<>max_row then {
for i=1 to cols {
swap a(diag, i), a(max_row, i)
}
swap b(diag), b(max_row)
}
invd= Inv(a(diag, diag))
if diag<=cols then {
for col=diag to cols {
mul a(diag, col), invd
}
}
mul b(diag), invd
for ro=1 to rows {
\\ work also d1=(a(ro,diag)(0), a(ro,diag)(1))
d1=copypointer(a(ro, diag))
if ro<>diag Then {
for col=diag to cols {subd a(ro, col), mul(d1, a(diag, col))}
subd b(ro), mul(d1, b(diag))
}
}
}
dim base 1, ans(len(b()))
for i=1 to cols {
ans(i)=res(b(i)) \\ : Print b(i) ' print pairs
}
=ans()
}
Function ArrayLines$(a(), leftmargin=6, maxwidth=8,decimals$="") {
\\ defualt no set decimals, can show any number
ex$={
}
const way$=", {0:"+decimals$+":-"+str$(maxwidth,"")+"}"
if dimension(a())=1 then {
m=each(a())
while m {ex$+=format$(way$,array(m))}
Insert 3, 2 ex$=string$(" ", leftmargin)
=ex$ : Break
}
for i=1 to dimension(a(),1) {
ex1$=""
for j=1 to dimension(a(),2 ) {
ex1$+=format$(way$,a(i,j))
}
Insert 1,2 ex1$=string$(" ", leftmargin)
ex$+=ex1$+{
}
}
=ex$
}
mm=GaussJordan(a(), b())
c=each(mm)
while c {
print array(c)
}
\\ check accuracy
link mm to r()
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
Print round(sum-b(i),26), b(i)
next i
\\ check accuracy
Document out$={Algorithm using pair of decimals as rational numbers
}+"Matrix A:"+ArrayLines$(a(),,,"2")+{
}+"Vector B:"+ArrayLines$(b(),,,"2")+{
}+"Solution: "+{
}
acc=26
for i=1 to dimension(a(),1)
sum=a(1,1)-a(1,1)
For j=1 to dimension(a(),2)
sum+=r(j)*a(i,j)
next j
p$=format$("Coef. {0::-2}, rounding to {1} decimal, compare {2:-5}, solution: {3}", i, acc, round(sum-b(i),acc)=0@, r(i) )
Print p$
Out$=p$+{
}
next i
Report out$
clipboard out$
}
Checkit2
Παρακάτω ακολουθούν αποτελέσματα με αλλαγές στο acc, και στην επιλογή τιμών στον πίνακα a().
Algorithm using decimals
Matrix A:
1,10, 0,12, 0,13, 0,12, 0,14, -0,12
1,21, 0,63, 0,39, 0,25, 0,16, 0,10
1,03, 1,26, 1,58, 1,98, 2,49, 3,13
1,06, 1,88, 3,55, 6,70, 12,62, 23,80
1,12, 2,51, 6,32, 15,88, 39,90, 100,28
1,16, 3,14, 9,87, 31,01, 97,41, 306,02
Vector B:
-0,01, 0,61, 0,91, 0,99, 0,60, 0,02
Solution:
Coef. 1, rounding to 26 decimal, compare True, solution: -0,0597391027501962649904316335
Coef. 2, rounding to 26 decimal, compare True, solution: 1,8501896672627829700670299288
Coef. 3, rounding to 26 decimal, compare True, solution: -1,9727833018116428175300387318
Coef. 4, rounding to 26 decimal, compare True, solution: 1,4697587750651240151384675034
Coef. 5, rounding to 26 decimal, compare True, solution: -0,5538741847821888403564152897
Coef. 6, rounding to 26 decimal, compare True, solution: 0,0723048745759411900531809852
Algorithm using pair of decimals as rational numbers
Matrix A:
1,10, 0,12, 0,13, 0,12, 0,14, -0,12
1,21, 0,63, 0,39, 0,25, 0,16, 0,10
1,03, 1,26, 1,58, 1,98, 2,49, 3,13
1,06, 1,88, 3,55, 6,70, 12,62, 23,80
1,12, 2,51, 6,32, 15,88, 39,90, 100,28
1,16, 3,14, 9,87, 31,01, 97,41, 306,02
Vector B:
-0,01, 0,61, 0,91, 0,99, 0,60, 0,02
Solution:
Coef. 1, rounding to 26 decimal, compare True, solution: -0,0597391027501962649904316335
Coef. 2, rounding to 26 decimal, compare True, solution: 1,8501896672627829700670299288
Coef. 3, rounding to 26 decimal, compare True, solution: -1,9727833018116428175300387317
Coef. 4, rounding to 26 decimal, compare True, solution: 1,4697587750651240151384675034
Coef. 5, rounding to 26 decimal, compare True, solution: -0,5538741847821888403564152897
Coef. 6, rounding to 26 decimal, compare True, solution: 0,0723048745759411900531809852
Algorithm using decimals
Matrix A:
1,00, 0,00, 0,00, 0,00, 0,00, 0,00
1,00, 0,63, 0,39, 0,25, 0,16, 0,10
1,00, 1,26, 1,58, 1,98, 2,49, 3,13
1,00, 1,88, 3,55, 6,70, 12,62, 23,80
1,00, 2,51, 6,32, 15,88, 39,90, 100,28
1,00, 3,14, 9,87, 31,01, 97,41, 306,02
Vector B:
-0,01, 0,61, 0,91, 0,99, 0,60, 0,02
Solution:
Coef. 1, rounding to 25 decimal, compare True, solution: -0,01
Coef. 2, rounding to 25 decimal, compare True, solution: 1,6027903945021139442641548525
Coef. 3, rounding to 25 decimal, compare True, solution: -1,6132030599055614189052834829
Coef. 4, rounding to 25 decimal, compare True, solution: 1,2454941213714367443882298102
Coef. 5, rounding to 25 decimal, compare True, solution: -0,4909897195846576129526569211
Coef. 6, rounding to 25 decimal, compare True, solution: 0,0657606961752320046201065486
Algorithm using pair of decimals as rational numbers
Matrix A:
1,00, 0,00, 0,00, 0,00, 0,00, 0,00
1,00, 0,63, 0,39, 0,25, 0,16, 0,10
1,00, 1,26, 1,58, 1,98, 2,49, 3,13
1,00, 1,88, 3,55, 6,70, 12,62, 23,80
1,00, 2,51, 6,32, 15,88, 39,90, 100,28
1,00, 3,14, 9,87, 31,01, 97,41, 306,02
Vector B:
-0,01, 0,61, 0,91, 0,99, 0,60, 0,02
Solution:
Coef. 1, rounding to 26 decimal, compare True, solution: -0,01
Coef. 2, rounding to 26 decimal, compare True, solution: 1,6027903945021139442641548522
Coef. 3, rounding to 26 decimal, compare True, solution: -1,6132030599055614189052834817
Coef. 4, rounding to 26 decimal, compare True, solution: 1,2454941213714367443882298085
Coef. 5, rounding to 26 decimal, compare True, solution: -0,4909897195846576129526569203
Coef. 6, rounding to 26 decimal, compare True, solution: 0,0657606961752320046201065485
Δεν υπάρχουν σχόλια:
Δημοσίευση σχολίου
You can feel free to write any suggestion, or idea on the subject.