# File bigdecimal/ludcmp.rb, line 5
 5:   def ludecomp(a,n,zero=0,one=1)
 6:     prec = BigDecimal.limit(nil)
 7:     ps     = []
 8:     scales = []
 9:     for i in 0...n do  # pick up largest(abs. val.) element in each row.
10:       ps <<= i
11:       nrmrow  = zero
12:       ixn = i*n
13:       for j in 0...n do
14:          biggst = a[ixn+j].abs
15:          nrmrow = biggst if biggst>nrmrow
16:       end
17:       if nrmrow>zero then
18:          scales <<= one.div(nrmrow,prec)
19:       else 
20:          raise "Singular matrix"
21:       end
22:     end
23:     n1          = n - 1
24:     for k in 0...n1 do # Gaussian elimination with partial pivoting.
25:       biggst  = zero;
26:       for i in k...n do
27:          size = a[ps[i]*n+k].abs*scales[ps[i]]
28:          if size>biggst then
29:             biggst = size
30:             pividx  = i
31:          end
32:       end
33:       raise "Singular matrix" if biggst<=zero
34:       if pividx!=k then
35:         j = ps[k]
36:         ps[k] = ps[pividx]
37:         ps[pividx] = j
38:       end
39:       pivot   = a[ps[k]*n+k]
40:       for i in (k+1)...n do
41:         psin = ps[i]*n
42:         a[psin+k] = mult = a[psin+k].div(pivot,prec)
43:         if mult!=zero then
44:            pskn = ps[k]*n
45:            for j in (k+1)...n do
46:              a[psin+j] -= mult.mult(a[pskn+j],prec)
47:            end
48:         end
49:       end
50:     end
51:     raise "Singular matrix" if a[ps[n1]*n+n1] == zero
52:     ps
53:   end
54: 
55:   def lusolve(a,b,ps,zero=0.0)
56:     prec = BigDecimal.limit(nil)
57:     n = ps.size
58:     x = []
59:     for i in 0...n do
60:       dot = zero
61:       psin = ps[i]*n
62:       for j in 0...i do
63:         dot = a[psin+j].mult(x[j],prec) + dot
64:       end
65:       x <<= b[ps[i]] - dot
66:     end
67:     (n-1).downto(0) do |i|
68:        dot = zero
69:        psin = ps[i]*n
70:        for j in (i+1)...n do
71:          dot = a[psin+j].mult(x[j],prec) + dot
72:        end
73:        x[i]  = (x[i]-dot).div(a[psin+i],prec)
74:     end
75:     x
76:   end
77: end

    # File bigdecimal/ludcmp.rb, line 55
55:   def lusolve(a,b,ps,zero=0.0)
56:     prec = BigDecimal.limit(nil)
57:     n = ps.size
58:     x = []
59:     for i in 0...n do
60:       dot = zero
61:       psin = ps[i]*n
62:       for j in 0...i do
63:         dot = a[psin+j].mult(x[j],prec) + dot
64:       end
65:       x <<= b[ps[i]] - dot
66:     end
67:     (n-1).downto(0) do |i|
68:        dot = zero
69:        psin = ps[i]*n
70:        for j in (i+1)...n do
71:          dot = a[psin+j].mult(x[j],prec) + dot
72:        end
73:        x[i]  = (x[i]-dot).div(a[psin+i],prec)
74:     end
75:     x
76:   end
77: end