1 #include <../src/sys/classes/random/randomimpl.h> 2 3 /* This shit is even more random than drand48(). It is random enough to run on Windows. Fuck you Microsoft, and your utter failure to implement POSIX. */ 4 #define RANDER48_SEED_0 (0x330e) 5 #define RANDER48_SEED_1 (0xabcd) 6 #define RANDER48_SEED_2 (0x1234) 7 #define RANDER48_MULT_0 (0xe66d) 8 #define RANDER48_MULT_1 (0xdeec) 9 #define RANDER48_MULT_2 (0x0005) 10 #define RANDER48_ADD (0x000b) 11 12 unsigned short _rander48_seed[3] = { 13 RANDER48_SEED_0, 14 RANDER48_SEED_1, 15 RANDER48_SEED_2 16 }; 17 unsigned short _rander48_mult[3] = { 18 RANDER48_MULT_0, 19 RANDER48_MULT_1, 20 RANDER48_MULT_2 21 }; 22 unsigned short _rander48_add = RANDER48_ADD; 23 24 void _dorander48(unsigned short xseed[3]) 25 { 26 unsigned long accu; 27 unsigned short temp[2]; 28 29 accu = (unsigned long) _rander48_mult[0] * (unsigned long) xseed[0] + (unsigned long)_rander48_add; 30 temp[0] = (unsigned short) accu; /* lower 16 bits */ 31 accu >>= sizeof(unsigned short) * 8; 32 accu += (unsigned long) _rander48_mult[0] * (unsigned long) xseed[1] + (unsigned long) _rander48_mult[1] * (unsigned long) xseed[0]; 33 temp[1] = (unsigned short)accu; /* middle 16 bits */ 34 accu >>= sizeof(unsigned short) * 8; 35 accu += _rander48_mult[0] * xseed[2] + _rander48_mult[1] * xseed[1] + _rander48_mult[2] * xseed[0]; 36 xseed[0] = temp[0]; 37 xseed[1] = temp[1]; 38 xseed[2] = (unsigned short) accu; 39 } 40 41 double erander48(unsigned short xseed[3]) 42 { 43 _dorander48(xseed); 44 return ldexp((double) xseed[0], -48) + ldexp((double) xseed[1], -32) + ldexp((double) xseed[2], -16); 45 } 46 47 double drander48() { 48 return erander48(_rander48_seed); 49 } 50 51 void srander48(long seed) { 52 _rander48_seed[0] = RANDER48_SEED_0; 53 _rander48_seed[1] = (unsigned short) seed; 54 _rander48_seed[2] = (unsigned short) (seed >> 16); 55 _rander48_mult[0] = RANDER48_MULT_0; 56 _rander48_mult[1] = RANDER48_MULT_1; 57 _rander48_mult[2] = RANDER48_MULT_2; 58 _rander48_add = RANDER48_ADD; 59 } 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "PetscRandomSeed_Rander48" 63 PetscErrorCode PetscRandomSeed_Rander48(PetscRandom r) 64 { 65 PetscFunctionBegin; 66 srander48(r->seed); 67 PetscFunctionReturn(0); 68 } 69 70 #undef __FUNCT__ 71 #define __FUNCT__ "PetscRandomGetValue_Rander48" 72 PetscErrorCode PetscRandomGetValue_Rander48(PetscRandom r, PetscScalar *val) 73 { 74 PetscFunctionBegin; 75 #if defined(PETSC_USE_COMPLEX) 76 if (r->iset) { 77 *val = PetscRealPart(r->low) + PetscImaginaryPart(r->low) * PETSC_i; 78 if (PetscRealPart(r->width)) { 79 *val += PetscRealPart(r->width)* drander48(); 80 } 81 if (PetscImaginaryPart(r->width)) { 82 *val += PetscImaginaryPart(r->width)* drander48() * PETSC_i; 83 } 84 } else { 85 *val = drander48() + drander48()*PETSC_i; 86 } 87 #else 88 if (r->iset) *val = r->width * drander48() + r->low; 89 else *val = drander48(); 90 #endif 91 PetscFunctionReturn(0); 92 } 93 94 #undef __FUNCT__ 95 #define __FUNCT__ "PetscRandomGetValueReal_Rander48" 96 PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val) 97 { 98 PetscFunctionBegin; 99 #if defined(PETSC_USE_COMPLEX) 100 if (r->iset) *val = PetscRealPart(r->width)*drander48() + PetscRealPart(r->low); 101 else *val = drander48(); 102 #else 103 if (r->iset) *val = r->width * drander48() + r->low; 104 else *val = drander48(); 105 #endif 106 PetscFunctionReturn(0); 107 } 108 109 static struct _PetscRandomOps PetscRandomOps_Values = { 110 /* 0 */ 111 PetscRandomSeed_Rander48, 112 PetscRandomGetValue_Rander48, 113 PetscRandomGetValueReal_Rander48, 114 0, 115 /* 5 */ 116 0 117 }; 118 119 /*MC 120 PETSCRANDER48 - simple reimplementation of basic Unix drand48() random number generator, because Microsoft 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