File: RA.c 10/15/1998, 5:32:14 PM
`
`1 /*
`2 RA.c
`3 (c) DJCM 98 09 28
`
` Repeat-accumulate code simulator
`
`45
`
`67
`
` read in code definition
`8 loop {
`9 encode source string
`10 add noise
`11 decode
`12 }
`13
`14 Code definition: (stored in "alist")
`15
`16 Use of alist allows arbitrary numbers of repetitions
`17 of each bit.
`18
`19 K source block length
`20 n_1 n_2 ... n_K number of repetitions of each source bit
`21 N = sum n_k
`22 alist defines permutation of N encoded bits
`23 note, an additional permutation of the N
` accumulated
`24 bits may be a good idea. (for non-memoryless channels)
`25
`26 transmitted bits are integral of encoded bits
`27
`28 Future plans:
`29 clump source bits into clumps. Have multiple parallel accumulated
` streams.
`30 Have little sub-matrices (like GF(q) ) defining response of
` accumulator to
`31 clumps.
`32
`33 */
`34
`35 #include "./ansi/r.h"
`36 #include "./ansi/rand2.h"
`37 #include "./ansi/mynr.h"
`38 #include "./ansi/cmatrix.h"
`39
`40 #include "./RA.h" /* this defines data_creation_param ; RA_control */
`41
`42 int RA_encode ( unsigned char * , RA_control * , unsigned char * ) ;
`43 static int t_to_b ( unsigned char * , RA_control * ) ;
`44 int RA_decode ( RA_control * ) ;
`45 int RA_horizontal_pass ( RA_control * ) ;
`46
`47 static void dc_defaults ( RA_control * ) ;
`48 static int process_command ( int , char ** , RA_control * ) ;
`49 static void print_usage ( char ** , FILE * , RA_control * );
`50 static int make_sense ( RA_control * ) ;
`51 static int make_space ( RA_control * ) ;
`52 static int score ( RA_control * ) ;
`53
`54 static void finalline ( FILE * , RA_control * , int ) ; /* int = 1 to
`get loads of info */
`
`Page: 1
`
`Hughes, Exh. 1019, p. 1
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`55
`56 static double bern ( int , int , double * , double * , double * ,
`double );
`57 static void histo ( FILE * , RA_control * ) ;
`58 static void snappyline ( RA_control * ) ;
`59 static void RA_free ( RA_control * ) ;
`60 static int check_alist_MN ( alist_matrix * , RA_control * ) ;
`61 static double h2 ( double ) ;
`62
`63 void main ( int , char ** ) ;
`64
`65 /*
`66 MAIN
`67 */
`68 void main ( int argc, char *argv[] )
`69 {
`70 FILE *fp ;
`71
`72 int k ;
`73 RA_control c ;
`74 dc_defaults ( &c ) ;
`75
`76 if ( process_command (argc, argv, &c ) < 0 ) exit (0) ;
`77 if ( read_allocate_alist ( &(c.a) , c.afile ) < 0 ) exit (0) ;
`78 if ( check_alist_MN ( &(c.a) , &c ) < 0 ) exit (0) ;
`79 if ( make_sense ( &c ) < 0 ) exit (0) ;
`80
`81 fprintf(stderr,"RA N=%d, K=%d, x=%6.3g xass=%6.3g fn=%6.3g
` fnass=%6.3g\n",
`82 c.N , c.K , c.gcx , c.gcxass , c.fn , c.fnass ) ;
`83 fflush(stderr);
`84
`85 if ( make_space ( &c ) < 0 ) exit (0) ;
`86
`87 if ( c.writelog ) {
`88 fp = fopen ( c.logfile , "w" ) ;
`89 if ( !fp ) {
`90 fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ;
`91 c.writelog = 0 ;
`92 } else fclose (fp ) ;
`93 }
`94
`95 if ( c.writelog ) {
`96 fp = fopen ( c.logfile , "w" ) ;
`97 if ( !fp ) {
`98 fprintf ( stderr , " couldn't open logfile %s\n" , c.logfile ) ;
`99 c.writelog = 0 ;
`100 } else fclose (fp ) ;
`101 }
`102
`103 if ( c.error_log ) {
`104 fp = fopen ( c.error_logfile , "w" ) ;
`105 if ( !fp ) {
`106 fprintf ( stderr , " couldn't open logfile %s\n" , c.error_logfile
` ) ;
`107 c.error_log = 0 ;
`108 } else fclose (fp ) ;
`109 }
`
`Page: 2
`
`Hughes, Exh. 1019, p. 2
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`110
`111 /*
`112 MAIN LOOP
`113 */
`114
`115 ran_seed ( c.vseed ) ;
`116 c.message = 1 ;
`117 for ( ; ( c.message <= c.MESSAGE ) &&
`118 ( ( c.failures==0 ) || ( c.failcount < c.failures ) ) ;
`119 c.message ++ ) {
`120
`121 snappyline( &c ) ;
`122 /* force parity bit at end */
`123 c.sourceweight = random_cvector ( c.s , c.fs , 1 , c.K ) ;
`124
`125 if ( c.verbose > 2) {
`126 printf ( "source vector:\n" ) ;
`127 for ( k = 1 ; k <= c.K ; k ++ ) {
`128 if ( c.s[k] ) printf ("1 ") ; else printf ( "0 " );
`129 }
`130 printf ( "\n" ) ;
`131 }
`132 RA_encode ( c.s , &c , c.t ) ;
`133 if ( c.verbose > 2) {
`134 printf ( "transmitted vector:\n" ) ;
`135 for ( k = 1 ; k <= c.N ; k ++ ) {
`136 if ( c.t[k] ) printf ("1 ") ; else printf ( "0 " );
`137 }
`138 printf ( "\n" ) ;
`139 }
`140 c.flipped = t_to_b ( c.t , &c ) ;
`141 if ( c.verbose > 2 ) {
`142 printf ( "received likelihoods:\n" ) ;
`143 for ( k = 1 ; k <= c.N ; k ++ ) {
`144 printf ("%1d " , (int) ( c.bias[k][1] * 10.0 ) ) ;
`145 }
`146 printf ( "\n" ) ;
`147 for ( k = 1 ; k <= c.N ; k ++ ) {
`148 printf ("%1d " , (int) ( c.bias[k][1] * 2.0 ) ) ;
`149 }
`150 printf ( "\n" ) ;
`151 }
`152 RA_decode ( &c ) ;
`153
`154 if ( score ( &c ) < 0 ) exit ( 0 ) ;
`155 if ( c.verbose > 0 ) finalline ( stdout , &c , 0 ) ;
`156 if ( c.printout ) { /* append */
`157 fp = fopen ( c.outfile , ((c.outappend)? "a":"w") ) ;
`158 if( !fp ) {
`159 fprintf( stderr, "No such file: %s\n", c.outfile ) ;
`160 finalline ( stderr , &c , 0 ) ;
`161 } else {
`162 finalline ( fp , &c , 0 ) ;
`163 fclose ( fp ) ;
`164 }
`165 }
`166 if ( (!( ( c.message+1 <= c.MESSAGE ) &&
`167 ( ( c.failures==0 )
`
`Page: 3
`
`Hughes, Exh. 1019, p. 3
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`168 || ( c.failcount < c.failures ) ) ) )
`169 || (!(c.message % c.big_write_period ) ) ) {
`170 if ( c.printtot ) { /* write */
`171 fp = fopen ( c.totoutfile , "w" ) ;
`172 if( !fp ) {
`173 fprintf( stderr, "No such file: %s\n", c.totoutfile ) ;
`174 finalline ( stderr , &c , 1 ) ; /* totalline */
`175 } else {
`176 finalline ( fp , &c , 1 ) ;
`177 fclose ( fp ) ;
`178 }
`179 }
`180 if ( c.printhisto && c.block_valid ) { /* update histogram file */
`181 fp = fopen ( c.histofile , "w" ) ;
`182 if( !fp ) {
`183 fprintf( stderr, "No such file: %s\n", c.histofile ) ;
`184 } else {
`185 histo ( fp , &c ) ;
`186 fclose ( fp ) ;
`187 }
`188 }
`189 }
`190 }
`191 snappyline( &c ) ; printf("\n") ;
`192 RA_free ( &c ) ;
`193 }
`194
`195 static void histo ( FILE *fp , RA_control *c ) {
`196 int l;
`197 double t , cum = 0.0 ;
`198 double tot = (double) c->block_valid ;
`199
`200 fprintf ( fp , "# total valid blocks %d\n" , c->block_valid ) ;
`201 for ( l = 1 ; l <= c->loops ; l ++ ) {
`202 t= (double) c->histo[l] ;
`203 cum += t ;
`204 fprintf ( fp , "%d\t%d\t%d\t%9.4g\t%9.4g\n" , l , (int)(t) ,
` (int)(cum) ,
`205 t/tot , cum/tot ) ;
`206 }
`207
`208 }
`209 static void snappyline ( RA_control *c ) {
`210 printf ( "%d:%du%dd%dl/%d\t" , c->block_errs , c->block_undet,
` c->block_det, c->block_detlw , c->message - 1 ) ; fflush ( stdout ) ;
`
`211
`212 }
`213
`214 /* <<<<N>>>>
`215 Encoding method:
`216 source bits d[1]..d[K] are mapped via an alist ^ 1 1 1
`217 into a pre-transmission vector K 1 1 1
`218 s[1]..s[N] . V 1 1 1
`219
`220 t[n] = t[n-1] ^ s[n]
`221
`222 s[n] s[n+1]
`223 | |
`
`Page: 4
`
`Hughes, Exh. 1019, p. 4
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`224 0 .. -+-> t[n] -+-> t[n+1] ..... t[N]
`225 | | |
`226 V V V
`227 y[n] y[n+1] y[N]
`228
`229 */
`230
`231 int RA_encode ( unsigned char *d , RA_control *c , unsigned char *t ) {
`232 int n , k ; int status = 0 ;
`233 alist_transpose_cvector_sparse_mod2 ( &c->a , d , t ) ; /* here 't'
` doubles as 's' */
`234 /* accumulate */
`235
`236 if ( c->verbose > 2) {
`237 printf ( "extended source vector:\n" ) ;
`238 for ( k = 1 ; k <= c->N ; k ++ ) {
`239 if ( t[k] ) printf ("1 ") ; else printf ( "0 " );
`240 }
`241 printf ( "\n" ) ;
`242 }
`243
`244 for ( n = 2 ; n <= c->N ; n ++ ) {
`245 t[n] = t[n]^t[n-1] ;
`246 }
`247
`248 if ( c->verbose > 2) {
`249 printf ( "accumulated transmission:\n" ) ;
`250 for ( k = 1 ; k <= c->N ; k ++ ) {
`251 if ( t[k] ) printf ("1 ") ; else printf ( "0 " );
`252 }
`253 printf ( "\n" ) ;
`254 }
`255 return status ;
`256 }
`257
`258 /*
`259 The channel outputs a normalized likelihood vector
`260 bias[n] = P( yn | tn = 1 )
`261
`262 The state of the decoder is q[1..N][0/1] and r[1..N][0/1]
`263
`264 q[n][s] = pseudoprior( s[n] = s ) s=0/1 initially 0.5
`265
`266 Use f/b algorithm to find:
`267 initial conditions:
`268 f[n][t] = P( y1...yn , tn=t ) f[0][0] = 1 ; f[0][1] = 0 ;
`269 b[n][t] = P( yn...yN | tn=t ) b[N+1][0] = 1 ; b[N+1][1] = 1 ;
`270
`271 Using
`272 f[n][t] = bias[n][t] * sum_{t': t'+s=t} ( f[n-1][t'] pi[n][s] )
`273 b[n][t] = bias[n][t] * sum_{t': t'+s=t} ( b[n+1][t'] pi[n+1][s] )
`274
`275 Find likelihood contribution at n:
`276 r[n][s] `=' P(y1..yN|s[n]=s) = sum_{s: t'+s=t} f[n-1][t] b[n][t']
`277
`278 Then in the vertical step we visit each incarnation of the source bit
`279 for( r = 1 .. repetitions[k] ) (number on mlist) n=nlist[r]
`280 d[r][s] = d[r-1][s] * r[n][s] ;
`
`Page: 5
`
`Hughes, Exh. 1019, p. 5
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`281 reverse,
`282 u[r][s] = u[r+1][s] * r[n][s] ;
`283 ->
`284 pi[n][s] `=' P( s[n] = s | y1..yN ) = prior0(omit) * d[r-1][s] *
` u[r+1][s]
`
`285
`286 If all l's have the same sign then the locally mp edges in the
` trellis
`287 form a valid codeword
`288
`289 */
`290
`291 int RA_decode ( RA_control *c ) {
`292 int u , n , N=c->N ;
`293
`294 for ( n = 1 ; n <= N ; n++ ) {
`295 for ( u = 0 ; u <= 1 ; u++ ) {
`296 c->q[n][u] = 1.0 ; /* assume all data strings equiprobable */
`297 }
`298 }
`299
`300 for ( c->loop = 1 , c->viols = 1 ;
`301 ( c->loop <= c->loops ) && c->viols ;
`302 c->loop ++ ) {
`303
`304 c->viols = RA_horizontal_pass ( c ) ;
`305 /* if ( c->writelog ) RA_fprint_state ( c ) ; */
`306
`307 }
`308 c->loop -- ;
`309 printf( "\tviols %3d its %2d\n",
`310 c->viols , c->loop ) ;
`311
`312 return ( c->viols ) ; /* zero if success */
`313
`314 }
`315
`316 /* for efficiency, the backward pass variables are not stored;
`317 they are immediately multiplied by the forward pass variables
`318 to give likelihood signals.
`319 The likelihood signals are (at present) written into the same
`320 array as the f's
`321 */
`322
`323 int RA_horizontal_pass ( RA_control *c ) {
`324 int u , umax ;
`325 int m , n , N=c->N ; /* ? */
`326 double **f = c->f ;
`327 /* f has been initialized with f[0][0] = 1 ; f[0][1] = 0 ; */
`328
`329 /* double **b=c->b ; b[N+1][0] = 1 ; b[N+1][1] = 1 ; */
`330 double **bias = c->bias ; /* likelihood */
`331 double **q = c->q ;
`332 double **r = c->r ;
`333 double p_t0 , p_t1 ;
`334 int *mlist ;
`335 int v0 , v1 , viols ; /* votes for 0 and 1 ; viols sees if they
` disagree */
`
`Page: 6
`
`Hughes, Exh. 1019, p. 6
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`336 /* int TERMINATED = 0 ; */
`337 double BOT = 0.99 ;
`338 double TOP = 1.01 ;
`339 alist_matrix *a = &(c->a) ;
`340 unsigned char *so = c->so ;
`341 unsigned char *to = c->to ;
`342 double f1 = 0.0 , f0 = 1.0 , b1 , b0 , p0 , p1 , q0 , q1 , s ;
`343 int clipwarn=0 ; /* whether we have spat a clip warning this time */
`344
`345 for ( n = 1 ; n <= N ; n++ ) { /* n=N is not really needed */
`346 q0 = q[n][0] ; q1 = q[n][1] ; /* latest prior probabilities */
`347 p0 = q0 * f0 + q1 * f1 ;
`348 p1 = q1 * f0 + q0 * f1 ;
`349 if ( c->verbose > 4) {
`350 printf ( "n = %d\n" , n ) ;
`351 printf ( "q: %10.4g %10.4g\n" , q0 , q1 ) ;
`352 printf ( "f: %10.4g %10.4g\n" , f0 , f1 ) ;
`353 printf ( "p: %10.4g %10.4g\n" , p0 , p1 ) ;
`354 printf ( "b: %10.4g %10.4g\n" , bias[n][0] , bias[n][1] ) ;
`355 printf ( "\n" ) ;
`356 pause_for_return() ;
`357 }
`358 f0 = p0 * bias[n][0] ; f1 = p1 * bias[n][1] ; s=f0+f1 ;
`359 if ( s < BOT || ( s > TOP ) ) { s = 1.0/s ; f0 *= s ; f1 *= s ; }
`360 f[n][0] = f0 ; f[n][1] = f1 ;
`361 /*
`362 f[n][1] = bias[n][1] * ( f[n-1][0] * q[n][1] + f[n-1][1] * q[n][0]
` );
`363 f[n][0] = bias[n][0] * ( f[n-1][0] * q[n][0] + f[n-1][1] * q[n][1]
` );
`
`364
`365 b[n][1] = bias[n][1] * ( b[n+1][0] * q[n+1][1] + b[n+1][1] *
` q[n+1][0] );
`366 b[n][0] = bias[n][0] * ( b[n+1][0] * q[n+1][0] + b[n+1][1] *
` q[n+1][1] );
`367 */
`368
`369 }
`370 if ( c->verbose > 3) {
`371 printf ( "state after forward pass:\n%10s %10s\n" , "f[n][0] " ,
` "f[n][1] " ) ;
`372 for ( n = 1 ; n <= N ; n ++ ) {
`373 printf ( "%-10.4g %-10.4g %2d\n" , f[n][0] , f[n][1] , n) ;
`374 }
`375 printf ( "\n" ) ;
`376 pause_for_return() ;
`377 }
`378 f0 = f[N][0] ; f1 = f[N][1] ;
`379 p0 = 1.0 ; p1 = 1.0 ;
`380 for ( n = N ; n >= 1 ; n-- ) {
`381 if ( n < N ) {
`382 q0 = q[n+1][0] ; q1 = q[n+1][1] ;
`383 p0 = q0 * b0 + q1 * b1 ; p1 = q1 * b0 + q0 * b1 ;
`384 }
`385 b0 = p0 * bias[n][0] ; b1 = p1 * bias[n][1] ; s=b0+b1 ;
`386 if ( s < BOT || ( s > TOP ) ) { s = 1.0/s ; b0 *= s ; b1 *= s ; }
`387
`388 /* we can now figure out the most probable state of trellis */
`
`Page: 7
`
`Hughes, Exh. 1019, p. 7
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`389 p_t0 = f0 * p0 ; p_t1 = f1 * p1 ;
`390 to[n] = ( p_t1 > p_t0 ) ;
`391
`392 /* now to the likelihoods */
`393 f0 = f[n-1][0] ; f1 = f[n-1][1] ;
`394
`395 /* b[n][0] = b0 ; b[n][1] = b1 ; */
`396 r[n][0] = f0*b0 + f1*b1 ;
`397 r[n][1] = f1*b0 + f0*b1 ; /* beware overflow here */
`398 }
`399
`400 if ( c->verbose > 3) {
`401 printf ( "likelihoods:\n%10s %10s\n" , "r[n][0] " , "r[n][1] " ) ;
`402 for ( n = 1 ; n <= N ; n ++ ) {
`403 printf ( "%-10.4g %-10.4g %2d\n" , r[n][0] , r[n][1] , n ) ;
`404 }
`405 printf ( "\n" ) ;
`406 pause_for_return() ;
`407 }
`408
`409 /* pseudo-vertical pass */
`410 viols = 0 ;
`411 for ( m = 1 ; m <= c->K ; m++ ) { /* run through source bits */
`412 umax = a->num_mlist[m] ; /* number of repetitions */
`413 mlist = a->mlist[m] ; /* list of bits */
`414 f0 = 1.0 ; f1 = 1.0 ; /* should insert prior here if
` source
`415 not dense */
`416 v0 = 0 ; v1 = 0 ; /* count votes for 0 and 1 */
`417 for ( u = 1 ; u <= umax ; u ++ ) {
`418 n = mlist[u] ;
`419 f0 *= r[n][0] ; f1 *= r[n][1] ; /* beware overflow here */
`420 if (to[n-1]^to[n]) v1 ++ ; else v0 ++ ;
`421 /* if ( r[n][0] > r[n][1] ) v0 ++ ; else v1 ++ ; */ /* detect
` inconsistent state */
`422 }
`423 if ( v0 * v1 ) { /* then this bit not yet consistent */
`424 viols ++ ;
`425 if ( c->verbose > 3) {
`426 printf ( "%dv" , m ) ;
`427 }
`428 }
`429 if ( f1 > f0 ) { so[m] = 1 ; }
`430 else { so[m] = 0 ; }
`431 if ( so[m] ^ ( v1 > v0 ) ) { /* disagreement between mp bit and
` votes */
`432 viols ++ ;
`433 if ( c->verbose > 5 ) {
`434 printf ( "so[m] = %d but v0 = %d and v1 = %d. (%f,%f)\n" ,
`435 so[m] , v0 , v1 , f0 , f1 ) ;
`436 }
`437 }
`438 for ( u = umax ; u >= 1 ; u -- ) {
`439 n = mlist[u] ;
`440 p0 = f0 / r[n][0] ; p1 = f1 / r[n][1] ; s = 1.0/(p0 + p1) ;
`441 q[n][0] = p0 * s ; q[n][1] = p1 * s ;
`442 if ( c->doclip ) { /* pragmatic clipping procedure to stop
` overflow */
`
`Page: 8
`
`Hughes, Exh. 1019, p. 8
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`443 if ( q[n][0] > c->clip ) {
`444 q[n][0] = c->clip ;
`445 q[n][1] = 1.0 - c->clip ;
`446 if ( !clipwarn ) { fprintf ( stderr , "c" ) ; fflush ( stderr ) ;
` clipwarn=1; }
`447 }
`448 else if ( q[n][1] > c->clip ) {
`449 q[n][1] = c->clip ;
`450 q[n][0] = 1.0 - c->clip ;
`451 if ( !clipwarn ) { fprintf ( stderr , "c" ) ; fflush ( stderr ) ;
` clipwarn=1; }
`452 }
`453 }
`454
`455 }
`456 }
`457 if ( c->verbose > 3) {
`458 printf ( "\nnew priors:\n%10s %10s\n" , "q[n][0] " , "q[n][1] " )
` ;
`459 for ( n = 1 ; n <= N ; n ++ ) {
`460 printf ( "%-10.4g %-10.4g %2d\n" , q[n][0] , q[n][1] , n ) ;
`461 }
`462 printf ( "viols = %d\n" , viols ) ;
`463 printf ( "\n" ) ;
`464 pause_for_return() ;
`465 }
`466 return ( viols ) ;
`467 }
`468
`469 static double bern ( int errors , int trials , double *p_point , double
`*p_upper
`470 , double *p_lower , double zscore ) {
`471 double factor ; double lowerbound = 1e-6 ;
`472 if ( errors == 0 ) {
`473 *p_point = lowerbound ;
`474 *p_lower = lowerbound ;
`475 *p_upper = 1.0 - exp( - zscore / (double)(trials) ) ;
`476 factor = 100.0 ;
`477 } else {
`478 *p_point = (double)(errors)/(double)(trials) ;
`479 /* factor corresponding z sd's in log p */
`480 factor = exp ( zscore * sqrt( (double)(trials - errors ) /
` ((double)(errors * trials ) ) )) ;
`481 *p_upper = *p_point * factor ;
`482 *p_lower = *p_point / factor ;
`483 }
`484 return factor ;
`485 }
`486
`487 static void finalline ( FILE *fp , RA_control *c , int totals )
`488 {
`489 double factor ;
`490 if ( totals ) {
`491 factor = bern ( c->block_errs , c->blocks , &(c->bep) , &(c->bep_u)
` , &(c->bep_l) , c->zscore ) ;
`492 c->bitep = (double) c->bit_errs / (double)( c->K * c->blocks ) ;
`493 c->bitep_u = c->bitep * factor ;
`494 c->bitep_l = c->bitep / factor ;
`
`Page: 9
`
`Hughes, Exh. 1019, p. 9
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`495 c->bep_det = (double) c->block_det / (double)( c->blocks ) ;
`496 c->bep_detlw = (double) c->block_detlw / (double)( c->blocks ) ;
`497 c->bep_undet = (double) c->block_undet / (double)( c->blocks ) ;
`498 }
`499
`500 if ( (!(c->outappend)) || (totals) ||
`501 ( c->pheading_period &&
`502 ! ( ( ++ c->pheading_count ) % c->pheading_period )) ) {
`503 if ( !totals ) {
`504 fprintf ( fp , "#es viols t_es ") ;
`505 }
`506 fprintf ( fp , "#blks es un det valid " ) ;
`507 fprintf ( fp , "Eb/No(dB) x xass ") ;
`508 if ( !totals ) {
`509 fprintf ( fp , " xact flipd ") ;
`510 }
`511 if ( totals ) {
`512 fprintf ( fp , " K N ") ;
`513 fprintf ( fp , "biters undet det ") ;
`514 fprintf ( fp , "block_ep u l ") ;
`515 fprintf ( fp , "bit_ep u l ") ;
`516 fprintf ( fp , "blep_det undet ") ;
`517 fprintf ( fp , "detlw blep_detlw " ) ;
`518 fprintf ( fp , "totloops lmax " ) ;
`519 }
`520 if ( !totals ) {
`521 fprintf ( fp , "loops ") ;
`522 }
`523 fprintf ( fp , "\n" ) ;
`524 }
`525
`526 if ( !totals ) {
`527 fprintf ( fp , "%4d " , c->count ) ; /* bits wrong this block */
`528 fprintf ( fp , "%4d " , c->viols ) ; /* number of detected problems
` this block */
`529 fprintf ( fp , "%4d " , c->tcount ) ; /* bits wrong this block */
`530 }
`531
`532 fprintf ( fp , "%4d " , c->blocks ) ; /* blocks transmitted */
`533 fprintf ( fp , "%3d " , c->block_errs ) ; /* block errors */
`534 fprintf ( fp , "%1d " , c->block_undet ) ; /* undetected */
`535 fprintf ( fp , "%3d " , c->block_det ) ; /* detected */
`536 fprintf ( fp , "%4d " , c->block_valid ) ; /* blocks that reached a
` valid decode state */
`
`537
`538 fprintf ( fp , "%7.4g " , c->ebno ) ; /* gaussian channel EbNo */
`539
`540 fprintf ( fp , "%7.4g " , c->gcx ) ; /* gaussian channel mean x */
`541 fprintf ( fp , "%7.4g " , c->gcxass ) ; /* assumed x */
`542 if ( !totals ) {
`543 fprintf ( fp , "%7.4g " , c->gcxact ) ; /* actual x this time
`544 (including fluctuation) */
`545 fprintf ( fp , "%4d " , c->flipped ) ; /* number of t bits "flipped"
` (a measure of noise level) */
`546 }
`547
`548 if ( totals ) {
`549
`
`Page: 10
`
`Hughes, Exh. 1019, p. 10
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`550 fprintf ( fp , "%4d " , c->K ) ; /* K source bits */
`551 fprintf ( fp , "%4d " , c->N ) ; /* N transmitted bits */
`552
`553 fprintf ( fp , "%4d " , c->bit_errs ) ; /* bit errors */
`554 fprintf ( fp , "%4d " , c->bit_undet ) ; /* undetected */
`555 fprintf ( fp , "%4d " , c->bit_det ) ; /* detected */
`556
`557 fprintf ( fp , "%8.4g " , c->bep ) ; /* block error probability */
`558 fprintf ( fp , "%7.2g " , c->bep_u ) ; /* block error probability
` (upper) */
`559 fprintf ( fp , "%7.2g " , c->bep_l ) ; /* block error probability
` (lower) */
`
`560
`561 fprintf ( fp , "%8.4g " , c->bitep ) ; /* bit error probability */
`562 fprintf ( fp , "%7.2g " , c->bitep_u ) ; /* error probability
` (upper) */
`563 fprintf ( fp , "%7.2g " , c->bitep_l ) ; /* error probability
` (lower) */
`
`564
`565 fprintf ( fp , "%9.3g " , c->bep_det ) ; /* detected block error
` probability */
`566 fprintf ( fp , "%9.3g " , c->bep_undet ) ; /* undetected block
` error probability */
`
`567
`568 fprintf ( fp , "%3d " , c->block_detlw ) ; /* detected errors
` involving low weight error */
`569 fprintf ( fp , "%9.3g " , c->bep_detlw ) ; /* detected block error
` probability involving low weight error words */
`570 fprintf ( fp , "%4d " , c->totloops ) ; /* totloops used in finding
` valid decodings */
`571 fprintf ( fp , "%4d " , c->loops ) ; /* when we would have stopped
` */
`572 } else {
`573 fprintf ( fp , "%4d " , c->loop ) ; /* iterations this time */
`574 }
`575
`576 fprintf ( fp , "\n" ) ;
`577
`578 }
`579
`580 static int make_sense ( RA_control *c )
`581 { /* correct silly RA_control parameters */
`582 int status = 0 ;
`583
`584 /* if no assumed noise level stated, set assumed equal to true */
`585 if ( c->gcxass < 0.0 ) { c->gcxass = c->gcx ; c->gcxact = c->gcx ; }
`586 if ( c->fnass < 0.0 ) { c->fnass = c->fn ; }
`587
`588 c->R = (double) (c->K) / (double) (c->N) ;
`589 c->ebno = 10.0 * log ( c->gcx * c->gcx / ( 2.0 * c->R ) ) / log (
` 10.0 ) ;
`
`590
`591 return status ;
`592 }
`593
`594 static int make_space ( RA_control *c ) {
`595 /* see also mnc_free */
`596 int status = 0 ;
`
`Page: 11
`
`Hughes, Exh. 1019, p. 11
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`597
`598 c->so = cvector ( 1 , c->K ) ;
`599 c->s = cvector ( 1 , c->K ) ;
`600 c->t = cvector ( 1 , c->N ) ;
`601 c->to = cvector ( 0 , c->N ) ;
`602 c->to[0] = 0 ;
`603
`604 c->histo = ivector ( 0 , c->loops+1 ) ; /* histogram of number of
` loops taken */
`605 set_ivector_const ( c->histo , 0 , c->loops+1 , 0 ) ;
`606
`607 /* vec->bias = dvector ( 1 , vec->N ) ; */
`608
`609 c->q = dmatrix ( 1 , c->N , 0 , 1 ) ;
`610 c->bias = dmatrix ( 1 , c->N , 0 , 1 ) ;
`611 c->f = dmatrix ( 0 , c->N , 0 , 1 ) ;
`612 c->r = c->f ; /* they can share space ! */
`613 c->f[0][0] = 1 ; c->f[0][1] = 0 ;
`614
`615 c->block_valid = 0 ;
`616 c->blocks=0 ;
`617 c->block_det = 0 ; c->block_undet = 0 ; c->block_errs = 0 ;
`618 c->bit_det = 0 ; c->bit_undet = 0 ; c->bit_errs = 0 ;
`619 c->block_detlw = 0 ;
`620 c->lowweight = c->N / 10 ;
`621 /* AD HOC definition of a 'low weight detected error event' */
`622
`623 return status ;
`624 }
`625
`626 int check_alist_MN ( alist_matrix *a , RA_control *c ) {
`627 int status = 0 ;
`628 if ( c->N != a->N
`629 || c->K != a->M ) {
`630 fprintf ( stderr , "setting N and K from alist: %d %d %d %d\n" ,
` c->N , a->N , c->K , a->M ) ;
`631 c->N = a->N ;
`632 c->K = a->M ;
`633 }
`634 return status ;
`635 }
`636
`637 static int score ( RA_control *c )
`638 /* compares so with s and works out quality of solution if wrong. */
`639 {
`640 int status = 0 ;
`641 int n , k , count = 0 ;
`642 FILE *fp ;
`643
`644 if ( c->verbose ) printf ( "Scoring answer --- \n" ) ;
`645 if ( c->verbose > 2) {
`646 printf ( "original source vector:\n" ) ;
`647 for ( k = 1 ; k <= c->K ; k ++ ) {
`648 if ( c->s[k] ) printf ("1 ") ; else printf ( "0 " );
`649 }
`650 printf ( "\n" ) ;
`651 }
`652 if ( c->verbose > 2) {
`
`Page: 12
`
`Hughes, Exh. 1019, p. 12
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`653 printf ( "reconstructed source vector:\n" ) ;
`654 for ( k = 1 ; k <= c->K ; k ++ ) {
`655 if ( c->so[k] ) printf ("1 ") ; else printf ( "0 " );
`656 }
`657 printf ( "\n" ) ;
`658 }
`659
`660 c->blocks ++ ;
`661 if ( !(c->viols) ) {
`662 c->block_valid ++ ;
`663 c->totloops += c->loop ;
`664 c->histo[c->loop] ++ ;
`665 }
`666 c->count = 0 ;
`667 for ( n = 1 ; n <= c->K ; n++ ) { /* this compares all bits */
`668 if ( c->so[n] != c->s[n] ) {
`669 c->count ++ ; /* count is the errors */
`670 }
`671 }
`672 c->tcount = 0 ;
`673 /* count the weight of the difference in codeword land */
`674 for ( n = 1 ; n <= c->N ; n++ ) {
`675 if ( c->to[n] != c->t[n] ) {
`676 c->tcount ++ ;
`677 }
`678 }
`679 if ( c->count > 0 ) {
`680 c->block_errs ++ ; c->failcount ++ ; /* why two registers!? */
`681 c->bit_errs += c->count ;
`682 if ( c->viols ) {
`683 c->block_det ++ ;
`684 c->bit_det += c->count ;
`685 /* check for near-codeword */
`686 if ( c->tcount < c->lowweight ) {
`687 c->block_detlw ++ ;
`688 }
`689 } else {
`690 printf ( "UNDETECTED ERROR\n" ) ;
`691 c->block_undet ++ ;
`692 c->bit_undet += c->count ;
`693 }
`694 if ( c->count && c->maxcount != 0 && c->error_log ) {
`695 fp = fopen ( c->error_logfile , "a" ) ;
`696 if ( !fp ) {
`697 fprintf ( stderr , " couldn't open logfile %s\n" , c->error_logfile )
` ;
`698 c->error_log = 0 ;
`699 } else {
`700 /* tell me more about this failure */ /* write to a file */
`701 finalline ( fp , c , 0 ) ;
`702 fprintf ( fp , "# " ) ;
`703 for ( n = 1 ; n <= c->K && ( ( count <= c->maxcount )
`704 || c->maxcount < 0 ) ; n++ ) {
`705 if ( c->so[n] != c->s[n] ) {
`706 count ++ ;
`707 if ( c->so[n] ) {
`708 fprintf ( fp , "+%4d " , n ) ; /* + indicates extra bit given in
` so */
`
`Page: 13
`
`Hughes, Exh. 1019, p. 13
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`709 } else {
`710 fprintf ( fp , "-%4d " , n ) ; /* - indicates bit lacking */
`711 }
`712 }
`713 }
`714 fprintf ( fp , "\n" ) ;
`715 fclose (fp ) ;
`716 }
`717 }
`718 } else {
`719 if ( c->verbose ) printf ( "SUCCESS\n" ) ;
`720 }
`721 return status ;
`722 }
`723
`724 static int t_to_b ( unsigned char *x , RA_control *c ) {
`725 int i , co = 0 ;
`726 int lo = 1 ; int hi = c->N ;
`727 double z , p ;
`728 double gcx = c->gcx ;
`729 double gcxass = c->gcxass ;
`730 double **b = c->bias ;
`731
`732 if ( c->gc ) {
`733 for ( i = lo ; i <= hi ; i ++ ) {
`734 /* make a random normal variate with s.d. 1.0 and mean +/- gcx */
`735 z = ( x[i] ? gcx : - gcx ) + rann() ;
`736 p = 1.0 / ( 1.0 + exp ( - 2.0 * z * gcxass ) ) ; /* likely to be
` large if x = 0 */
`737 /* was b[i] = 1.0 - p ; */
`738 b[i][1] = p ; b[i][0] = 1.0 - p ;
`739 if ( !(( p < 0.5 )^(x[i])) ) co++ ; /* counts as a flipped bit */
`740 }
`741
`742 } else { /* BSC
`743 co = random_cvector ( c->x , c->fn , lo, hi ) ;
`744 for ( i = lo ; i <= hi ; i ++) {
`745 if ( (c->x[i]) ^ (t[i]) ) {
`746 c->bias[i] = 1.0 - c->fnass ;
`747 } else {
`748 c->bias[i] = c->fnass ;
`749 }
`750 } */
`751 }
`752
`753 return co ;
`754 }
`755
`756 static void dc_defaults ( RA_control *c )
`757 {
`758 #include "RAdc_var_def.c"
`759 }
`760
`761 static int process_command ( int argc , char **argv , RA_control *c ) {
`762 int p_usage = 0 ;
`763 int status = 0 ;
`764 int cs , i ;
`765
`
`Page: 14
`
`Hughes, Exh. 1019, p. 14
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`766 if ( argc < 1 ) {
`767 p_usage = 1 ;
`768 status -- ;
`769 }
`770
`771 #define ERROR1 fprintf ( stderr , "arg to `%s' missing\n" , argv[i] ) ;
`\
`772 status --
`773 #define ERROR2 fprintf ( stderr , "args to `%s' missing\n" , argv[i] ) ;
`\
`774 status --
`775 for (i = 1 ; i < argc; i++) {
`776 cs = 1 ;
`777 if ( strcmp (argv[i], "-V") == 0 ) {
`778 c->verbose = 1;
`779 }
`780 else if ( strcmp (argv[i], "-VV") == 0 ) {
`781 c->verbose = 2;
`782 }
`783 #include "RAdc_var_clr.c"
`784 else {
`785 fprintf ( stderr , "arg `%s' not recognised\n" , argv[i] ) ;
`786 p_usage = 1 ;
`787 status -- ;
`788 }
`789 if ( cs == 0 ) {
`790 fprintf ( stderr , "arg at or before `%s' has incorrect format\n"
` ,
`791 argv[i] ) ;
`792 p_usage = 1 ;
`793 status -- ;
`794 }
`795 }
`796 if ( p_usage ) print_usage ( argv , stderr , c ) ;
`797 return ( status ) ;
`798 }
`799 #undef ERROR1
`800 #undef ERROR2
`801
`802 #define DNT fprintf( fp, "\n ")
`803 #define NLNE fprintf( fp, "\n")
`804
`805 static void print_usage ( char **argv , FILE * fp , RA_control *c )
`806 {
`807 fprintf( fp, "Usage: %s ",argv[0]);
`808 fprintf( fp, " [optional arguments]");
`809
`810 DNT; fprintf( fp, "-V | -VV (verbose or very
` verbose)");
`811 NLNE; fprintf( fp, " Data creation:" ) ;
`812 #include "RAdc_var_usg.c"
`813 pause_for_return();
`814 fprintf( fp, "\n");
`815 return ;
`816 }
`817
`818 #undef DNT
`819 #undef NLNE
`
`Page: 15
`
`Hughes, Exh. 1019, p. 15
`
`

`

`File: RA.c 10/15/1998, 5:32:14 PM
`
`820
`821 static void RA_free ( RA_control *vec ) {
`822
`823 free_cvector ( vec->so, 1 , vec->K ) ;
`824 free_cvector ( vec->t , 1 , vec->N ) ;
`825
`826 free_cvector ( vec->s , 1 , vec->K ) ;
`827
`828 free_dmatrix ( vec->bias , 1 , vec->N , 0 , 1 ) ;
`829 free_dmatrix ( vec->q , 1 , vec->N , 0 , 1 ) ;
`830 free_dmatrix ( vec->f , 0 , vec->N , 0 , 1 ) ;
`831
`832 }
`833
`834 double h2 ( double x ) {
`835 double tmp ;
`836 if ( x <= 0.0 || x>= 1.0 ) tmp = 0.0 ;
`837 tmp = x * log ( x ) + ( 1.0 - x ) * log ( 1.0 - x ) ;
`838
`839 return - tmp / log ( 2.0) ;
`840 }
`841
`842 /*
`843 <!-- hhmts start -->
`844 <!-- hhmts end -->
`845 */
`
`Page: 16
`
`Hughes, Exh. 1019, p. 16
`
`

`

`File: r.h 8/15/2002, 11:46:22 AM
`
`1 /* r.h header for assorted minor subroutines library
`release 1.1
`
` Copyright (c) 2002 David J.C. MacKay
`
`23
`
`45
`
` This library is free software; you can redistribute it and/or
`6 modify it under the terms of the GNU Lesser General Public
`7 License as published by the Free Software Foundation; either
`8 version 2.1 of the License, or (at your option) any later version.
`9
`10 This library is distributed in the hope that it will be useful,
`11 but WITHOUT ANY WARRANTY; without even the implied warranty of
`12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
`13 Lesser General Public License for more details.
`14
`15 You should have received a copy of the GNU Lesser General Public
`16 License along with this library; if not, write to the Free Software
`17 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
` USA
`
`18
`19 GNU licenses are here :
`20 http://www.gnu.org/licenses/licenses.html
`21
`22 Author contact details are here :
`23 http://www.inference.phy.cam.ac.uk/mackay/
` mackay@mrao.cam.ac.uk
`24 */
`25 /* ANSI Version 1 */
`26 /* 9 6 92 */
`27
`28 #include <stdio.h>
`29 #include <math.h>
`30 #include <stdlib.h>
`31 #include <string.h>
`32 #include <memory.h>
`33 #include "nrutil.h"
`34
`35 #ifndef INT
`36 #define INT( a ) ( ( ( a ) > 0.0 ) ? ( (int) ( (a)+0.5 ) ) : ( (int) (
`(a)-0.5 ) ) )
`37 #endif
`38 #ifndef MIN
`39 #define MIN( a , b ) ( ( a ) < ( b ) ? ( a ) : ( b ) )
`40 #endif
`41 #ifndef MAX
`42 #define MAX( a , b ) ( ( a ) > ( b ) ? ( a ) : ( b ) )
`43 #endif
`44 #ifndef PI
`45 #define PI 3.1415926535
`46 #endif
`47 #ifndef STPI
`48 #define STPI 2.50663 /* sqrt ( 2*PI ) */
`49 #endif
`50 #ifndef LTPI
`51 #define LTPI 1.83788 /* log ( 2*PI ) */
`52 #endif
`53 #ifndef FLUSH
`54 #define FLUSH fflush( stdout ) ;
`
`Page: 1
`
`Hughes, Exh. 1019, p. 17
`
`

`

`File: r.h 8/15/2002, 11:46:22 AM
`
`55 #endif
`56 #define fnewline fprintf( fp , "\n" )
`57 #define junkstring fscanf( fp , "%s" , junk )
`58 #define FNEWLINE fprintf( fp , "\n" )
`59 #define JUNKSTRING fscanf( fp , "%s" , junk )
`60
`61 /* list of functions that should be acce