-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.go
97 lines (91 loc) · 1.82 KB
/
main.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
// Package infinity provides functions for working with infinite numbers.
package infinity
import (
"math"
"github.com/ready-steady/linear/matrix"
)
var (
Plus = math.Inf(1.0) // +∞
Minus = math.Inf(-1.0) // -∞
)
// Linear multiplies an m-by-n matrix by an n-element vector. The elements of
// the vector can potentially be infinite.
func Linear(A, x []float64, m, n uint) []float64 {
y := make([]float64, m)
ok, s := inspect(x, n)
if ok {
matrix.Multiply(A, x, y, m, n, 1)
return y
}
for i := uint(0); i < m; i++ {
fin, inf := 0.0, 0.0
for j := uint(0); j < n; j++ {
a := A[j*m+i]
if a == 0.0 {
continue
}
if s[j] == 0.0 {
fin += a * x[j]
} else {
inf += a * s[j]
}
}
if inf != 0.0 {
y[i] = inf * Plus
} else {
y[i] = fin
}
}
return y
}
// Quadratic multiplies an m-by-m matrix by an m-element vector from both sides.
// The elements of the vector can potentially be infinite.
func Quadratic(A, x []float64, m uint) float64 {
ok, s := inspect(x, m)
if ok {
y := make([]float64, m)
matrix.Multiply(A, x, y, m, m, 1)
return matrix.Dot(x, y, m)
}
Fin, Inf, INF := 0.0, 0.0, 0.0
for i := uint(0); i < m; i++ {
fin, inf := 0.0, 0.0
for j := uint(0); j < m; j++ {
a := A[j*m+i]
if a == 0.0 {
continue
}
if s[j] == 0.0 {
fin += a * x[j]
} else {
inf += a * s[j]
}
}
if s[i] == 0.0 {
Fin += x[i] * fin
Inf += x[i] * inf
} else {
Inf += s[i] * fin
INF += s[i] * inf
}
}
if INF != 0.0 {
return INF * Plus
} else if Inf != 0.0 {
return Inf * Plus
} else {
return Fin
}
}
func inspect(x []float64, m uint) (bool, []float64) {
ok, signs := true, make([]float64, m)
for i := uint(0); i < m; i++ {
switch x[i] {
case Plus:
ok, signs[i] = false, +1.0
case Minus:
ok, signs[i] = false, -1.0
}
}
return ok, signs
}