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)) {
57: *val += PetscRealPart(r->width)* _dorander48(r48);
58: }
59: if (PetscImaginaryPart(r->width)) {
60: *val += PetscImaginaryPart(r->width)* _dorander48(r48) * PETSC_i;
61: }
62: } else {
63: *val = _dorander48(r48) + _dorander48(r48)*PETSC_i;
64: }
65: #else
66: if (r->iset) *val = r->width * _dorander48(r48) + r->low;
67: else *val = _dorander48(r48);
68: #endif
69: return 0;
70: }
72: static PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val)
73: {
74: PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data;
76: #if defined(PETSC_USE_COMPLEX)
77: if (r->iset) *val = PetscRealPart(r->width)*_dorander48(r48) + PetscRealPart(r->low);
78: else *val = _dorander48(r48);
79: #else
80: if (r->iset) *val = r->width * _dorander48(r48) + r->low;
81: else *val = _dorander48(r48);
82: #endif
83: return 0;
84: }
86: static PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r)
87: {
88: PetscFree(r->data);
89: return 0;
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: };
101: /*MC
102: PETSCRANDER48 - simple portable reimplementation of basic Unix drand48() random number generator that should generate the
103: exact same random numbers on any system.
105: Options Database Keys:
106: . -random_type <rand,rand48,rander48,sprng>
108: Notes:
109: This is the default random number generate provided by PetscRandomCreate() if you do not set a particular implementation.
111: Each PetscRandom object created with this type has its own seed and its own history, so multiple PetscRandom objects of this type
112: will not interfer with random numbers generated by other objects. Each PETSc object of this type will produce the exact same set of
113: random numbers so if you wish different PetscObjects of this type set different seeds for each one after you create them with
114: PetscRandomSetSeed() followed by PetscRandomSeed().
116: Level: beginner
118: .seealso: PetscRandomCreate(), PetscRandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG, PetscRandomSetSeed(), PetscRandomSeed()
119: M*/
121: PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r)
122: {
123: PetscRandom_Rander48 *r48;
125: PetscNewLog(r,&r48);
126: /* r48 does not need to be initialized because PetscRandomSeed() is always called before use and sets the needed values */
127: r->data = r48;
128: PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));
129: PetscObjectChangeTypeName((PetscObject) r, PETSCRANDER48);
130: return 0;
131: }