package main

import (
	"fmt"
	"math"
)

const (
	alpha = 2.6  // change
	beta
	eps = 0.01
)

func getF(A [][]float64) [][]float64 {
	var F [][]float64
	for i := 0; i < len(A); i++ {
		var temp []float64
		for j := 0; j < len(A[i]); j++ {
			if i == j {
				temp = append(temp, 0.0)
			} else {
				temp = append(temp, (-1 * A[i][j] / A[i][i]))
			}
		}
		F = append(F, temp)
	}
	return F
}
func getc(b []float64) []float64 {
	var c []float64
	for i := 0; i < len(b); i++ {
		c = append(c, b[i]/alpha)
	}
	return c
}
func norm(F [][]float64) float64 {
	var max float64
	for i := 0; i < len(F); i++ {
		temp := 0.0
		for j := 0; j < len(F[i]); j++ {
			temp += math.Abs(F[i][j])
		}
		if i == 0 && temp < 1 {
			max = temp
		}
		if temp > max && temp < 1 {
			max = temp
		}
	}
	return max
}
func subvect(a, b []float64) []float64 {
	var res []float64
	for i := 0; i < len(a); i++ {
		res = append(res, a[i]-b[i])
	}
	return res
}
func sumvect(a, b []float64) []float64 {
	var res []float64
	for i := 0; i < len(a); i++ {
		res = append(res, a[i]+b[i])
	}
	return res
}
func normvect(in []float64) float64 {
	max := math.Abs(in[0])
	for i := 1; i < len(in); i++ {
		if math.Abs(in[i]) > max {
			max = math.Abs(in[i])
		}
	}
	return math.Abs(max)
}
func absolut(F [][]float64, xk, xk1 []float64) float64 {
	normF := norm(F)
	normX := normvect(subvect(xk, xk1))
	res := normF / (1 - normF) * normX
	return res
}
func otnos(absoluth float64, xk []float64) float64 {
	return absoluth / normvect(xk)
}
func pogr(x1, x2 []float64) float64 {
	normX := normvect(subvect(x1, x2))
	return normX
}
func mul_m(A [][]float64, x []float64) []float64 {
	var out []float64
	for i := 0; i < len(A); i++ {
		out = append(out, 0)
		for j := 0; j < len(A[i]); j++ {
			out[i] += A[i][j] * x[j]
		}
	}
	return out
}
func main() {

	A := [][]float64{{10.0 + alpha, -1.0, 0.2, 2.0},
		{1.0, 12.0 - alpha, -2.0, 0.1},
		{0.3, -4.0, 12.0 - alpha, 1.0},
		{0.2, -0.3, -0.5, 8.0 - alpha}}
	b := []float64{1.0 + beta, 2.0 - beta, 3.0, 1.0}
	F := getF(A)
	c := getc(b)
	var x [][]float64
	x = append(x, c)

	fmt.Println("F", F)
	fmt.Println("c", c)
	fmt.Println("|F|", norm(F))
	z:=0
	for {
		var tempX []float64
		for i := 0; i < 4; i++ {
			sumf := 0.0
			for j := 0; j < 4; j++ {
				sumf += F[i][j] * x[len(x)-1][j]
			}
			tempX = append(tempX, sumf+c[i])
		}
		x = append(x, tempX)
		abspogr := absolut(F, x[len(x)-1], x[len(x)-2])
		othospogr := otnos(abspogr, x[len(x)-1])
		fmt.Println("step",z)
		fmt.Println("absolute error", abspogr)
		fmt.Println("relative error", othospogr)
		if math.Abs(othospogr) < eps {
			break
		}
		z++
		//	fmt.Println(x)
	}


	fmt.Println("Iterations", x[len(x)-1])
	//Z
	var xDop [][]float64
	xDop = append(xDop, c)
	for {
		var tempX []float64
		for i := 0; i < 4; i++ {
			sumf1 := 0.0
			for j := 0; j <= i-1; j++ {
				sumf1 += F[i][j] * tempX[j]
			}

			sumf2 := 0.0
			for j := i; j < 4; j++ {
				sumf2 += F[i][j] * xDop[len(xDop)-1][j]
			}
			tempX = append(tempX, sumf1+sumf2+c[i])
		}
		xDop = append(xDop, tempX)
		//abspogr := absolut(F, xDop[len(xDop)-1], xDop[len(xDop)-2])
		//othospogr := otnos(abspogr, xDop[len(xDop)-1])
		//fmt.Println("absolut", abspogr)
		//fmt.Println("othos", othospogr)
		if pogr(xDop[len(xDop)-1], xDop[len(xDop)-2]) < eps {
			break
		}
		//	fmt.Println(x)
	}

	fmt.Println("Zeidel", xDop[len(xDop)-1])

	fmt.Println(sumvect(mul_m(F, x[len(x)-1]), c))
	fmt.Println(sumvect(mul_m(F, xDop[len(xDop)-1]), c))
}