1 #include <petsc/private/randomimpl.h>
2
3 typedef struct {
4 unsigned short seed[3];
5 unsigned short mult[3];
6 unsigned short add;
7 } PetscRandom_Rander48;
8
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)
16
_dorander48(PetscRandom_Rander48 * r48)17 static double _dorander48(PetscRandom_Rander48 *r48)
18 {
19 unsigned long accu;
20 unsigned short temp[2];
21
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 }
34
PetscRandomSeed_Rander48(PetscRandom r)35 static PetscErrorCode PetscRandomSeed_Rander48(PetscRandom r)
36 {
37 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;
38
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 }
49
PetscRandomGetValue_Rander48(PetscRandom r,PetscScalar * val)50 static PetscErrorCode PetscRandomGetValue_Rander48(PetscRandom r, PetscScalar *val)
51 {
52 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;
53
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 }
69
PetscRandomGetValueReal_Rander48(PetscRandom r,PetscReal * val)70 static PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val)
71 {
72 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48 *)r->data;
73
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 }
84
PetscRandomDestroy_Rander48(PetscRandom r)85 static PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r)
86 {
87 PetscFunctionBegin;
88 PetscCall(PetscFree(r->data));
89 PetscFunctionReturn(PETSC_SUCCESS);
90 }
91
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 };
101
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.
105
106 Options Database Key:
107 . -random_type <rand,rand48,rander48,sprng> - select the random number generator at runtime
108
109 Level: beginner
110
111 Notes:
112 This is the default random number generate provided by `PetscRandomCreate()` if you do not set a particular implementation.
113
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()`.
118
119 .seealso: `PetscRandomCreate()`, `PetscRandomSetType()`, `PETSCRAND`, `PETSCRAND48`, `PETSCRANDER48`, `PETSCSPRNG`, `PetscRandomSetSeed()`,
120 `PetscRandomSeed()`, `PetscRandomSetFromOptions()`
121 M*/
122
PetscRandomCreate_Rander48(PetscRandom r)123 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r)
124 {
125 PetscRandom_Rander48 *r48;
126
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 }
135