Karim BELABAS on Thu, 1 Oct 1998 18:34:13 +0200 (MET DST)

[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

 limit number of random

```[Kiyoshi Ohgishi writes:]
> I want get large random number.
>
> But version 2.0.4.alpha,
>
> ? random(2^100)
>   ***   invalid bound in random.
>
> I think "random" work at version 2.0.alpha

random() never produced anything bigger than 2^31-1. But I agree this could
be useful. The following patch does what you want (and improves a bit the
previous version even if the bound is < 2^31). Apply it to 2.0.11, or wait
for 2.0.12 to be released.

Cheers,

Karim.

*** src/basemath/bibli2.c.orig	Mon Sep 28 14:29:45 1998
--- src/basemath/bibli2.c	Thu Oct  1 18:27:47 1998
***************
*** 652,680 ****
long
mymyrand()
{
! #ifdef LONG_IS_64BIT
!   if (BITS_IN_RANDOM==64)
!     pari_randseed = 1000000000000654397*pari_randseed + 12347;
!   else
!     pari_randseed = (1000276549*pari_randseed + 12347) % 2147483648;
#else
!   pari_randseed = 1000276549*pari_randseed + 12347;
!   if (pari_randseed<0) pari_randseed &= (~HIGHBIT);
#endif
return pari_randseed;
}

GEN
genrand(GEN N)
{
!   long av,tetpil;

if (!N) return stoi(mymyrand());
!   if (typ(N)!=t_INT || signe(N)<=0 || is_bigint(N))
!     err(talker,"invalid bound in random");

!   av = avma; N = mulis(N, mymyrand()); tetpil = avma;
!   return gerepile(av,tetpil, shifti(N, -(BITS_IN_RANDOM-1)));
}

GEN
--- 652,699 ----
long
mymyrand()
{
! #if BITS_IN_RANDOM == 64
!   pari_randseed = (1000000000000654397*pari_randseed + 12347) & ~HIGHBIT;
#else
!   pari_randseed = (1000276549*pari_randseed + 12347) & 0x7fffffff;
#endif
return pari_randseed;
}

+ GEN muluu(ulong x, ulong y);
+
+ #define GLUE(hi, lo) (((hi) << BITS_IN_HALFULONG) | (lo))
+ static ulong
+ gp_rand()
+ {
+ }
+
GEN
genrand(GEN N)
{
!   long lx,i;
!   GEN x;

if (!N) return stoi(mymyrand());
!   if (typ(N)!=t_INT || signe(N)<=0) err(talker,"invalid bound in random");

!   lx = lgefint(N); x = new_chunk(lx);
!   for (i=2; i<lx; i++)
!   {
!     long av = avma, n = N[i];
!     ulong r = gp_rand();
!     if (i < lx-1) n++; else if (!n) r = 0;
!     if (n) { GEN p1 = muluu(n,r); r = (lgefint(p1)<=3)? 0: p1[2]; }
!     avma = av; x[i] = r;
!     if (r < (ulong)N[i]) break;
!   }
!   for (i++; i<lx; i++) x[i] = gp_rand();
!   i=2; while (i<lx && !x[i]) i++;
!   i -= 2; x += i; lx -= i;
!   x[1] = evalsigne(lx>2) | evallgefint(lx);
!   x[0] = evaltyp(t_INT) | evallg(lx);
!   avma = (long)x; return x;
}

GEN
*** src/kernel/none/mp.c.orig	Wed Sep 30 19:03:00 1998
--- src/kernel/none/mp.c	Thu Oct  1 17:45:41 1998
***************
*** 648,653 ****
--- 648,671 ----
z[2]=p1; return z;
}

+ GEN
+ muluu(ulong x, ulong y)
+ {
+   long p1;
+   GEN z;
+   LOCAL_HIREMAINDER;
+
+   if (!x || !y) return gzero;
+   p1 = mulll(x,y);
+   if (hiremainder)
+   {
+     z=cgeti(4); z[1] = evalsigne(1) | evallgefint(4);
+     z[2]=hiremainder; z[3]=p1; return z;
+   }
+   z=cgeti(3); z[1] = evalsigne(1) | evallgefint(3);
+   z[2]=p1; return z;
+ }
+
/* assume ny > 0 */
#ifdef INLINE
INLINE
--
Karim Belabas                    email: Karim.Belabas@math.u-psud.fr
Dep. de Mathematiques, Bat. 425
Universite Paris-Sud             Tel: (00 33) 1 69 15 57 48
F-91405 Orsay (France)           Fax: (00 33) 1 69 15 60 19
--