Mercurial > hg > octave-lojdl
changeset 13212:e81b93284605
Various minor stylistic improvements to sprandsym.m
author | Jordi Gutiérrez Hermoso <jordigh@octave.org> |
---|---|
date | Sun, 25 Sep 2011 01:42:37 -0500 |
parents | 78744376463a |
children | 544304a09e42 |
files | scripts/sparse/sprandsym.m |
diffstat | 1 files changed, 17 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/sparse/sprandsym.m +++ b/scripts/sparse/sprandsym.m @@ -65,9 +65,8 @@ ondiag = randperm (n, r); offdiag = randperm (n*(n - 1)/2, m); - ## solve with quadratic formula n^2 - n - 2*offdiag = 0 to get the row index - i = floor((1 + sqrt(1 + 8*offdiag))/2); - i(i.^2 - i - 2*offdiag != 0) += 1; + ## Row index + i = lookup (cumsum (1:n), offdiag) + 2; ## Column index j = offdiag - (i - 1).*(i - 2)/2; @@ -106,27 +105,31 @@ ## Q = ------------------------------- ## (D + 2)*(D + 1)*(A - M + 1) ## - ## Then the cumprod of these quotients is + ## Then, after prepending 1, the cumprod of these quotients is ## - ## C = [ T(2)/T(1), T(3)/T(1), ..., T(N)/T(1) ] + ## C = [ T(1)/T(1), T(2)/T(1), T(3)/T(1), ..., T(N)/T(1) ] ## - ## Their sum + 1 is thus S = sum (T)/T(1), and then C(i)/S is the - ## desired probability P(i) for i = 2:N. The first P(1) can be - ## obtained by the condition that sum (P) = 1, and the cumsum will - ## give the distribution function for computing the random number of - ## entries on the diagonal R. + ## Their sum is thus S = sum (T)/T(1), and then C(i)/S is the desired + ## probability P(i) for i=1:N. The cumsum will finally give the + ## distribution function for computing the random number of entries on + ## the diagonal R. + + ## Degenerate case + if k == 1 + r = 1; + return + endif ## Compute the stuff described above a = n*(n - 1)/2; d = [mod(k,2):2:min(n,k)-2]; m = (k - d)/2; q = (n - d).*(n - d - 1).*m ./ (d + 2)./(d + 1)./(a - m + 1); - c = cumprod (q); - s = sum (c) + 1; + c = [1 cumprod (q)]; + s = sum (c); p = c/s; - ## Add missing entries - p = [1 - sum(p), p]; + ## Add final d d(end+1) = d(end) + 2; ## Pick a random r using this distribution