vppinfra/vppinfra/random_isaac.c

   1 /*
   2  * Copyright (c) 2015 Cisco and/or its affiliates.
   3  * Licensed under the Apache License, Version 2.0 (the "License");
   4  * you may not use this file except in compliance with the License.
   5  * You may obtain a copy of the License at:
   6  *
   7  *     http://www.apache.org/licenses/LICENSE-2.0
   8  *
   9  * Unless required by applicable law or agreed to in writing, software
  10  * distributed under the License is distributed on an "AS IS" BASIS,
  11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12  * See the License for the specific language governing permissions and
  13  * limitations under the License.
  14  */
  15 /*
  16   ------------------------------------------------------------------------------
  17   By Bob Jenkins, 1996, Public Domain
  18   MODIFIED:
  19   960327: Creation (addition of randinit, really)
  20   970719: use context, not global variables, for internal state
  21   980324: renamed seed to flag
  22   980605: recommend ISAAC_LOG2_SIZE=4 for noncryptography.
  23   010626: note this is public domain
  24   ------------------------------------------------------------------------------
  25
  26   Modified for CLIB by Eliot Dresselhaus.
  27   Dear Bob, Thanks for all the great work. - Eliot
  28
  29   modifications copyright (c) 2003 Eliot Dresselhaus
  30
  31   Permission is hereby granted, free of charge, to any person obtaining
  32   a copy of this software and associated documentation files (the
  33   "Software"), to deal in the Software without restriction, including
  34   without limitation the rights to use, copy, modify, merge, publish,
  35   distribute, sublicense, and/or sell copies of the Software, and to
  36   permit persons to whom the Software is furnished to do so, subject to
  37   the following conditions:
  38
  39   The above copyright notice and this permission notice shall be
  40   included in all copies or substantial portions of the Software.
  41
  42   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  43   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  44   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  45   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
  46   LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
  47   OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
  48   WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  49 */
  50
  51 /* ISAAC is Bob Jenkins' random number generator.
  52    http://burtleburtle.net/bob/rand/isaacafa.html */
  53
  54 #include <vppinfra/random_isaac.h>
  55
  56 #if uword_bits != 32 && uword_bits != 64
  57 #error "isaac only works for 32 or 64 bit words"
  58 #endif
  59
  60 #if uword_bits == 32
  61
  62 #define ind32(mm,x)  (*(u32 *)((u8 *)(mm) + ((x) & ((ISAAC_SIZE-1)<<2))))
  63 #define rngstep32(mix,a,b,mm,m,m2,r,x,y)                \
  64 {                                                       \
  65   x = *m;                                               \
  66   a = (a^(mix)) + *(m2++);                              \
  67   *(m++) = y = ind32(mm,x) + a + b;                     \
  68   *(r++) = b = ind32(mm,y>>ISAAC_LOG2_SIZE) + x;        \
  69 }
  70
  71 void isaac (isaac_t * ctx, uword * results)
  72 {
  73   u32 a, b, c, x, y, * m, * mm, * m2, * r, * mend;
  74
  75   mm = ctx->memory;
  76   r = results;
  77   a = ctx->a;
  78   b = ctx->b;
  79   c = ctx->c;
  80
  81   b += ++c;
  82   mend = m2 = mm + ARRAY_LEN (ctx->memory) / 2;
  83   m = mm;
  84   while (m < mend)
  85     {
  86       rngstep32 (a<<13, a, b, mm, m, m2, r, x, y);
  87       rngstep32 (a>>6 , a, b, mm, m, m2, r, x, y);
  88       rngstep32 (a<<2 , a, b, mm, m, m2, r, x, y);
  89       rngstep32 (a>>16, a, b, mm, m, m2, r, x, y);
  90     }
  91
  92   m2 = mm;
  93   while (m2 < mend)
  94     {
  95       rngstep32 (a<<13, a, b, mm, m, m2, r, x, y);
  96       rngstep32 (a>>6 , a, b, mm, m, m2, r, x, y);
  97       rngstep32 (a<<2 , a, b, mm, m, m2, r, x, y);
  98       rngstep32 (a>>16, a, b, mm, m, m2, r, x, y);
  99     }
 100
 101   ctx->a = a;
 102   ctx->b = b;
 103   ctx->c = c;
 104 }
 105
 106 /* Perform 2 isaac runs with different contexts simultaneously. */
 107 void isaac2 (isaac_t * ctx, uword * results)
 108 {
 109 #define _(n) \
 110   u32 a##n, b##n, c##n, x##n, y##n, * m##n, * mm##n, * m2##n, * r##n, * mend##n
 111
 112   _ (0); _ (1);
 113
 114 #undef _
 115
 116 #define _(n)                                                    \
 117 do {                                                            \
 118   mm##n = ctx[(n)].memory;                                      \
 119   r##n = results + (n) * ISAAC_SIZE;                            \
 120   a##n = ctx[(n)].a;                                            \
 121   b##n = ctx[(n)].b;                                            \
 122   c##n = ctx[(n)].c;                                            \
 123   b##n += ++c##n;                                               \
 124   mend##n = m2##n = mm##n + ARRAY_LEN (ctx[(n)].memory) / 2;    \
 125   m##n = mm##n;                                                 \
 126 } while (0)
 127
 128   _ (0); _ (1);
 129
 130 #undef _
 131
 132   while (m0 < mend0)
 133     {
 134       rngstep32 (a0<<13, a0, b0, mm0, m0, m20, r0, x0, y0);
 135       rngstep32 (a1<<13, a1, b1, mm1, m1, m21, r1, x1, y1);
 136       rngstep32 (a0>>6 , a0, b0, mm0, m0, m20, r0, x0, y0);
 137       rngstep32 (a1>>6 , a1, b1, mm1, m1, m21, r1, x1, y1);
 138       rngstep32 (a0<<2 , a0, b0, mm0, m0, m20, r0, x0, y0);
 139       rngstep32 (a1<<2 , a1, b1, mm1, m1, m21, r1, x1, y1);
 140       rngstep32 (a0>>16, a0, b0, mm0, m0, m20, r0, x0, y0);
 141       rngstep32 (a1>>16, a1, b1, mm1, m1, m21, r1, x1, y1);
 142     }
 143
 144   m20 = mm0;
 145   m21 = mm1;
 146   while (m20 < mend0)
 147     {
 148       rngstep32 (a0<<13, a0, b0, mm0, m0, m20, r0, x0, y0);
 149       rngstep32 (a1<<13, a1, b1, mm1, m1, m21, r1, x1, y1);
 150       rngstep32 (a0>>6 , a0, b0, mm0, m0, m20, r0, x0, y0);
 151       rngstep32 (a1>>6 , a1, b1, mm1, m1, m21, r1, x1, y1);
 152       rngstep32 (a0<<2 , a0, b0, mm0, m0, m20, r0, x0, y0);
 153       rngstep32 (a1<<2 , a1, b1, mm1, m1, m21, r1, x1, y1);
 154       rngstep32 (a0>>16, a0, b0, mm0, m0, m20, r0, x0, y0);
 155       rngstep32 (a1>>16, a1, b1, mm1, m1, m21, r1, x1, y1);
 156     }
 157
 158   ctx[0].a = a0;
 159   ctx[0].b = b0;
 160   ctx[0].c = c0;
 161   ctx[1].a = a1;
 162   ctx[1].b = b1;
 163   ctx[1].c = c1;
 164 }
 165
 166 #define mix32(a,b,c,d,e,f,g,h)                  \
 167 {                                               \
 168    a^=b<<11; d+=a; b+=c;                        \
 169    b^=c>>2;  e+=b; c+=d;                        \
 170    c^=d<<8;  f+=c; d+=e;                        \
 171    d^=e>>16; g+=d; e+=f;                        \
 172    e^=f<<10; h+=e; f+=g;                        \
 173    f^=g>>4;  a+=f; g+=h;                        \
 174    g^=h<<8;  b+=g; h+=a;                        \
 175    h^=a>>9;  c+=h; a+=b;                        \
 176 }
 177
 178 void isaac_init (isaac_t * ctx, uword * seeds)
 179 {
 180    word i;
 181    u32 a, b, c, d, e, f, g, h, * m, * r;
 182
 183    ctx->a = ctx->b = ctx->c = 0;
 184    m = ctx->memory;
 185    r = seeds;
 186
 187    a = b = c = d = e = f = g = h = 0x9e3779b9;  /* the golden ratio */
 188
 189    for (i = 0; i < 4; ++i)          /* scramble it */
 190      mix32(a,b,c,d,e,f,g,h);
 191
 192    /* initialize using the contents of r[] as the seed */
 193    for (i=0; i<ISAAC_SIZE; i+=8)
 194      {
 195        a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
 196        e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];
 197        mix32(a,b,c,d,e,f,g,h);
 198        m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
 199        m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
 200      }
 201
 202    /* do a second pass to make all of the seed affect all of m */
 203    for (i=0; i<ISAAC_SIZE; i+=8)
 204      {
 205        a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
 206        e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];
 207        mix32(a,b,c,d,e,f,g,h);
 208        m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
 209        m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
 210      }
 211 }
 212 #endif /* uword_bits == 32 */
 213
 214 #if uword_bits == 64
 215
 216 #define ind64(mm,x)  (*(u64 *)((u8 *)(mm) + ((x) & ((ISAAC_SIZE-1)<<3))))
 217 #define rngstep64(mix,a,b,mm,m,m2,r,x,y)                \
 218 {                                                       \
 219   x = *m;                                               \
 220   a = (mix) + *(m2++);                                  \
 221   *(m++) = y = ind64(mm,x) + a + b;                     \
 222   *(r++) = b = ind64(mm,y>>ISAAC_LOG2_SIZE) + x;        \
 223 }
 224
 225 void isaac (isaac_t * ctx, uword * results)
 226 {
 227   u64 a, b, c, x, y, * m, * mm, * m2, * r, * mend;
 228
 229   mm = ctx->memory;
 230   r = results;
 231   a = ctx->a;
 232   b = ctx->b;
 233   c = ctx->c;
 234
 235   b += ++c;
 236   mend = m2 = mm + ARRAY_LEN (ctx->memory) / 2;
 237   m = mm;
 238   while (m < mend)
 239     {
 240       rngstep64 (~(a^(a<<21)), a, b, mm, m, m2, r, x, y);
 241       rngstep64 (  a^(a>>5)  , a, b, mm, m, m2, r, x, y);
 242       rngstep64 (  a^(a<<12) , a, b, mm, m, m2, r, x, y);
 243       rngstep64 (  a^(a>>33) , a, b, mm, m, m2, r, x, y);
 244     }
 245
 246   m2 = mm;
 247   while (m2 < mend)
 248     {
 249       rngstep64 (~(a^(a<<21)), a, b, mm, m, m2, r, x, y);
 250       rngstep64 (  a^(a>>5)  , a, b, mm, m, m2, r, x, y);
 251       rngstep64 (  a^(a<<12) , a, b, mm, m, m2, r, x, y);
 252       rngstep64 (  a^(a>>33) , a, b, mm, m, m2, r, x, y);
 253     }
 254
 255   ctx->a = a;
 256   ctx->b = b;
 257   ctx->c = c;
 258 }
 259
 260 /* Perform 2 isaac runs with different contexts simultaneously. */
 261 void isaac2 (isaac_t * ctx, uword * results)
 262 {
 263 #define _(n) \
 264   u64 a##n, b##n, c##n, x##n, y##n, * m##n, * mm##n, * m2##n, * r##n, * mend##n
 265
 266   _ (0); _ (1);
 267
 268 #undef _
 269
 270 #define _(n)                                                    \
 271 do {                                                            \
 272   mm##n = ctx[(n)].memory;                                      \
 273   r##n = results + (n) * ISAAC_SIZE;                            \
 274   a##n = ctx[(n)].a;                                            \
 275   b##n = ctx[(n)].b;                                            \
 276   c##n = ctx[(n)].c;                                            \
 277   b##n += ++c##n;                                               \
 278   mend##n = m2##n = mm##n + ARRAY_LEN (ctx[(n)].memory) / 2;    \
 279   m##n = mm##n;                                                 \
 280 } while (0)
 281
 282   _ (0); _ (1);
 283
 284 #undef _
 285
 286   (void) mend1;                 /* compiler warning */
 287
 288   while (m0 < mend0)
 289     {
 290       rngstep64 (~(a0^(a0<<21)), a0, b0, mm0, m0, m20, r0, x0, y0);
 291       rngstep64 (~(a1^(a1<<21)), a1, b1, mm1, m1, m21, r1, x1, y1);
 292       rngstep64 (  a0^(a0>>5)  , a0, b0, mm0, m0, m20, r0, x0, y0);
 293       rngstep64 (  a1^(a1>>5)  , a1, b1, mm1, m1, m21, r1, x1, y1);
 294       rngstep64 (  a0^(a0<<12) , a0, b0, mm0, m0, m20, r0, x0, y0);
 295       rngstep64 (  a1^(a1<<12) , a1, b1, mm1, m1, m21, r1, x1, y1);
 296       rngstep64 (  a0^(a0>>33) , a0, b0, mm0, m0, m20, r0, x0, y0);
 297       rngstep64 (  a1^(a1>>33) , a1, b1, mm1, m1, m21, r1, x1, y1);
 298     }
 299
 300   m20 = mm0;
 301   m21 = mm1;
 302   while (m20 < mend0)
 303     {
 304       rngstep64 (~(a0^(a0<<21)), a0, b0, mm0, m0, m20, r0, x0, y0);
 305       rngstep64 (~(a1^(a1<<21)), a1, b1, mm1, m1, m21, r1, x1, y1);
 306       rngstep64 (  a0^(a0>>5)  , a0, b0, mm0, m0, m20, r0, x0, y0);
 307       rngstep64 (  a1^(a1>>5)  , a1, b1, mm1, m1, m21, r1, x1, y1);
 308       rngstep64 (  a0^(a0<<12) , a0, b0, mm0, m0, m20, r0, x0, y0);
 309       rngstep64 (  a1^(a1<<12) , a1, b1, mm1, m1, m21, r1, x1, y1);
 310       rngstep64 (  a0^(a0>>33) , a0, b0, mm0, m0, m20, r0, x0, y0);
 311       rngstep64 (  a1^(a1>>33) , a1, b1, mm1, m1, m21, r1, x1, y1);
 312     }
 313
 314   ctx[0].a = a0;
 315   ctx[0].b = b0;
 316   ctx[0].c = c0;
 317   ctx[1].a = a1;
 318   ctx[1].b = b1;
 319   ctx[1].c = c1;
 320 }
 321
 322 #define mix64(a,b,c,d,e,f,g,h)                  \
 323 {                                               \
 324    a-=e; f^=h>>9;  h+=a;                        \
 325    b-=f; g^=a<<9;  a+=b;                        \
 326    c-=g; h^=b>>23; b+=c;                        \
 327    d-=h; a^=c<<15; c+=d;                        \
 328    e-=a; b^=d>>14; d+=e;                        \
 329    f-=b; c^=e<<20; e+=f;                        \
 330    g-=c; d^=f>>17; f+=g;                        \
 331    h-=d; e^=g<<14; g+=h;                        \
 332 }
 333
 334 void isaac_init (isaac_t * ctx, uword * seeds)
 335 {
 336   word i;
 337   u64 a, b, c, d, e, f, g, h, * m, * r;
 338
 339   ctx->a = ctx->b = ctx->c = 0;
 340   m = ctx->memory;
 341   r = seeds;
 342
 343   a = b = c = d = e = f = g = h = 0x9e3779b97f4a7c13LL;  /* the golden ratio */
 344
 345   for (i=0; i<4; ++i)                    /* scramble it */
 346     mix64(a,b,c,d,e,f,g,h);
 347
 348   for (i=0; i<ISAAC_SIZE; i+=8)   /* fill in mm[] with messy stuff */
 349     {
 350       a+=r[i  ]; b+=r[i+1]; c+=r[i+2]; d+=r[i+3];
 351       e+=r[i+4]; f+=r[i+5]; g+=r[i+6]; h+=r[i+7];
 352       mix64(a,b,c,d,e,f,g,h);
 353       m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
 354       m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
 355     }
 356
 357   /* do a second pass to make all of the seed affect all of mm */
 358   for (i=0; i<ISAAC_SIZE; i+=8)
 359     {
 360       a+=m[i  ]; b+=m[i+1]; c+=m[i+2]; d+=m[i+3];
 361       e+=m[i+4]; f+=m[i+5]; g+=m[i+6]; h+=m[i+7];
 362       mix64(a,b,c,d,e,f,g,h);
 363       m[i  ]=a; m[i+1]=b; m[i+2]=c; m[i+3]=d;
 364       m[i+4]=e; m[i+5]=f; m[i+6]=g; m[i+7]=h;
 365     }
 366 }
 367 #endif /* uword_bits == 64 */
 368