Robust/src/Benchmarks/oooJava/moldyn/JGFMolDynBench.java

   1 /**************************************************************************
   2  * * Java Grande Forum Benchmark Suite - Thread Version 1.0 * * produced by * *
   3  * Java Grande Benchmarking Project * * at * * Edinburgh Parallel Computing
   4  * Centre * * email: epcc-javagrande@epcc.ed.ac.uk * * * This version copyright
   5  * (c) The University of Edinburgh, 2001. * All rights reserved. * *
   6  **************************************************************************/
   7 public class JGFMolDynBench {
   8   public int ITERS;
   9   public double LENGTH;
  10   public double m;
  11   public double mu;
  12   public double kb;
  13   public double TSIM;
  14   public double deltat;
  15
  16   public int PARTSIZE;
  17
  18   public double[] epot;
  19   public double[] vir;
  20   public double[] ek;
  21
  22   int size, mm;
  23   int[] datasizes;
  24
  25   public int interactions;
  26   public int[] interacts;
  27
  28   public int nthreads;
  29   public int workload;
  30
  31   public JGFMolDynBench(int nthreads,int workload) {
  32     this.nthreads = nthreads;
  33     this.workload=workload;
  34   }
  35
  36   public void JGFsetsize(int size) {
  37     this.size = size;
  38   }
  39
  40   public void JGFinitialise() {
  41     interactions = 0;
  42     datasizes = new int[3];
  43     datasizes[0] = 8;
  44     datasizes[1] = 13;
  45     datasizes[2] = 11;
  46
  47     mm = datasizes[size];
  48     PARTSIZE = mm * mm * mm * 4;
  49     ITERS = 100;
  50     LENGTH = 50e-10;
  51     m = 4.0026;
  52     mu = 1.66056e-27;
  53     kb = 1.38066e-23;
  54     TSIM = 50;
  55     deltat = 5e-16;
  56   }
  57
  58   public static void JGFapplication(JGFMolDynBench mold) {
  59     double sh_force[][];
  60     double sh_force2[][][];
  61     int partsize, numthreads;
  62     partsize = mold.PARTSIZE;
  63     numthreads = mold.nthreads;
  64
  65     sh_force = new double[3][partsize];
  66     sh_force2 = new double[3][numthreads][partsize];
  67     mold.epot = new double[numthreads];
  68     mold.vir = new double[numthreads];
  69     mold.ek = new double[numthreads];
  70     mold.interacts = new int[numthreads];
  71     // for(int i=0;i<numthreads;i++) {
  72     // mold.epot[i]=new double();
  73     // mold.vir[i]=new double();
  74     // mold.ek[i]=new double();
  75     // mold.interacts[i]=new IntWrapper();
  76     // }
  77
  78     // spawn threads
  79     MDWrap[] thobjects = new MDWrap[numthreads];
  80
  81     for (int i = 0; i < numthreads; i++) {
  82       thobjects[i] = new MDWrap(new mdRunner(i, mold.mm, sh_force, sh_force2, mold.nthreads, mold,mold.workload));
  83     }
  84
  85     /*
  86      * boolean waitfordone=true; while(waitfordone) { if (mybarr.done)
  87      * waitfordone=false; }
  88      */
  89     long start=System.currentTimeMillis();
  90     for (int i = 0; i < numthreads; i++) {
  91       // thobjects[i].md.start(mid[i]);
  92       thobjects[i].md.run();
  93     }
  94     long end=System.currentTimeMillis();
  95 //    System.out.println("Total="+(end-start));
  96   }
  97
  98   public void JGFvalidate() {
  99     double[] refval = new double[2];
 100     refval[0] = 1731.4306625334357;
 101     refval[1] = 7397.392307839352;
 102     double dev = Math.fabs(ek[0] - refval[size]);
 103     if (dev > 1.0e-10) {
 104       // System.printString("Validation failed\n");
 105       // System.printString("Kinetic Energy = " + (long)ek[0] + "  " + (long)dev
 106       // + "  " + size + "\n");
 107     }
 108   }
 109 }
 110
 111 class mdRunner {
 112
 113   double count;
 114   int id, i, j, k, lg, mm;
 115   double l, rcoff, rcoffs, side, sideh, hsq, hsq2, vel, velt;
 116   double a, r, sum, tscale, sc, ekin, ts, sp;
 117   double den;
 118   double tref;
 119   double h;
 120   double vaver, vaverh, rand;
 121   double etot, temp, pres, rp;
 122   double u1, u2, s, xx, yy, zz;
 123   double xvelocity, yvelocity, zvelocity;
 124
 125   double[][] sh_force;
 126   double[][][] sh_force2;
 127
 128   int ijk, npartm, iseed, tint;
 129   int irep;
 130   int istop;
 131   int iprint;
 132
 133   JGFMolDynBench mymd;
 134   int nthreads;
 135   int workload;
 136
 137   public mdRunner(int id, int mm, double[][] sh_force, double[][][] sh_force2, int nthreads,
 138       JGFMolDynBench mymd, int workload) {
 139     this.id = id;
 140     this.mm = mm;
 141     this.sh_force = sh_force;
 142     this.sh_force2 = sh_force2;
 143     this.nthreads = nthreads;
 144     this.mymd = mymd;
 145     count = 0.0;
 146     den = 0.83134;
 147     tref = 0.722;
 148     h = 0.064;
 149     irep = 10;
 150     istop = 19;
 151     iprint = 10;
 152     this.workload=workload;
 153   }
 154
 155   public void init(particle[] one, int mdsize) {
 156     int id = this.id;
 157     for (int lg = 0; lg <= 1; lg++) {
 158       for (int i = 0; i < mm; i++) {
 159         for (int j = 0; j < mm; j++) {
 160           for (int k = 0; k < mm; k++) {
 161             one[ijk] =
 162                 new particle((i * a + lg * a * 0.5), (j * a + lg * a * 0.5), (k * a), xvelocity,
 163                     yvelocity, zvelocity, sh_force, sh_force2, id, one);
 164             ijk = ijk + 1;
 165           }
 166         }
 167       }
 168     }
 169
 170     for (int lg = 1; lg <= 2; lg++) {
 171       for (int i = 0; i < mm; i++) {
 172         for (int j = 0; j < mm; j++) {
 173           for (int k = 0; k < mm; k++) {
 174             one[ijk] =
 175                 new particle((i * a + (2 - lg) * a * 0.5), (j * a + (lg - 1) * a * 0.5),
 176                     (k * a + a * 0.5), xvelocity, yvelocity, zvelocity, sh_force, sh_force2, id,
 177                     one);
 178             ijk = ijk + 1;
 179           }
 180         }
 181       }
 182     }
 183
 184     /* Initialise velocities */
 185
 186     iseed = 0;
 187     double v1 = 0.0;
 188     double v2 = 0.0;
 189     random randnum = new random(iseed, v1, v2);
 190
 191     for (int i = 0; i < mdsize; i += 2) {
 192       r = randnum.seed();
 193       one[i].xvelocity = r * randnum.v1;
 194       one[i + 1].xvelocity = r * randnum.v2;
 195     }
 196
 197     for (int i = 0; i < mdsize; i += 2) {
 198       r = randnum.seed();
 199       one[i].yvelocity = r * randnum.v1;
 200       one[i + 1].yvelocity = r * randnum.v2;
 201     }
 202
 203     for (int i = 0; i < mdsize; i += 2) {
 204       r = randnum.seed();
 205       one[i].zvelocity = r * randnum.v1;
 206       one[i + 1].zvelocity = r * randnum.v2;
 207     }
 208
 209     /* velocity scaling */
 210
 211     ekin = 0.0;
 212     sp = 0.0;
 213
 214     for (int i = 0; i < mdsize; i++) {
 215       sp = sp + one[i].xvelocity;
 216     }
 217     sp = sp / mdsize;
 218
 219     for (int i = 0; i < mdsize; i++) {
 220       one[i].xvelocity = one[i].xvelocity - sp;
 221       ekin = ekin + one[i].xvelocity * one[i].xvelocity;
 222     }
 223
 224     sp = 0.0;
 225     for (int i = 0; i < mdsize; i++) {
 226       sp = sp + one[i].yvelocity;
 227     }
 228     sp = sp / mdsize;
 229
 230     for (int i = 0; i < mdsize; i++) {
 231       one[i].yvelocity = one[i].yvelocity - sp;
 232       ekin = ekin + one[i].yvelocity * one[i].yvelocity;
 233     }
 234
 235     sp = 0.0;
 236     for (int i = 0; i < mdsize; i++) {
 237       sp = sp + one[i].zvelocity;
 238     }
 239     sp = sp / mdsize;
 240
 241     for (int i = 0; i < mdsize; i++) {
 242       one[i].zvelocity = one[i].zvelocity - sp;
 243       ekin = ekin + one[i].zvelocity * one[i].zvelocity;
 244     }
 245
 246     ts = tscale * ekin;
 247     sc = h * Math.sqrt(tref / ts);
 248
 249     for (int i = 0; i < mdsize; i++) {
 250
 251       one[i].xvelocity = one[i].xvelocity * sc;
 252       one[i].yvelocity = one[i].yvelocity * sc;
 253       one[i].zvelocity = one[i].zvelocity * sc;
 254
 255     }
 256
 257   }
 258
 259   public void doinit(int mdsize) {
 260     for (int j = 0; j < 3; j++) {
 261       double[] sh = sh_force[j];
 262       for (int i = 0; i < sh.length; i++) {
 263         sh[i] = 0.0;
 264       }
 265     }
 266   }
 267
 268   public void doinit2(int mdsize) {
 269     for (int k = 0; k < 3; k++) {
 270       double[] sh = sh_force[k];
 271       double[][] sha = sh_force2[k];
 272       for (int j = 0; j < nthreads; j++) {
 273         double[] sha2 = sha[j];
 274         for (int i = 0; i < mdsize; i++) {
 275           sh[i] += sha2[i];
 276         }
 277       }
 278     }
 279
 280     for (int k = 0; k < 3; k++) {
 281       double[][] sh1 = sh_force2[k];
 282       for (int j = 0; j < nthreads; j++) {
 283         double[] sh2 = sh1[j];
 284         for (int i = 0; i < mdsize; i++) {
 285
 286           sh2[i] = 0.0;
 287         }
 288       }
 289     }
 290
 291     for (int j = 1; j < nthreads; j++) {
 292       mymd.epot[0] += mymd.epot[j];
 293       mymd.vir[0] += mymd.vir[j];
 294     }
 295     for (int j = 1; j < nthreads; j++) {
 296       mymd.epot[j] = mymd.epot[0];
 297       mymd.vir[j] = mymd.vir[0];
 298     }
 299     for (int j = 0; j < nthreads; j++) {
 300       mymd.interactions += mymd.interacts[j];
 301     }
 302
 303     for (int j = 0; j < 3; j++) {
 304       double sh[] = sh_force[j];
 305       for (int i = 0; i < mdsize; i++) {
 306         sh[i] = sh[i] * hsq2;
 307       }
 308     }
 309   }
 310
 311   public void run() {
 312     /* Parameter determination */
 313
 314     int mdsize;
 315     double tmpden;
 316     int movemx = 50;
 317     particle[] one;
 318     int id;
 319     id = this.id;
 320     mdsize = mymd.PARTSIZE;
 321     one = new particle[mdsize];
 322     l = mymd.LENGTH;
 323     tmpden = den;
 324     side = Math.pow((mdsize / tmpden), 0.3333333);
 325     rcoff = mm / 4.0;
 326
 327     a = side / mm;
 328     sideh = side * 0.5;
 329     hsq = h * h;
 330     hsq2 = hsq * 0.5;
 331     npartm = mdsize - 1;
 332     rcoffs = rcoff * rcoff;
 333     tscale = 16.0 / (1.0 * mdsize - 1.0);
 334     vaver = 1.13 * Math.sqrt(tref / 24.0);
 335     vaverh = vaver * h;
 336
 337     /* Particle Generation */
 338
 339     xvelocity = 0.0;
 340     yvelocity = 0.0;
 341     zvelocity = 0.0;
 342     ijk = 0;
 343     init(one, mdsize);
 344
 345     /* Synchronise threads and start timer before MD simulation */
 346
 347     /* MD simulation */
 348
 349     JGFMolDynBench l_mymd=mymd;
 350
 351     int numP= (mdsize / workload)+1;
 352
 353     double scratchpad[][][];
 354     scratchpad=new double[numP][3][mdsize];
 355
 356     long par_time=0;
 357
 358     for (int move = 0; move < movemx; move++) {
 359       /* move the particles and update velocities */
 360       for (int i = 0; i < one.length; i++) {
 361         one[i].domove(side, i);
 362       }
 363
 364       if (id == 0) {
 365         doinit(mdsize);
 366       }
 367
 368       mymd.epot[id] = 0.0;
 369       mymd.vir[id] = 0.0;
 370       mymd.interacts[id] = 0;
 371
 372       /* compute forces */
 373       int numThread = nthreads;
 374       int lworkload = workload;
 375       // for (int i=0+id;i<mdsize;i+=numThread) {
 376       int scratch_idx=0;
 377
 378       double l_epot=0.0;
 379       double l_vir=0.0;
 380       int l_interacts=0;
 381
 382       for (int i = 0 ; i < mdsize; i += lworkload) {
 383
 384         int ilow = i;
 385         int iupper = i + lworkload;
 386         if (iupper > mdsize) {
 387           iupper = mdsize;
 388         }
 389         int l_size = iupper - ilow;
 390
 391         double workingpad[][]=scratchpad[scratch_idx++];
 392         long par_start=System.currentTimeMillis();
 393         sese parallel{
 394           for(int j=0;j<3;j++){
 395             for(int l=0;l<mdsize;l++){
 396               workingpad[j][l]=0;
 397             }
 398           }
 399           MDStore store = new MDStore();
 400           for(int idx=ilow;idx<iupper;idx++){
 401 //            one[i].force(side, rcoff, mdsize, i, xx, yy, zz, mymd, worker);
 402             one[idx].force(side, rcoff, mdsize, idx, xx, yy, zz, mymd, store,workingpad);
 403           }
 404         }
 405         long par_end=System.currentTimeMillis();
 406         par_time+=(par_end-par_start);
 407
 408         sese serial{
 409           for (int k = 0; k < 3; k++) {
 410             for (int j = 0; j < mdsize; j++) {
 411 //              sh_force[k][j] += worker.sh_force2[k][j];
 412               sh_force[k][j] += workingpad[k][j];
 413             }
 414           }
 415           l_epot+=store.epot;
 416           l_vir+=store.vir;
 417           l_interacts+=store.interacts;
 418 //          mymd.epot[0] += store.epot;
 419 //          mymd.vir[0] += store.vir;
 420 //          mymd.interactions += store.interacts;
 421         }
 422
 423       }
 424
 425       mymd.epot[0] =l_epot;
 426 //      System.out.println("SEQ_START");
 427       mymd.vir[0] = l_vir;
 428       mymd.interactions = l_interacts;
 429
 430       for (int k = 0; k < 3; k++) {
 431         for (int j = 0; j <  sh_force[k].length; j++) {
 432           sh_force[k][j] = sh_force[k][j] * hsq2;
 433         }
 434       }
 435
 436       /* update force arrays */
 437       // if(id == 0) {
 438       // doinit2(mdsize);
 439       // }
 440
 441       /* scale forces, update velocities */
 442       double l_sum=0.0;
 443       int maxIdx=one.length;
 444       for (int i = 0; i < maxIdx; i++) {
 445         l_sum = l_sum + one[i].mkekin(hsq2, i);
 446       }
 447       sum=l_sum;
 448
 449       ekin = sum / hsq;
 450
 451       vel = 0.0;
 452       count = 0.0;
 453
 454       /* average velocity */
 455
 456       for (int i = 0; i < mdsize; i++) {
 457         velt = one[i].velavg(vaverh, h);
 458         if (velt > vaverh) {
 459           count = count + 1.0;
 460         }
 461         vel = vel + velt;
 462       }
 463
 464       vel = vel / h;
 465
 466       /* temperature scale if required */
 467
 468       if ((move < istop) && (((move + 1) % irep) == 0)) {
 469         sc = Math.sqrt(tref / (tscale * ekin));
 470         for (int i = 0; i < mdsize; i++) {
 471           one[i].dscal(sc, 1);
 472         }
 473         ekin = tref / tscale;
 474       }
 475
 476       /* sum to get full potential energy and virial */
 477
 478       if (((move + 1) % iprint) == 0) {
 479         mymd.ek[id] = 24.0 * ekin;
 480         mymd.epot[id] = 4.0 * mymd.epot[id];
 481         etot = mymd.ek[id] + mymd.epot[id];
 482         temp = tscale * ekin;
 483         pres = tmpden * 16.0 * (ekin - mymd.vir[id]) / mdsize;
 484         vel = vel / mdsize;
 485         rp = (count / mdsize) * 100.0;
 486       }
 487
 488       // if (id == 0) JGFInstrumentor.stopTimer("Section3:MolDyn:Run",
 489       // instr.timers);
 490     }
 491 //    System.out.println("par time="+par_time);
 492   }
 493 }
 494
 495 class particle {
 496
 497   public double xcoord, ycoord, zcoord;
 498   public double xvelocity, yvelocity, zvelocity;
 499   int part_id;
 500   int id;
 501   double[][] sh_force;
 502   double[][][] sh_force2;
 503   particle[] one;
 504
 505   public particle(double xcoord, double ycoord, double zcoord, double xvelocity, double yvelocity,
 506       double zvelocity, double[][] sh_force, double[][][] sh_force2, int id, particle[] one) {
 507
 508     this.xcoord = xcoord;
 509     this.ycoord = ycoord;
 510     this.zcoord = zcoord;
 511     this.xvelocity = xvelocity;
 512     this.yvelocity = yvelocity;
 513     this.zvelocity = zvelocity;
 514     this.sh_force = sh_force;
 515     this.sh_force2 = sh_force2;
 516     this.id = id;
 517     this.one = one;
 518   }
 519
 520   public void domove(double side, int part_id) {
 521
 522     xcoord = xcoord + xvelocity + sh_force[0][part_id];
 523     ycoord = ycoord + yvelocity + sh_force[1][part_id];
 524     zcoord = zcoord + zvelocity + sh_force[2][part_id];
 525
 526     if (xcoord < 0) {
 527       xcoord = xcoord + side;
 528     }
 529     if (xcoord > side) {
 530       xcoord = xcoord - side;
 531     }
 532     if (ycoord < 0) {
 533       ycoord = ycoord + side;
 534     }
 535     if (ycoord > side) {
 536       ycoord = ycoord - side;
 537     }
 538     if (zcoord < 0) {
 539       zcoord = zcoord + side;
 540     }
 541     if (zcoord > side) {
 542       zcoord = zcoord - side;
 543     }
 544
 545     xvelocity = xvelocity + sh_force[0][part_id];
 546     yvelocity = yvelocity + sh_force[1][part_id];
 547     zvelocity = zvelocity + sh_force[2][part_id];
 548
 549   }
 550
 551   // public void force(double side, double rcoff,int mdsize,int x, double xx,
 552   // double yy, double zz, JGFMolDynBench mymd) {
 553   public void force(double side, double rcoff, int mdsize, int x, double xx, double yy, double zz,
 554       JGFMolDynBench mymd, MDStore store,double workingpad[][]) {
 555
 556     double sideh;
 557     double rcoffs;
 558
 559     double fxi, fyi, fzi;
 560     double rd, rrd, rrd2, rrd3, rrd4, rrd6, rrd7, r148;
 561     double forcex, forcey, forcez;
 562     int id = this.id;
 563     sideh = 0.5 * side;
 564     rcoffs = rcoff * rcoff;
 565
 566     fxi = 0.0;
 567     fyi = 0.0;
 568     fzi = 0.0;
 569
 570     for (int i = x + 1; i < mdsize; i++) {
 571       xx = this.xcoord - one[i].xcoord;
 572       yy = this.ycoord - one[i].ycoord;
 573       zz = this.zcoord - one[i].zcoord;
 574
 575       if (xx < (-sideh)) {
 576         xx = xx + side;
 577       }
 578       if (xx > (sideh)) {
 579         xx = xx - side;
 580       }
 581       if (yy < (-sideh)) {
 582         yy = yy + side;
 583       }
 584       if (yy > (sideh)) {
 585         yy = yy - side;
 586       }
 587       if (zz < (-sideh)) {
 588         zz = zz + side;
 589       }
 590       if (zz > (sideh)) {
 591         zz = zz - side;
 592       }
 593
 594       rd = xx * xx + yy * yy + zz * zz;
 595
 596       if (rd <= rcoffs) {
 597         rrd = 1.0 / rd;
 598         rrd2 = rrd * rrd;
 599         rrd3 = rrd2 * rrd;
 600         rrd4 = rrd2 * rrd2;
 601         rrd6 = rrd2 * rrd4;
 602         rrd7 = rrd6 * rrd;
 603         // mymd.epot[id] += (rrd6 - rrd3);
 604         store.epot += (rrd6 - rrd3);
 605         r148 = rrd7 - 0.5 * rrd4;
 606         // mymd.vir[id] += - rd*r148;
 607         store.vir += -rd * r148;
 608         forcex = xx * r148;
 609         fxi = fxi + forcex;
 610
 611         // sh_force2[0][id][i] = sh_force2[0][id][i] - forcex;
 612 //        worker.sh_force2[0][i] = worker.sh_force2[0][i] - forcex;
 613         workingpad[0][i] = workingpad[0][i] - forcex;
 614
 615         forcey = yy * r148;
 616         fyi = fyi + forcey;
 617
 618         // sh_force2[1][id][i] = sh_force2[1][id][i] - forcey;
 619 //        worker.sh_force2[1][i] = worker.sh_force2[1][i] - forcey;
 620         workingpad[1][i] = workingpad[1][i] - forcey;
 621
 622         forcez = zz * r148;
 623         fzi = fzi + forcez;
 624
 625         // sh_force2[2][id][i] = sh_force2[2][id][i] - forcez;
 626 //        worker.sh_force2[2][i] = worker.sh_force2[2][i] - forcez;
 627         workingpad[2][i] = workingpad[2][i] - forcez;
 628
 629         // mymd.interacts[id]++;
 630         store.interacts++;
 631       }
 632
 633     }
 634
 635     // sh_force2[0][id][x] = sh_force2[0][id][x] + fxi;
 636     // sh_force2[1][id][x] = sh_force2[1][id][x] + fyi;
 637     // sh_force2[2][id][x] = sh_force2[2][id][x] + fzi;
 638
 639 //    worker.sh_force2[0][x] = worker.sh_force2[0][x] + fxi;
 640 //    worker.sh_force2[1][x] = worker.sh_force2[1][x] + fyi;
 641 //    worker.sh_force2[2][x] = worker.sh_force2[2][x] + fzi;
 642
 643     workingpad[0][x] = workingpad[0][x] + fxi;
 644     workingpad[1][x] = workingpad[1][x] + fyi;
 645     workingpad[2][x] = workingpad[2][x] + fzi;
 646
 647   }
 648
 649   public double mkekin(double hsq2, int part_id) {
 650
 651     double sumt = 0.0;
 652
 653     xvelocity = xvelocity + sh_force[0][part_id];
 654     yvelocity = yvelocity + sh_force[1][part_id];
 655     zvelocity = zvelocity + sh_force[2][part_id];
 656
 657     sumt = (xvelocity * xvelocity) + (yvelocity * yvelocity) + (zvelocity * zvelocity);
 658     return sumt;
 659   }
 660
 661   public double velavg(double vaverh, double h) {
 662
 663     double velt;
 664     double sq;
 665
 666     sq = Math.sqrt(xvelocity * xvelocity + yvelocity * yvelocity + zvelocity * zvelocity);
 667
 668     velt = sq;
 669     return velt;
 670   }
 671
 672   public void dscal(double sc, int incx) {
 673     xvelocity = xvelocity * sc;
 674     yvelocity = yvelocity * sc;
 675     zvelocity = zvelocity * sc;
 676   }
 677 }
 678
 679 class random {
 680
 681   public int iseed;
 682   public double v1, v2;
 683
 684   public random(int iseed, double v1, double v2) {
 685     this.iseed = iseed;
 686     this.v1 = v1;
 687     this.v2 = v2;
 688   }
 689
 690   public double update() {
 691
 692     double rand;
 693     double scale = 4.656612875e-10;
 694
 695     int is1, is2, iss2;
 696     int imult = 16807;
 697     int imod = 2147483647;
 698
 699     if (iseed <= 0) {
 700       iseed = 1;
 701     }
 702
 703     is2 = iseed % 32768;
 704     is1 = (iseed - is2) / 32768;
 705     iss2 = is2 * imult;
 706     is2 = iss2 % 32768;
 707     is1 = (is1 * imult + (iss2 - is2) / 32768) % (65536);
 708
 709     iseed = (is1 * 32768 + is2) % imod;
 710
 711     rand = scale * iseed;
 712
 713     return rand;
 714
 715   }
 716
 717   public double seed() {
 718
 719     double s, u1, u2, r;
 720     s = 1.0;
 721     do {
 722       u1 = update();
 723       u2 = update();
 724
 725       v1 = 2.0 * u1 - 1.0;
 726       v2 = 2.0 * u2 - 1.0;
 727       s = v1 * v1 + v2 * v2;
 728
 729     } while (s >= 1.0);
 730
 731     r = Math.sqrt(-2.0 * Math.log(s) / s);
 732
 733     return r;
 734
 735   }
 736 }