1: #include <../src/sys/classes/random/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;

 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:   return(0);
 48: }

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

 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)) {
 59:       *val += PetscRealPart(r->width)* _dorander48(r48);
 60:     }
 61:     if (PetscImaginaryPart(r->width)) {
 62:       *val += PetscImaginaryPart(r->width)* _dorander48(r48) * PETSC_i;
 63:     }
 64:   } else {
 65:     *val = _dorander48(r48) +  _dorander48(r48)*PETSC_i;
 66:   }
 67: #else
 68:   if (r->iset) *val = r->width * _dorander48(r48) + r->low;
 69:   else         *val = _dorander48(r48);
 70: #endif
 71:   return(0);
 72: }

 74: static PetscErrorCode  PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val)
 75: {
 76:   PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data;

 79: #if defined(PETSC_USE_COMPLEX)
 80:   if (r->iset) *val = PetscRealPart(r->width)*_dorander48(r48) + PetscRealPart(r->low);
 81:   else         *val = _dorander48(r48);
 82: #else
 83:   if (r->iset) *val = r->width * _dorander48(r48) + r->low;
 84:   else         *val = _dorander48(r48);
 85: #endif
 86:   return(0);
 87: }

 89: static PetscErrorCode  PetscRandomDestroy_Rander48(PetscRandom r)
 90: {

 94:   PetscFree(r->data);
 95:   return(0);
 96: }

 98: static struct _PetscRandomOps PetscRandomOps_Values = {
 99:   /* 0 */
100:   PetscRandomSeed_Rander48,
101:   PetscRandomGetValue_Rander48,
102:   PetscRandomGetValueReal_Rander48,
103:   PetscRandomDestroy_Rander48,
104:   /* 5 */
105:   0
106: };

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

112:    Options Database Keys:
113: . -random_type <rand,rand48,rander48,sprng>

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

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

123:   Level: beginner

125: .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG, PetscRandomSetSeed(), PetscRandomSeed()
126: M*/

128: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r)
129: {
130:   PetscErrorCode       ierr;
131:   PetscRandom_Rander48 *r48;

134:   PetscNewLog(r,&r48);
135:   /* r48 does not need to be initialized because PetscRandomSeed() is always called before use and sets the needed values */
136:   r->data = r48;
137:   PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));
138:   PetscObjectChangeTypeName((PetscObject) r, PETSCRANDER48);
139:   return(0);
140: }