aparler/RosettaCodeData
1
0
Fork 0
Code Issues 12 Pull requests 3 Projects Releases Packages Wiki Activity Actions

Update all new Tasks

Browse source
This commit is contained in:
Ingy döt Net 2015-02-20 09:02:09 -05:00
parent 00a190b0a6
commit 91df62d461
5697 changed files with 93386 additions and 804 deletions
Download patch file Download diff file Expand all files Collapse all files

1
Task/Gaussian-elimination/00DESCRIPTION Normal file
View file

@ -0,0 +1 @@
Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.

3
Task/Gaussian-elimination/00META.yaml Normal file
View file

@ -0,0 +1,3 @@
---
category:
- Matrices

26
Task/Gaussian-elimination/ALGOL-68/gaussian-elimination-1.alg Normal file
View file

@ -0,0 +1,26 @@
# -*- coding: utf-8 -*- #
COMMENT PROVIDES
MODE FIXED; INT fixed exception, unfixed exception;
PROC (STRING message) FIXED raise, raise value error
END COMMENT
# Note: ℵ indicates attribute is "private", and
should not be used outside of this prelude #
MODE FIXED = BOOL; # if an exception is detected, can it be fixed "on-site"? #
FIXED fixed exception = TRUE, unfixed exception = FALSE;
MODE #ℵ#SIMPLEOUTV = [0]UNION(CHAR, STRING, INT, REAL, BOOL, BITS);
MODE #ℵ#SIMPLEOUTM = [0]#ℵ#SIMPLEOUTV;
MODE #ℵ#SIMPLEOUTT = [0]#ℵ#SIMPLEOUTM;
MODE SIMPLEOUT = [0]#ℵ#SIMPLEOUTT;
PROC raise = (#ℵ#SIMPLEOUT message)FIXED: (
putf(stand error, ($"Exception:"$, $xg$, message, $l$));
stop
);
PROC raise value error = (#ℵ#SIMPLEOUT message)FIXED:
IF raise(message) NE fixed exception THEN exception value error; FALSE FI;
SKIP

53
Task/Gaussian-elimination/ALGOL-68/gaussian-elimination-2.alg Normal file
View file

@ -0,0 +1,53 @@
# -*- coding: utf-8 -*- #
COMMENT PRELUDE REQUIRES
MODE SCAL = REAL;
FORMAT scal repr = real repr
# and various SCAL OPerators #
END COMMENT
COMMENT PRELUDE PROIVIDES
MODE VEC, MAT;
OP :=:, -:=, +:=, *:=, /:=;
FORMAT sub, sep, bus;
FORMAT vec repr, mat repr
END COMMENT
# Note: ℵ indicates attribute is "private", and
should not be used outside of this prelude #
INT #ℵ#lwb vec := 1, #ℵ#upb vec := 0;
INT #ℵ#lwb mat := 1, #ℵ#upb mat := 0;
MODE VEC = [lwb vec:upb vec]SCAL,
MAT = [lwb mat:upb mat,lwb vec:upb vec]SCAL;
FORMAT sub := $"( "$, sep := $", "$, bus := $")"$, nl:=$lxx$;
FORMAT vec repr := $f(sub)n(upb vec - lwb vec)(f(scal repr)f(sep))f(scal repr)f(bus)$;
FORMAT mat repr := $f(sub)n(upb mat - lwb mat)(f( vec repr)f(nl))f( vec repr)f(bus)$;
# OPerators to swap the contents of two VECtors #
PRIO =:= = 1;
OP =:= = (REF VEC u, v)VOID:
FOR i TO UPB u DO SCAL scal=u[i]; u[i]:=v[i]; v[i]:=scal OD;
OP +:= = (REF VEC lhs, VEC rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] +:= rhs[i] OD;
lhs
);
OP -:= = (REF VEC lhs, VEC rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] -:= rhs[i] OD;
lhs
);
OP *:= = (REF VEC lhs, SCAL rhs)REF VEC: (
FOR i TO UPB lhs DO lhs[i] *:= rhs OD;
lhs
);
OP /:= = (REF VEC lhs, SCAL rhs)REF VEC: (
SCAL inv = 1 / rhs; # multiplication is faster #
FOR i TO UPB lhs DO lhs[i] *:= inv OD;
lhs
);
SKIP

78
Task/Gaussian-elimination/ALGOL-68/gaussian-elimination-3.alg Normal file
View file

@ -0,0 +1,78 @@
# -*- coding: utf-8 -*- #
COMMENT PRELUDE REQUIRES
MODE SCAL = REAL,
REAL near min scal = min real ** 0.99,
MODE VEC = []REAL,
MODE MAT = [,]REAL,
FORMAT scal repr = real repr,
and various OPerators of MAT and VEC
END COMMENT
COMMENT PRELUDE PROVIDES
PROC(MAT a, b)MAT gaussian elimination;
PROC(REF MAT a, b)REF MAT in situ gaussian elimination
END COMMENT
####################################################
# using Gaussian elimination, find x where A*x = b #
####################################################
PROC in situ gaussian elimination = (REF MAT a, b)REF MAT: (
# Note: a and b are modified "in situ", and b is returned as x #
FOR diag TO UPB a-1 DO
INT pivot row := diag; SCAL pivot factor := ABS a[diag,diag];
FOR row FROM diag + 1 TO UPB a DO # Full pivoting #
SCAL abs a diag = ABS a[row,diag];
IF abs a diag>=pivot factor THEN
pivot row := row; pivot factor := abs a diag FI
OD;
# now we have the "best" diag to full pivot, do the actual pivot #
IF diag NE pivot row THEN
# a[pivot row,] =:= a[diag,]; XXX: unoptimised # #DB#
a[pivot row,diag:] =:= a[diag,diag:]; # XXX: optimised #
b[pivot row,] =:= b[diag,] # swap/pivot the diags of a & b #
FI;
IF ABS a[diag,diag] <= near min scal THEN
raise value error("singular matrix") FI;
SCAL a diag reciprocal := 1 / a[diag, diag];
FOR row FROM diag+1 TO UPB a DO
SCAL factor = a[row,diag] * a diag reciprocal;
# a[row,] -:= factor * a[diag,] XXX: "unoptimised" # #DB#
a[row,diag+1:] -:= factor * a[diag,diag+1:];# XXX: "optimised" #
b[row,] -:= factor * b[diag,]
OD
OD;
# We have a triangular matrix, at this point we can traverse backwards
up the diagonal calculating b\A Converting it initial to a diagonal
matrix, then to the identity. #
FOR diag FROM UPB a BY -1 TO 1+LWB a DO
IF ABS a[diag,diag] <= near min scal THEN
raise value error("Zero pivot encountered?") FI;
SCAL a diag reciprocal = 1 / a[diag,diag];
FOR row TO diag-1 DO
SCAL factor = a[row,diag] * a diag reciprocal;
# a[row,diag] -:= factor * a[diag,diag]; XXX: "unoptimised" so remove # #DB#
b[row,] -:= factor * b[diag,]
OD;
# Now we have only diagonal elements we can simply divide b
by the values along the diagonal of A. #
b[diag,] *:= a diag reciprocal
OD;
b # EXIT #
);
PROC gaussian elimination = (MAT in a, in b)MAT: (
# Note: a and b are cloned and not modified "in situ" #
[UPB in a, 2 UPB in a]SCAL a := in a;
[UPB in b, 2 UPB in b]SCAL b := in b;
in situ gaussian elimination(a,b)
);
SKIP

9
Task/Gaussian-elimination/ALGOL-68/gaussian-elimination-4.alg Normal file
View file

@ -0,0 +1,9 @@
# -*- coding: utf-8 -*- #
COMMENT POSTLUDE PROIVIDES
PROC VOID exception too many iterations, exception value error;
END COMMENT
SKIP EXIT
exception too many iterations:
exception value error:
stop

29
Task/Gaussian-elimination/ALGOL-68/gaussian-elimination-5.alg Normal file
View file

@ -0,0 +1,29 @@
#!/usr/bin/algol68g-full --script #
# -*- coding: utf-8 -*- #
PR READ "prelude_exception.a68" PR;
# define the attributes of the scalar field being used #
MODE SCAL = REAL;
FORMAT scal repr = $g(-0,real width)$;
# create "near min scal" as is scales better then small real #
SCAL near min scal = min real ** 0.99;
PR READ "prelude_mat_lib.a68" PR;
PR READ "prelude_gaussian_elimination.a68" PR;
MAT 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));
VEC b = (-0.01, 0.61, 0.91, 0.99, 0.60, 0.02);
[UPB b,1]SCAL col b; col b[,1]:= b;
upb vec := 2 UPB a;
printf((vec repr, gaussian elimination(a,col b)));
PR READ "postlude_exception.a68" PR

73
Task/Gaussian-elimination/C/gaussian-elimination.c Normal file
View file

@ -0,0 +1,73 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define mat_elem(a, y, x, n) (a + ((y) * (n) + (x)))
void swap_row(double *a, double *b, int r1, int r2, int n)
{
double tmp, *p1, *p2;
int i;
if (r1 == r2) return;
for (i = 0; i < n; i++) {
p1 = mat_elem(a, r1, i, n);
p2 = mat_elem(a, r2, i, n);
tmp = *p1, *p1 = *p2, *p2 = tmp;
}
tmp = b[r1], b[r1] = b[r2], b[r2] = tmp;
}
void gauss_eliminate(double *a, double *b, double *x, int n)
{
#define A(y, x) (*mat_elem(a, y, x, n))
int i, j, col, row, max_row,dia;
double max, tmp;
for (dia = 0; dia < n; dia++) {
max_row = dia, max = A(dia, dia);
for (row = dia + 1; row < n; row++)
if ((tmp = fabs(A(row, dia))) > max)
max_row = row, max = tmp;
swap_row(a, b, dia, max_row, n);
for (row = dia + 1; row < n; row++) {
tmp = A(row, dia) / A(dia, dia);
for (col = dia+1; col < n; col++)
A(row, col) -= tmp * A(dia, col);
A(row, dia) = 0;
b[row] -= tmp * b[dia];
}
}
for (row = n - 1; row >= 0; row--) {
tmp = b[row];
for (j = n - 1; j > row; j--)
tmp -= x[j] * A(row, j);
x[row] = tmp / A(row, row);
}
#undef A
}
int main(void)
{
double 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
};
double b[] = { -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 };
double x[6];
int i;
gauss_eliminate(a, b, x, 6);
for (i = 0; i < 6; i++)
printf("%g\n", x[i]);
return 0;
}

69
Task/Gaussian-elimination/D/gaussian-elimination.d Normal file
View file

@ -0,0 +1,69 @@
import std.stdio, std.math, std.algorithm, std.range, std.numeric,
std.typecons;
Tuple!(double[],"x", string,"err")
gaussPartial(in double[][] a0, in double[] b0) pure /*nothrow*/
in {
assert(a0.length == a0[0].length);
assert(a0.length == b0.length);
assert(a0.all!(row => row.length == a0[0].length));
} body {
enum eps = 1e-6;
immutable m = b0.length;
// Make augmented matrix.
//auto a = a0.zip(b0).map!(c => c[0] ~ c[1]).array; // Not mutable.
auto a = a0.zip(b0).map!(c => [] ~ c[0] ~ c[1]).array;
// Wikipedia algorithm from Gaussian elimination page,
// produces row-eschelon form.
foreach (immutable k; 0 .. a.length) {
// Find pivot for column k and swap.
a[k .. m].minPos!((x, y) => x[k] > y[k]).front.swap(a[k]);
if (a[k][k].abs < eps)
return typeof(return)(null, "singular");
// Do for all rows below pivot.
foreach (immutable i; k + 1 .. m) {
// Do for all remaining elements in current row.
a[i][k+1 .. m+1] -= a[k][k+1 .. m+1] * (a[i][k] / a[k][k]);
a[i][k] = 0; // Fill lower triangular matrix with zeros.
}
}
// End of WP algorithm. Now back substitute to get result.
auto x = new double[m];
foreach_reverse (immutable i; 0 .. m)
x[i] = (a[i][m] - a[i][i+1 .. m].dotProduct(x[i+1 .. m])) / a[i][i];
return typeof(return)(x, null);
}
void main() {
// The test case result is correct to this tolerance.
enum eps = 1e-13;
// Common RC example. Result computed with rational arithmetic
// then converted to double, and so should be about as close to
// correct as double represention allows.
immutable 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]];
immutable b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02];
immutable r = gaussPartial(a, b);
if (!r.err.empty)
return writeln("Error: ", r.err);
r.x.writeln;
immutable result = [-0.01, 1.602790394502114,
-1.6132030599055613, 1.2454941213714368,
-0.4909897195846576, 0.065760696175232];
foreach (immutable i, immutable xi; result)
if (abs(r.x[i] - xi) > eps)
return writeln("Out of tolerance: ", r.x[i], " ", xi);
}

43
Task/Gaussian-elimination/Fortran/gaussian-elimination.f Normal file
View file

@ -0,0 +1,43 @@
program ge
real, allocatable :: a(:,:),b(:)
a = reshape( &
[1.0, 1.00, 1.00, 1.00, 1.00, 1.00, &
0.0, 0.63, 1.26, 1.88, 2.51, 3.14, &
0.0, 0.39, 1.58, 3.55, 6.32, 9.87, &
0.0, 0.25, 1.98, 6.70, 15.88, 31.01, &
0.0, 0.16, 2.49, 12.62, 39.90, 97.41, &
0.0, 0.10, 3.13, 23.80, 100.28, 306.02], [6,6] )
b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]
print'(f15.7)',solve_wbs(ge_wpp(a,b))
contains
function solve_wbs(u) result(x) ! solve with backward substitution
real :: u(:,:)
integer :: i,n
real , allocatable :: x(:)
n = size(u,1)
allocate(x(n))
forall (i=n:1:-1) x(i) = ( u(i,n+1) - sum(u(i,i+1:n)*x(i+1:n)) ) / u(i,i)
end function
function ge_wpp(a,b) result(u) ! gaussian eliminate with partial pivoting
real :: a(:,:),b(:),upi
integer :: i,j,n,p
real , allocatable :: u(:,:)
n = size(a,1)
u = reshape( [a,b], [n,n+1] )
do j=1,n
p = maxloc(abs(u(j:n,j)),1) + j-1 ! maxloc returns indices between (1,n-j+1)
if (p /= j) u([p,j],j) = u([j,p],j)
u(j+1:,j) = u(j+1:,j)/u(j,j)
do i=j+1,n+1
upi = u(p,i)
if (p /= j) u([p,j],i) = u([j,p],i)
u(j+1:n,i) = u(j+1:n,i) - upi*u(j+1:n,j)
end do
end do
end function
end program

97
Task/Gaussian-elimination/Go/gaussian-elimination.go Normal file
View file

@ -0,0 +1,97 @@
package main
import (
"errors"
"fmt"
"log"
"math"
)
type testCase struct {
a [][]float64
b []float64
x []float64
}
var tc = testCase{
// common RC example. Result x computed with rational arithmetic then
// converted to float64, and so should be about as close to correct as
// float64 represention allows.
a: [][]float64{
{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}},
b: []float64{-0.01, 0.61, 0.91, 0.99, 0.60, 0.02},
x: []float64{-0.01, 1.602790394502114, -1.6132030599055613,
1.2454941213714368, -0.4909897195846576, 0.065760696175232},
}
// result from above test case turns out to be correct to this tolerance.
const ε = 1e-13
func main() {
x, err := GaussPartial(tc.a, tc.b)
if err != nil {
log.Fatal(err)
}
fmt.Println(x)
for i, xi := range x {
if math.Abs(tc.x[i]-xi) > ε {
log.Println("out of tolerance")
log.Fatal("expected", tc.x)
}
}
}
func GaussPartial(a0 [][]float64, b0 []float64) ([]float64, error) {
// make augmented matrix
m := len(b0)
a := make([][]float64, m)
for i, ai := range a0 {
row := make([]float64, m+1)
copy(row, ai)
row[m] = b0[i]
a[i] = row
}
// WP algorithm from Gaussian elimination page
// produces row-eschelon form
for k := range a {
// Find pivot for column k:
iMax := k
max := math.Abs(a[k][k])
for i := k + 1; i < m; i++ {
if abs := math.Abs(a[i][k]); abs > max {
iMax = i
max = abs
}
}
if a[iMax][k] == 0 {
return nil, errors.New("singular")
}
// swap rows(k, i_max)
a[k], a[iMax] = a[iMax], a[k]
// Do for all rows below pivot:
for i := k + 1; i < m; i++ {
// Do for all remaining elements in current row:
for j := k + 1; j <= m; j++ {
a[i][j] -= a[k][j] * (a[i][k] / a[k][k])
}
// Fill lower triangular matrix with zeros:
a[i][k] = 0
}
}
// end of WP algorithm.
// now back substitute to get result.
x := make([]float64, m)
for i := m - 1; i >= 0; i-- {
x[i] = a[i][m]
for j := i + 1; j < m; j++ {
x[i] -= a[i][j] * x[j]
}
x[i] /= a[i][i]
}
return x, nil
}

57
Task/Gaussian-elimination/J/gaussian-elimination.j Normal file
View file

@ -0,0 +1,57 @@
f=: 6j2&": NB. formatting verb
sin=: 1&o. NB. verb to evaluate circle function 1, the sine
add_noise=: ] + (* (_0.5 + 0 ?@:#~ #)) NB. AMPLITUDE add_noise SIGNAL
f RADIANS=: o.@:(%~ i.@:>:)5 NB. monadic circle function is pi times
0.00 0.63 1.26 1.88 2.51 3.14
f SINES=: sin RADIANS
0.00 0.59 0.95 0.95 0.59 0.00
f NOISY_SINES=: 0.1 add_noise SINES
_0.01 0.61 0.91 0.99 0.60 0.02
A=: (^/ i.@:#) RADIANS NB. A is the quintic coefficient matrix
NB. display the equation to solve
(f A) ; 'x' ; '=' ; f@:,. NOISY_SINES
┌────────────────────────────────────┬─┬─┬──────┐
│ 1.00 0.00 0.00 0.00 0.00 0.00│x│=│ _0.01│
│ 1.00 0.63 0.39 0.25 0.16 0.10│ │ │ 0.61│
│ 1.00 1.26 1.58 1.98 2.49 3.13│ │ │ 0.91│
│ 1.00 1.88 3.55 6.70 12.62 23.80│ │ │ 0.99│
│ 1.00 2.51 6.32 15.88 39.90100.28│ │ │ 0.60│
│ 1.00 3.14 9.87 31.01 97.41306.02│ │ │ 0.02│
└────────────────────────────────────┴─┴─┴──────┘
f QUINTIC_COEFFICIENTS=: NOISY_SINES %. A NB. %. solves the linear system
_0.01 1.71 _1.88 1.48 _0.58 0.08
quintic=: QUINTIC_COEFFICIENTS&p. NB. verb to evaluate the polynomial
NB. %. also solves the least squares fit for overdetermined system
quadratic=: (NOISY_SINES %. (^/ i.@:3:) RADIANS)&p. NB. verb to evaluate quadratic.
quadratic
_0.0200630695393961729 1.26066877804926536 _0.398275112136019516&p.
NB. The quintic is agrees with the noisy data, as it should
f@:(NOISY_SINES ,. sin ,. quadratic ,. quintic) RADIANS
_0.01 0.00 _0.02 _0.01
0.61 0.59 0.61 0.61
0.91 0.95 0.94 0.91
0.99 0.95 0.94 0.99
0.60 0.59 0.63 0.60
0.02 0.00 0.01 0.02
f MID_POINTS=: (+ -:@:(-/@:(2&{.)))RADIANS
_0.31 0.31 0.94 1.57 2.20 2.83
f@:(sin ,. quadratic ,. quintic) MID_POINTS
_0.31 _0.46 _0.79
0.31 0.34 0.38
0.81 0.81 0.77
1.00 0.98 1.00
0.81 0.83 0.86
0.31 0.36 0.27

1
Task/Gaussian-elimination/Julia/gaussian-elimination.julia Normal file
View file

@ -0,0 +1 @@
x = A \ b

70
Task/Gaussian-elimination/MATLAB/gaussian-elimination.m Normal file
View file

@ -0,0 +1,70 @@
function [ x ] = GaussElim( A, b)
% Ensures A is n by n
sz = size(A);
if sz(1)~=sz(2)
fprintf('A is not n by n\n');
clear x;
return;
end
n = sz(1);
% Ensures b is n x 1.
if n~=sz(1)
fprintf('b is not 1 by n.\n');
return
end
x = zeros(n,1);
aug = [A b];
tempmatrix = aug;
for i=2:sz(1)
% Find maximum of row and divide by the maximum
tempmatrix(1,:) = tempmatrix(1,:)/max(tempmatrix(1,:));
% Finds the maximum in column
temp = find(abs(tempmatrix) - max(abs(tempmatrix(:,1))));
if length(temp)>2
for j=1:length(temp)-1
if j~=temp(j)
maxi = j; %maxi = column number of maximum
break;
end
end
else % length(temp)==2
maxi=1;
end
% Row swap if maxi is not 1
if maxi~=1
temp = tempmatrix(maxi,:);
tempmatrix(maxi,:) = tempmatrix(1,:);
tempmatrix(1,:) = temp;
end
% Row reducing
for j=2:length(tempmatrix)-1
tempmatrix(j,:) = tempmatrix(j,:)-tempmatrix(j,1)/tempmatrix(1,1)*tempmatrix(1,:);
if tempmatrix(j,j)==0 || isnan(tempmatrix(j,j)) || abs(tempmatrix(j,j))==Inf
fprintf('Error: Matrix is singular.\n');
clear x;
return
end
end
aug(i-1:end,i-1:end) = tempmatrix;
% Decrease matrix size
tempmatrix = tempmatrix(2:end,2:end);
end
% Backwards Substitution
x(end) = aug(end,end)/aug(end,end-1);
for i=n-1:-1:1
x(i) = (aug(i,end)-dot(aug(i,1:end-1),x))/aug(i,i);
end
end

1
Task/Gaussian-elimination/Mathematica/gaussian-elimination.math Normal file
View file

@ -0,0 +1 @@
GaussianElimination[A_?MatrixQ, b_?VectorQ] := Last /@ RowReduce[Flatten /@ Transpose[{A, b}]]

28
Task/Gaussian-elimination/OCaml/gaussian-elimination-1.ocaml Normal file
View file

@ -0,0 +1,28 @@
module Array = struct
include Array
(* Computes: f a.(0) + f a.(1) + ... where + is 'g'. *)
let foldmap g f a =
let n = Array.length a in
let rec aux acc i =
if i >= n then acc else aux (g acc (f a.(i))) (succ i)
in aux (f a.(0)) 1
(* like the stdlib fold_left, but also provides index to f *)
let foldi_left f x a =
let r = ref x in
for i = 0 to length a - 1 do
r := f i !r (unsafe_get a i)
done;
!r
end
let foldmap_range g f (a,b) =
let rec aux acc n =
let n = succ n in
if n > b then acc else aux (g acc (f n)) n
in aux (f a) a
let fold_range f init (a,b) =
let rec aux acc n =
if n > b then acc else aux (f acc n) (succ n)
in aux init a

38
Task/Gaussian-elimination/OCaml/gaussian-elimination-2.ocaml Normal file
View file

@ -0,0 +1,38 @@
(* Some less-general support functions for 'solve'. *)
let swap_elem m i j = let x = m.(i) in m.(i) <- m.(j); m.(j) <- x
let maxtup a b = if (snd a) > (snd b) then a else b
let augmented_matrix m b =
Array.(init (length m) ( fun i -> append m.(i) [|b.(i)|] ))
(* Solve Ax=b for x, using gaussian elimination with scaled partial pivot,
* and then back-substitution of the resulting row-echelon matrix. *)
let solve m b =
let n = Array.length m in
let n' = pred n in (* last index = n-1 *)
let s = Array.(map (foldmap max abs_float) m) in (* scaling vector *)
let a = augmented_matrix m b in
for k = 0 to pred n' do
(* Scaled partial pivot, to preserve precision *)
let pair i = (i, abs_float a.(i).(k) /. s.(i)) in
let i_max,v = foldmap_range maxtup pair (k,n') in
if v < epsilon_float then failwith "Matrix is singular.";
swap_elem a k i_max;
swap_elem s k i_max;
(* Eliminate one column *)
for i = succ k to n' do
let tmp = a.(i).(k) /. a.(k).(k) in
for j = succ k to n do
a.(i).(j) <- a.(i).(j) -. tmp *. a.(k).(j);
done
done
done;
(* Backward substitution; 'b' is in the 'nth' column of 'a' *)
let x = Array.copy b in (* just a fresh array of the right size and type *)
for i = n' downto 0 do
let minus_dprod t j = t -. x.(j) *. a.(i).(j) in
x.(i) <- fold_range minus_dprod a.(i).(n) (i+1,n') /. a.(i).(i);
done;
x

8
Task/Gaussian-elimination/OCaml/gaussian-elimination-3.ocaml Normal file
View file

@ -0,0 +1,8 @@
let 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 |] |]
let b = [| -0.01; 0.61; 0.91; 0.99; 0.60; 0.02 |]

4
Task/Gaussian-elimination/OCaml/gaussian-elimination-4.ocaml Normal file
View file

@ -0,0 +1,4 @@
# let x = solve a b;;
val x : float array =
[|-0.0100000000000000991; 1.60279039450210536; -1.61320305990553226;
1.24549412137140547; -0.490989719584644546; 0.0657606961752301433|]

8
Task/Gaussian-elimination/OCaml/gaussian-elimination-5.ocaml Normal file
View file

@ -0,0 +1,8 @@
let mul m v =
Array.mapi (fun i u ->
Array.foldi_left (fun j sum uj ->
sum +. uj *. v.(j)
) 0. u
) m
let sub u v = Array.mapi (fun i e -> e -. v.(i)) u

4
Task/Gaussian-elimination/OCaml/gaussian-elimination-6.ocaml Normal file
View file

@ -0,0 +1,4 @@
# let residual = sub b (mul a x);;
val residual : float array =
[|9.8879238130678e-17; 1.11022302462515654e-16; 2.22044604925031308e-16;
8.88178419700125232e-16; -5.5511151231257827e-16; 4.26741975090294545e-16|]

1
Task/Gaussian-elimination/PARI-GP/gaussian-elimination.pari Normal file
View file

@ -0,0 +1 @@
matsolve(A,B)

75
Task/Gaussian-elimination/PHP/gaussian-elimination.php Normal file
View file

@ -0,0 +1,75 @@
function swap_rows(&$a, &$b, $r1, $r2)
{
if ($r1 == $r2) return;
$tmp = $a[$r1];
$a[$r1] = $a[$r2];
$a[$r2] = $tmp;
$tmp = $b[$r1];
$b[$r1] = $b[$r2];
$b[$r2] = $tmp;
}
function gauss_eliminate($A, $b, $N)
{
for ($col = 0; $col < $N; $col++)
{
$j = $col;
$max = $A[$j][$j];
for ($i = $col + 1; $i < $N; $i++)
{
$tmp = abs($A[$i][$col]);
if ($tmp > $max)
{
$j = $i;
$max = $tmp;
}
}
swap_rows($A, $b, $col, $j);
for ($i = $col + 1; $i < $N; $i++)
{
$tmp = $A[$i][$col] / $A[$col][$col];
for ($j = $col + 1; $j < $N; $j++)
{
$A[$i][$j] -= $tmp * $A[$col][$j];
}
$A[$i][$col] = 0;
$b[$i] -= $tmp * $b[$col];
}
}
$x = array();
for ($col = $N - 1; $col >= 0; $col--)
{
$tmp = $b[$col];
for ($j = $N - 1; $j > $col; $j--)
{
$tmp -= $x[$j] * $A[$col][$j];
}
$x[$col] = $tmp / $A[$col][$col];
}
return $x;
}
function test_gauss()
{
$a = array(
array(1.00, 0.00, 0.00, 0.00, 0.00, 0.00),
array(1.00, 0.63, 0.39, 0.25, 0.16, 0.10),
array(1.00, 1.26, 1.58, 1.98, 2.49, 3.13),
array(1.00, 1.88, 3.55, 6.70, 12.62, 23.80),
array(1.00, 2.51, 6.32, 15.88, 39.90, 100.28),
array(1.00, 3.14, 9.87, 31.01, 97.41, 306.02)
);
$b = array( -0.01, 0.61, 0.91, 0.99, 0.60, 0.02 );
$x = gauss_eliminate($a, $b, 6);
ksort($x);
print_r($x);
}
test_gauss();

71
Task/Gaussian-elimination/PL-I/gaussian-elimination.pli Normal file
View file

@ -0,0 +1,71 @@
Solve: procedure options (main); /* 11 January 2014 */
declare n fixed binary;
put ('Program to solve n simultaneous equations of the form Ax = b. Please type n:' );
get (n);
begin;
declare (A(n, n), b(n), x(n)) float(18);
declare (SA(n,n), Sb(n)) float (18);
declare i fixed binary;
put skip list ('Please type A:');
get (a);
put skip list ('Please type the right-hand sides, b:');
get (b);
SA = A; Sb = b;
put skip list ('The equations are:');
do i = 1 to n;
put skip edit (A(i,*), b(i)) (f(5), x(1));
end;
call Gauss_elimination (A, b);
call Backward_substitution (A, b, x);
put skip list ('Solutions:'); put skip data (x);
/* Check solutions: */
put skip list ('Residuals:');
do i = 1 to n;
put skip list (sum(SA(i,*) * x(*)) - Sb(i));
end;
end;
Gauss_elimination: procedure (A, b) options (reorder); /* Triangularise */
declare (A(*,*), b(*)) float(18);
declare n fixed binary initial (hbound(A, 1));
declare (i, j, k) fixed binary;
declare t float(18);
do j = 1 to n;
do i = j+1 to n; /* For each of the rows beneath the current (pivot) row. */
t = A(j,j) / A(i,j);
do k = j+1 to n; /* Subtract a multiple of row i from row j. */
A(i,k) = A(j,k) - t*A(i,k);
end;
b(i) = b(j) - t*b(i); /* ... and the right-hand side. */
end;
end;
end Gauss_elimination;
Backward_substitution: procedure (A, b, x) options (reorder);
declare (A(*,*), b(*), x(*)) float(18);
declare t float(18);
declare n fixed binary initial (hbound(A, 1));
declare (i, j) fixed binary;
x(n) = b(n) / a(n,n);
do j = n-1 to 1 by -1;
t = 0;
do i = j+1 to n;
t = t + a(j,i)*x(i);
end;
x(j) = (b(j) - t) / a(j,j);
end;
end Backward_substitution;
end Solve;

54
Task/Gaussian-elimination/Perl-6/gaussian-elimination.pl6 Normal file
View file

@ -0,0 +1,54 @@
sub mat_elem(@a, $y, $x, $n) is rw { @a[ $y * $n + $x ] }
sub swap_row(@a, @b, $r1, $r2, $n) {
return if $r1 == $r2;
for ^$n -> $i {
(
mat_elem(@a, $r1, $i, $n),
mat_elem(@a, $r2, $i, $n)
).=reverse;
}
@b[$r1, $r2].=reverse;
}
sub gauss_eliminate(@a, @b, $n) {
sub A($y, $x) is rw { mat_elem(@a, $y, $x, $n) }
my ($i, $j, $col, $row, $max_row, $dia);
my ($max, $tmp);
for ^$n -> $dia {
for $dia ^..^ $n -> $row {
swap_row @a, @b, $dia,
max(:by({ abs(A($_, $dia)) }), $dia ^..^ $n),
$n;
$tmp = A($row, $dia) / A($dia, $dia);
for $dia ^..^ $n -> $col {
A($row, $col) -= $tmp * A($dia, $col);
}
A($row, $dia) = 0;
@b[$row] -= $tmp * @b[$dia];
}
}
my @x;
for $n - 1, $n - 2 ... 0 -> $row {
$tmp = @b[$row];
for $n - 1, $n - 2 ...^ $row -> $j {
$tmp -= @x[$j] * A($row, $j);
}
@x[$row] = $tmp / A($row, $row);
}
return @x;
}
sub MAIN {
my @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
>;
my @b = <
-0.01 0.61 0.91 0.99 0.60 0.02
>;
.say for gauss_eliminate(@a, @b, 6);
}

9
Task/Gaussian-elimination/Perl/gaussian-elimination.pl Normal file
View file

@ -0,0 +1,9 @@
use Math::Matrix;
my $a = Math::Matrix->new([0,1,0],
[0,0,1],
[2,0,1]);
my $b = Math::Matrix->new([1],
[2],
[4]);
my $x = $a->concat($b)->solve;
print $x;

96
Task/Gaussian-elimination/Python/gaussian-elimination.py Normal file
View file

@ -0,0 +1,96 @@
# The 'gauss' function takes two matrices, 'a' and 'b', with 'a' square, and it return the determinant of 'a' and a matrix 'x' such that a*x = b.
# If 'b' is the identity, then 'x' is the inverse of 'a'.
import copy
from fractions import Fraction
def gauss(a, b):
a = copy.deepcopy(a)
b = copy.deepcopy(b)
n = len(a)
p = len(b[0])
det = 1
for i in range(n - 1):
k = i
for j in range(i + 1, n):
if abs(a[j][i]) > abs(a[k][i]):
k = j
if k != i:
a[i], a[k] = a[k], a[i]
b[i], b[k] = b[k], b[i]
det = -det
for j in range(i + 1, n):
t = a[j][i]/a[i][i]
for k in range(i + 1, n):
a[j][k] -= t*a[i][k]
for k in range(p):
b[j][k] -= t*b[i][k]
for i in range(n - 1, -1, -1):
for j in range(i + 1, n):
t = a[i][j]
for k in range(p):
b[i][k] -= t*b[j][k]
t = 1/a[i][i]
det *= a[i][i]
for j in range(p):
b[i][j] *= t
return det, b
def zeromat(p, q):
return [[0]*q for i in range(p)]
def matmul(a, b):
n, p = len(a), len(a[0])
p1, q = len(b), len(b[0])
if p != p1:
raise ValueError("Incompatible dimensions")
c = zeromat(n, q)
for i in range(n):
for j in range(q):
c[i][j] = sum(a[i][k]*b[k][j] for k in range(p))
return c
def mapmat(f, a):
return [list(map(f, v)) for v in a]
def ratmat(a):
return mapmat(Fraction, a)
# As an example, compute the determinant and inverse of 3x3 magic square
a = [[2, 9, 4], [7, 5, 3], [6, 1, 8]]
b = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
det, c = gauss(a, b)
det
-360.0
c
[[-0.10277777777777776, 0.18888888888888888, -0.019444444444444438],
[0.10555555555555554, 0.02222222222222223, -0.061111111111111116],
[0.0638888888888889, -0.14444444444444446, 0.14722222222222223]]
# Check product
matmul(a, c)
[[1.0, 0.0, 0.0], [5.551115123125783e-17, 1.0, 0.0],
[1.1102230246251565e-16, -2.220446049250313e-16, 1.0]]
# Same with fractions, so the result is exact
det, c = gauss(ratmat(a), ratmat(b))
det
Fraction(-360, 1)
c
[[Fraction(-37, 360), Fraction(17, 90), Fraction(-7, 360)],
[Fraction(19, 180), Fraction(1, 45), Fraction(-11, 180)],
[Fraction(23, 360), Fraction(-13, 90), Fraction(53, 360)]]
matmul(a, c)
[[Fraction(1, 1), Fraction(0, 1), Fraction(0, 1)],
[Fraction(0, 1), Fraction(1, 1), Fraction(0, 1)],
[Fraction(0, 1), Fraction(0, 1), Fraction(1, 1)]]

1
Task/Gaussian-elimination/README Normal file
View file

@ -0,0 +1 @@
Data source: http://rosettacode.org/wiki/Gaussian_elimination

79
Task/Gaussian-elimination/REXX/gaussian-elimination-1.rexx Normal file
View file

@ -0,0 +1,79 @@
/* REXX ---------------------------------------------------------------
* 07.08.2014 Walter Pachl translated from PL/I)
*--------------------------------------------------------------------*/
Numeric Digits 20
Parse Arg t
n=3
Parse Value '1 2 3 14' With a.1.1 a.1.2 a.1.3 b.1
Parse Value '2 1 3 13' With a.2.1 a.2.2 a.2.3 b.2
Parse Value '3 -2 -1 -4' With a.3.1 a.3.2 a.3.3 b.3
If t=6 Then Do
n=6
Parse Value '1.00 0.00 0.00 0.00 0.00 0.00 ' With a.1.1 a.1.2 a.1.3 a.1.4 a.1.5 a.1.6 .
Parse Value '1.00 0.63 0.39 0.25 0.16 0.10 ' With a.2.1 a.2.2 a.2.3 a.2.4 a.2.5 a.2.6 .
Parse Value '1.00 1.26 1.58 1.98 2.49 3.13 ' With a.3.1 a.3.2 a.3.3 a.3.4 a.3.5 a.3.6 .
Parse Value '1.00 1.88 3.55 6.70 12.62 23.80 ' With a.4.1 a.4.2 a.4.3 a.4.4 a.4.5 a.4.6 .
Parse Value '1.00 2.51 6.32 15.88 39.90 100.28' With a.5.1 a.5.2 a.5.3 a.5.4 a.5.5 a.5.6 .
Parse Value '1.00 3.14 9.87 31.01 97.41 306.02' With a.6.1 a.6.2 a.6.3 a.6.4 a.6.5 a.6.6 .
Parse Value '-0.01 0.61 0.91 0.99 0.60 0.02' With b.1 b.2 b.3 b.4 b.5 b.6 .
End
Do i=1 To n
Do j=1 To n
sa.i.j=a.i.j
End
sb.i=b.i
End
Say 'The equations are:'
do i = 1 to n;
ol=''
Do j=1 To n
ol=ol format(a.i.j,4,4)
End
ol=ol' 'format(b.i,4,4)
Say ol
end
call Gauss_elimination
call Backward_substitution
Say 'Solutions:'
Do i=1 To n
Say 'x('i')='||x.i
End
/* Check solutions: */
Say 'Residuals:'
do i = 1 to n
res=0
Do j=1 To n
res=res+(sa.i.j*x.j)
End
res=res-sb.i
Say 'res('i')='res
End
Exit
Gauss_elimination:
do j = 1 to n
do i = j+1 to n /* For each of the rows beneath the current (pivot) row. */
t = a.j.j / a.i.j
do k = j+1 to n /* Subtract a multiple of row i from row j. */
a.i.k = a.j.k - t*a.i.k
end
b.i = b.j - t*b.i /* ... and the right-hand side. */
end
end
Return
Backward_substitution:
x.n = b.n / a.n.n
do j = n-1 to 1 by -1
t = 0
do i = j+1 to n
t = t + a.j.i*x.i
end
x.j = (b.j - t) / a.j.j
end
Return

35
Task/Gaussian-elimination/REXX/gaussian-elimination-2.rexx Normal file
View file

@ -0,0 +1,35 @@
/*REXX pgm solves Ax=b with Gaussian elimination &backwards substitution*/
parse arg iFID .; if iFID=='' then iFID='GAUSS_E.DAT'
do files=1 while lines(iFID)\==0; #=0 /*read equation sets*/
numeric digits 200 /*reduces rounding. */
do $=1 while lines(iFID)\==0 /*process equation. */
z=linein(iFID); if z='' then leave /*Blank line? e─o─d.*/
if $==1 then do; say; say center(' equations ',75,'▒'); say; end
say z /*show an equation. */
if left(space(z),1)=='*' then iterate /*ignore comments. */
#=#+1; n=words(z)-1 /*assign equation #s*/
do e=1 for n; a.#.e=word(z,e); end /*e*/ /*process the A #s.*/
b.#=word(z,n+1) /* " " B #.*/
end /*$*/
if #==0 then iterate /*extra blank line? */
call Gauss_elimination /*invoke Gauss elim.*/
say; say center('solution',75,'═'); say /* [↓] show solution*/
numeric digits 8; pad=left('',30) /*only show 8 digits*/
do s=1 for n; say pad 'x['s"] = " left('', x.s>=0) x.s/1; end
end /*files*/
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────GAUSS_ELIMINATION subroutine────────*/
Gauss_elimination: do j=1 for n; jp=j+1
do i=jp to n; _=a.j.j/a.i.j
do k=jp to n; a.i.k=a.j.k-_*a.i.k; end /*k*/
b.i=b.j-_*b.i
end /*i*/
end /*j*/
x.n=b.n/a.n.n
do j=n-1 to 1 by -1; _=0
do i=j+1 to n; _=_+a.j.i*x.i; end /*i*/
x.j=(b.j-_)/a.j.j
end /*j*/ /* [↑] backwards substitution.*/
return

13
Task/Gaussian-elimination/Racket/gaussian-elimination-1.rkt Normal file
View file

@ -0,0 +1,13 @@
#lang racket
(require math/matrix)
(define A
(matrix [[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]]))
(define b (col-matrix [-0.01 0.61 0.91 0.99 0.60 0.02]))
(matrix-solve A b)

8
Task/Gaussian-elimination/Racket/gaussian-elimination-2.rkt Normal file
View file

@ -0,0 +1,8 @@
#<array
'#(6 1)
#[-0.01
1.602790394502109
-1.613203059905556
1.2454941213714346
-0.4909897195846582
0.06576069617523222]>

18
Task/Gaussian-elimination/Ruby/gaussian-elimination.rb Normal file
View file

@ -0,0 +1,18 @@
require 'bigdecimal/ludcmp'
include LUSolve
BigDecimal::limit(30)
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].map{|i|BigDecimal(i,16)}
b = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02].map{|i|BigDecimal(i,16)}
n = 6
zero = BigDecimal("0.0")
one = BigDecimal("1.0")
lusolve(a, b, ludecomp(a, n, zero,one), zero).each{|v| puts v.to_s('F')[0..20]}

12
Task/Gaussian-elimination/Tcl/gaussian-elimination.tcl Normal file
View file

@ -0,0 +1,12 @@
package require math::linearalgebra
set 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}
}
set b {-0.01 0.61 0.91 0.99 0.60 0.02}
puts -nonewline [math::linearalgebra::show [math::linearalgebra::solveGauss $A $b] "%.2f"]