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->width)*(PetscReal) drander48() + PetscRealPart(r->low) + (PetscImaginaryPart(r->width)*(PetscReal) drander48() + PetscImaginaryPart(r->low)) * PETSC_i; 78 } else { 79 *val = (PetscReal) drander48() + ((PetscReal) drander48())*PETSC_i; 80 } 81 #else 82 if (r->iset) *val = r->width * drander48() + r->low; 83 else *val = drander48(); 84 #endif 85 PetscFunctionReturn(0); 86 } 87 88 #undef __FUNCT__ 89 #define __FUNCT__ "PetscRandomGetValueReal_Rander48" 90 PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val) 91 { 92 PetscFunctionBegin; 93 #if defined(PETSC_USE_COMPLEX) 94 if (r->iset) *val = PetscRealPart(r->width)*drander48() + PetscRealPart(r->low); 95 else *val = drander48(); 96 #else 97 if (r->iset) *val = r->width * drander48() + r->low; 98 else *val = drander48(); 99 #endif 100 PetscFunctionReturn(0); 101 } 102 103 static struct _PetscRandomOps PetscRandomOps_Values = { 104 /* 0 */ 105 PetscRandomSeed_Rander48, 106 PetscRandomGetValue_Rander48, 107 PetscRandomGetValueReal_Rander48, 108 0, 109 /* 5 */ 110 0 111 }; 112 113 /*MC 114 PETSCRANDER48 - simple reimplementation of basic Unix drand48() random number generator, because Microsoft 115 116 Options Database Keys: 117 . -random_type <rand,rand48,rander48,sprng> 118 119 Level: beginner 120 121 .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG 122 M*/ 123 124 #undef __FUNCT__ 125 #define __FUNCT__ "PetscRandomCreate_Rander48" 126 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r) 127 { 128 PetscErrorCode ierr; 129 130 PetscFunctionBegin; 131 ierr = PetscMemcpy(r->ops, &PetscRandomOps_Values, sizeof(PetscRandomOps_Values));CHKERRQ(ierr); 132 ierr = PetscObjectChangeTypeName((PetscObject) r, PETSCRANDER48);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135