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 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 += r48->mult[0] * r48->seed[2] + r48->mult[1] * r48->seed[1] + 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 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(0); 48 } 49 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)) { 59 *val += PetscRealPart(r->width)* _dorander48(r48); 60 } 61 if (PetscImaginaryPart(r->width)) { 62 *val += PetscImaginaryPart(r->width)* _dorander48(r48) * PETSC_i; 63 } 64 } else { 65 *val = _dorander48(r48) + _dorander48(r48)*PETSC_i; 66 } 67 #else 68 if (r->iset) *val = r->width * _dorander48(r48) + r->low; 69 else *val = _dorander48(r48); 70 #endif 71 PetscFunctionReturn(0); 72 } 73 74 static PetscErrorCode PetscRandomGetValueReal_Rander48(PetscRandom r, PetscReal *val) 75 { 76 PetscRandom_Rander48 *r48 = (PetscRandom_Rander48*)r->data; 77 78 PetscFunctionBegin; 79 #if defined(PETSC_USE_COMPLEX) 80 if (r->iset) *val = PetscRealPart(r->width)*_dorander48(r48) + PetscRealPart(r->low); 81 else *val = _dorander48(r48); 82 #else 83 if (r->iset) *val = r->width * _dorander48(r48) + r->low; 84 else *val = _dorander48(r48); 85 #endif 86 PetscFunctionReturn(0); 87 } 88 89 static PetscErrorCode PetscRandomDestroy_Rander48(PetscRandom r) 90 { 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 ierr = PetscFree(r->data);CHKERRQ(ierr); 95 PetscFunctionReturn(0); 96 } 97 98 static struct _PetscRandomOps PetscRandomOps_Values = { 99 /* 0 */ 100 PetscRandomSeed_Rander48, 101 PetscRandomGetValue_Rander48, 102 PetscRandomGetValueReal_Rander48, 103 PetscRandomDestroy_Rander48, 104 /* 5 */ 105 NULL 106 }; 107 108 /*MC 109 PETSCRANDER48 - simple portable reimplementation of basic Unix drand48() random number generator that should generate the 110 exact same random numbers on any system. 111 112 Options Database Keys: 113 . -random_type <rand,rand48,rander48,sprng> 114 115 Notes: 116 This is the default random number generate provided by PetscRandomCreate() if you do not set a particular implementation. 117 118 Each PetscRandom object created with this type has its own seed and its own history, so multiple PetscRandom objects of this type 119 will not interfer with random numbers generated by other objects. Each PETSc object of this type will produce the exact same set of 120 random numbers so if you wish different PetscObjects of this type set different seeds for each one after you create them with 121 PetscRandomSetSeed() followed by PetscRandomSeed(). 122 123 Level: beginner 124 125 .seealso: RandomCreate(), RandomSetType(), PETSCRAND, PETSCRAND48, PETSCRANDER48, PETSCSPRNG, PetscRandomSetSeed(), PetscRandomSeed() 126 M*/ 127 128 PETSC_EXTERN PetscErrorCode PetscRandomCreate_Rander48(PetscRandom r) 129 { 130 PetscErrorCode ierr; 131 PetscRandom_Rander48 *r48; 132 133 PetscFunctionBegin; 134 ierr = PetscNewLog(r,&r48);CHKERRQ(ierr); 135 /* r48 does not need to be initialized because PetscRandomSeed() is always called before use and sets the needed values */ 136 r->data = r48; 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