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 += (unsigned long)r48->mult[0] * r48->seed[2] + (unsigned long)r48->mult[1] * r48->seed[1] + (unsigned long)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:   PetscFunctionBegin;
 40:   r48->seed[0] = RANDER48_SEED_0;
 41:   r48->seed[1] = (unsigned short)r->seed;
 42:   r48->seed[2] = (unsigned short)(r->seed >> 16);
 43:   r48->mult[0] = RANDER48_MULT_0;
 44:   r48->mult[1] = RANDER48_MULT_1;
 45:   r48->mult[2] = RANDER48_MULT_2;
 46:   r48->add     = RANDER48_ADD;
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

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

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

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

 74:   PetscFunctionBegin;
 75: #if defined(PETSC_USE_COMPLEX)
 76:   if (r->iset) *val = PetscRealPart(r->width) * _dorander48(r48) + PetscRealPart(r->low);
 77:   else *val = _dorander48(r48);
 78: #else
 79:   if (r->iset) *val = r->width * _dorander48(r48) + r->low;
 80:   else *val = _dorander48(r48);
 81: #endif
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: static PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r)
 86: {
 87:   PetscFunctionBegin;
 88:   PetscCall(PetscFree(r->data));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static struct _PetscRandomOps PetscRandomOps_Values = {
 93:   PetscDesignatedInitializer(seed, PetscRandomSeed_Rander48),
 94:   PetscDesignatedInitializer(getvalue, PetscRandomGetValue_Rander48),
 95:   PetscDesignatedInitializer(getvaluereal, PetscRandomGetValueReal_Rander48),
 96:   PetscDesignatedInitializer(getvalues, NULL),
 97:   PetscDesignatedInitializer(getvaluesreal, NULL),
 98:   PetscDesignatedInitializer(destroy, PetscRandomDestroy_Rander48),
 99:   PetscDesignatedInitializer(setfromoptions, NULL),
100: };

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

106:    Options Database Key:
107: . -random_type <rand,rand48,rander48,sprng> - select the random number generator at runtime

109:   Level: beginner

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

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

119: .seealso: `PetscRandomCreate()`, `PetscRandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCRANDER48`, `PETSCSPRNG`, `PetscRandomSetSeed()`,
120:           `PetscRandomSeed()`, `PetscRandomSetFromOptions()`
121: M*/

123: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r)
124: {
125:   PetscRandom_Rander48 *r48;

127:   PetscFunctionBegin;
128:   PetscCall(PetscNew(&r48));
129:   /* r48 does not need to be initialized because PetscRandomSeed() is always called before use and sets the needed values */
130:   r->data   = r48;
131:   r->ops[0] = PetscRandomOps_Values;
132:   PetscCall(PetscObjectChangeTypeName((PetscObject)r, PETSCRANDER48));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }