Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - random.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 157 157 100.0 %
Date: 2022-07-03 07:33:15 Functions: 17 17 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : /********************************************************************/
      15             : /*                                                                  */
      16             : /*                      PSEUDO-RANDOM INTEGERS                      */
      17             : /*                                                                  */
      18             : /********************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : /********************************************************************/
      22             : /*                    XORGEN (Richard P. Brent)                     */
      23             : /*          http://wwwmaths.anu.edu.au/~brent/random.html           */
      24             : /*        (initial adaptation to PARI/GP by Randall Rathbun)        */
      25             : /********************************************************************/
      26             : /* Adapted from xorgens.c version 3.04, Richard P. Brent, 20060628 (GPL).
      27             :  * 32-bit or 64-bit integer random number generator with period at
      28             :  * least 2**4096-1. It is assumed that "ulong" is a 32-bit or 64-bit integer */
      29             : 
      30             : #ifdef LONG_IS_64BIT
      31             :   typedef ulong u64;
      32             : #else
      33             :   typedef unsigned long long u64;
      34             : static u64
      35      448306 : _32to64(ulong a, ulong b) { u64 v = a; return (v<<32)|b; }
      36             : static void
      37     7317752 : _64to32(u64 v, ulong *a, ulong *b) { *a = v>>32; *b = v&0xFFFFFFFF; }
      38             : #endif
      39             : static THREAD u64 state[64];
      40             : static THREAD u64 xorgen_w;
      41             : static THREAD int xorgen_i;
      42             : /* weyl = odd approximation to 2^64*(sqrt(5)-1)/2. */
      43             : static const u64 weyl = (((u64)0x61c88646U)<<32)|((u64)0x80b583ebU);
      44             : 
      45             : static u64
      46   122268912 : block(void)
      47             : {
      48   122268912 :   const int r = 64;
      49   122268912 :   const int a = 33, b = 26, c = 27, d = 29, s = 53;
      50             :   u64 t, v, w;
      51   122268912 :   xorgen_i = (xorgen_i+1)&(r-1);
      52   122268912 :   t = state[xorgen_i];
      53   122268912 :   v = state[(xorgen_i+(r-s))&(r-1)];   /* index is (i-s) mod r */
      54   122268912 :   t ^= t<<a; t ^= t>>b;                   /* (I + L^a)(I + R^b) */
      55   122268912 :   v ^= v<<c; v ^= v>>d;                   /* (I + L^c)(I + R^d) */
      56   122268912 :   w = t^v;
      57   122268912 :   return state[xorgen_i] = w;               /* update circular array */
      58             : }
      59             : 
      60             : /* v > 0 */
      61             : static void
      62      284670 : init_xor4096i(u64 v)
      63             : {
      64      284670 :   const int r = 64;
      65             :   int k;
      66             : 
      67    17814736 :   for (k = r; k > 0; k--) {/* avoid correlations for close seeds */
      68    17530066 :     v ^= v<<10; v ^= v>>15; /* recurrence has period 2**64-1 */
      69    17530066 :     v ^= v<<4;  v ^= v>>13;
      70             :   }
      71    17620118 :   for (xorgen_w = v, k = 0; k < r; k++) { /* initialise circular array */
      72    17335448 :     v ^= v<<10; v ^= v>>15;
      73    17335448 :     v ^= v<<4;  v ^= v>>13;
      74    17335448 :     state[k] = v + (xorgen_w+=weyl);
      75             :   }
      76             :   /* discard first 4*r results */
      77    55826402 :   for (xorgen_i = r-1, k = 4*r; k > 0; k--) (void)block();
      78      278396 : }
      79             : 
      80      230288 : void pari_init_rand(void) { init_xor4096i(1UL); }
      81             : 
      82             : static u64
      83    66706850 : rand64(void)
      84             : {
      85    66706850 :   u64 v = block();
      86    66709756 :   xorgen_w += weyl; /* update Weyl generator */
      87    66709756 :   return v + (xorgen_w ^ (xorgen_w>>27));
      88             : }
      89             : 
      90             : /* One random number uniformly distributed in [0..2**BIL) is returned, where
      91             :  * BIL = 8*sizeof(ulong) = 32 or 64. */
      92             : ulong
      93       78717 : pari_rand(void) { return rand64(); }
      94             : 
      95             : void
      96      102844 : setrand(GEN x)
      97             : {
      98      102844 :   const int r2 = numberof(state);
      99             :   long i, lx;
     100             :   u64 v;
     101             :   GEN xp;
     102      102844 :   if (typ(x)!=t_INT) pari_err_TYPE("setrand",x);
     103      102844 :   if (signe(x) <= 0) pari_err_DOMAIN("setrand","n", "<=", gen_0, x);
     104      102837 :   lx = lgefint(x);
     105      102837 :   if (lx == 3) { v = x[2]; init_xor4096i(v); return; }
     106             : #ifndef LONG_IS_64BIT
     107        6899 :   if (lx == 4)
     108             :   {
     109           1 :     v = _32to64(*int_W(x,1),*int_W(x,0));
     110           1 :     init_xor4096i(v); return;
     111             :   }
     112             : #endif
     113       48286 :   xp = int_LSW(x);
     114             : #ifdef LONG_IS_64BIT
     115       41388 :   if (lx != 2 + r2+2)
     116           6 :     pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
     117     2689830 :   for (i = 0; i < r2; i++, xp = int_nextW(xp)) state[i] = *xp;
     118       41382 :   xorgen_w = *xp; xp = int_nextW(xp);
     119             : #else
     120        6898 :   if (lx != 2 + 2*r2+3)
     121           1 :     pari_err_DOMAIN("setrand", "n", "!=", strtoGENstr("getrand()"), x);
     122      448305 :   for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
     123      441408 :     state[i] = _32to64(*int_nextW(xp), *xp);
     124        6897 :   xorgen_w = _32to64(*int_nextW(xp), *xp); xp = int_nextW(int_nextW(xp));
     125             : #endif
     126       48279 :   xorgen_i =  (*xp) & 63;
     127             : }
     128             : 
     129             : GEN
     130      772989 : getrand(void)
     131             : {
     132      772989 :   const int r2 = numberof(state);
     133             :   GEN x;
     134             :   ulong *xp;
     135             :   long i;
     136      772989 :   if (xorgen_i < 0) init_xor4096i(1UL);
     137             : 
     138             : #ifdef LONG_IS_64BIT
     139      662287 :   x = cgetipos(2+r2+2); xp = (ulong *) int_LSW(x);
     140    43048568 :   for (i = 0; i < r2; i++, xp = int_nextW(xp)) *xp = state[i];
     141      662287 :   *xp = xorgen_w; xp = int_nextW(xp);
     142             : #else
     143      110702 :   x = cgetipos(2+2*r2+3); xp = (ulong *) int_LSW(x);
     144     7195630 :   for (i = 0; i < r2; i++, xp = int_nextW(int_nextW(xp)))
     145     7084928 :     _64to32(state[i], int_nextW(xp), xp);
     146      110702 :   _64to32(xorgen_w, int_nextW(xp), xp); xp = int_nextW(int_nextW(xp));
     147             : #endif
     148      772989 :   *xp = xorgen_i? xorgen_i: 64; return x;
     149             : }
     150             : 
     151             : /* assume 0 <= k <= BITS_IN_LONG. Return uniform random 0 <= x < (1<<k) */
     152             : long
     153    17023111 : random_bits(long k) { return rand64() >> (64-k); }
     154             : 
     155             : /********************************************************************/
     156             : /*                                                                  */
     157             : /*                         GENERIC ROUTINES                         */
     158             : /*                                                                  */
     159             : /********************************************************************/
     160             : 
     161             : /* assume n > 0 */
     162             : ulong
     163    34294952 : random_Fl(ulong n)
     164             : {
     165             :   ulong d;
     166             :   int shift;
     167             : #ifdef LONG_IS_64BIT
     168    28808795 :   int SHIFT = 0;
     169             : #else
     170     5486157 :   int SHIFT = 32;
     171             : #endif
     172             : 
     173    34294952 :   if (n == 1) return 0;
     174             : 
     175    34024290 :   shift = bfffo(n); /* 2^(BIL-shift) > n >= 2^(BIL-shift-1)*/
     176             :   /* if N a power of 2, increment shift. No reject */
     177    34024290 :   if ((n << shift) == HIGHBIT) return rand64() >> (SHIFT+shift+1);
     178             :   for (;;) {
     179    47071304 :     d = rand64() >> (SHIFT+shift); /* d < 2^(64-shift) uniformly distributed */
     180             :     /* reject strategy: proba success = n 2^(shift-64), in [1/2, 1[ */
     181    47072158 :     if (d < n) return d;
     182             :   }
     183             : }
     184             : 
     185             : /* assume N > 0, see random_Fl() for algorithm. Make sure that 32-bit and
     186             :  * 64-bit architectures produce the same integers (consuming random bits
     187             :  * by packets of 64) */
     188             : GEN
     189     5962275 : randomi(GEN N)
     190             : {
     191     5962275 :   long lx = lgefint(N);
     192             :   GEN x, d;
     193             :   int shift;
     194             : 
     195     5962275 :   if (lx == 3) return utoi( random_Fl(N[2]) );
     196             : 
     197      161076 :   shift = bfffo(*int_MSW(N));
     198             :   /* if N a power of 2, increment shift */
     199      161076 :   if (Z_ispow2(N) && ++shift == BITS_IN_LONG) { shift = 0; lx--; }
     200      161076 :   x = cgetipos(lx);
     201      117930 :   for (;;) {
     202      279006 :     GEN y, MSW = int_MSW(x), STOP = MSW;
     203             : #ifdef LONG_IS_64BIT
     204      679587 :     for (d = int_LSW(x); d != STOP; d = int_nextW(d)) *d = rand64();
     205      226099 :     *d = rand64() >> shift;
     206             : #else
     207       52907 :     if (!odd(lx)) STOP = int_precW(STOP);
     208             :     /* STOP points to where MSW would in 64-bit */
     209      122122 :     for (d = int_LSW(x); d != STOP; d = int_nextW(d))
     210             :     {
     211       69215 :       ulong a, b; _64to32(rand64(), &a,&b);
     212       69215 :       *d = b; d = int_nextW(d);
     213       69215 :       *d = a;
     214             :     }
     215             :     {
     216       52907 :       ulong a, b; _64to32(rand64() >> shift, &a,&b);
     217       52907 :       if (d == MSW) /* 32 bits needed */
     218       29312 :         *d = a;
     219             :       else
     220             :       { /* 64 bits needed */
     221       23595 :         *d = b; d = int_nextW(d);
     222       23595 :         *d = a;
     223             :       }
     224             :     }
     225             : #endif
     226      279006 :     y = int_normalize(x, 0);
     227      279006 :     if (abscmpii(y, N) < 0) return y;
     228             :   }
     229             : }
     230             : 
     231             : GEN
     232      544855 : random_F2x(long d, long vs)
     233             : {
     234      544855 :   ulong db, dl = dvmduBIL(d,&db);
     235      544855 :   long i, l = 2 + dl + !!db;
     236      544855 :   GEN y = cgetg(l,t_VECSMALL); y[1] = vs;
     237             : #ifdef LONG_IS_64BIT
     238      934630 :   for (i=2; i<l; i++) uel(y,i) = rand64();
     239             : #else
     240       77978 :   for (i=2; i<l-1; i+=2)
     241             :   {
     242         135 :     u64 v = rand64();
     243         135 :     uel(y,i)   = (ulong) v;
     244         135 :     uel(y,i+1) = (ulong) (v>>32);
     245             :   }
     246       77843 :   if (i<l) uel(y,i) = (ulong) rand64();
     247             : #endif
     248      544855 :   if (db) uel(y,l-1) &= ((1UL<<db)-1UL);
     249      544855 :   return F2x_renormalize(y,l);
     250             : }
     251             : 
     252             : GEN
     253          56 : random_zv(long n)
     254             : {
     255          56 :   GEN y = cgetg(n+1, t_VECSMALL);
     256             :   long i;
     257       78711 :   for (i=1; i<=n; i++) uel(y,i) = pari_rand();
     258          56 :   return y;
     259             : }
     260             : 
     261             : GEN
     262           7 : randomr(long prec)
     263             : {
     264             :   pari_sp av;
     265             :   long b;
     266             :   GEN x, y;
     267           7 :   if (prec <= 2) return real_0_bit(0);
     268           7 :   x = cgetr(prec); av = avma;
     269           7 :   b = prec2nbits(prec);
     270           7 :   y = randomi(int2n(b));
     271           7 :   if (!signe(y)) return real_0_bit(b);
     272           7 :   affir(y, x); shiftr_inplace(x, - b);
     273           7 :   set_avma(av); return x;
     274             : }
     275             : 
     276             : static GEN
     277      156296 : polrandom(GEN N) /* assume N!=0 */
     278             : {
     279      156296 :   long i, d = lg(N);
     280      156296 :   GEN z = leading_coeff(N);
     281      156296 :   GEN y = cgetg(d,t_POL);
     282      156296 :   y[1] = evalsigne(1) | evalvarn(varn(N));
     283      653268 :   for (i=2; i<d; i++) gel(y,i) = genrand(z);
     284      156296 :   return normalizepol_lg(y,d);
     285             : }
     286             : 
     287             : GEN
     288     5520599 : genrand(GEN N)
     289             : {
     290             :   GEN z;
     291     5520599 :   if (!N) return utoi( random_bits(31) );
     292     5519920 :   switch(typ(N))
     293             :   {
     294      736995 :     case t_INT:
     295      736995 :       switch(signe(N))
     296             :       {
     297             :         pari_sp av;
     298             :         GEN d;
     299      282674 :         case 1:
     300      282674 :           return randomi(N);
     301      454314 :         case -1:
     302      454314 :           av = avma; N = addiu(N, 1); d = subui(1, shifti(N, 1));
     303      454314 :           return gerepileuptoint(av, addii(N, randomi(d)));
     304           7 :         default: pari_err_DOMAIN("random","N","=",gen_0,gen_0);
     305             :       }
     306           7 :     case t_REAL:
     307           7 :       return randomr(realprec(N));
     308      619339 :     case t_INTMOD:
     309      619339 :       z = cgetg(3, t_INTMOD);
     310      619339 :       gel(z,1) = icopy(gel(N,1));
     311      619339 :       gel(z,2) = randomi(gel(N,1)); return z;
     312      187474 :     case t_FFELT:
     313      187474 :       return ffrandom(N);
     314      156296 :     case t_POL:
     315      156296 :       if (signe(N)==0) return pol_0(varn(N));
     316      156296 :       return polrandom(N);
     317     3819802 :     case t_VEC:
     318     3819802 :       if (lg(N) == 3)
     319             :       {
     320     3575824 :         pari_sp av = avma;
     321     3575824 :         GEN a = gel(N,1), b = gel(N,2), d;
     322     3575824 :         if (typ(a) != t_INT) a = gceil(a);
     323     3575824 :         if (typ(b) != t_INT) b = gfloor(b);
     324     3575824 :         if (typ(a) != t_INT || typ(b) != t_INT) pari_err_TYPE("random", N);
     325     3575824 :         d = subii(b,a);
     326     3575824 :         if (signe(d) < 0) pari_err_TYPE("random([a,b]) (a > b)", N);
     327     3575824 :         return gerepileuptoint(av, addii(a, randomi(addiu(d,1))));
     328             :       }
     329      243978 :       return ellrandom(N);
     330           7 :     default:
     331           7 :       pari_err_TYPE("genrand",N);
     332             :       return NULL;/*LCOV_EXCL_LINE*/
     333             :   }
     334             : }

Generated by: LCOV version 1.13