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