limit number of random

[Karim BELABAS <Karim.Belabas@math.u-psud.fr>]

[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()
+ {
+   return GLUE(((mymyrand()>>12)&LOWMASK), ((mymyrand()>>12)&LOWMASK));
+ }
+ 
  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
--
PARI/GP Home Page: http://pari.home.ml.org