Actual source code: rander48.c

  1: #include <petsc/private/randomimpl.h>

  3: typedef struct {
  4:   unsigned short seed[3];
  5:   unsigned short mult[3];
  6:   unsigned short add;
  7: } PetscRandom_Rander48;

  9: #define RANDER48_SEED_0 (0x330e)
 10: #define RANDER48_SEED_1 (0xabcd)
 11: #define RANDER48_SEED_2 (0x1234)
 12: #define RANDER48_MULT_0 (0xe66d)
 13: #define RANDER48_MULT_1 (0xdeec)
 14: #define RANDER48_MULT_2 (0x0005)
 15: #define RANDER48_ADD    (0x000b)

 17: static double _dorander48(PetscRandom_Rander48 *r48)
 18: {
 19:   unsigned long  accu;
 20:   unsigned short temp[2];

 22:   accu    = (unsigned long)r48->mult[0] * (unsigned long)r48->seed[0] + (unsigned long)r48->add;
 23:   temp[0] = (unsigned short)accu; /* lower 16 bits */
 24:   accu >>= sizeof(unsigned short) * 8;
 25:   accu += (unsigned long)r48->mult[0] * (unsigned long)r48->seed[1] + (unsigned long)r48->mult[1] * (unsigned long)r48->seed[0];
 26:   temp[1] = (unsigned short)accu; /* middle 16 bits */
 27:   accu >>= sizeof(unsigned short) * 8;
 28:   accu += r48->mult[0] * r48->seed[2] + r48->mult[1] * r48->seed[1] + r48->mult[2] * r48->seed[0];
 29:   r48->seed[0] = temp[0];
 30:   r48->seed[1] = temp[1];
 31:   r48->seed[2] = (unsigned short)accu;
 32:   return ldexp((double)r48->seed[0], -48) + ldexp((double)r48->seed[1], -32) + ldexp((double)r48->seed[2], -16);
 33: }

 35: static PetscErrorCode PetscRandomSeed_Rander48(PetscRandom r)
 36: {
 37:   PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;

 39:   r48->seed[0] = RANDER48_SEED_0;
 40:   r48->seed[1] = (unsigned short)r->seed;
 41:   r48->seed[2] = (unsigned short)(r->seed >> 16);
 42:   r48->mult[0] = RANDER48_MULT_0;
 43:   r48->mult[1] = RANDER48_MULT_1;
 44:   r48->mult[2] = RANDER48_MULT_2;
 45:   r48->add     = RANDER48_ADD;
 46:   return 0;
 47: }

 49: static PetscErrorCode PetscRandomGetValue_Rander48(PetscRandom r, PetscScalar *val)
 50: {
 51:   PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;

 53: #if defined(PETSC_USE_COMPLEX)
 54:   if (r->iset) {
 55:     *val = PetscRealPart(r->low) + PetscImaginaryPart(r->low) * PETSC_i;
 56:     if (PetscRealPart(r->width)) *val += PetscRealPart(r->width) * _dorander48(r48);
 57:     if (PetscImaginaryPart(r->width)) *val += PetscImaginaryPart(r->width) * _dorander48(r48) * PETSC_i;
 58:   } else {
 59:     *val = _dorander48(r48) + _dorander48(r48) * PETSC_i;
 60:   }
 61: #else
 62:   if (r->iset) *val = r->width * _dorander48(r48) + r->low;
 63:   else *val = _dorander48(r48);
 64: #endif
 65:   return 0;
 66: }

 68: static PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val)
 69: {
 70:   PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;

 72: #if defined(PETSC_USE_COMPLEX)
 73:   if (r->iset) *val = PetscRealPart(r->width) * _dorander48(r48) + PetscRealPart(r->low);
 74:   else *val = _dorander48(r48);
 75: #else
 76:   if (r->iset) *val = r->width * _dorander48(r48) + r->low;
 77:   else *val = _dorander48(r48);
 78: #endif
 79:   return 0;
 80: }

 82: static PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r)
 83: {
 84:   PetscFree(r->data);
 85:   return 0;
 86: }

 88: static struct _PetscRandomOps PetscRandomOps_Values = {
 89:   PetscDesignatedInitializer(seed, PetscRandomSeed_Rander48),
 90:   PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Rander48),
 91:   PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Rander48),
 92:   PetscDesignatedInitializer(getvalues, NULL),
 93:   PetscDesignatedInitializer(getvaluesreal, NULL),
 94:   PetscDesignatedInitializer(destroy, PetscRandomDestroy_Rander48),
 95: };

 97: /*MC
 98:    PETSCRANDER48 - simple portable reimplementation of basic Unix drand48() random number generator that should generate the
 99:         exact same random numbers on any system.

101:    Options Database Keys:
102: . -random_type <rand,rand48,rander48,sprng> - select the random number generator at runtime

104:   Notes:
105:     This is the default random number generate provided by `PetscRandomCreate()` if you do not set a particular implementation.

107:   Each `PetscRandom` object created with this type has its own seed and its own history, so multiple `PetscRandom` objects of this type
108:   will not interfere with random numbers generated by other objects. Each PETSc object of this type will produce the exact same set of
109:   random numbers so if you wish different `PetscRandom` objects of this type set different seeds for each one after you create them with
110:   `PetscRandomSetSeed()` followed by `PetscRandomSeed()`.

112:   Level: beginner

114: .seealso: `PetscRandomCreate()`, `PetscRandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCRANDER48`, `PETSCSPRNG`, `PetscRandomSetSeed()`,
115:           `PetscRandomSeed()`, `PetscRandomSetFromOptions()`
116: M*/

118: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r)
119: {
120:   PetscRandom_Rander48 *r48;

122:   PetscNew(&r48);
123:   /* r48 does not need to be initialized because PetscRandomSeed() is always called before use and sets the needed values */
124:   r->data = r48;
125:   PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));
126:   PetscObjectChangeTypeName((PetscObject)r, PETSCRANDER48);
127:   return 0;
128: }