1: .text 2: .globl _ranm 3: _ranm_: 4: 5: 6: .text 7: _ranm: 8: jsr r5,csv 9: mov r5,sr5 10: jsr pc,ranm 11: mov r0,t0 12: mov r1,t1 13: mov r2,t2 14: mov r3,t3 15: movf t0, fr0 16: mov sr5,r5 17: jmp cret 18: 19: 20: 21: .bss 22: .even 23: t0: .=.+2 24: t1: .=.+2 25: t2: .=.+2 26: t3: .=.+2 27: sr5: .=.+2 28: 29: // rany: random number generator 30: //created by d. lehmer and dave hutchinson 31: //period is about 2 billion 32: //even the low order bits are wonderfully random 33: //coded up by d. w. krumme, june 1977 34: 35: //entry iran 36: // initialize generators 37: // the six integer arguments are used as the starting numbers 38: // 39: .globl _iran 40: .text 41: _iran: 42: jsr r5,csv 43: mov 4(r5),r0 44: mov 6(r5),r1 45: mov 10(r5),r2 46: jsr pc,2f 47: mov r0,xno1 48: mov r1,rno1 49: mov r2,rno1+2 50: mov 12(r5),r0 51: mov 14(r5),r1 52: mov 16(r5),r2 53: jsr pc,2f 54: mov r0,xno2 55: mov r1,rno2 56: mov r2,rno2+2 57: inc initfg 58: 59: 2: 60: bic $100000,r0 61: bic $100000,r1 62: bis $1,r0 63: bis $1,r2 64: rts pc 65: 66: // filtab 67: // fill up the mixing table 68: 69: filtab: 70: mov $mixtab,r5 71: 1: /alternate enerators in pairs 72: mov $rno1,r4 73: jsr pc,tabone 74: jsr pc,tabone 75: mov $rno2,r4 76: jsr pc,tabone 77: jsr pc,tabone 78: cmp r5,$etable 79: blo 1b 80: clr initfg 81: rts pc 82: tabone: 83: mov r5,-(sp) 84: mov r4,r5 85: add $4,r5 86: jsr pc,modmul 87: mov (sp)+,r5 88: mov r0,(r5)+ 89: mov r1,(r5)+ 90: rts pc 91: 92: // entry ranm 93: // random numbers with mixed up sequence 94: 95: .globl ranm 96: ranm: 97: tst initfg 98: beq 1f 99: jsr pc,filtab /on the first call, fill up the table 100: 1: mov $xmul1,r4 101: jsr pc,mixint /get high 31 bits 102: mov r3,-(sp) /and save 103: mov r2,-(sp) 104: mov $xmul2,r4 105: jsr pc,mixint /get low 31 bits 106: mov (sp)+,r0 107: mov (sp)+,r1 108: ashc $1,r2 /concatenate high and low 109: /now normalize 110: mov $200,r5 /exponent zero 111: norm: ashc $1,r2 112: rol r1 113: rol r0 114: bmi 1f 115: sob r5,norm 116: /here all would be zero 117: /now form floating point 118: 1: 119: movb r1,r4 120: ashc $-8.,r0 121: ashc $-8.,r2 122: swab r2 123: clrb r2 124: bisb r4,r2 125: swab r2 126: bic $177600,r0 127: swab r5 /clears c-bit 128: ror r5 129: bis r5,r0 130: rts pc 131: 132: // mixint 133: // get a 31-bit integer by way of the mixing generator and table 134: 135: mixint: 136: mov (r4)+,r0 137: mul (r4),r0 /run the mixing generator 138: bic $100000,r1 139: mov r1,(r4)+ 140: mov r4,r5 /r4: rno 141: add $4,r5 /r5: rmul 142: jsr pc,modmul /get next integer in r0,r1 143: movb -1(r4),r5 /get 7 bits from the mixing genertor 144: ash $2,r5 145: add $mixtab,r5 / run integer through table 146: mov (r5),r2 147: mov r0,(r5)+ 148: mov (r5),r3 149: mov r1,(r5)+ 150: rts pc /return result in r2,r3 151: .data 152: .even 153: statev: / 1 + 12 = 13, + 256 = 269 words 154: 155: initfg: 1 156: 157: mixtab: .=.+512. / 256 words 158: etable: 159: 160: xmul1: 12555 161: xno1: 12555 162: rno1: 24223;046343 163: rmul1: 24223;046343 /680742115 = 7**23 (mod 2**31 - 1) 164: 165: xmul2: 13265 166: xno2: 13265 167: rno2: 1450;012656 168: rmul2: 1450;012656 /52958638 = 7**17 (mod 2**31 -1) 169: 170: // since the low words of both multipliers have bit 15 clear, 171: // the corection factor code for them has been commented out below 172: // if the multiliers are changed, be careful 173: 174: 175: 176: // multiply two 31-bit integers modulo 2**31 - 1 177: // this procedure generates random integers exhausting [1,2**31-2] 178: // when the right multipliers are used 179: 180: .text 181: modmul: 182: 183: // the following code multiplies two 31-bit (positive) integers 184: // producing a 62-bit result in r0,r1,r2,r3 185: // the numbers are wa + b and wc + d, where w=2**16 186: // the full 32 bit unsigned multiplication were desired, more correction 187: // factors would have to be included, so that each multiplication 188: // would be done as bd is done here 189: // correction factors are necessary because the pdp-11 multiplier interprrets 190: // negative numbers. mul x,y gives the 32-bit results: xy/ -x(w-y)/ 191: // -(w-x)y/ (w-x)(w-y)-w**2 0<x,y/ y<0<x/ x<0<y/ x,y<0 respectively 192: // a is at (r4), b at 2(r4) 193: // c is at (r5), d at 2(r5) 194: 195: mov (r4),r0 /ac 196: mul (r5)+,r0 /no correction 197: mov (r4),r2 /ad 198: mul (r5),r2 199: // tst (r5) /d is known to be ok! 200: // bpl 1f 201: // add (r4),r2 /correct d<0 202: 1: 203: add r2,r1 204: adc r0 205: mov r3,-(sp) /keep this slot on the stack 206: 207: tst (r4)+ /bc 208: mov (r4),r2 209: mul -(r5),r2 210: tst (r4) 211: bpl 2f 212: add (r5),r2 /correct b<0 213: 2: 214: add r2,r1 215: adc r0 216: add r3,(sp) 217: adc r1 218: adc r0 219: 220: mov (r4),r2 /bd 221: tst (r5)+ 222: mul (r5),r2 223: tst (r4) 224: bpl 3f 225: add (r5),r2 /correct b<0 226: 3:/ tst (r5) /d is known to be ok! 227: // bpl 4f 228: // add (r4),r2 /correct d<0 229: 4: 230: add (sp)+,r2 231: adc r1 232: adc r0 /62-bit result in r0-r3 233: 234: // the following code takes modulo 2**31 - 1 by casting out (2**31-1)'s in base 2**31 235: // r0 must have top two bits zero 236: // the result is left in r0,r1 237: 238: rol r2 /align to base 2**31 239: rol r1 240: rol r0 241: ror r2 /c-bit was clear here 242: 243: add r3,r1 244: adc r0 245: bpl 8f 246: add $100000,r0 247: adc r1 /can't carry 248: 249: 8: 250: add r2,r0 251: bpl 9f 252: add $100000,r0 253: adc r1 254: adc r0 /can't carry into bit 255: 9: /we should reduce 2**31-1 to 0 here 256: /but it can't happen 257: 258: // all done 259: // put the result at (r4) 260: // and restore r4 261: 262: mov r1,(r4) 263: mov r0,-(r4) 264: rts pc 265: 266: ///********************************************************************* 267: /// new added section. 268: // 3/6/80 j reeds. save and restore state of generators 269: // 270: // integer bigtab(269), smtab(12) 271: // 272: // call ransav(bigtab) 273: // call ranres(bigtab) ! save, retore complete state 274: // 275: // these two not 276: // implemented yet: 277: // call ranmup(smtab) 278: // call ranmdn(smtab) ! save and restore seeds, etc but not 279: // !mixing table 280: // 281: .text 282: blkmov: 283: dec r2 284: blt 2f /while (--r2 >= 0 BEGIN 285: cmp r3,$0 286: beq 1f 287: mov (r0)+,(r1)+ /if (r3 == 0) *r0++ = *r1++ 288: br blkmov 289: 1: 290: mov (r1)+,(r0)+ /else *r1++ = *r0++ 291: br blkmov 292: 2: 293: rts pc / END elihw 294: .globl _ransav 295: _ransav: 296: _ransav_: 297: 298: jsr r5,csv 299: mov $1,r3 /r3 = 1 so params = r0 -> r1 = table 300: br ranrs2 301: .globl _ranres 302: _ranres_: 303: 304: _ranres: 305: jsr r5,csv 306: clr r3 /r3 = 0 so table = r1 -> r0 = params 307: ranrs2: 308: / cmp (r5)+,$1 / arg count == 1? if not, call is no-op 309: / beq 1b 310: / rts pc 311: /1: 312: mov 4(r5),r1 / r1 is start of user array 313: mov $statev,r0 / r0 is rng's statevector 314: mov $269.,r2 / which is 269 words long 315: 316: jsr pc,blkmov 317: 318: jmp cret 319: 320: //.globl ranmup 321: //ranmup: mov $1,r3 322: // br fred 323: //.globl ranmdn 324: //ranmdn: clr r3 325: //fred: 326: // tst (r5)+ 327: // mov (r5),r1 328: // mov $statev,r0 329: // mov $12.,r2 330: // 331: // jsr pc,blkmov 332: // rts pc 333: / .end 334: rts pc 335: 336: time = 13. 337: .text 338: rndmze: 339: _randomi: 340: _rndmze_: 341: sys time 342: mov r0,-(sp) 343: mov r1,-(sp) 344: xor r0,r1 345: mov r1,-(sp) 346: ror r0 347: rol r1 348: mov r0,-(sp) 349: mov r1,-(sp) 350: xor r0,r1 351: mov r1,-(sp) 352: 353: jsr pc,_iran 354: 355: add $12.,sp 356: rts pc
|
Last modified: 1981-05-01
Generated: 2016-12-26 |
Generated by src2html V0.67
page hit count: 554 |
|