1 #include <../src/sys/classes/random/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 17 unsigned short _rander48_seed[3] = { 18 RANDER48_SEED_0, 19 RANDER48_SEED_1, 20 RANDER48_SEED_2 21 }; 22 unsigned short _rander48_mult[3] = { 23 RANDER48_MULT_0, 24 RANDER48_MULT_1, 25 RANDER48_MULT_2 26 }; 27 unsigned short _rander48_add = RANDER48_ADD; 28 29 void _dorander48(unsigned short xseed[3]) 30 { 31 unsigned long accu; 32 unsigned short temp[2]; 33 34 accu = (unsigned long) _rander48_mult[0] * (unsigned long) xseed[0] + (unsigned long)_rander48_add; 35 temp[0] = (unsigned short) accu; /* lower 16 bits */ 36 accu >>= sizeof(unsigned short) * 8; 37 accu += (unsigned long) _rander48_mult[0] * (unsigned long) xseed[1] + (unsigned long) _rander48_mult[1] * (unsigned long) xseed[0]; 38 temp[1] = (unsigned short)accu; /* middle 16 bits */ 39 accu >>= sizeof(unsigned short) * 8; 40 accu += _rander48_mult[0] * xseed[2] + _rander48_mult[1] * xseed[1] + _rander48_mult[2] * xseed[0]; 41 xseed[0] = temp[0]; 42 xseed[1] = temp[1]; 43 xseed[2] = (unsigned short) accu; 44 } 45 46 double erander48(unsigned short xseed[3]) 47 { 48 _dorander48(xseed); 49 return ldexp((double) xseed[0], -48) + ldexp((double) xseed[1], -32) + ldexp((double) xseed[2], -16); 50 } 51 52 double drander48() { 53 return erander48(_rander48_seed); 54 } 55 56 void srander48(long seed) { 57 _rander48_seed[0] = RANDER48_SEED_0; 58 _rander48_seed[1] = (unsigned short) seed; 59 _rander48_seed[2] = (unsigned short) (seed >> 16); 60 _rander48_mult[0] = RANDER48_MULT_0; 61 _rander48_mult[1] = RANDER48_MULT_1; 62 _rander48_mult[2] = RANDER48_MULT_2; 63 _rander48_add = RANDER48_ADD; 64 } 65 66 #undef __FUNCT__ 67 #define __FUNCT__ "PetscRandomSeed_Rander48" 68 PetscErrorCode PetscRandomSeed_Rander48(PetscRandom r) 69 { 70 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data; 71 72 PetscFunctionBegin; 73 r48->seed[0] = RANDER48_SEED_0; 74 r48->seed[1] = (unsigned short) r->seed; 75 r48->seed[2] = (unsigned short) (r->seed >> 16); 76 r48->mult[0] = RANDER48_MULT_0; 77 r48->mult[1] = RANDER48_MULT_1; 78 r48->mult[2] = RANDER48_MULT_2; 79 r48->add = RANDER48_ADD; 80 srander48(r->seed); 81 PetscFunctionReturn(0); 82 } 83 84 #undef __FUNCT__ 85 #define __FUNCT__ "PetscRandomGetValue_Rander48" 86 PetscErrorCode PetscRandomGetValue_Rander48(PetscRandom r, PetscScalar *val) 87 { 88 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data; 89 90 PetscFunctionBegin; 91 #if defined(PETSC_USE_COMPLEX) 92 if (r->iset) { 93 *val = PetscRealPart(r->low) + PetscImaginaryPart(r->low) * PETSC_i; 94 if (PetscRealPart(r->width)) { 95 *val += PetscRealPart(r->width)* erander48(r48->seed); 96 } 97 if (PetscImaginaryPart(r->width)) { 98 *val += PetscImaginaryPart(r->width)* erander48(r48->seed) * PETSC_i; 99 } 100 } else { 101 *val = erander48(r48->seed) + erander48(r48->seed)*PETSC_i; 102 } 103 #else 104 if (r->iset) *val = r->width * erander48(r48->seed) + r->low; 105 else *val = erander48(r48->seed); 106 #endif 107 PetscFunctionReturn(0); 108 } 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "PetscRandomGetValueReal_Rander48" 112 PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val) 113 { 114 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data; 115 116 PetscFunctionBegin; 117 #if defined(PETSC_USE_COMPLEX) 118 if (r->iset) *val = PetscRealPart(r->width)*erander48(r48->seed) + PetscRealPart(r->low); 119 else *val = erander48(r48->seed); 120 #else 121 if (r->iset) *val = r->width * erander48(r48->seed) + r->low; 122 else *val = erander48(r48->seed); 123 #endif 124 PetscFunctionReturn(0); 125 } 126 127 #undef __FUNCT__ 128 #define __FUNCT__ "PetscRandomDestroy_Rander48" 129 PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r) 130 { 131 PetscErrorCode ierr; 132 133 PetscFunctionBegin; 134 ierr = PetscFree(r->data);CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 static struct _PetscRandomOps PetscRandomOps_Values = { 139 /* 0 */ 140 PetscRandomSeed_Rander48, 141 PetscRandomGetValue_Rander48, 142 PetscRandomGetValueReal_Rander48, 143 PetscRandomDestroy_Rander48, 144 /* 5 */ 145 0 146 }; 147 148 /*MC 149 PETSCRANDER48 - simple portable reimplementation of basic Unix drand48() random number generator that should generate the 150 exact same random numbers on any system. 151 152 Options Database Keys: 153 . -random_type <rand,rand48,rander48,sprng> 154 155 Notes: This is the default random number generate provided by PetscRandomCreate() if you do not set a particular implementation. 156 157 Level: beginner 158 159 .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG 160 M*/ 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "PetscRandomCreate_Rander48" 164 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r) 165 { 166 PetscErrorCode ierr; 167 PetscRandom_Rander48 *r48; 168 169 PetscFunctionBegin; 170 ierr = PetscNewLog(r,&r48);CHKERRQ(ierr); 171 r->data = r48; 172 ierr = PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));CHKERRQ(ierr); 173 ierr = PetscObjectChangeTypeName((PetscObject) r, PETSCRANDER48);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176