1 #include <../src/sys/classes/random/randomimpl.h> 2 3 #define RANDER48_SEED_0 (0x330e) 4 #define RANDER48_SEED_1 (0xabcd) 5 #define RANDER48_SEED_2 (0x1234) 6 #define RANDER48_MULT_0 (0xe66d) 7 #define RANDER48_MULT_1 (0xdeec) 8 #define RANDER48_MULT_2 (0x0005) 9 #define RANDER48_ADD (0x000b) 10 11 unsigned short _rander48_seed[3] = { 12 RANDER48_SEED_0, 13 RANDER48_SEED_1, 14 RANDER48_SEED_2 15 }; 16 unsigned short _rander48_mult[3] = { 17 RANDER48_MULT_0, 18 RANDER48_MULT_1, 19 RANDER48_MULT_2 20 }; 21 unsigned short _rander48_add = RANDER48_ADD; 22 23 void _dorander48(unsigned short xseed[3]) 24 { 25 unsigned long accu; 26 unsigned short temp[2]; 27 28 accu = (unsigned long) _rander48_mult[0] * (unsigned long) xseed[0] + (unsigned long)_rander48_add; 29 temp[0] = (unsigned short) accu; /* lower 16 bits */ 30 accu >>= sizeof(unsigned short) * 8; 31 accu += (unsigned long) _rander48_mult[0] * (unsigned long) xseed[1] + (unsigned long) _rander48_mult[1] * (unsigned long) xseed[0]; 32 temp[1] = (unsigned short)accu; /* middle 16 bits */ 33 accu >>= sizeof(unsigned short) * 8; 34 accu += _rander48_mult[0] * xseed[2] + _rander48_mult[1] * xseed[1] + _rander48_mult[2] * xseed[0]; 35 xseed[0] = temp[0]; 36 xseed[1] = temp[1]; 37 xseed[2] = (unsigned short) accu; 38 } 39 40 double erander48(unsigned short xseed[3]) 41 { 42 _dorander48(xseed); 43 return ldexp((double) xseed[0], -48) + ldexp((double) xseed[1], -32) + ldexp((double) xseed[2], -16); 44 } 45 46 double drander48() { 47 return erander48(_rander48_seed); 48 } 49 50 void srander48(long seed) { 51 _rander48_seed[0] = RANDER48_SEED_0; 52 _rander48_seed[1] = (unsigned short) seed; 53 _rander48_seed[2] = (unsigned short) (seed >> 16); 54 _rander48_mult[0] = RANDER48_MULT_0; 55 _rander48_mult[1] = RANDER48_MULT_1; 56 _rander48_mult[2] = RANDER48_MULT_2; 57 _rander48_add = RANDER48_ADD; 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "PetscRandomSeed_Rander48" 62 PetscErrorCode PetscRandomSeed_Rander48(PetscRandom r) 63 { 64 PetscFunctionBegin; 65 srander48(r->seed); 66 PetscFunctionReturn(0); 67 } 68 69 #undef __FUNCT__ 70 #define __FUNCT__ "PetscRandomGetValue_Rander48" 71 PetscErrorCode PetscRandomGetValue_Rander48(PetscRandom r, PetscScalar *val) 72 { 73 PetscFunctionBegin; 74 #if defined(PETSC_USE_COMPLEX) 75 if (r->iset) { 76 *val = PetscRealPart(r->low) + PetscImaginaryPart(r->low) * PETSC_i; 77 if (PetscRealPart(r->width)) { 78 *val += PetscRealPart(r->width)* drander48(); 79 } 80 if (PetscImaginaryPart(r->width)) { 81 *val += PetscImaginaryPart(r->width)* drander48() * PETSC_i; 82 } 83 } else { 84 *val = drander48() + drander48()*PETSC_i; 85 } 86 #else 87 if (r->iset) *val = r->width * drander48() + r->low; 88 else *val = drander48(); 89 #endif 90 PetscFunctionReturn(0); 91 } 92 93 #undef __FUNCT__ 94 #define __FUNCT__ "PetscRandomGetValueReal_Rander48" 95 PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val) 96 { 97 PetscFunctionBegin; 98 #if defined(PETSC_USE_COMPLEX) 99 if (r->iset) *val = PetscRealPart(r->width)*drander48() + PetscRealPart(r->low); 100 else *val = drander48(); 101 #else 102 if (r->iset) *val = r->width * drander48() + r->low; 103 else *val = drander48(); 104 #endif 105 PetscFunctionReturn(0); 106 } 107 108 static struct _PetscRandomOps PetscRandomOps_Values = { 109 /* 0 */ 110 PetscRandomSeed_Rander48, 111 PetscRandomGetValue_Rander48, 112 PetscRandomGetValueReal_Rander48, 113 0, 114 /* 5 */ 115 0 116 }; 117 118 /*MC 119 PETSCRANDER48 - simple portable reimplementation of basic Unix drand48() random number generator that should generate the 120 exact same random numbers on any system. 121 122 Options Database Keys: 123 . -random_type <rand,rand48,rander48,sprng> 124 125 Level: beginner 126 127 .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG 128 M*/ 129 130 #undef __FUNCT__ 131 #define __FUNCT__ "PetscRandomCreate_Rander48" 132 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r) 133 { 134 PetscErrorCode ierr; 135 136 PetscFunctionBegin; 137 ierr = PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));CHKERRQ(ierr); 138 ierr = PetscObjectChangeTypeName((PetscObject) r, PETSCRANDER48);CHKERRQ(ierr); 139 PetscFunctionReturn(0); 140 } 141