From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: (qmail 15806 invoked by alias); 5 Aug 2002 19:10:09 -0000 Mailing-List: contact gsl-discuss-help@sources.redhat.com; run by ezmlm Precedence: bulk List-Subscribe: List-Archive: List-Post: List-Help: , Sender: gsl-discuss-owner@sources.redhat.com Received: (qmail 15798 invoked from network); 5 Aug 2002 19:10:07 -0000 Received: from unknown (HELO thunder6.cwihosting.com) (64.39.15.209) by sources.redhat.com with SMTP; 5 Aug 2002 19:10:07 -0000 Received: from three by thunder6.cwihosting.com with local (Exim 3.35 #1) id 17bnF0-0004P7-00; Mon, 05 Aug 2002 14:10:06 -0500 From: "Gavin Crooks" To: gsl-discuss@sources.redhat.com Subject: Dirchlet RNG X-IPAddress: 128.32.236.51 MIME-Version: 1.0 Content-Type: multipart/mixed; boundary="----=NEOMAIL_ATT_0.44137657254851" Message-Id: Date: Mon, 05 Aug 2002 12:10:00 -0000 X-AntiAbuse: This header was added to track abuse, please include it with any abuse report X-AntiAbuse: Primary Hostname - thunder6.cwihosting.com X-AntiAbuse: Original Domain - sources.redhat.com X-AntiAbuse: Originator/Caller UID/GID - [588 591] / [588 591] X-AntiAbuse: Sender Address Domain - thunder6.cwihosting.com X-SW-Source: 2002-q3/txt/msg00100.txt.bz2 This is a multi-part message in MIME format. ------=NEOMAIL_ATT_0.44137657254851 Content-Type: text/plain; charset=iso-8859-1 Content-length: 2725 Hi, I have written a Dirchlet random number generator in the style of gsl, which I thought might be of some use. Gavin Crooks gec@compbio.berkeley.edu http://threeplusone.com/ Computational Biology Research Group University of California, Berkeley -- /* randist/dirchlet.c * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or (at * your option) any later version. * * This program is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ /*#include Is this needed?*/ #include #include #include #include #include /* The Dirchlet probability distribution of order K-1 is p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K = (1/Z) \prod_i=1,K \theta_i^(alpha_i - 1) \delta(1 -\sum_i=1,K \theta_i) The normalization factor Z can be expressed in terms of gamma functions: Z = \prod_i=1,K \gamma(\alpha_i) / \gamma( \sum_i=1,K \alpha_i) The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters, \theta_1,...,\theta_K are nonnegative and sum to 1. This code generates a Dirchlet by sampling K values from gamma distributions with parameters \alpha_i, and renormalizing. See A.M. Law, W.D. Kelton, 1991, Simulation Modeling and Analysis Gavin E. Crooks (2002) */ void gsl_ran_dirchlet (const gsl_rng * r, const size_t K, const double alpha[], double theta[]) { size_t i; double norm=0.0; for(i = 0; i < K; i++) theta[i] = gsl_ran_gamma(r, alpha[i], 1.0); for(i = 0; i < K; i++) norm += theta[i]; for(i = 0; i < K; i++) theta[i] /= norm; } double gsl_ran_dirchlet_pdf (const size_t K, const double alpha[], const double theta[]) { /*We calculate the log of the pdf to minimize the possibility of overflow */ size_t i; double log_p = 0.0; double sum_alpha = 0.0; double log_Z =0.0; for(i = 0; i < K; i++) log_p += (alpha[i]-1.0) * log(theta[i]); for(i = 0; i < K; i++) sum_alpha += alpha[i]; log_p += gsl_sf_lngamma ( sum_alpha ); for(i = 0; i < K; i++) log_p -= gsl_sf_lngamma ( alpha[i] ); return exp(log_p); } ------=NEOMAIL_ATT_0.44137657254851 Content-Type: application/octet-stream; name="dirchlet.c" Content-Transfer-Encoding: base64 Content-length: 3339 LyogcmFuZGlzdC9kaXJjaGxldC5jCiAqIAogKiBUaGlzIHByb2dyYW0gaXMg ZnJlZSBzb2Z0d2FyZTsgeW91IGNhbiByZWRpc3RyaWJ1dGUgaXQgYW5kL29y IG1vZGlmeQogKiBpdCB1bmRlciB0aGUgdGVybXMgb2YgdGhlIEdOVSBHZW5l cmFsIFB1YmxpYyBMaWNlbnNlIGFzIHB1Ymxpc2hlZCBieQogKiB0aGUgRnJl ZSBTb2Z0d2FyZSBGb3VuZGF0aW9uOyBlaXRoZXIgdmVyc2lvbiAyIG9mIHRo ZSBMaWNlbnNlLCBvciAoYXQKICogeW91ciBvcHRpb24pIGFueSBsYXRlciB2 ZXJzaW9uLgogKiAKICogVGhpcyBwcm9ncmFtIGlzIGRpc3RyaWJ1dGVkIGlu IHRoZSBob3BlIHRoYXQgaXQgd2lsbCBiZSB1c2VmdWwsIGJ1dAogKiBXSVRI T1VUIEFOWSBXQVJSQU5UWTsgd2l0aG91dCBldmVuIHRoZSBpbXBsaWVkIHdh cnJhbnR5IG9mCiAqIE1FUkNIQU5UQUJJTElUWSBvciBGSVRORVNTIEZPUiBB IFBBUlRJQ1VMQVIgUFVSUE9TRS4gIFNlZSB0aGUgR05VCiAqIEdlbmVyYWwg UHVibGljIExpY2Vuc2UgZm9yIG1vcmUgZGV0YWlscy4KICogCiAqIFlvdSBz aG91bGQgaGF2ZSByZWNlaXZlZCBhIGNvcHkgb2YgdGhlIEdOVSBHZW5lcmFs IFB1YmxpYyBMaWNlbnNlCiAqIGFsb25nIHdpdGggdGhpcyBwcm9ncmFtOyBp ZiBub3QsIHdyaXRlIHRvIHRoZSBGcmVlIFNvZnR3YXJlCiAqIEZvdW5kYXRp b24sIEluYy4sIDY3NSBNYXNzIEF2ZSwgQ2FtYnJpZGdlLCBNQSAwMjEzOSwg VVNBLgogKi8KCi8qI2luY2x1ZGUgPGNvbmZpZy5oPiBJcyB0aGlzIG5lZWRl ZD8qLwojaW5jbHVkZSA8bWF0aC5oPgojaW5jbHVkZSA8Z3NsL2dzbF9tYXRo Lmg+CiNpbmNsdWRlIDxnc2wvZ3NsX3JuZy5oPgojaW5jbHVkZSA8Z3NsL2dz bF9yYW5kaXN0Lmg+CiNpbmNsdWRlIDxnc2wvZ3NsX3NmLmg+CgovKiBUaGUg RGlyY2hsZXQgcHJvYmFiaWxpdHkgZGlzdHJpYnV0aW9uIG9mIG9yZGVyIEst MSBpcyAKCiAgICAgcChcdGhldGFfMSwuLi4sXHRoZXRhX0spIGRcdGhldGFf MSAuLi4gZFx0aGV0YV9LID0gCiAgICAgICAgKDEvWikgXHByb2RfaT0xLEsg XHRoZXRhX2leKGFscGhhX2kgLSAxKSBcZGVsdGEoMSAtXHN1bV9pPTEsSyBc dGhldGFfaSkKCiAgIFRoZSBub3JtYWxpemF0aW9uIGZhY3RvciBaIGNhbiBi ZSBleHByZXNzZWQgaW4gdGVybXMgb2YgZ2FtbWEgZnVuY3Rpb25zOgoKICAg ICAgWiA9IFxwcm9kX2k9MSxLIFxnYW1tYShcYWxwaGFfaSkgLyBcZ2FtbWEo IFxzdW1faT0xLEsgXGFscGhhX2kpICAKCiAgIFRoZSBLIGNvbnN0YW50cywg XGFscGhhXzEsLi4uLFxhbHBoYV9LLCBtdXN0IGJlIHBvc2l0aXZlLiBUaGUg SyBwYXJhbWV0ZXJzLCAKICAgXHRoZXRhXzEsLi4uLFx0aGV0YV9LIGFyZSBu b25uZWdhdGl2ZSBhbmQgc3VtIHRvIDEuCgogICBUaGlzIGNvZGUgZ2VuZXJh dGVzIGEgRGlyY2hsZXQgYnkgc2FtcGxpbmcgSyB2YWx1ZXMgZnJvbSBnYW1t YQogICBkaXN0cmlidXRpb25zIHdpdGggcGFyYW1ldGVycyBcYWxwaGFfaSwg YW5kIHJlbm9ybWFsaXppbmcuIAogICBTZWUgQS5NLiBMYXcsIFcuRC4gS2Vs dG9uLCAxOTkxLCBTaW11bGF0aW9uIE1vZGVsaW5nIGFuZCBBbmFseXNpcwoK ICAgR2F2aW4gRS4gQ3Jvb2tzIDxnZWNAY29tcGJpby5iZXJrZWxleS5lZHU+ ICgyMDAyKQoqLwoKdm9pZApnc2xfcmFuX2RpcmNobGV0IChjb25zdCBnc2xf cm5nICogciwgY29uc3Qgc2l6ZV90IEssIGNvbnN0IGRvdWJsZSBhbHBoYVtd LCBkb3VibGUgdGhldGFbXSkgCnsKICBzaXplX3QgaTsKICBkb3VibGUgbm9y bT0wLjA7CgogIGZvcihpID0gMDsgaSA8IEs7IGkrKykgCiAgICB0aGV0YVtp XSA9IGdzbF9yYW5fZ2FtbWEociwgYWxwaGFbaV0sIDEuMCk7CgogIGZvcihp ID0gMDsgaSA8IEs7IGkrKykgCiAgICBub3JtICs9IHRoZXRhW2ldOwogIAog IGZvcihpID0gMDsgaSA8IEs7IGkrKykgCiAgICB0aGV0YVtpXSAvPSBub3Jt Owp9CgoKZG91YmxlCmdzbF9yYW5fZGlyY2hsZXRfcGRmIChjb25zdCBzaXpl X3QgSywgY29uc3QgZG91YmxlIGFscGhhW10sIGNvbnN0IGRvdWJsZSB0aGV0 YVtdKSAKewogIC8qV2UgY2FsY3VsYXRlIHRoZSBsb2cgb2YgdGhlIHBkZiB0 byBtaW5pbWl6ZSB0aGUgcG9zc2liaWxpdHkgb2Ygb3ZlcmZsb3cgKi8gCiAg c2l6ZV90IGk7CiAgZG91YmxlIGxvZ19wID0gMC4wOwogIGRvdWJsZSBzdW1f YWxwaGEgPSAwLjA7CiAgZG91YmxlIGxvZ19aICA9MC4wOwoKICBmb3IoaSA9 IDA7IGkgPCBLOyBpKyspIAogICAgbG9nX3AgKz0gKGFscGhhW2ldLTEuMCkg KiBsb2codGhldGFbaV0pOyAKCiAgZm9yKGkgPSAwOyBpIDwgSzsgaSsrKSAK ICAgIHN1bV9hbHBoYSArPSBhbHBoYVtpXTsKCiAgbG9nX3AgKz0gZ3NsX3Nm X2xuZ2FtbWEgKCBzdW1fYWxwaGEgKTsgCgoKICBmb3IoaSA9IDA7IGkgPCBL OyBpKyspIAogICAgbG9nX3AgLT0gZ3NsX3NmX2xuZ2FtbWEgKCBhbHBoYVtp XSApOyAKCgogIHJldHVybiBleHAobG9nX3ApOwp9CgoK ------=NEOMAIL_ATT_0.44137657254851 Content-Type: application/octet-stream; name="test_dirchlet.c" Content-Transfer-Encoding: base64 Content-length: 6523 LyogbW9kaWZpZWQgcmFuZGlzdC90ZXN0LmMKICogCiAqIENvcHlyaWdodCAo QykgMTk5NiwgMTk5NywgMTk5OCwgMTk5OSwgMjAwMCBKYW1lcyBUaGVpbGVy LCBCcmlhbiBHb3VnaAogKiAKICogVGhpcyBwcm9ncmFtIGlzIGZyZWUgc29m dHdhcmU7IHlvdSBjYW4gcmVkaXN0cmlidXRlIGl0IGFuZC9vciBtb2RpZnkK ICogaXQgdW5kZXIgdGhlIHRlcm1zIG9mIHRoZSBHTlUgR2VuZXJhbCBQdWJs aWMgTGljZW5zZSBhcyBwdWJsaXNoZWQgYnkKICogdGhlIEZyZWUgU29mdHdh cmUgRm91bmRhdGlvbjsgZWl0aGVyIHZlcnNpb24gMiBvZiB0aGUgTGljZW5z ZSwgb3IgKGF0CiAqIHlvdXIgb3B0aW9uKSBhbnkgbGF0ZXIgdmVyc2lvbi4K ICogCiAqIFRoaXMgcHJvZ3JhbSBpcyBkaXN0cmlidXRlZCBpbiB0aGUgaG9w ZSB0aGF0IGl0IHdpbGwgYmUgdXNlZnVsLCBidXQKICogV0lUSE9VVCBBTlkg V0FSUkFOVFk7IHdpdGhvdXQgZXZlbiB0aGUgaW1wbGllZCB3YXJyYW50eSBv ZgogKiBNRVJDSEFOVEFCSUxJVFkgb3IgRklUTkVTUyBGT1IgQSBQQVJUSUNV TEFSIFBVUlBPU0UuICBTZWUgdGhlIEdOVQogKiBHZW5lcmFsIFB1YmxpYyBM aWNlbnNlIGZvciBtb3JlIGRldGFpbHMuCiAqIAogKiBZb3Ugc2hvdWxkIGhh dmUgcmVjZWl2ZWQgYSBjb3B5IG9mIHRoZSBHTlUgR2VuZXJhbCBQdWJsaWMg TGljZW5zZQogKiBhbG9uZyB3aXRoIHRoaXMgcHJvZ3JhbTsgaWYgbm90LCB3 cml0ZSB0byB0aGUgRnJlZSBTb2Z0d2FyZQogKiBGb3VuZGF0aW9uLCBJbmMu LCA2NzUgTWFzcyBBdmUsIENhbWJyaWRnZSwgTUEgMDIxMzksIFVTQS4KICov CgovKiAjaW5jbHVkZSA8Y29uZmlnLmg+IElzIHRoaXMgbmVlZGVkPyAqLwoj aW5jbHVkZSA8c3RkaW8uaD4KI2luY2x1ZGUgPHN0ZGxpYi5oPgojaW5jbHVk ZSA8bWF0aC5oPgojaW5jbHVkZSA8Z3NsL2dzbF9tYXRoLmg+CiNpbmNsdWRl IDxnc2wvZ3NsX3JhbmRpc3QuaD4KI2luY2x1ZGUgPGdzbC9nc2xfcm5nLmg+ CiNpbmNsdWRlIDxnc2wvZ3NsX3Rlc3QuaD4KI2luY2x1ZGUgPGdzbC9nc2xf aWVlZV91dGlscy5oPgojaW5jbHVkZSA8Z3NsL2dzbF9zZi5oPgoKI2RlZmlu ZSBOIDEwMDAwMAoKCnZvaWQgZ3NsX3Jhbl9kaXJjaGxldCAoY29uc3QgZ3Ns X3JuZyAqIHIsIGNvbnN0IGludCBLLCBjb25zdCBkb3VibGUgYWxwaGFbXSwg ZG91YmxlIHRoZXRhW10pOwpkb3VibGUgZ3NsX3Jhbl9kaXJjaGxldF9wZGYg KGNvbnN0IGludCBLLCBkb3VibGUgYWxwaGFbXSwgZG91YmxlIHRoZXRhW10p OwoKdm9pZCB0ZXN0X2RpcmNobGV0X21vbWVudHModm9pZCk7Cgpkb3VibGUg dGVzdF9kaXJjaGxldCh2b2lkKTsKZG91YmxlIHRlc3RfZGlyY2hsZXRfcGRm KGRvdWJsZSB4KTsKdm9pZCB0ZXN0UERGIChkb3VibGUgKCpmKSAodm9pZCks IGRvdWJsZSAoKnBkZikoZG91YmxlKSwgY29uc3QgY2hhciAqbmFtZSk7CgoK Z3NsX3JuZyAqcl9nbG9iYWw7CgppbnQKbWFpbiAodm9pZCkKewogIGdzbF9p ZWVlX2Vudl9zZXR1cCAoKTsKCiAgZ3NsX3JuZ19lbnZfc2V0dXAoKSA7CiAg cl9nbG9iYWwgPSBnc2xfcm5nX2FsbG9jIChnc2xfcm5nX2RlZmF1bHQpOwoK I2RlZmluZSBGVU5DKHgpICB0ZXN0XyAjIyB4LCAgICAgICAgICAgICAgICAg ICAgICJ0ZXN0IGdzbF9yYW5fIiAjeAojZGVmaW5lIEZVTkMyKHgpIHRlc3Rf ICMjIHgsIHRlc3RfICMjIHggIyMgX3BkZiwgInRlc3QgZ3NsX3Jhbl8iICN4 CgogIHRlc3RQREYgKEZVTkMyKGRpcmNobGV0KSk7CiAgdGVzdF9kaXJjaGxl dF9tb21lbnRzKCk7CgogIGV4aXQgKGdzbF90ZXN0X3N1bW1hcnkoKSk7Cn0K CgoKLyogQ2hlY2sgdGhhdCB0aGUgb2JzZXJ2ZWQgbWVhbnMgb2YgdGhlIGRp cmNobGV0IHZhcmlhYmxlcyBhcmUKICAgd2l0aGluIHJlYXNvbmFibGUgc3Rh dGlzdGljYWwgZXJyb3JzIG9mIHRoZWlyIGNvcnJlY3QgdmFsdWVzLgoqLwoK I2RlZmluZSBESVJDSExFVF9LIDEwCnZvaWQKdGVzdF9kaXJjaGxldF9tb21l bnRzKHZvaWQpIAp7CiAgZG91YmxlIGFscGhhW0RJUkNITEVUX0tdOwogIGRv dWJsZSB0aGV0YVtESVJDSExFVF9LXTsKICBkb3VibGUgdGhldGFfc3VtW0RJ UkNITEVUX0tdOwogIAogIGRvdWJsZSBhbHBoYV9zdW0gPTAuMDsKICBkb3Vi bGUgbWVhbiwgb2JzX21lYW4sIHNkLCBzaWdtYTsKICBpbnQgc3RhdHVzOyAg CiAgaW50IGssbjsKCgogIGZvcihrPTA7IGs8RElSQ0hMRVRfSzsgaysrKQog ICAgewogICAgICBhbHBoYVtrXSA9IGdzbF9yYW5fZXhwb25lbnRpYWwocl9n bG9iYWwsIDAuMSk7CiAgICAgIGFscGhhX3N1bSArPSBhbHBoYVtrXTsKICAg ICAgdGhldGFfc3VtW2tdID0gMC4wOwogICAgfQoKICBmb3Iobj0wOyBuPE47 IG4rKykgCiAgICB7CiAgICAgIGdzbF9yYW5fZGlyY2hsZXQocl9nbG9iYWws IERJUkNITEVUX0ssIGFscGhhLCB0aGV0YSk7CiAgICAgIGZvcihrPTA7azxE SVJDSExFVF9LOyBrKyspIAoJdGhldGFfc3VtW2tdICs9IHRoZXRhW2tdOwog ICAgfQogIAogIGZvcihrPTA7IGs8RElSQ0hMRVRfSzsgaysrKQogICAgewog ICAgICBtZWFuID0gYWxwaGFba10vIGFscGhhX3N1bTsKICAgICAgc2QgICA9 IHNxcnQoIChhbHBoYVtrXSAqICgxLiAtIGFscGhhW2tdL2FscGhhX3N1bSkp LyhhbHBoYV9zdW0qKGFscGhhX3N1bSsxLikpKTsKICAgICAgb2JzX21lYW4g PSB0aGV0YV9zdW1ba10gLyBOOwogICAgICBzaWdtYSA9IHNxcnQoTikqZmFi cyhtZWFuLW9ic19tZWFuKS9zZDsKCiAgICAgIHN0YXR1cyA9IChzaWdtYT4z Lik7CgogICAgICBnc2xfdGVzdCAoc3RhdHVzLCAidGVzdCBnc2xfcmFuX2Rp cmNobGV0OiBtZWFuICVnIGV4cGVjdGVkLCAlZyBvYnNlcnZlZCIsCgkgICAg bWVhbiwgb2JzX21lYW4pOwogICAgfQp9CgoKCmRvdWJsZQp0ZXN0X2RpcmNo bGV0ICh2b2lkKQp7CiAgLyogVGhpcyBpcyBhIGJpdCBvZiBhIGxhbWUgdGVz dCwgc2luY2Ugd2hlbiBLPTIsIHRoZSBkaXJjaGxldCBkaXN0cmlidXRpb24K ICAgICBiZWNvbWVzIGEgYmV0YSBkaXN0cmlidXRpb24gKi8KICAKICBpbnQg aTsKICBpbnQgSyA9MjsKICBkb3VibGUgYWxwaGFbXSA9IHsyLjUsIDUuMH07 CiAgZG91YmxlIHRoZXRhWzJdOwoKICBnc2xfcmFuX2RpcmNobGV0ICggcl9n bG9iYWwsIEssIGFscGhhLCB0aGV0YSk7CgogIHJldHVybiB0aGV0YVswXTsK fQoKZG91YmxlCnRlc3RfZGlyY2hsZXRfcGRmIChkb3VibGUgeCkKewogIGlu dCBpOwogIGludCBLID0yOwogIGRvdWJsZSBhbHBoYVtdID0gezIuNSwgNS4w fTsKICBkb3VibGUgdGhldGFbMl07CgogIGlmKCB4PD0wLjAgfHwgeD49IDEu MCkgcmV0dXJuIDAuMDsgLyogT3V0IG9mIHJhbmdlICovCgogIHRoZXRhWzBd PXg7CiAgdGhldGFbMV0gPSAxLjAgLXg7CgogIHJldHVybiAgZ3NsX3Jhbl9k aXJjaGxldF9wZGYgKEssIGFscGhhLCB0aGV0YSk7Cn0KCgoKCgojZGVmaW5l IEJJTlMgMTAwCgp2b2lkCnRlc3RQREYgKGRvdWJsZSAoKmYpICh2b2lkKSwg ZG91YmxlICgqcGRmKShkb3VibGUpLCBjb25zdCBjaGFyICpuYW1lKQp7CiAg ZG91YmxlIGNvdW50W0JJTlNdLCBwW0JJTlNdOwogIGRvdWJsZSBhID0gLTUu MCwgYiA9ICs1LjAgOwogIGRvdWJsZSBkeCA9IChiIC0gYSkgLyBCSU5TIDsK ICBpbnQgaSxqLHN0YXR1cyA9IDAsIHN0YXR1c19pID0wIDsKCiAgZm9yIChp ID0gMDsgaSA8IEJJTlM7IGkrKykKICAgIGNvdW50W2ldID0gMCA7CgogIGZv ciAoaSA9IDA7IGkgPCBOOyBpKyspCiAgICB7CiAgICAgIGRvdWJsZSByID0g ZiAoKTsKICAgICAgaWYgKHIgPCBiICYmIHIgPiBhKQoJeyAKCSAgaiA9ICAo aW50KSgociAtIGEpL2R4KSA7CgkgIGNvdW50W2pdKys7Cgl9CiAgICB9CiAg CiAgZm9yIChpID0gMDsgaSA8IEJJTlM7IGkrKykKICAgIHsKICAgICAgLyog Q29tcHV0ZSBhbiBhcHByb3hpbWF0aW9uIHRvIHRoZSBpbnRlZ3JhbCBvZiBw KHgpIGZyb20geCB0bwogICAgICAgICB4K2R4IHVzaW5nIFNpbXBzb24ncyBy dWxlICovCgogICAgICBkb3VibGUgeCA9IGEgKyBpICogZHggOwojZGVmaW5l IFNURVBTIDEwMAogICAgICBkb3VibGUgc3VtID0gMCA7CiAgICAgIAogICAg ICBpZiAoZmFicyh4KSA8IDFlLTEwKSAvKiBoaXQgdGhlIG9yaWdpbiBleGFj dGx5ICovCgl4ID0gMC4wIDsgCiAgICAgIAogICAgICBmb3IgKGogPSAxOyBq IDwgU1RFUFM7IGorKykKCXN1bSArPSBwZGYoeCArIGogKiBkeCAvIFNURVBT KSA7CgogICAgICBwW2ldID0gIDAuNSAqIChwZGYoeCkgKyAyKnN1bSArIHBk Zih4ICsgZHggLSAxZS03KSkgKiBkeCAvIFNURVBTIDsKICAgIH0KCiAgZm9y IChpID0gMDsgaSA8IEJJTlM7IGkrKykKICAgIHsKICAgICAgZG91YmxlIHgg PSBhICsgaSAqIGR4IDsKICAgICAgZG91YmxlIGQgPSBmYWJzKGNvdW50W2ld IC0gTipwW2ldKSA7CiAgICAgIGlmIChwW2ldICE9IDApCgl7CgkgIGRvdWJs ZSBzID0gZCAvIHNxcnQoTipwW2ldKSA7CgkgIHN0YXR1c19pID0gKHMgPiA1 KSAmJiAoZCA+IDEpIDsKCX0KICAgICAgZWxzZQoJewoJICBzdGF0dXNfaSA9 IChjb3VudFtpXSAhPSAwKSA7Cgl9CiAgICAgIHN0YXR1cyB8PSBzdGF0dXNf aSA7CiAgICAgIGlmIChzdGF0dXNfaSkgCglnc2xfdGVzdCAoc3RhdHVzX2ks ICIlcyBbJWcsJWcpICglZy8lZD0lZyBvYnNlcnZlZCB2cyAlZyBleHBlY3Rl ZCkiLCAKCQkgIG5hbWUsIHgsIHgrZHgsIGNvdW50W2ldLE4sY291bnRbaV0v TiwgcFtpXSkgOwogICAgfQoKICBpZiAoc3RhdHVzID09IDApCiAgICBnc2xf dGVzdCAoc3RhdHVzLCAiJXMsIHNhbXBsaW5nIGFnYWluc3QgcGRmIG92ZXIg cmFuZ2UgWyVnLCVnKSAiLCAKCSAgICAgIG5hbWUsIGEsIGIpIDsKfQo= ------=NEOMAIL_ATT_0.44137657254851--