مولد الأرقام العشوائية ليمر

مولد الأرقام العشوائية ليمر ^{[ 1 ]} (نسبةً إلى دي إتش ليمر )، والذي يُشار إليه أحيانًا باسم مولد الأرقام العشوائية بارك-ميلر (نسبةً إلى ستيفن ك.  بارك وكيث دبليو. ميلر )، هو نوع من مولدات التوافق الخطي (LCG) التي تعمل في مجموعة ضربية من الأعداد الصحيحة بتردد n . الصيغة العامة هي

${\displaystyle X_{k+1}=a\cdot X_{k}{\bmod {m}},}$

حيث يكون المعامل m عددًا أوليًا أو قوة لعدد أولي ، ويكون المضاعف a عنصرًا من رتبة ضربية عالية modulo m (على سبيل المثال، جذر بدائي modulo n )، ويكون البذرة X ₀ عددًا أوليًا نسبيًا مع m .

ومن الأسماء الأخرى مولد التوافق الخطي المضاعف (MLCG) ^{[ 2 ]} ومولد التوافق المضاعف (MCG) .

في عام 1988، اقترح بارك وميلر ^{[ 3 ]} مولد أرقام عشوائية من نوع ليمر بمعاملات محددة: m = 2^ ³¹ - 1 = 2,147,483,647 ( عدد أولي ميرسين M^ ₃₁ ) و a = 7 ^{^5} = 16,807 (جذر أولي modulo M^ ₃₁ )، والمعروف الآن باسم MINSTD . على الرغم من أن MINSTD تعرض لانتقادات لاحقة من قبل مارساجليا وسوليفان (1993) ^{[ 4 ] [ 5 ]،} إلا أنه لا يزال مستخدمًا حتى اليوم (خاصةً في مكتبة CarbonLib ولغة C++11 ) minstd_rand0. رد بارك وميلر وستوكمير على هذه الانتقادات (1993) ^{[ 6 ]} قائلين:

نظراً لطبيعة هذا المجال الديناميكية، يصعب على غير المتخصصين اتخاذ قرارات بشأن المولد المناسب. "أعطني شيئاً أستطيع فهمه وتطبيقه ونقله... ليس بالضرورة أن يكون متطوراً للغاية، يكفي أن يكون جيداً وفعالاً بشكل معقول." كانت مقالتنا ومولد المعايير الأدنى المرتبط بها محاولةً للاستجابة لهذا الطلب. بعد خمس سنوات، لا نرى حاجة لتغيير ردنا سوى اقتراح استخدام المضاعف a  =  48271 بدلاً من 16807.

يتم استخدام هذا الثابت المعدل في مولد الأرقام العشوائية الخاص بـ C++11 .minstd_rand

يستخدم جهاز Sinclair ZX81 والإصدارات اللاحقة منه مولد الأرقام العشوائية Lehmer RNG بمعاملات m  =  2^ ¹⁶  +  1 = 65537 ( عدد أولي من نوع Fermat F^ ₄ ) و a  =  75 (جذر أولي modulo F _^4 ). ^{[ 7 ] [ 8 ]} مولد الأرقام العشوائية RANF من CRAY هو مولد أرقام عشوائية Lehmer RNG بمعامل m = 2^ ⁴⁸ و a = 44485709377909. ^{[ 9 ]} تتضمن مكتبة GNU العلمية العديد من مولدات الأرقام العشوائية من نوع Lehmer، بما في ذلك MINSTD وRANF ومولد الأرقام العشوائية RANDU سيئ السمعة من IBM . ^{[ 9 ]}    

في أغلب الأحيان، يُختار المعامل كعدد أولي، مما يجعل اختيار بذرة أولية فيما بينها أمرًا بسيطًا (أي عدد بين 0 و ₀ ). ينتج عن ذلك أفضل جودة للمخرجات، ولكنه يُضيف بعض التعقيد إلى التنفيذ، ومن غير المرجح أن يتطابق نطاق المخرجات مع التطبيق المطلوب؛ إذ يتطلب التحويل إلى النطاق المطلوب عملية ضرب إضافية.

استخدام معامل m، وهو قوة للعدد اثنين، يُسهّل عملية التنفيذ الحاسوبي بشكل خاص، ولكنه يأتي بتكلفة: فالدورة لا تتجاوز m /4، ودورات البتات الأدنى أقصر من ذلك. والسبب في ذلك هو أن البتات k الأدنى تُشكّل مولدًا بمعامل 2k ^{بمفردها ؛ فالبتات الأعلى لا تؤثر أبدًا} _على البتات الأدنى. ^{[ 10 ]} القيم Xi دائمًا فردية (البت 0 لا يتغير أبدًا)، ويتناوب البتّان 2 و1 (تتكرر البتات الثلاثة الأدنى بدورة 2)، وتتكرر البتات الأربعة الأدنى بدورة 4، وهكذا. لذلك، يجب على التطبيق الذي يستخدم هذه الأرقام العشوائية استخدام البتات الأكثر أهمية؛ إذ أن تقليص النطاق باستخدام عملية بمعامل زوجي سيؤدي إلى نتائج كارثية. ^{[ 11 ]}

لتحقيق هذه الفترة، يجب أن يحقق المضاعف a  ≡  ±3 (mod  8)، ^{[ 12 ]} ويجب أن تكون البذرة X ₀ فردية.

يُمكن استخدام مُعامل مُركّب، ولكن يجب تهيئة المُولّد بقيمة أولية نسبية مع m ، وإلا سيتقلص طول الدورة بشكل كبير. على سبيل المثال، قد يبدو مُعامل F₅ = ^2³²  +  1 جذابًا، حيث يُمكن بسهولة ربط المُخرجات بكلمة 32 بت 0 ≤ Xᵢ _- 1 < 2³² . مع ذلك، _فإن تهيئة قيمة أولية X₀ = _6700417 (التي تقسم ^2³² + 1) أو أي مُضاعف  لها ستؤدي إلى مُخرج بدورة 640 فقط ^.     

مولد آخر ذو معامل مركب هو الذي أوصى به ناكازاوا وناكازاوا: ^{[ 13 ]}

• م =134 265 023 ×134 475 827 =18 055 400 005 099 021 ≈ 2 ⁵⁴
• أ =7 759 097 958 782 935 (أي من ± a ^±1 (mod m ) سيفي بالغرض أيضًا)

بما أن كلا عاملي المعامل أقل من 2^ ³² ، فمن الممكن الحفاظ على الحالة بتردد كل عامل، وبناء قيمة الخرج باستخدام نظرية الباقي الصينية ، باستخدام عمليات حسابية وسيطة لا تتجاوز 64 بت . ^{[ 13 ] : 70}

يُعدّ مولد التوافق الخطي المُدمج أحد أكثر التطبيقات شيوعًا للفترات الطويلة ؛ إذ يُكافئ دمج عدة مولدات (مثلاً بجمع مخرجاتها) مخرج مولد واحد يكون معامله هو حاصل ضرب معاملات المولدات المكونة له. ^{[ 14 ]} وتكون دورته هي المضاعف المشترك الأصغر لدورات المولدات المكونة. على الرغم من أن الدورات تشترك في قاسم مشترك هو 2، إلا أنه يمكن اختيار المعاملات بحيث يكون هو القاسم المشترك الوحيد، وتكون الدورة الناتجة هي ( _m1 - ₁ ⁾ ( m2  - 1)···( mk _- 1 )/2k ^- 1. ^{[ 2 ] : 744.} ومن الأمثلة على ذلك مولد ويشمان-هيل .     

على الرغم من إمكانية اعتبار مولد الأرقام العشوائية ليمر حالة خاصة من مولد التوافق الخطي عندما يكون c = 0 _، إلا أنها حالة خاصة تستلزم قيودًا وخصائص معينة. فعلى وجه الخصوص، في مولد الأرقام العشوائية ليمر، يجب أن يكون البذرة الأولية X₀ أوليًا نسبيًا مع المعامل m ، وهو شرط غير مطلوب في مولدات التوافق الخطي عمومًا. كما أن اختيار المعامل m والمضاعف a أكثر تقييدًا في مولد الأرقام العشوائية ليمر. وعلى عكس مولد التوافق الخطي، فإن أقصى دورة لمولد الأرقام العشوائية ليمر تساوي m - 1، وتكون كذلك عندما يكون m عددًا أوليًا ويكون a جذرًا أوليًا بتردد m .  

من ناحية أخرى، فإن اللوغاريتمات المنفصلة (للأساس a أو أي جذر أولي modulo m ) لـ X _k في${\displaystyle \mathbb {Z} _{m}}$تمثل تسلسلًا متطابقًا خطيًا modulo totient أويلر${\displaystyle \varphi (m)}$.

يتطلب حساب المعامل الأولي حساب حاصل ضرب مزدوج العرض وخطوة اختزال صريحة. إذا استُخدم معامل أقل بقليل من قوة العدد 2 ( الأعداد الأولية لمرسين ^2³¹ -  1  و2⁶¹ - 1 شائعة، وكذلك 2³² ^- 5 و2⁶⁴ ^- 59 ) ^، فإن الاختزال بمعامل m = 2e ^- d يمكن تنفيذه بتكلفة أقل من القسمة العامة مزدوجة العرض باستخدام المتطابقة 2e ≡ d (mod m ⁾ .      

تقسم خطوة الاختزال الأساسية الناتج إلى جزأين، كل منهما مكون من e بت، ثم تضرب الجزء الأعلى في d ، وتجمعهما: ( ax mod 2e ⁾ + d ⌊ ax / ^2e ⌋ . يمكن بعد ذلك طرح m حتى يصبح الناتج ضمن النطاق المحدد. يقتصر عدد عمليات الطرح على ad / m ، ويمكن تقليصه بسهولة إلى عملية واحدة إذا كانت d صغيرة وتم اختيار a < m / d . (يضمن هذا الشرط أيضًا أن يكون d ⌊ ax /2e ^⌋ ناتج ضرب أحادي العرض؛ إذا لم يتحقق هذا الشرط، يجب حساب ناتج ضرب ثنائي العرض).

عندما يكون المعامل عددًا أوليًا من أعداد ميرسين ( d  =  1)، تكون العملية بسيطة للغاية. فعملية الضرب في d ليست بديهية فحسب، بل يمكن استبدال الطرح المشروط بإزاحة وجمع غير مشروطين. ولتوضيح ذلك، لاحظ أن الخوارزمية تضمن أن x ≢ 0 (mod m ) ، مما يعني أن x  =  0 و x  = m كلاهما مستحيل. وهذا يُغني عن الحاجة إلى النظر في تمثيلات e -bit المكافئة للحالة؛ إذ لا تحتاج إلى اختزال إلا القيم التي تكون فيها البتات العليا غير صفرية. 

لا يمكن أن تمثل البتات المنخفضة (e) في حاصل ضرب ax قيمة أكبر من m ، ولن تحمل البتات العليا قيمة أكبر من a  −  1  ≤  m  −  2. وبالتالي، تُنتج خطوة الاختزال الأولى قيمة لا تتجاوز m  + a − 1 ≤ 2 m − 2 = 2 ^{e + 1} − 4. هذا عدد مكون من ( e + 1) بت، والذي قد يكون أكبر من m (أي قد تكون البتة e فيه مُفعّلة)، لكن النصف الأعلى لا يتجاوز 1، وإذا كان كذلك، فستكون البتات المنخفضة (e) أقل من m . لذا، سواء كانت البتة العليا 1 أو 0، فإن خطوة الاختزال الثانية (جمع النصفين) لن تتجاوز عدد البتات e ، وسيكون المجموع هو القيمة المطلوبة.         

إذا كانت قيمة d أكبر من 1، يمكن تجنب الطرح الشرطي، لكن العملية تصبح أكثر تعقيدًا. يكمن التحدي الأساسي في حساب معامل مثل 2^ ³² - 5 في ضمان إنتاج تمثيل واحد فقط لقيم مثل 1 ≡ 2^ ³² - 4. الحل هو إضافة d مؤقتًا ، بحيث يكون نطاق القيم الممكنة من d إلى 2 ^{^e} - 1، ثم تقليل القيم الأكبر من e بت بطريقة لا تُنتج تمثيلات أقل من d . أخيرًا، ينتج عن طرح الإزاحة المؤقتة القيمة المطلوبة.          

لنفترض مبدئيًا أن لدينا قيمة y مختزلة جزئيًا ومحدودة بحيث 0  ≤ y < 2m = 2e + ^{1 −} 2d . في هذه الحالة، ستنتج خطوة طرح واحدة إزاحة 0 ≤ y ′ = (( y + d ) mod 2e ⁾ + d⌊ ( y + d )/ ^2e⌋ − d < m . ولتوضيح ذلك، لننظر في حالتين :            

0 ≤ y < m = 2 ^e − d

في هذه الحالة، y + d < 2 ^e و y ′  = y < m ، كما هو مطلوب.   

م ≤ ص < ٢ م

في هذه الحالة، 2e ^≤ y + d < 2e + ^{1 هو} عدد مكون من ( e + 1) بت، و ⌊ ( y + d )/ ^2e⌋ = 1. بالتالي، y ′ = ( y + d ) − 2e + d ⁻ d = y − 2e + ^d = y − m < m ، كما هو مطلوب. ولأن الجزء العلوي المضروب هو d ، فإن المجموع يساوي d على الأقل ، وطرح الإزاحة لا يُسبب أبدًا حدوث تدفق سفلي.                              

(في حالة مولد ليمر تحديدًا، لن تحدث حالة الصفر أو صورتها y  = m أبدًا، لذا فإن إزاحة d − 1 ستعمل بنفس الطريقة، إذا كان ذلك أكثر ملاءمة. هذا يقلل الإزاحة إلى 0 في حالة عدد ميرسين الأولي، عندما d = 1.)     

يمكن تقليل ناتج أكبر ax إلى أقل من 2 m  =  2 ^{e +1}  −  2 d عن طريق خطوة اختزال واحدة أو أكثر بدون إزاحة.

إذا كان ad  ≤ m ، فإن خطوة اختزال إضافية واحدة تكفي. بما أن x < m ، فإن ax < am ≤ ( a − 1)2e ^{، وخطوة اختزال واحدة تحول هذا إلى 2e} − 1 + ( a − 1) d = m + ad − 1 على الأكثر. وهذا ضمن حد 2m إذا كان ad − 1 < m ، وهو الافتراض الأولي.                        

إذا كان ad > m ، فمن الممكن أن تُنتج خطوة الاختزال الأولى مجموعًا أكبر من 2m = 2e + ^{1 −} 2d ، وهو كبير جدًا بالنسبة لخطوة الاختزال النهائية. (كما يتطلب الأمر الضرب في d لإنتاج ناتج أكبر من e بت، كما ذُكر أعلاه). ومع ذلك ، طالما ^أن d² < 2e ^،        سيؤدي التخفيض الأول إلى قيمة في النطاق المطلوب لتطبيق خطوتين تخفيض في الحالة السابقة.

إذا لم يكن الضرب ذو العرض المزدوج متاحًا، فيمكن استخدام طريقة شراج ، ^{[ 15 ] [ 16 ]} والتي تسمى أيضًا طريقة التحليل التقريبي، ^{[ 17 ]} لحساب ax mod m ، ولكن هذا يأتي على حساب:

• يجب أن يكون المعامل قابلاً للتمثيل في عدد صحيح موقّع ؛ يجب أن تسمح العمليات الحسابية بنطاق ± m .
• إن اختيار المضاعف a محدود. نشترط أن يكون m mod a ≤ ⌊ m / a ⌋ ، ويتحقق ذلك عادةً باختيار a  ≤ √ m . 
• يلزم إجراء عملية قسمة واحدة (مع الباقي) لكل تكرار.

على الرغم من شيوع هذه التقنية في التطبيقات المحمولة بلغات البرمجة عالية المستوى التي تفتقر إلى عمليات الضرب المضاعف، ^{[ 2 ] : 744، إلا} أنه في الحواسيب الحديثة، يُنفذ القسمة على ثابت عادةً باستخدام الضرب المضاعف، لذا يُنصح بتجنب هذه التقنية إذا كانت الكفاءة مهمة. حتى في لغات البرمجة عالية المستوى، إذا كان المضاعف a محدودًا بـ √m ، فيمكن حساب حاصل الضرب المضاعف ax باستخدام عمليتي ضرب أحاديتي العرض، ثم اختزاله باستخدام التقنيات المذكورة أعلاه.

لاستخدام طريقة شراج، قم أولاً بتحليل المعادلة m = qa + r ، أي احسب مسبقًا الثوابت المساعدة r = m mod a و q = ⌊ m / a ⌋ = ( m − r )/ a . ثم، في كل تكرار، احسب ax ≡ a ( x mod q ) − r ⌊ x / q ⌋ (mod m ) .

هذه المساواة قائمة لأن

${\displaystyle {\begin{aligned}aq&=a\cdot (mr)/a\\&=(a/a)\cdot (mr)\\&=(mr)\\&\equiv -r{\pmod {m}}\\\end{aligned}}}$

إذا قمنا بتحليل x = ( x mod q ) + q ⌊ x / q ⌋ ، فسنحصل على:

${\displaystyle {\begin{align}ax&=a(x{\bmod {q}})+aq\lfloor x/q\rfloor \\&=a(x{\bmod {q}})+(mr)\lfloor x/q\rfloor \\&\equiv a(x{\bmod {q}})-r\lfloor x/q\rfloor {\pmod {م}}\\\النهاية{محاذاة}}}$

السبب في عدم حدوث تجاوز هو أن كلا الحدين أصغر من m . بما أن x  mod q < q ≤ m / a ، فإن الحد الأول أصغر تمامًا من am / a = m ويمكن حسابه باستخدام ضرب أحادي العرض.   

إذا تم اختيار قيمة a بحيث يكون r  ≤ q (وبالتالي r / q ≤ 1)، فإن الحد الثاني يكون أيضًا أقل من m : r ⌊ x / q ⌋ ≤ rx / q = x ( r / q ) ≤ x (1) = x < m . وبالتالي، يقع الفرق في النطاق [1− m , m −1] ويمكن اختزاله إلى [0, m −1] بجمع شرطي واحد. ^{[ 18 ]}     

يمكن توسيع هذه التقنية للسماح بقيمة r سالبة (− q  ≤ r < 0)، مما يؤدي إلى تغيير الاختزال النهائي إلى طرح مشروط.   

يمكن توسيع هذه التقنية للسماح بقيم أكبر لـ a بتطبيقها بشكل متكرر. ^{[ 17 ] : 102} من بين الحدين المطروحين للحصول على النتيجة النهائية، فإن الحد الثاني فقط ( r ⌊ x / q ⌋ ) معرض لخطر تجاوز السعة. ولكن هذا في حد ذاته ضرب معياري بثابت وقت الترجمة r ، ويمكن تنفيذه بنفس التقنية. ولأن كل خطوة، في المتوسط، تقسم حجم المضاعف إلى النصف (0 ≤ r < a ، القيمة المتوسطة ( a − 1) / 2)، يبدو أن هذا يتطلب خطوة واحدة لكل بت، وهو أمر غير فعال بشكل ملحوظ. ومع ذلك، تقسم كل خطوة أيضًا x على خارج قسمة متزايد باستمرار q = ⌊ m / a ⌋ ، وسرعان ما يتم الوصول إلى نقطة يكون فيها الوسيط 0، ويمكن إنهاء الاستدعاء المتكرر.    

باستخدام لغة البرمجة C ، يمكن كتابة مولد الأرقام العشوائية Park-Miller على النحو التالي:

uint32_t lcg_parkmiller ( uint32_t * state ) { return * state = ( uint64_t ) * state * 48271 % 0x7fffffff ; }

يمكن استدعاء هذه الدالة بشكل متكرر لتوليد أرقام شبه عشوائية، شريطة أن يحرص المستدعي على تهيئة الحالة بأي عدد أكبر من الصفر وأصغر من المعامل. في هذا التطبيق، يلزم استخدام حسابات 64 بت؛ وإلا فقد يحدث تجاوز في حاصل ضرب عددين صحيحين 32 بت.

لتجنب القسمة على 64 بت، قم بإجراء عملية الاختزال يدويًا:

uint32_t lcg_parkmiller ( uint32_t * state ) { uint64_t product = ( uint64_t ) * state * 48271 ; uint32_t x = ( product & 0x7fffffff ) + ( product >> 31 );x = ( x & 0x7fffffff ) + ( x >> 31 ); return * state = x ; }

لاستخدام العمليات الحسابية ذات 32 بت فقط، استخدم طريقة شراج:

uint32_t lcg_parkmiller ( uint32_t * state ) { // معلمات مُحسوبة مسبقًا لطريقة شراج const uint32_t M = 0x7fffffff ; const uint32_t A = 48271 ; const uint32_t Q = M / A ; // 44488 const uint32_t R = M % A ; // 3399uint32_t div = * state / Q ; // الحد الأقصى: M / Q = A = 48,271 uint32_t rem = * state % Q ; // الحد الأقصى: Q - 1 = 44,487int32_t s = rem * A ; // الحد الأقصى: 44,487 * 48,271 = 2,147,431,977 = 0x7fff3629 int32_t t = div * R ; // الحد الأقصى: 48,271 * 3,399 = 164,073,129 int32_t result = s - t ;إذا كانت ( النتيجة < 0 ) النتيجة += M ؛return * state = result ; }

أو استخدم عمليتي ضرب 16×16 بت:

uint32_t lcg_parkmiller ( uint32_t * state ) { const uint32_t A = 48271 ;uint32_t low = ( * state & 0x7fff ) * A ; // الحد الأقصى: 32,767 * 48,271 = 1,581,695,857 = 0x5e46c371 uint32_t high = ( * state >> 15 ) * A ; // الحد الأقصى: 65,535 * 48,271 = 3,163,439,985 = 0xbc8e4371 uint32_t x = low + (( high & 0xffff ) << 15 ) + ( high >> 16 ); // الحد الأقصى: 0x5e46c371 + 0x7fff8000 + 0xbc8e = 0xde46ffffx = ( x & 0x7fffffff ) + ( x >> 31 ); return * state = x ; }

يستخدم مولد ليمر شائع آخر المعامل الأولي 2 ³² −5:

uint32_t lcg_rand ( uint32_t * state ) { return * state = ( uint64_t ) * state * 279470273u % 0xfffffffb ; }

ويمكن كتابة هذا أيضًا بدون قسمة 64 بت:

uint32_t lcg_rand ( uint32_t * state ) { uint64_t product = ( uint64_t ) * state * 279470273u ; uint32_t x ;// غير مطلوب لأن 5 * 279470273 = 0x5349e3c5 يناسب 32 بت. // الناتج = (الناتج & 0xffffffff) + 5 * (الناتج >> 32)؛ // المضاعف الأكبر من 0x33333333 = 858,993,459 سيحتاج إليه.// نتيجة الضرب تتسع في 32 بت، لكن المجموع قد يكون 33 بت. المنتج = ( المنتج & 0xffffffff ) + 5 * ( uint32_t )( المنتج >> 32 );product += 4 ; // هذا المجموع مضمون أن يكون 32 بت. x = ( uint32_t ) product + 5 * ( uint32_t )( product >> 32 ); return * state = x - 4 ; }

تتمتع العديد من مولدات ليمر الأخرى بخصائص جيدة. يتطلب مولد ليمر التالي ذو المقياس 2^ ^{128 دعمًا من} المترجم لـ 128 بت ، ويستخدم مضاعفًا محسوبًا بواسطة ليكويير. ^{[ 19 ]} دورته 2^ ¹²⁶ .

static unsigned __int128 state ;/* يجب تهيئة الحالة بقيمة فردية. */ void seed ( unsigned __int128 seed ) { state = seed << 1 | 1 ; }uint64_t next ( void ) { // لا يمكن لـ GCC كتابة قيم حرفية 128 بت، لذلك نستخدم تعبيرًا const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; state *= mult ; return state >> 64 ; }

يقوم المولد بحساب قيمة فردية مكونة من 128 بت ويعيد أعلى 64 بت منها.

يجتاز هذا المولد اختبار BigCrush من TestU01 ، لكنه يفشل في اختبار TMFn من PractRand . صُمم هذا الاختبار خصيصًا لاكتشاف عيب هذا النوع من المولدات: بما أن المعامل هو قوة للعدد 2، فإن دورة البت الأدنى في الناتج هي 2⁶² فقط ^، وليست ^2¹²⁶ . وتتشابه المولدات الخطية التوافقية ذات المعامل من قوة العدد 2 في سلوكها.

تعمل الروتينية الأساسية التالية على تحسين سرعة الكود أعلاه لأحمال العمل الصحيحة (إذا سمح المترجم بتحسين تعريف الثابت خارج حلقة الحساب):

uint64_t next ( void ) { uint64_t result = state >> 64 ; // لا يمكن لـ GCC كتابة قيم حرفية 128 بت، لذلك نستخدم تعبيرًا const unsigned __int128 mult = ( unsigned __int128 ) 0x12e15e35b500f16e << 64 | 0x2e714eb2b37916a5 ; state *= mult ; return result ; }

ومع ذلك، ولأن عملية الضرب مؤجلة، فإنها غير مناسبة للتجزئة، حيث أن الاستدعاء الأول يعيد ببساطة أعلى 64 بت من حالة البذرة.

• ليمر، د. هـ. (1949). "الأساليب الرياضية في وحدات الحوسبة واسعة النطاق". وقائع الندوة الثانية حول آلات الحوسبة الرقمية واسعة النطاق . الصفحات 141-146 . MR 0044899 .  (نسخة المجلة: حوليات مختبر الحوسبة بجامعة هارفارد ، المجلد 26 (1951)).
• ستيف بارك، مولدات الأرقام العشوائية

• قد تكون الأعداد الأولية التي تقل قليلاً عن قوة العدد اثنين مفيدة لاختيار المعاملات. جزء من صفحات الأعداد الأولية .