Actual source code: rander48.c
petsc-3.13.6 2020-09-29
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: NULL
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: }