1: //
   2: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
   3: // integration for neuronal models.
   4: //
   5: 
   6: //////////////////////////////////////////////////////////////////////////////
   7: //'
   8: //' Heccer : testbed C implementation
   9: //'
  10: //' Copyright (C) 2006-2008 Hugo Cornelis
  11: //'
  12: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
  13: //'
  14: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
  15: //'
  16: //////////////////////////////////////////////////////////////////////////////
  17: 
  18: 
  19: #include <stdio.h>
  20: 
  21: #include "heccer/random.h"
  22: 
  23: 
  24: /*
  25: ** Random number routines from "Numerical Recipes in C", Chapter 7
  26: ** Entered by Michael D. Speight, 10th September 1991
  27: ** All references to '*idum' changed to 'idum'
  28: */
  29: #define M1 259200
  30: #define IA1 7141
  31: #define IC1 54773
  32: #define RM1 (1.0)/M1
  33: #define M2 134456
  34: #define IA2 8121
  35: #define IC2 28411
  36: #define RM2 (1.0)/M2
  37: #define M3 243000
  38: #define IA3 4561
  39: #define IC3 51349
  40: 
  41: float ran1(time_t idum)
  42: {
  43:   static long ix1,ix2,ix3;
  44:   static float r[98];
  45:   float temp;
  46:   static int iff=0;
  47:   int j;
  48:   if (idum < 0 || iff == 0) {
  49:     iff=1;
  50:     ix1=(IC1-(idum)) % M1;
  51:     ix1=(IA1*ix1+IC1) % M1;
  52:     ix2=ix1 % M2;
  53:     ix1=(IA1*ix1+IC1) % M1;
  54:     ix3=ix1 % M3;
  55:     for (j=1;j<=97;j++) {
  56:       ix1=(IA1*ix1+IC1) % M1;
  57:       ix2=(IA2*ix2+IC2) % M2;
  58:       r[j]=(ix1+ix2*RM2)*RM1;
  59:     }
  60:     idum=1;
  61:   }
  62:   ix1=(IA1*ix1+IC1) % M1;
  63:   ix2=(IA2*ix2+IC2) % M2;
  64:   ix3=(IA3*ix3+IC3) % M3;
  65:   j=1 + ((97*ix3)/M3);
  66:   if (j>97||j<1) printf("RAN1: This cannot happen.\n");
  67:   temp=r[j];
  68:   r[j]=(ix1+ix2*RM2)*RM1;
  69:   return temp;
  70: }
  71: 
  72: