SCM

SCM Repository

[rmetrics] Diff of /pkg/randtoolbox/src/randtoolbox.c
ViewVC logotype

Diff of /pkg/randtoolbox/src/randtoolbox.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 6135, Thu Nov 29 11:07:56 2018 UTC revision 6136, Tue Apr 30 12:56:05 2019 UTC
# Line 6  Line 6 
6   * @author Petr Savicky   * @author Petr Savicky
7   *   *
8   *   *
9   * Copyright (C) 2009, Christophe Dutang,   * Copyright (C) 2019, Christophe Dutang,
10   * Petr Savicky, Academy of Sciences of the Czech Republic.   * Petr Savicky, Academy of Sciences of the Czech Republic.
11     * Christophe Dutang, see http://dutangc.free.fr
12   * All rights reserved.   * All rights reserved.
13   *   *
14   * The new BSD License is applied to this software.   * The new BSD License is applied to this software.
15   * Copyright (c) 2009 Christophe Dutang, Petr Savicky.   * Copyright (c) 2019 Christophe Dutang, Petr Savicky.
16   * All rights reserved.   * All rights reserved.
17   *   *
18   *      Redistribution and use in source and binary forms, with or without   *      Redistribution and use in source and binary forms, with or without
# Line 165  Line 166 
166      //init the seed of SF Mersenne Twister algo      //init the seed of SF Mersenne Twister algo
167      SFMT_init_gen_rand(seed);      SFMT_init_gen_rand(seed);
168    
   
     //Rprintf("state %lu\n", seed);  
   
169      //compute      //compute
170      for(j = 0; j < dim; j++)      for(j = 0; j < dim; j++)
171      {      {
172        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
173        {        {
174          state = SFMT_gen_rand32();          state = SFMT_gen_rand32();
         //Rprintf("state %lu\n", state);  
   
175          u[i + j * nb] = fracPart( state * sqrt( prime[j] ) ) ;          u[i + j * nb] = fracPart( state * sqrt( prime[j] ) ) ;
176        }        }
177      }      }
# Line 189  Line 185 
185        state = ((unsigned int) seed >> 16);        state = ((unsigned int) seed >> 16);
186      }else      }else
187        state  = offset;        state  = offset;
   
   
     //Rprintf("state %u %lu\n", state, state);  
188      //compute      //compute
189      for(j = 0; j < dim; j++)      for(j = 0; j < dim; j++)
190        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
# Line 263  Line 256 
256        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
257        {        {
258          state = SFMT_gen_rand32();          state = SFMT_gen_rand32();
         Rprintf("state %u %lu\n", state, state);  
   
259          u[i + j * nb] = HALTONREC( j, state ) ;          u[i + j * nb] = HALTONREC( j, state ) ;
260        }        }
261      }      }
# Line 279  Line 270 
270      }else      }else
271        state  = offset;        state  = offset;
272    
     //Rprintf("state %u %lu\n", state, state);  
273      for(j = 0; j < dim; j++)      for(j = 0; j < dim; j++)
274        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
275          u[i + j * nb] = HALTONREC( j, state + i ) ;          u[i + j * nb] = HALTONREC( j, state + i ) ;
# Line 324  Line 314 
314    
315    
316    
317  //compute the vector sequence of the Halton algorithm  //compute the vector sequence of the Sobol algorithm
318  void sobol_c(double *u, int nb, int dim, int offset, int ismixed, int usetime, int mexp)  void sobol_c(double *u, int nb, int dim, int offset, int ismixed, int usetime, int mexp)
319  {  {
320    //temporary working variables    //temporary working variables
# Line 390  Line 380 
380    double incrtemp = asReal( increment ) ; //increment as a double numeric    double incrtemp = asReal( increment ) ; //increment as a double numeric
381    uint64_t mod, mult, incr, mask; //modulus, multiplier, increment, mask    uint64_t mod, mult, incr, mask; //modulus, multiplier, increment, mask
382    
   Rprintf("mod %f\n", modultemp);  
   
383    //define a mask as in congruRand.c (function put_state_congru())    //define a mask as in congruRand.c (function put_state_congru())
384    if (modultemp < two_64_d) //modulus lesser than 2^64 => mask=2^d-1 if mod = 2^d    if (modultemp < two_64_d) //modulus lesser than 2^64 => mask=2^d-1 if mod = 2^d
385    {    {
# Line 414  Line 402 
402    else    else
403      error(_("increment greater than 2^64-1"));      error(_("increment greater than 2^64-1"));
404    
   Rprintf("mask: %llu\n", mask);  
   Rprintf("modulus: %llu\n", mod);  
   Rprintf("multiplier: %llu\n", mult);  
   Rprintf("increment: %llu\n", incr);  
   
405    //result    //result
406    double *u ; //result in C    double *u ; //result in C
407    SEXP resultinR; //result in R    SEXP resultinR; //result in R
# Line 543  Line 526 
526    //int rest = nb % blocksize;    //int rest = nb % blocksize;
527    
528    
   
   //Rprintf("zog\n");  
   
   //    PROTECT(array = allocMatrix(INTSXP, nb, dim)); //allocate a n x d matrix  
   
   
529    /*    /*
530     * @param size the number of 32-bit pseudorandom integers to be     * @param size the number of 32-bit pseudorandom integers to be
531     * generated.  size must be a multiple of 4, and greater than or equal     * generated.  size must be a multiple of 4, and greater than or equal
# Line 557  Line 534 
534  #if defined(HAVE_SSE2)  #if defined(HAVE_SSE2)
535    if(nb * dim >= blocksize)    if(nb * dim >= blocksize)
536    {    {
   
     //unsigned int *array;  
     //        static __m128i array1[BLOCK_SIZE / 4];  
     //        static __m128i array2[10000 / 4];  
   
537      __m128i array1[nb * dim*4];      __m128i array1[nb * dim*4];
538    
539      uint32_t *array32 = (uint32_t *)array1;      uint32_t *array32 = (uint32_t *)array1;
540    
     //        array = (unsigned int *) R_alloc(nb * dim, sizeof(__m128i));  
   
541      fill_array32(array32, nb*dim);      fill_array32(array32, nb*dim);
542    
543      // compute u_ij      // compute u_ij
# Line 575  Line 545 
545        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
546          u[i + j * nb] = to_real3( array32[i + j * nb] ); // real on ]0,1[ interval          u[i + j * nb] = to_real3( array32[i + j * nb] ); // real on ]0,1[ interval
547    
   
     //Free(array);  
548    }    }
549    else    else
550    {    {
551  #endif  #endif
552    
   
553      // compute u_ij      // compute u_ij
554      for(j = 0; j < dim; j++)      for(j = 0; j < dim; j++)
555        for(i = 0; i < nb; i++)        for(i = 0; i < nb; i++)
# Line 669  Line 636 
636    //init TAOCP RNG    //init TAOCP RNG
637    ranf_start( seed );    ranf_start( seed );
638    
   //ranf_arr_cycle();  
   
639    // compute u_ij's    // compute u_ij's
640    // declare an array a little bit longer than KK (100) long lag if too short    // declare an array a little bit longer than KK (100) long lag if too short
641    // see Knuth's file for details    // see Knuth's file for details
# Line 688  Line 653 
653    else    else
654      ranf_array( u, nb*dim );      ranf_array( u, nb*dim );
655    
   //Rprintf("1st term %.20f --- seed  %u\n", u[0], seed);  
656    isInit = 0;    isInit = 0;
657  }  }
658    

Legend:
Removed from v.6135  
changed lines
  Added in v.6136

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge