//
// Math.seedrandom('yipee'); Sets Math.random to a function that is
// initialized using the given explicit seed.
//
// Math.seedrandom(); Sets Math.random to a function that is
// seeded using the current time, dom state,
// and other accumulated local entropy.
// The generated seed string is returned.
//
// Math.seedrandom('yowza', true);
// Seeds using the given explicit seed mixed
// together with accumulated entropy.
//
//
// Seeds using physical random bits downloaded
// from random.org.
//
// Seeds using urandom bits from call.jsonlib.com,
// which is faster than random.org.
//
// Examples:
//
// Math.seedrandom("hello"); // Use "hello" as the seed.
// document.write(Math.random()); // Always 0.5463663768140734
// document.write(Math.random()); // Always 0.43973793770592234
// var rng1 = Math.random; // Remember the current prng.
//
// var autoseed = Math.seedrandom(); // New prng with an automatic seed.
// document.write(Math.random()); // Pretty much unpredictable.
//
// Math.random = rng1; // Continue "hello" prng sequence.
// document.write(Math.random()); // Always 0.554769432473455
//
// Math.seedrandom(autoseed); // Restart at the previous seed.
// document.write(Math.random()); // Repeat the 'unpredictable' value.
//
// Notes:
//
// Each time seedrandom('arg') is called, entropy from the passed seed
// is accumulated in a pool to help generate future seeds for the
// zero-argument form of Math.seedrandom, so entropy can be injected over
// time by calling seedrandom with explicit data repeatedly.
//
// On speed - This javascript implementation of Math.random() is about
// 3-10x slower than the built-in Math.random() because it is not native
// code, but this is typically fast enough anyway. Seeding is more expensive,
// especially if you use auto-seeding. Some details (timings on Chrome 4):
//
// Our Math.random() - avg less than 0.002 milliseconds per call
// seedrandom('explicit') - avg less than 0.5 milliseconds per call
// seedrandom('explicit', true) - avg less than 2 milliseconds per call
// seedrandom() - avg about 38 milliseconds per call
//
// LICENSE (BSD):
//
// Copyright 2010 David Bau, all rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of this module nor the names of its contributors may
// be used to endorse or promote products derived from this software
// without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
/**
* All code is in an anonymous closure to keep the global namespace clean.
*
* @param {number=} overflow
* @param {number=} startdenom
*/
// Patched by Seb so that seedrandom.js does not pollute the Math object.
// My tests suggest that doing Math.trouble = 1 makes Math lookups about 5%
// slower.
numeric.seedrandom = { pow:Math.pow, random:Math.random };
(function (pool, math, width, chunks, significance, overflow, startdenom) {
//
// seedrandom()
// This is the seedrandom function described above.
//
math['seedrandom'] = function seedrandom(seed, use_entropy) {
var key = [];
var arc4;
// Flatten the seed string or build one from local entropy if needed.
seed = mixkey(flatten(
use_entropy ? [seed, pool] :
arguments.length ? seed :
[new Date().getTime(), pool, window], 3), key);
// Use the seed to initialize an ARC4 generator.
arc4 = new ARC4(key);
// Mix the randomness into accumulated entropy.
mixkey(arc4.S, pool);
// Override Math.random
// This function returns a random double in [0, 1) that contains
// randomness in every bit of the mantissa of the IEEE 754 value.
math['random'] = function random() { // Closure to return a random double:
var n = arc4.g(chunks); // Start with a numerator n < 2 ^ 48
var d = startdenom; // and denominator d = 2 ^ 48.
var x = 0; // and no 'extra last byte'.
while (n < significance) { // Fill up all significant digits by
n = (n + x) * width; // shifting numerator and
d *= width; // denominator and generating a
x = arc4.g(1); // new least-significant-byte.
}
while (n >= overflow) { // To avoid rounding up, before adding
n /= 2; // last byte, shift everything
d /= 2; // right using integer math until
x >>>= 1; // we have exactly the desired bits.
}
return (n + x) / d; // Form the number within [0, 1).
};
// Return the seed that was used
return seed;
};
//
// ARC4
//
// An ARC4 implementation. The constructor takes a key in the form of
// an array of at most (width) integers that should be 0 <= x < (width).
//
// The g(count) method returns a pseudorandom integer that concatenates
// the next (count) outputs from ARC4. Its return value is a number x
// that is in the range 0 <= x < (width ^ count).
//
/** @constructor */
function ARC4(key) {
var t, u, me = this, keylen = key.length;
var i = 0, j = me.i = me.j = me.m = 0;
me.S = [];
me.c = [];
// The empty key [] is treated as [0].
if (!keylen) { key = [keylen++]; }
// Set up S using the standard key scheduling algorithm.
while (i < width) { me.S[i] = i++; }
for (i = 0; i < width; i++) {
t = me.S[i];
j = lowbits(j + t + key[i % keylen]);
u = me.S[j];
me.S[i] = u;
me.S[j] = t;
}
// The "g" method returns the next (count) outputs as one number.
me.g = function getnext(count) {
var s = me.S;
var i = lowbits(me.i + 1); var t = s[i];
var j = lowbits(me.j + t); var u = s[j];
s[i] = u;
s[j] = t;
var r = s[lowbits(t + u)];
while (--count) {
i = lowbits(i + 1); t = s[i];
j = lowbits(j + t); u = s[j];
s[i] = u;
s[j] = t;
r = r * width + s[lowbits(t + u)];
}
me.i = i;
me.j = j;
return r;
};
// For robust unpredictability discard an initial batch of values.
// See http://www.rsa.com/rsalabs/node.asp?id=2009
me.g(width);
}
//
// flatten()
// Converts an object tree to nested arrays of strings.
//
/** @param {Object=} result
* @param {string=} prop
* @param {string=} typ */
function flatten(obj, depth, result, prop, typ) {
result = [];
typ = typeof(obj);
if (depth && typ == 'object') {
for (prop in obj) {
if (prop.indexOf('S') < 5) { // Avoid FF3 bug (local/sessionStorage)
try { result.push(flatten(obj[prop], depth - 1)); } catch (e) {}
}
}
}
return (result.length ? result : obj + (typ != 'string' ? '\0' : ''));
}
//
// mixkey()
// Mixes a string seed into a key that is an array of integers, and
// returns a shortened string seed that is equivalent to the result key.
//
/** @param {number=} smear
* @param {number=} j */
function mixkey(seed, key, smear, j) {
seed += ''; // Ensure the seed is a string
smear = 0;
for (j = 0; j < seed.length; j++) {
key[lowbits(j)] =
lowbits((smear ^= key[lowbits(j)] * 19) + seed.charCodeAt(j));
}
seed = '';
for (j in key) { seed += String.fromCharCode(key[j]); }
return seed;
}
//
// lowbits()
// A quick "n mod width" for width a power of 2.
//
function lowbits(n) { return n & (width - 1); }
//
// The following constants are related to IEEE 754 limits.
//
startdenom = math.pow(width, chunks);
significance = math.pow(2, significance);
overflow = significance * 2;
//
// When seedrandom.js is loaded, we immediately mix a few bits
// from the built-in RNG into the entropy pool. Because we do
// not want to intefere with determinstic PRNG state later,
// seedrandom will not call math.random on its own again after
// initialization.
//
mixkey(math.random(), pool);
// End anonymous scope, and pass initial values.
}(
[], // pool: entropy pool starts empty
numeric.seedrandom, // math: package containing random, pow, and seedrandom
256, // width: each RC4 output is 0 <= x < 256
6, // chunks: at least six RC4 outputs for each double
52 // significance: there are 52 significant digits in a double
));
/* This file is a slightly modified version of quadprog.js from Alberto Santini.
* It has been slightly modified by Sébastien Loisel to make sure that it handles
* 0-based Arrays instead of 1-based Arrays.
* License is in resources/LICENSE.quadprog */
(function(exports) {
function base0to1(A) {
if(typeof A !== "object") { return A; }
var ret = [], i,n=A.length;
for(i=0;i meq) {
work[l] = sum;
} else {
work[l] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][i] = -amat[j][i];
}
bvec[i] = -bvec[i];
}
}
}
for (i = 1; i <= nact; i = i + 1) {
work[iwsv + iact[i]] = 0;
}
nvl = 0;
temp = 0;
for (i = 1; i <= q; i = i + 1) {
if (work[iwsv + i] < temp * work[iwnbv + i]) {
nvl = i;
temp = work[iwsv + i] / work[iwnbv + i];
}
}
if (nvl === 0) {
return 999;
}
return 0;
}
function fn_goto_55() {
for (i = 1; i <= n; i = i + 1) {
sum = 0;
for (j = 1; j <= n; j = j + 1) {
sum = sum + dmat[j][i] * amat[j][nvl];
}
work[i] = sum;
}
l1 = iwzv;
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = 0;
}
for (j = nact + 1; j <= n; j = j + 1) {
for (i = 1; i <= n; i = i + 1) {
work[l1 + i] = work[l1 + i] + dmat[i][j] * work[j];
}
}
t1inf = true;
for (i = nact; i >= 1; i = i - 1) {
sum = work[i];
l = iwrm + (i * (i + 3)) / 2;
l1 = l - i;
for (j = i + 1; j <= nact; j = j + 1) {
sum = sum - work[l] * work[iwrv + j];
l = l + j;
}
sum = sum / work[l1];
work[iwrv + i] = sum;
if (iact[i] < meq) {
// continue;
break;
}
if (sum < 0) {
// continue;
break;
}
t1inf = false;
it1 = i;
}
if (!t1inf) {
t1 = work[iwuv + it1] / work[iwrv + it1];
for (i = 1; i <= nact; i = i + 1) {
if (iact[i] < meq) {
// continue;
break;
}
if (work[iwrv + i] < 0) {
// continue;
break;
}
temp = work[iwuv + i] / work[iwrv + i];
if (temp < t1) {
t1 = temp;
it1 = i;
}
}
}
sum = 0;
for (i = iwzv + 1; i <= iwzv + n; i = i + 1) {
sum = sum + work[i] * work[i];
}
if (Math.abs(sum) <= vsmall) {
if (t1inf) {
ierr[1] = 1;
// GOTO 999
return 999;
} else {
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - t1 * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + t1;
// GOTO 700
return 700;
}
} else {
sum = 0;
for (i = 1; i <= n; i = i + 1) {
sum = sum + work[iwzv + i] * amat[i][nvl];
}
tt = -work[iwsv + nvl] / sum;
t2min = true;
if (!t1inf) {
if (t1 < tt) {
tt = t1;
t2min = false;
}
}
for (i = 1; i <= n; i = i + 1) {
sol[i] = sol[i] + tt * work[iwzv + i];
if (Math.abs(sol[i]) < vsmall) {
sol[i] = 0;
}
}
crval[1] = crval[1] + tt * sum * (tt / 2 + work[iwuv + nact + 1]);
for (i = 1; i <= nact; i = i + 1) {
work[iwuv + i] = work[iwuv + i] - tt * work[iwrv + i];
}
work[iwuv + nact + 1] = work[iwuv + nact + 1] + tt;
if (t2min) {
nact = nact + 1;
iact[nact] = nvl;
l = iwrm + ((nact - 1) * nact) / 2 + 1;
for (i = 1; i <= nact - 1; i = i + 1) {
work[l] = work[i];
l = l + 1;
}
if (nact === n) {
work[l] = work[n];
} else {
for (i = n; i >= nact + 1; i = i - 1) {
if (work[i] === 0) {
// continue;
break;
}
gc = Math.max(Math.abs(work[i - 1]), Math.abs(work[i]));
gs = Math.min(Math.abs(work[i - 1]), Math.abs(work[i]));
if (work[i - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[i - 1] / temp;
gs = work[i] / temp;
if (gc === 1) {
// continue;
break;
}
if (gc === 0) {
work[i - 1] = gs * temp;
for (j = 1; j <= n; j = j + 1) {
temp = dmat[j][i - 1];
dmat[j][i - 1] = dmat[j][i];
dmat[j][i] = temp;
}
} else {
work[i - 1] = temp;
nu = gs / (1 + gc);
for (j = 1; j <= n; j = j + 1) {
temp = gc * dmat[j][i - 1] + gs * dmat[j][i];
dmat[j][i] = nu * (dmat[j][i - 1] + temp) - dmat[j][i];
dmat[j][i - 1] = temp;
}
}
}
work[l] = work[nact];
}
} else {
sum = -bvec[nvl];
for (j = 1; j <= n; j = j + 1) {
sum = sum + sol[j] * amat[j][nvl];
}
if (nvl > meq) {
work[iwsv + nvl] = sum;
} else {
work[iwsv + nvl] = -Math.abs(sum);
if (sum > 0) {
for (j = 1; j <= n; j = j + 1) {
amat[j][nvl] = -amat[j][nvl];
}
bvec[nvl] = -bvec[nvl];
}
}
// GOTO 700
return 700;
}
}
return 0;
}
function fn_goto_797() {
l = iwrm + (it1 * (it1 + 1)) / 2 + 1;
l1 = l + it1;
if (work[l1] === 0) {
// GOTO 798
return 798;
}
gc = Math.max(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
gs = Math.min(Math.abs(work[l1 - 1]), Math.abs(work[l1]));
if (work[l1 - 1] >= 0) {
temp = Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
} else {
temp = -Math.abs(gc * Math.sqrt(1 + gs * gs / (gc * gc)));
}
gc = work[l1 - 1] / temp;
gs = work[l1] / temp;
if (gc === 1) {
// GOTO 798
return 798;
}
if (gc === 0) {
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = work[l1 - 1];
work[l1 - 1] = work[l1];
work[l1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = dmat[i][it1];
dmat[i][it1] = dmat[i][it1 + 1];
dmat[i][it1 + 1] = temp;
}
} else {
nu = gs / (1 + gc);
for (i = it1 + 1; i <= nact; i = i + 1) {
temp = gc * work[l1 - 1] + gs * work[l1];
work[l1] = nu * (work[l1 - 1] + temp) - work[l1];
work[l1 - 1] = temp;
l1 = l1 + i;
}
for (i = 1; i <= n; i = i + 1) {
temp = gc * dmat[i][it1] + gs * dmat[i][it1 + 1];
dmat[i][it1 + 1] = nu * (dmat[i][it1] + temp) - dmat[i][it1 + 1];
dmat[i][it1] = temp;
}
}
return 0;
}
function fn_goto_798() {
l1 = l - it1;
for (i = 1; i <= it1; i = i + 1) {
work[l1] = work[l];
l = l + 1;
l1 = l1 + 1;
}
work[iwuv + it1] = work[iwuv + it1 + 1];
iact[it1] = iact[it1 + 1];
it1 = it1 + 1;
if (it1 < nact) {
// GOTO 797
return 797;
}
return 0;
}
function fn_goto_799() {
work[iwuv + nact] = work[iwuv + nact + 1];
work[iwuv + nact + 1] = 0;
iact[nact] = 0;
nact = nact - 1;
iter[2] = iter[2] + 1;
return 0;
}
go = 0;
while (true) {
go = fn_goto_50();
if (go === 999) {
return;
}
while (true) {
go = fn_goto_55();
if (go === 0) {
break;
}
if (go === 999) {
return;
}
if (go === 700) {
if (it1 === nact) {
fn_goto_799();
} else {
while (true) {
fn_goto_797();
go = fn_goto_798();
if (go !== 797) {
break;
}
}
fn_goto_799();
}
}
}
}
}
function solveQP(Dmat, dvec, Amat, bvec, meq, factorized) {
Dmat = base0to1(Dmat);
dvec = base0to1(dvec);
Amat = base0to1(Amat);
var i, n, q,
nact, r,
crval = [], iact = [], sol = [], work = [], iter = [],
message;
meq = meq || 0;
factorized = factorized ? base0to1(factorized) : [undefined, 0];
bvec = bvec ? base0to1(bvec) : [];
// In Fortran the array index starts from 1
n = Dmat.length - 1;
q = Amat[1].length - 1;
if (!bvec) {
for (i = 1; i <= q; i = i + 1) {
bvec[i] = 0;
}
}
for (i = 1; i <= q; i = i + 1) {
iact[i] = 0;
}
nact = 0;
r = Math.min(n, q);
for (i = 1; i <= n; i = i + 1) {
sol[i] = 0;
}
crval[1] = 0;
for (i = 1; i <= (2 * n + (r * (r + 5)) / 2 + 2 * q + 1); i = i + 1) {
work[i] = 0;
}
for (i = 1; i <= 2; i = i + 1) {
iter[i] = 0;
}
qpgen2(Dmat, dvec, n, n, sol, crval, Amat,
bvec, n, q, meq, iact, nact, iter, work, factorized);
message = "";
if (factorized[1] === 1) {
message = "constraints are inconsistent, no solution!";
}
if (factorized[1] === 2) {
message = "matrix D in quadratic function is not positive definite!";
}
return {
solution: base1to0(sol),
value: base1to0(crval),
unconstrained_solution: base1to0(dvec),
iterations: base1to0(iter),
iact: base1to0(iact),
message: message
};
}
exports.solveQP = solveQP;
}(numeric));
/*
Shanti Rao sent me this routine by private email. I had to modify it
slightly to work on Arrays instead of using a Matrix object.
It is apparently translated from http://stitchpanorama.sourceforge.net/Python/svd.py
*/
numeric.svd= function svd(A) {
var temp;
//Compute the thin SVD from G. H. Golub and C. Reinsch, Numer. Math. 14, 403-420 (1970)
var prec= numeric.epsilon; //Math.pow(2,-52) // assumes double prec
var tolerance= 1.e-64/prec;
var itmax= 50;
var c=0;
var i=0;
var j=0;
var k=0;
var l=0;
var u= numeric.clone(A);
var m= u.length;
var n= u[0].length;
if (m < n) throw "Need more rows than columns"
var e = new Array(n);
var q = new Array(n);
for (i=0; i b)
return a*Math.sqrt(1.0+(b*b/a/a))
else if (b == 0.0)
return a
return b*Math.sqrt(1.0+(a*a/b/b))
}
//Householder's reduction to bidiagonal form
var f= 0.0;
var g= 0.0;
var h= 0.0;
var x= 0.0;
var y= 0.0;
var z= 0.0;
var s= 0.0;
for (i=0; i < n; i++)
{
e[i]= g;
s= 0.0;
l= i+1;
for (j=i; j < m; j++)
s += (u[j][i]*u[j][i]);
if (s <= tolerance)
g= 0.0;
else
{
f= u[i][i];
g= Math.sqrt(s);
if (f >= 0.0) g= -g;
h= f*g-s
u[i][i]=f-g;
for (j=l; j < n; j++)
{
s= 0.0
for (k=i; k < m; k++)
s += u[k][i]*u[k][j]
f= s/h
for (k=i; k < m; k++)
u[k][j]+=f*u[k][i]
}
}
q[i]= g
s= 0.0
for (j=l; j < n; j++)
s= s + u[i][j]*u[i][j]
if (s <= tolerance)
g= 0.0
else
{
f= u[i][i+1]
g= Math.sqrt(s)
if (f >= 0.0) g= -g
h= f*g - s
u[i][i+1] = f-g;
for (j=l; j < n; j++) e[j]= u[i][j]/h
for (j=l; j < m; j++)
{
s=0.0
for (k=l; k < n; k++)
s += (u[j][k]*u[i][k])
for (k=l; k < n; k++)
u[j][k]+=s*e[k]
}
}
y= Math.abs(q[i])+Math.abs(e[i])
if (y>x)
x=y
}
// accumulation of right hand gtransformations
for (i=n-1; i != -1; i+= -1)
{
if (g != 0.0)
{
h= g*u[i][i+1]
for (j=l; j < n; j++)
v[j][i]=u[i][j]/h
for (j=l; j < n; j++)
{
s=0.0
for (k=l; k < n; k++)
s += u[i][k]*v[k][j]
for (k=l; k < n; k++)
v[k][j]+=(s*v[k][i])
}
}
for (j=l; j < n; j++)
{
v[i][j] = 0;
v[j][i] = 0;
}
v[i][i] = 1;
g= e[i]
l= i
}
// accumulation of left hand transformations
for (i=n-1; i != -1; i+= -1)
{
l= i+1
g= q[i]
for (j=l; j < n; j++)
u[i][j] = 0;
if (g != 0.0)
{
h= u[i][i]*g
for (j=l; j < n; j++)
{
s=0.0
for (k=l; k < m; k++) s += u[k][i]*u[k][j];
f= s/h
for (k=i; k < m; k++) u[k][j]+=f*u[k][i];
}
for (j=i; j < m; j++) u[j][i] = u[j][i]/g;
}
else
for (j=i; j < m; j++) u[j][i] = 0;
u[i][i] += 1;
}
// diagonalization of the bidiagonal form
prec= prec*x
for (k=n-1; k != -1; k+= -1)
{
for (var iteration=0; iteration < itmax; iteration++)
{ // test f splitting
var test_convergence = false
for (l=k; l != -1; l+= -1)
{
if (Math.abs(e[l]) <= prec)
{ test_convergence= true
break
}
if (Math.abs(q[l-1]) <= prec)
break
}
if (!test_convergence)
{ // cancellation of e[l] if l>0
c= 0.0
s= 1.0
var l1= l-1
for (i =l; i= itmax-1)
throw 'Error: no convergence.'
// shift from bottom 2x2 minor
x= q[l]
y= q[k-1]
g= e[k-1]
h= e[k]
f= ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y)
g= pythag(f,1.0)
if (f < 0.0)
f= ((x-z)*(x+z)+h*(y/(f-g)-h))/x
else
f= ((x-z)*(x+z)+h*(y/(f+g)-h))/x
// next QR transformation
c= 1.0
s= 1.0
for (i=l+1; i< k+1; i++)
{
g= e[i]
y= q[i]
h= s*g
g= c*g
z= pythag(f,h)
e[i-1]= z
c= f/z
s= h/z
f= x*c+g*s
g= -x*s+g*c
h= y*s
y= y*c
for (j=0; j < n; j++)
{
x= v[j][i-1]
z= v[j][i]
v[j][i-1] = x*c+z*s
v[j][i] = -x*s+z*c
}
z= pythag(f,h)
q[i-1]= z
c= f/z
s= h/z
f= c*g+s*y
x= -s*g+c*y
for (j=0; j < m; j++)
{
y= u[j][i-1]
z= u[j][i]
u[j][i-1] = y*c+z*s
u[j][i] = -y*s+z*c
}
}
e[l]= 0.0
e[k]= f
q[k]= x
}
}
//vt= transpose(v)
//return (u,q,vt)
for (i=0;i= 0; j--)
{
if (q[j] < q[i])
{
// writeln(i,'-',j)
c = q[j]
q[j] = q[i]
q[i] = c
for(k=0;k