#define TRUE (0<1)
#define FALSE (1<0)
#define MSWINDOWS FALSE
#if MSWINDOWS
  #include "stdafx.h"  /*needed if Microsoft Windows OS, unwanted if LINUX; will be more of that*/
#endif
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <ctype.h>
#include <sys/types.h>
#if MSWINDOWS
  #include <time.h>
  #include <sys/timeb.h>
  #include <process.h>
#else 
  #include <sys/time.h>
  #include <unistd.h>
#endif


/* #define NDEBUG    uncomment this line if want to turn off asserts for speed */
#define VERSION 3.24
#define VERSIONYEAR 2007
#define VERSIONMONTH 2

/**IEVS ***** Warren D. Smiths's 2007 infinitely extendible voting system comparison engine ****
Copyright (C): Warren D. Smith 2007.  Anybody can use this non-commercially 
provided they acknowledge me, notify me (and make available for my unrestricted use)
re any & all code improvements created or bugs found.  If you want commercial use,
then we'll negotiate.  

On Unix, Mac OSX, and Linux: Compilation:  gcc IEVS.c -Wall -lm -O6 -o IEVS
    for extra speed (but less safety) add the #define NDEBUG line at top.
On MSWINDOWS:   change #define MSWINDOWS to TRUE.  Should work, see notes below by Cary -
    thanks much to David Cary for the MSWINDOWS port!!  (Except that he introduced a bunch
    of bugs with careless typecasts, which hopefully I now have fixed...)
What election methods are available?:    fgrep EMETH IEVS.c
What utility-generators now available?:  fgrep UTGEN IEVS.c

This program should also work (with minor modifications) under microsoft windows OS.
David Cary successfully ported it and his notes on how he did it, are below.
As a side effect, Cary also noticed how to speed up Yee-picture-generation significantly
and fixed a bug in OutputCompressedBarray(). :)
If you do such porting, please report your experience, how you compiled it, 
etc.   warren.wds AT gmail.com.

Software Architecture:
I. voting methods.
II. voting strategies.
III. ignorance generators
IV. utility generators

Anybody can add new voting methods, new strategies, or new utility generators.
The idea is to build a "Chinese menu" system which can investigate A*B*C*D kinds of scenarios
BUT the effort to write it is only A+B+C+D.

The information-flow-direction is  IV-->III-->II-->I-->winners-->regrets.

I am initially writing this with not very many of each.  It is a "skeleton"
system which can be fleshed out later by adding more of each of I,II,III,IV.
Please send your contributed routines to warren.wds AT gmail.com .
You should be able to make an implementation of a new method pretty easily by
imitating some already-implemented methods similar to yours.
BUGGY: Rouse??, Copeland??
As-yet unimplemented voting methods include:
ICA, Maxtree, range+runoff(based on range ballots with ties discarded),
RRV+runoff, Sarvo-methods, Banks set, Boehm, Dodgson.
Asset & Candidate-withdrawal-option-IRV (both don't really fit into our model).
To add a new voting method you need to:
add your new EMETH subroutine, and you need to change 3 more lines:
  #define NumMethods PutCorrect#Methodshere
  and the case lines in GimmeWinner() and PrintMethName().
(If you add a new "core" method also must change NumCoreMethods and 
perhaps renumber the methods...)

Currently only rating-based, strict-ranking-based, and approval-based methods (& combinations)
are supported - equalities are not allowed in rankings.  ("3-slot" methods
also supported.)  To extend the code to allow
equalities, I suggest adding a boolean vector
with v[x]=TRUE meaning x is really tied with, not below,
the candidate immediately preceding.  That is a later goal.

HonestyStrat() now supports an arbitrary mixture of honest and (one kind of) strategic voters.  
More strategies can/will be added later.

3 utility generators now implemented (some parameterized), plus 
another "reality-based" utility generator based on
Tideman-et-al's real-world election collection.

Future plans/ideas:  translate from C into D?
(Templates & automatic array bounds checking would really help.)
auto-searcher for property violations?  
L1 utils should be skewed distn.
Allow a user-selected voting-method subset.
Targeted BR finder (new user interface to BR finder).
Other strat generators.
Use of 2-candidate elections is knd of silly.
Need to make better controls on summarizers so you can ask it to summarize only
SUBSETS of the data (such as, honest voters only, large #s of candidates only, etc)...
More realistic strategies for plur+runoff, MCA, IRNR, Benham2AppRunoff?
Build good hybrid voting methods?
Add feature to input votes from user and tally election.
Two-humped Gaussian distributions?
RandomWinner regret mean & variance should be computed more exactly.
DMC & the like - do the approvals and rankings correspond with strategy? Should they?
Summarizer: Also look at dominance relations, "approval".
***************************

Now notes by David Cary, 2007-02-14 re modifications to port IEVS 2.58 from LINUX to 
Windows XP using Microsoft Visual Studio .Net 2003.  

David Cary's Changes (not listing ones WDS did anyhow) include:
 1) Fix compile errors due to passing an int[] to a function that expects a uint[],
 2) Fix compile errors due to passing a uint[] to a function that expects an int[],
 3) Fix runtime error caused by opening the .bmp output file in text mode, to avoid
      LF bytes being written as CR LF combinations.
 4) Add an include for "stdafx.h".

 This program was successfully compiled and run as a C++ program after making the following
 adjustments to the default Visual Studio Project:
 1)  Increase the stack size to 4 MB (under Linker/System)
        (the default is 1 MB, but sizeof(edata) is about 2.2 MB and lives on the stack).

 All reported warnings were eliminated by:
 1) Updating the Visual Studio Project to not report 64-bit issues
      (under "C/C++" / General)
 2) Demoting warning messages #4018 & 4244 to warning level 4 by adding command line options
      /w44018 /w44244 under "C/C++" / Command Line.
      These messages warn about implicit int/uint conversions and conversions from __int64
      to unsigned int.
****************************/

#define until(x) while(!(x))
#define uint unsigned int
#define bool char
#define uint32 unsigned int
#define uint64 unsigned long long
#define uchar unsigned char
#define schar signed char
#define real double

#define PI 3.14159265358979323844

#define MaxNumCands 32
#define MaxNumVoters 2048
#define MaxNumIssues 8
#define MaxScenarios 17424
#define NumMethods 68  /*EMETH the # in this line needs to change if add new voting method*/
#define NumSlowMethods 8
#define NumFastMethods (NumMethods-NumSlowMethods)
#define NumUtilGens 11  /*UTGEN the # in this line needs to change if add new util gen*/
#define NumCoreMethods 12
#define CWSPEEDUP FALSE   /* TRUE causes it to be faster, but FALSE better for diagnosing bugs */

#define HTMLMODE 1  /*output adding "HTML table" formatting commands*/
#define TEXMODE 2   /*output adding "TeX table" formatting commands*/
#define NORMALIZEREGRETS 4   /*Output "normalized" Regrets so RandomWinner=1, SociallyBest=0. */
#define SORTMODE 8           /*Output with voting methods in sorted order by increasing regrets.*/
#define SHENTRUPVSR 16   /*Output "Shentrup normalized" Regrets so RandomWinner=0%, SociallyBest=100%. */
#define OMITERRORBARS 32   
#define VBCONDMODE 64  /*vote-based condorcet winner agreement counts; versus true utility based CWs*/
#define DOAGREETABLES 128
#define ALLMETHS 256
#define TOP10METHS 512
uint BROutputMode=0;

/******************** GENERAL PURPOSE routines having nothing to do with voting: ******/

int EulerPrimePoly(int n){ return(n*n+n+41); } /*prime for n=0..39*/
const int Pow2Primes[] = {2, 3, 7, 13, 31, 61, 127, 251, 509, 1021, 2039, 4093, 8191, 16381};
/** Greatest prime <=2^n. **/

/******** Fns to deal with sets represented by bit-vectors in 1 machine word: ******/
bool SingletonSet(uint x){ return ((x&(x-1))==0); } /*assumes non-empty*/
bool StrictSuperset(uint x, uint y){ return ((x&y)==y && x!=y); }
bool EmptySet(uint x){ return (x==0); }
int CardinalitySet(uint x){ 
   int ct=0;
   while( x ){ ct++; x &= x-1; }
   return(ct);
}

/****** convenient fns: *******/
real SquareReal(real x){  return x*x; }
int SquareInt(int x){  return x*x; }
int AbsInt(int x){  if(x<0) return -x; else return x; }
real PosReal(real x){ if(x>0.0) return x; else return 0.0; }
uint PosInt(int x){ if(x>0) return x; else return 0; }
int SignInt(int x){ if(x>0) return 1; else if(x<0) return -1; else return 0; }
int SignReal(real x){ if(x>0) return 1; else if(x<0) return -1; else return 0; }
int MaxInt(int a, int b){ return (((a)>(b)) ? (a):(b)); }
real MaxReal(real a, real b){ return (((a)>(b)) ? (a):(b)); }
int MinInt(int a, int b){ return (((a)<(b)) ? (a):(b)); }
real MinReal(real a, real b){ return (((a)<(b)) ? (a):(b)); }

uint GCD(uint a, uint b){ /*Euclid alg*/
  if(a>b){ a %= b; if(a==0) return b; }
  for(;;){
    b %= a;
    if(b==0) return a;
    a %= b;
    if(a==0) return b;
  }
}

uint LCMfact[32];

void BuildLCMfact(){
  uint j,x;
  printf("\nComputing sequence LCM(1,2,...N) for N=1..22:\n");
  LCMfact[0] = 1;
   for(j=1; j<23; j++){ /*LCMfact[23]=5354228880 won't fit in 32 bit word*/
     LCMfact[j] = LCMfact[j-1];
     x = LCMfact[j]%j;
     if(x) LCMfact[j] *= j/GCD(x,j);
   }
   printf("LCMfact[%d]=%u\n", j, LCMfact[j]);
}


/***************************************** other "Artin primes" are 
3, 5, 11, 13, 19, 29, 37, 53, 59, 61, 67, 83, 101, 107, 131, 139, 149, 163, 173, 179, 181, 
197, 211, 227, 269, 293, 317, 347, 349, 373, 379, 389, 419, 421, 443, 461, 467, 491, 509, 
523, 541, 547, 557, 563, 587, 613, 619, 653, 659, 661, 677, 701, 709, 757, 773, 787, 797
i.e. these are primes which have 2 as a primitive root.
This routine finds the greatest Artin prime P with P<=x.  (Not intended to be fast.)
  *******************************************/
uint ARTINPRIME;
uint FindArtinPrime(uint x){
  uint j,p,k;
  p = x;
  if(p%2==0) p--;  /* make p be odd */
  for(; p>2; p -= 2){
    for(j=4, k=3; j!=2; j=(j*2)%p, k++){
      if(k >= p){ return(p); }
    }
  }
  printf("FindArtinPrime failed - terminating\n");
  exit(EXIT_FAILURE);
}

void PrintNSpaces(int N){
  int i;
  for(i=0; i<N; i++) putchar(' ');
}
/****** convenient constants: *******/
#define BIGINT 0x7FFFFFFF
#define MAXUINT ((uint)((255<<48)|(255<<40)|(255<<32)|(255<<24)|(255<<16)|(255<<8)|(255)))
/* defn works on 8,16,32, and 64-bit machines */


uint32 BLC32x[60];  /* 32*60=1920 bits of state. Must be nonzero mod P. */
int BLC32NumLeft;

/********************************************************
Warren D. Smith 2001
**********************************************************
Linear congruential pseudo-random number generator mod P,
where P is the enormous prime (578 base-10 digits long; 
60 words long if each word is 32 bits)
  P = [2^(48*32) - 1] * 2^(12*32) + 1.
This prime can yield PRNGs suitable for use on 
computers with w-bit words, w=8,16,32,48,64,96,128.
The following code is intended for w=32.
The fact that 2^(60*32) = 2^(12*32) - 1 (mod P)
makes modular arithmetic mod P easy, and is the
reason this particular P was chosen.
The period of our generator is P-1.
***********************************************************
Although it is usually easy to detect the nonrandomness of
linear congruential generators because they generate d-tuples
lying on a lattice, in the present case the "grain size" of the
lattice is invisibly small (far less than a least significant bit),
for 32-bit words, if 1<=d<=180. If you want 64-bit words, then need
to concatenate two 32-bit words, and then grain size invisibly small
if 1<=d<=90. These bounds are conservative; in fact I suspect
one is ok, even for 64-bit words, even in up to 1000 dimensions.
***********************************************************
Some even more enormous primes which we have not used are:
[2^{59*32} - 1] * 2^{8 *32} + 1,
[2^{63*32} - 1] * 2^{24*32} + 1,
[2^{69*32} - 1] * 2^{14*32} + 1,
[2^{95*32} - 1] * 2^{67*32} + 1,
[2^{99*32} - 1] * 2^{35*32} + 1;
these would also be suitable for (8,16,32)-bit computers,
and the second of them would also be good for (48,96)-bit computers.
Unfortunately the factorization of P-1 is not known for the last 3 
I've listed here, preventing you from being certain you've found a
primitive root mod that P. A still more enormous prime is
  [2^4253 - 1] * 2^4580 + 1    [2659 digits long!]
(where note 2^4253-1 is also prime so that factorizing P-1 is
trivial) but doing arithmetic mod this P is (although still fairly
easy) less pleasant because bit-shifting is required.
*************************************************************/
uint32 BigLinCong32(){
   uint32 y[120];
   int i;
   uint64 u;

   if(BLC32NumLeft==0){
      /* Need to refill BLC32x[0..59] with 60 new random numbers: */

 /****************************************************************
 * If BLC32x[0..59] is the digits, LS..MS, of a number in base 2^w,
 * then the following code fragment puts A times that number 
 * in y[0..119].  Here
 *  A = 1284507170 * 2^(w*3) + 847441413 * 2^(w*44) + 650134147 * 2^(w*59)
 * is a "good" primitive root mod P, if w=32.
 *****************************************************************/
#define lohalf(x) (uint32)(x)
#define A1 (uint64)1284507170
#define A2 (uint64)847441413
#define A3 (uint64)650134147
      for(i=0; i<3; i++){
	 y[i] = 0;
      }
      u=0;
      for(/*i=3*/; i<44; i++){
	 u += A1 * BLC32x[i-3];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=44*/; i<59; i++){
	 u += A1 * BLC32x[i-3]; 
	 u += A2 * BLC32x[i-44];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=59*/; i<60+3; i++){
	 u += A1 * BLC32x[i-3]; 
	 u += A2 * BLC32x[i-44]; 
	 u += A3 * BLC32x[i-59];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=60+3*/; i<60+44; i++){
	 u += A2 * BLC32x[i-44]; 
	 u += A3 * BLC32x[i-59];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=60+44*/; i<60+59; i++){
	 u += A3 * BLC32x[i-59];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      /*i=60+59=119*/
      y[i] = lohalf(u);
#undef A1 
#undef A2 
#undef A3 
 /*************************************************************
 * If y[0..119] is the digits, LS..MS, of a number in base 2^w,
 * then the following code fragment replaces that number with
 * its remainder mod P in y[0..59]  (conceivably the result will
 * be >P, but this does not matter; it will never be >=2^(w*60)).
 **************************************************************/
      u=1; /*borrow*/
#define AllF 0xffffffff
      /* Step 1: y[0..72] = y[0..59] + y[60..119]shift12 - y[60..119]: */
      for(i=0; i<12; i++){
	 u += y[i]; 
	 u += (uint64)~y[60+i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=12*/; i<60; i++){
	 u += y[i]; 
	 u += y[48+i]; 
	 u += (uint64)~y[60+i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=60*/; i<72; i++){
	 u += AllF; 
	 u += y[48+i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      assert(u>0);
      y[72] = (uint32)(u-1); /*unborrow*/

      /*  Step 2: y[0..60] = y[0..59] + y[60..72]shift12  - y[60..72]: */
      u=1; /*borrow*/
      for(i=0; i<12; i++){
	 u += y[i]; 
	 u += (uint64)~y[60+i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      /*i=12*/
      u += y[i] + y[48+i]; 
      u += (uint64)~y[60+i];
      y[i] = lohalf(u);
      u = u>>32;
      i++;
      for(/*i=13*/; i<25; i++){
	 u += AllF; 
	 u += y[i]; 
	 u += y[48+i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      for(/*i=25*/; i<60; i++){
	 u += AllF; 
	 u += y[i];
	 y[i] = lohalf(u);
	 u = u>>32;
      }
      /*i=60*/
      assert(u>0);
      y[i] = (uint32)(u-1); /*unborrow*/

     /*It is rare that any iterations of this loop are needed:*/
      while(y[60]!=0){ 
         printf("rare loop\n");
	 /*Step 3+:  y[0..60] = y[0..59] + y[60]shift12 - y[60]:*/
	 u=1; /*borrow*/
	 u += y[0]; 
	 u += (uint64)~y[60];
	 y[0] = lohalf(u);
	 u = u>>32;
	 for(i=1; i<12; i++){
	    u += AllF; 
	    u += y[i];
	    y[i] = lohalf(u);
	    u = u>>32;
	 }
	 /*i=12*/
	 u += AllF; 
	 u += y[i]; 
	 u += y[60];
	 y[i] = lohalf(u);
	 u = u>>32;
	 i++;
	 for(/*i=13*/; i<60; i++){
	    u += AllF; 
	    u += y[i];
	    y[i] = lohalf(u);
	    u = u>>32;
	 }
	 /*i=60*/
	 assert(u>0);
	 y[i] = (uint32)(u-1); /*unborrow*/
      }
#undef AllF 
#undef lohalf

      /* Copy y[0..59] into BLC32x[0..59]: */
      for(i=0; i<60; i++){ 
	 BLC32x[i] = y[i]; 
      }
      /*printf("%u\n", BLC32x[44]);*/
      BLC32NumLeft=60;
   }
   /* (Else) We have random numbers left, so return one: */
   BLC32NumLeft--;
   return BLC32x[BLC32NumLeft];
}

void testbiglincong(){
   int i;
   /* lexically minimal permissible seed: */
   for(i=0; i<60; i++){ BLC32x[i]=0; }
   BLC32x[0] = 1; 
   BLC32NumLeft = 0;

   for(i=0; i<599; i++){ BigLinCong32(); }
   printf("%u %u %u %u %u\n", 
       BigLinCong32(),BigLinCong32(),
       BigLinCong32(),BigLinCong32(),BigLinCong32());

   for(i=0; i<12; i++){ 
      printf("%8x %8x %8x %8x %8x\n",
       BigLinCong32(),BigLinCong32(),
       BigLinCong32(),BigLinCong32(),BigLinCong32());
   }
}

/********************************
MAPLE script to check it works:

w := 32;
A := 1284507170 * 2^(w*3) + 847441413 * 2^(w*44) + 650134147 * 2^(w*59);
P := (2^(48*32) - 1) * 2^(12*32) + 1;
for i from 1 to 11 do
   print( floor((A &^ i mod P) / 2^(44*w)) mod (2^w) );
od;
ap := A &^ 11 mod P;
ap mod (2^w);
quit;
********************
MAPLE output:
847441413
4038410930
102374915
470100141
3896743552
243412576
1911259311
1640083353
4014446395
2679952782
4087228849
and
2475219032

C output:
847441413
4038410930
102374915
470100141
3896743552
243412576
oops
2990118053
2614294516
3539391572
1589778147
1758817216
2847725135 1364008005 3563722108 2020382641 1091616930
*************************/

uint32 RandUint(){ /* returns random uint32 */
  return BigLinCong32();
}

real Rand01(){ /* returns random uniform in [0,1] */
  return ((BigLinCong32()+0.5)/(1.0 + MAXUINT) + BigLinCong32())/(1.0 + MAXUINT);
}

void InitRand(uint seed){ /* initializes the randgen */
   int i;
   int seed_sec=0, processId=0;
   uint seed_usec=0;
#if     MSWINDOWS
   tm* locTime;
   _timeb currTime;
   time_t now;
#else
   struct timeval tv;
#endif
   if(seed==0){
     printf("using time of day and PID to generate seed");
#if MSWINDOWS
     now = time(NULL);
     locTime = localtime(&now);
     _ftime(&currTime);
     seed_usec = currTime.millitm*1001;
     seed_sec = locTime->tm_sec + 60*(locTime->tm_min + 60*locTime->tm_hour);
     processId = _getpid();
#else
     gettimeofday(&tv,0);
     seed_sec = tv.tv_sec;
     seed_usec = tv.tv_usec;
     processId = getpid();
#endif
     seed = 1000000*(uint)seed_sec + (uint)seed_usec + 
       (((uint)processId)<<20) + (((uint)processId)<<10);
     printf("=%u\n", seed);
   }
   for(i=0; i<60; i++){ BLC32x[i]=0; }
   BLC32x[0] = seed; 
   BLC32NumLeft = 0;
   for(i=0; i<599; i++){ BigLinCong32(); }
   printf("Random generator initialized with seed=%u:\n", seed);
   for(i=0; i<7; i++){ 
     printf("%.6f ", Rand01());
   }
   printf("\n");
}

void TestRand01(){ 
  int i,y, ct[10];
   real x,s,mx,mn,v;
   s=0.0; v=0.0;
   mn = HUGE;
   mx = -HUGE;
   for(i=0; i<10; i++) ct[i]=0; 
   printf("Performing 100000 randgen calls to test that randgen[0,1] behaving ok:\n");
   for(i=0; i<100000; i++){ 
     x = Rand01();
     s += x;
     if(mx<x) mx=x;
     if(x<mn) mn=x;
     v += SquareReal(x-0.5);
     y = (int)(x*10.0);
     if(x>=0 && y<10) ct[y]++;
   }
   printf("mean=%g(should be 1/2)  min=%f  max=%g   variance=%g(should be 1/12=%g)\n", 
      s/100000.0, mn, mx, v/100000.0, 1/12.0);
   printf("counts in 10 bins 0-0.1, 0.1-0.2, etc: ");
   for(i=0; i<10; i++) printf(" %d", ct[i]);
   printf("\n");
   fflush(stdout);
}

bool RandBool(){ /* returns random boolean */
  if( Rand01() > 0.5 ) return TRUE;
  return FALSE;
}

real RandExpl(){ /* returns standard exponential (density=exp(-x) for x>0) random deviate */
  real x;
  do{ 
    x = Rand01();
  }while( x==0.0 );
  return -log(x);
}

void TestRandExpl(){ 
  int i, y, ct[10];
  real x,s,mx,mn,v;
  s=0.0; v=0.0;
  mn = HUGE;
  mx = -HUGE;
  for(i=0; i<10; i++) ct[i]=0; 
  printf("Performing 100000 randgen calls to test that expl randgen behaving ok:\n");
  for(i=0; i<100000; i++){ 
     x = RandExpl();
     s += x;
     if(mx<x) mx=x;
     if(x<mn) mn=x;
     v += x*x;
     y = (int)(x*10);
     if(x>=0 && y<10) ct[y]++;
  }
  printf("mean=%g(should be 1)  min=%f  max=%g   variance=%g(should be 2?)\n", 
	 s/100000.0, mn, mx, v/100000.0);
  printf("counts in 10 bins 0-0.1, 0.1-0.2, etc: ");
  for(i=0; i<10; i++) printf(" %d", ct[i]);
  printf("\n");
  fflush(stdout);
}

real RandNormal(){ /* returns standard Normal (gaussian variance=1 mean=0) deviate */
  real w, x1;
  static real x2;
  static bool ready = FALSE;
  if(ready){
    ready = FALSE;
    return x2;
  }
  do{
    x1 = 2*Rand01() - 1.0;
    x2 = 2*Rand01() - 1.0;
    w = x1*x1 + x2*x2;
  }while ( w > 1.0 || w==0.0 );
  w = sqrt( (-2.0*log(w)) / w );
  x1 *= w;
  x2 *= w;  /* Now x1 and x2 are two indep normals (Box-Muller polar method) */
  ready = TRUE;
  return x1;
}

void TestNormalRand(){ 
  int i, y, ct[10];
  real x,s,mx,mn,v;
  s=0.0; v=0.0;
  mn = HUGE;
  mx = -HUGE;
  for(i=0; i<10; i++) ct[i]=0; 
  printf("Performing 100000 randgen calls to test that normal randgen behaving ok:\n");
  for(i=0; i<100000; i++){ 
     x = RandNormal();
     s += x;
     if(mx<x) mx=x;
     if(x<mn) mn=x;
     v += x*x;
     y = (int)(x*10);
     if(x>=0 && y<10) ct[y]++;
  }
  printf("mean=%g(should be 0)  min=%f  max=%g   variance=%g(should be 1)\n", 
	 s/100000.0, mn, mx, v/100000.0);
  printf("counts in 10 bins 0-0.1, 0.1-0.2, etc: ");
  for(i=0; i<10; i++) printf(" %d", ct[i]);
  printf("\n");
  fflush(stdout);
}

/* If a 2D normal [x & y coords of which are i.i.d. standard normals]
 * is known to lie on a ray thru center; this returns the distance along the ray. 
 * Equivalently, sin(2*pi*T)*R and cos(2*pi*T)*R are iid standard normals if T is
 * random uniform in [0,1] and if R is the independent output of RandRadialNormal().  ****/
real RandRadialNormal(){ 
  real w;
  do{ w = Rand01(); }while(w==0.0);
  w = sqrt( -2.0*log(w) );
  return w;
}

void TestRadialNormalRand(){ 
  int i, y, ct[10];
  real x,s,mx,mn,v;
  s=0.0; v=0.0;
  mn = HUGE;
  mx = -HUGE;
  for(i=0; i<10; i++) ct[i]=0; 
  printf("Performing 100000 randgen calls to test that radial normal randgen behaving ok:\n");
  for(i=0; i<100000; i++){ 
     x = RandRadialNormal();
     s += x;
     if(mx<x) mx=x;
     if(x<mn) mn=x;
     v += x*x;
     y = (int)(x*10);
     if(x>=0 && y<10) ct[y]++;
  }
  printf("mean=%g(should be ~1.25)  min=%f  max=%g   meansq=%g(should be 2)\n", 
	 s/100000.0, mn, mx, v/100000.0);
  printf("counts in 10 bins 0-0.1, 0.1-0.2, etc: ");
  for(i=0; i<10; i++) printf(" %d", ct[i]);
  printf("\n");
  fflush(stdout);
}

void TestRadialNormalRand2(){ 
  int i, y, ct[10];
  real x,s,mx,mn,v,w;
  s=0.0; v=0.0;
  mn = HUGE;
  mx = -HUGE;
  for(i=0; i<10; i++) ct[i]=0; 
  printf("Performing 100000 randgen calls to test that normal randgen behaving ok radially:\n");
  for(i=0; i<100000; i++){ 
     x = RandNormal();
     w = RandNormal();
     x = sqrt(x*x+w*w);
     s += x;
     if(mx<x) mx=x;
     if(x<mn) mn=x;
     v += x*x;
     y = (int)(x*10);
     if(x>=0 && y<10) ct[y]++;
  }
  printf("mean=%g(should be ~1.25)  min=%f  max=%g   meansq=%g(should be 2)\n", 
	 s/100000.0, mn, mx, v/100000.0);
  printf("counts in 10 bins 0-0.1, 0.1-0.2, etc: ");
  for(i=0; i<10; i++) printf(" %d", ct[i]);
  printf("\n");
  fflush(stdout);
}

#define RECIPRTPI 0.564189583547756286948079451560772585844050   /* 1/sqrt(pi) */
real RandSkew(){ /* returns mean=0 skewed deviate; the fact the mean is really 0, has been 
		  * confirmed by Monte Carlo to +-0.0001 */
  real x,y;
  x = RandNormal();
  y = RandNormal();
  if(x<y) x=y;
  return( 1.21129 * (x-RECIPRTPI) );  /*1.21129 chosen so that variance is 1*/
}

void GenRandNormalArr(int N, real Arr[]){ /* returns Arr[0..N-1] of standard normal randoms */
  int i;
  for(i=0; i<N; i++){
    Arr[i] = RandNormal();
  }
}

void GenRandSkewArr(int N, real Arr[]){ /* returns Arr[0..N-1] of skew randoms */
  int i;
  for(i=0; i<N; i++){ Arr[i] = RandSkew(); }
}

uint RandInt(uint N){ /* returns random integer in {0,1, ..., N-1} */
    return (int)(Rand01()*N);
}

real wkc(int a, int b){
  int si;
  assert(b>0);
  si = (a%2) ? 1 : -1; /*si=pseudorandom sign*/
  return  si*cos( (2*a+1)*PI / (4*b) );
}

real wks(int a, int b){
  assert(b>0);
  return sin( (2*a+1)*PI / (4*b) );
}

/* This is a goofy probability density on N-vectors.  Its mean is the zero vector.
 * The variance of each component is 1.  Its support set is all of N-space.
 * It is NOT centrosymmetric.    It has 3 or fewer modes.
 * Centrosymmetric densities (e.g. normal) have the disadvantage that they cause every Condorcet
 * voting method to be identical (and cause Condorcet winner always to exist) 
 * with distance-based utilities...  which is not realistic.  Hence this.
 **********/
#define RT85 1.264911064067351732799557417773  /* sqrt(8/5) */
#define RT25 0.632455532033675866399778708886  /* sqrt(2/5) */
void GenRandWackyArr(int N, real Arr[]){ /* returns Arr[0..N-1] of skew randoms */
  int k, which;
  which = (int)RandInt(3);
  switch(which){
  case(0) : for(k=0; k<N; k++){ Arr[k] = RandSkew(); } break;
  case(1) : for(k=0; k<N; k++){ Arr[k] =  wkc(k, N) + RT85*wks(k, N) * RandNormal(); } break;
  case(2) : for(k=0; k<N; k++){ Arr[k] = -wkc(k, N) + RT25*wks(k, N) * RandNormal(); } break;
  default : 
    printf("impossible case in GenRandWackyArr\n"); 
    fflush(stdout); 
    exit(EXIT_FAILURE);
  }
}

void OldSortedUpRand01Arr(int N, real Arr[]){ /*makes Arr[0..N-1] of uniform01 randoms, SORTED increasing*/
  int i;
  real k,x,p;
  k = 1.0;
  for(i=N; i>0; i--){
    p = 1.0/i;
    do{
      x = Rand01();
      x = pow(x, p);
    }while(x<=0.0 || x>=1.0);
    k = k*x;
    /*printf("i=%d k=%f\n", i,k);*/
    assert(k<=1.0);
    assert(0.0<=k);
    Arr[i-1] = k;
  }
}

/* J.Bentley & J.Saxe: Generating Sorted Lists of Random Numbers, 
 * ACM Trans. Math'l. Software 6,3 (1980) 359-364. */
void SortedUpRand01Arr(int N, real Arr[]){ /*makes Arr[0..N-1] of uniform01 randoms, SORTED increasing*/
  int k;
  real s;
  s = 0.0;
  for(k=0; k<N; k++){
    s -= log(Rand01());
    Arr[k] = s;
  }
  s -= log(Rand01());
  for(k=0; k<N; k++){
    Arr[k] /= s;
  }
}


void TestsOfRand(){
  TestRand01(); TestNormalRand(); TestRandExpl(); 
  TestRadialNormalRand();
  TestRadialNormalRand2();
}

bool IsPerm( uint N, uint Perm[] ){ /* true if is a perm of [0..N-1] */
  int i;
  int ct[MaxNumCands];
  assert(N<MaxNumCands);
  for(i=0; i<(int)N; i++){ ct[i] = 0; }
  for(i=0; i<(int)N; i++){ ct[ Perm[i] ]++; }
  for(i=0; i<(int)N; i++){ if(ct[i]!=1) return FALSE; }
  return TRUE;
}

void MakeIdentityPerm( uint N, uint Perm[] ){ 
  int i;
  for(i=0; i<(int)N; i++){ Perm[i] = i; }
}

void RandomlyPermute( uint N, uint RandPerm[] ){ /* randomly permutes RandPerm[0..N-1] */
  int i, j, t;
  for(i=N-1; i>0; i--){
    j = (int)RandInt((uint)i);
    t = RandPerm[j];
    RandPerm[j] = RandPerm[i];
    RandPerm[i] = t;
  }
  assert(IsPerm(N,RandPerm));
}

/******* vector handling: **********/
void CopyIntArray(uint N, int src[], int dest[] ){ 
  int i;
  for(i=N-1; i>=0; i--) dest[i] = src[i];
}

void CopyRealArray(uint N, real src[], real dest[] ){ 
  int i;
  for(i=N-1; i>=0; i--) dest[i] = src[i];
}

void FillBoolArray(uint N, bool Arr[], bool Filler ){ /* sets Arr[0..N-1] = Filler */
  int i;
  for(i=N-1; i>=0; i--) Arr[i] = Filler;
}

void FillRealArray(uint N, real Arr[], real Filler ){ /* sets Arr[0..N-1] = Filler */
  int i;
  for(i=N-1; i>=0; i--) Arr[i] = Filler;
}

void FillIntArray(uint N, int Arr[], int Filler ){ /* sets Arr[0..N-1] = Filler */
  int i;
  for(i=N-1; i>=0; i--) Arr[i] = Filler;
}

void ZeroIntArray(uint N, int Arr[] ){ /* sets Arr[0..N-1] = 0 */
   FillIntArray( N, Arr, 0 );
}

void ZeroRealArray(uint N, real Arr[] ){ /* sets Arr[0..N-1] = 0 */
   FillRealArray( N, Arr, 0.0 );
}

/* Assumes RandPerm[0..N-1] contains perm. Returns index of random min entry of Arr[0..N-1]. */
int ArgMinIntArr(uint N, int Arr[], int RandPerm[] ){ 
  int minc, i, r, winner;
  winner = -1;
  minc = BIGINT;
  RandomlyPermute( N, (uint*)RandPerm );
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]<minc){
      minc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  assert( Arr[winner] <= Arr[0] );
  assert( Arr[winner] <= Arr[N-1] );
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains perm. Returns index of random max entry of Arr[0..N-1]. */
int ArgMaxIntArr(uint N, int Arr[], int RandPerm[] ){
  int maxc, i, r, winner;
  winner = -1;
  maxc = -BIGINT;
  RandomlyPermute( N, (uint*)RandPerm );
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r] > maxc){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains perm. Returns index of random max entry of Arr[0..N-1]. */
int ArgMaxUIntArr(uint N, uint Arr[], int RandPerm[] ){
  uint maxc;
  int i, r, winner;
  winner = -1;
  maxc = 0;
  RandomlyPermute( N, (uint*)RandPerm );
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]>=maxc){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains perm. Returns index of random min entry of Arr[0..N-1]. */
int ArgMinRealArr(uint N, real Arr[], int RandPerm[] ){ 
  real minc;
  int i,r,winner; 
  winner = -1;
  minc = HUGE;
  RandomlyPermute( N, (uint*)RandPerm );
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]<minc){
      minc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains perm. Returns index of random max entry of Arr[0..N-1]. */
int ArgMaxRealArr(uint N, real Arr[], int RandPerm[] ){
  real maxc;
  int i,r,winner;
  winner = -1;
  maxc = -HUGE;
  RandomlyPermute( N, (uint*)RandPerm );
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r] > maxc){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains random perm and MinInd is index for Arr[] yielding min value.
 * Returns index of second-min. */
int Arg2MinIntArr(uint N, int Arr[], int RandPerm[], int MinInd ){
  int minc, i, r, winner;
  winner = -1;
  minc = BIGINT;
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]<minc && r!=MinInd){
      minc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains random perm and MaxInd is index for Arr[] yielding max.
 * Returns index of second-max. */
int Arg2MaxIntArr(uint N, int Arr[], int RandPerm[], int MaxInd ){
  int maxc, i, r, winner;
  winner = -1;
  maxc = -BIGINT;
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]>maxc && r!=MaxInd){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains random perm and MaxInd is index for Arr[] yielding max.
 * Returns index of second-max. */
int Arg2MaxUIntArr(uint N, uint Arr[], int RandPerm[], int MaxInd ){
  uint maxc;
  int i, r, winner;
  winner = -1;
  maxc = 0;
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]>=maxc && r!=MaxInd){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

/* Assumes RandPerm[0..N-1] contains random perm and MaxInd is index for Arr[] yielding max.
 * Returns index of second-max. */
int Arg2MaxRealArr(uint N, real Arr[], int RandPerm[], int MaxInd ){
  real maxc;
  int i, r, winner;
  winner = -1;
  maxc = -HUGE;
  for(i=0; i<(int)N; i++){
    r = RandPerm[i];
    if(Arr[r]>maxc && r!=MaxInd){
      maxc=Arr[r]; 
      winner=r; 
    }
  }
  assert(winner>=0);
  return(winner);
}

void ScaleRealVec(uint N, real a[], real scalefac){
  int i;
  for(i=0; i<(int)N; i++) a[i] *= scalefac;
}

void ScaleIntVec(uint N, int a[], int scalefac){
  int i;
  for(i=0; i<(int)N; i++) a[i] *= scalefac;
}

/*Donald Knuth's "The Art of Computer Programming, Volume 2: Seminumerical Algorithms", section 4.2.2
describes how to compute mean and standard deviation using a recurrence relation, like this:
M(1) = x(1),   M(k) = M(k-1) + (x(k) - M(k-1)) / k
S(1) = 0,      S(k) = S(k-1) + (x(k) - M(k-1)) * (x(k) - M(k))
for 2 <= k <= n, then
sigma = sqrt(S(n) / (n - 1))
Attributes this method to B.P. Welford, Technometrics 4 (1962) 419-420.
*******/
inline void WelfordUpdateMeanSD(real NewDatum, int *Count, real *M, real *S){
  real OldMean;
  OldMean = *M;
  (*Count)++;
  *M += (NewDatum - OldMean)/(*Count);
  *S += (NewDatum - OldMean)*(NewDatum - *M);
  return;
}

real DistanceSquared(uint N, real a[], real b[] ){
 real d = 0.0;
  int i;
  for(i=0; i<(int)N; i++) d += SquareReal( a[i]-b[i] );
  return d;
}

real L1Distance(uint N, real a[], real b[] ){
  real d = 0.0;
  int i;
  for(i=0; i<(int)N; i++) d += fabs( a[i]-b[i] );
  return d;
}

real LpDistanceSquared(uint N, real a[], real b[], real Lp ){
 real d = 0.0;
  int i;
  assert(Lp >= 1.0);
  if(Lp==1.0) return SquareReal(L1Distance(N,a,b));
  if(Lp==2.0) return DistanceSquared(N,a,b);
  for(i=0; i<(int)N; i++) d += pow( fabs( a[i]-b[i] ), Lp );
  return pow( d, 2.0/Lp );
}

real LpDistance(uint N, real a[], real b[], real Lp ){
 real d = 0.0;
  int i;
  assert(Lp >= 1.0);
  if(Lp==1.0) return L1Distance(N,a,b);
  if(Lp==2.0) return sqrt( DistanceSquared(N,a,b) );
  for(i=0; i<(int)N; i++) d += pow( fabs( a[i]-b[i] ), Lp );
  return pow( d, 1.0/Lp );
}

real  SumRealArray( uint N, real a[] ){
  real s = 0.0;
  int i;
  for(i=0; i<(int)N; i++) s += a[i];
  return s;
}

real DotProd(uint N, real a[], real b[] ){
 real d = 0.0;
  int i;
  for(i=0; i<(int)N; i++) d += a[i]*b[i];
  return d;
}

/******* sorting: **********/
int SortedReal( uint N, real Arr[] ){ /* +1 if Arr[0..N-1] increasing, -1 if decreasing, 2 if unsorted, 0 if all-eq */
  int i,s;
  s = SignReal( Arr[N-1]-Arr[0] );
  for(i=1; i<(int)N; i++){
    if( s*SignReal( Arr[i] - Arr[i-1] ) < 0 ) return(2);
  }
  return s;
}

int SortedInt( uint N, int Arr[] ){ /* +1 if Arr[0..N-1] increasing, -1 if decreasing, 2 if unsorted, 0 if all-eq */
  int i,s;
  s = SignInt( Arr[N-1]-Arr[0] );
  for(i=1; i<(int)N; i++){
    if( s*SignInt( Arr[i] - Arr[i-1] ) < 0 ) return(2);
  }
  return s;
}

int SortedRealKey( uint N, int Arr[], real Key[] ){
  int i,s;
  s = SignReal( Key[Arr[N-1]]-Key[Arr[0]] );
  for(i=1; i<(int)N; i++){
    if( s*SignReal( Key[Arr[i]] - Key[Arr[i-1]] ) < 0 ) return(2);
  }
  return s;
}

int SortedIntKey( uint N, int Arr[], int Key[] ){
  int i,s;
  s = SignInt( Key[Arr[N-1]]-Key[Arr[0]] );
  for(i=1; i<(int)N; i++){
    if( s*SignInt( Key[Arr[i]] - Key[Arr[i-1]] ) < 0 ) return(2);
  }
  return s;
}

const int ShellIncs[] = {1750, 701, 301, 132, 57, 23, 10, 4, 1, 0};
/* Marcin Ciura: Best Increments for the Average Case of Shellsort,
   13th International Symposium on Fundamentals of Computation Theory,
   Riga, Latvia, 22-24 August 2001;
   Springer Lecture Notes in Computer Science #2138, pp.106-117. 
Here 1750 is unsure and how the sequence continues past 1750 is unknown.
 ***/

void RealShellSortUp( uint N, real Arr[] ){ /* Sorts Arr[0..N-1] into increasing order */
  int h,i,j,k;
  real x;
  for(k=0; (h=ShellIncs[k])>0; k++){
      for(i=h; i<(int)N; i++){
	x=Arr[i];
	for(j=i-h; j>=0 && Arr[j]>x; j -= h){ Arr[j+h]=Arr[j]; }
	Arr[j+h]=x;
      }
  }
  assert((SortedReal(N,Arr)&(~1))==0);
}

void IntShellSortUp( uint N, int Arr[] ){ /* Sorts Arr[0..N-1] into increasing order */
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
  
      for(i=h; i<(int)N; i++){
	x=Arr[i];
	for(j=i-h; j>=0 && Arr[j]>x; j -= h){ Arr[j+h]=Arr[j]; }
	Arr[j+h]=x;
      }
  }
  assert((SortedInt(N,Arr)&(~1))==0);
}

void RealShellSortDown( uint N, real Arr[] ){ /* Sorts Arr[0..N-1] into decreasing order */
  int h,i,j,k;
  real x;
  for(k=0; (h=ShellIncs[k])>0; k++){
      for(i=h; i<(int)N; i++){
	x=Arr[i];
	for(j=i-h; j>=0 && Arr[j]<x; j -= h){ Arr[j+h]=Arr[j]; }
	Arr[j+h]=x;
      }
  }
  assert(SortedReal(N,Arr)<=0);
}

void IntShellSortDown( uint N, int Arr[] ){ /* Sorts Arr[0..N-1] into decreasing order */
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
      for(i=h; i<(int)N; i++){
	x=Arr[i];
	for(j=i-h; j>=0 && Arr[j]<x; j -= h){ Arr[j+h]=Arr[j]; }
	Arr[j+h]=x;
      }
  }
  assert(SortedInt(N,Arr)<=0);
}

/* Rearranges Perm[0..N-1] so Key[Perm[0..N-1]] is in increasing order: */
void IntPermShellSortUp(  uint N, int Perm[], int Key[] ){
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
    for(i=h; i<(int)N; i++){
      x=Perm[i];
      for(j=i-h; j>=0 && Key[Perm[j]]>Key[x]; j -= h){ Perm[j+h]=Perm[j]; }
      Perm[j+h]=x;
    }
  }
  assert((SortedIntKey(N,Perm,Key)&(~1))==0);
}

/* Rearranges Perm[0..N-1] so Key[Perm[0..N-1]] is in increasing order: */
void RealPermShellSortUp(  uint N, int Perm[], real Key[] ){
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
    for(i=h; i<(int)N; i++){
      x=Perm[i];
      for(j=i-h; j>=0 && Key[Perm[j]]>Key[x]; j -= h){ Perm[j+h]=Perm[j]; }
      Perm[j+h]=x;
    }
  }
  assert((SortedRealKey(N,Perm,Key)&(~1))==0);
}

/* Rearranges Perm[0..N-1] so Key[Perm[0..N-1]] is in decreasing order.  Key[] not altered: */
void IntPermShellSortDown(  uint N, int Perm[], int Key[] ){
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
    for(i=h; i<(int)N; i++){
      x=Perm[i];
      for(j=i-h; j>=0 && Key[Perm[j]]<Key[x]; j -= h){ Perm[j+h]=Perm[j]; }
      Perm[j+h]=x;
    }
  }
  assert(SortedIntKey(N,Perm,Key)<=0);
}

/* Rearranges Perm[0..N-1] so Key[Perm[0..N-1]] is in decreasing order.  Key[] not altered: */
void ScharPermShellSortDown(  uint N, schar Perm[], uint Key[] ){
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
    for(i=h; i<(int)N; i++){
      x=Perm[i];
      for(j=i-h; j>=0 && Key[Perm[j]]<Key[x]; j -= h){ Perm[j+h]=Perm[j]; }
      Perm[j+h]=x;
    }
  }
}

/* Rearranges Perm[0..N-1] so Key[Perm[0..N-1]] is in decreasing order: */
void RealPermShellSortDown(  uint N, int Perm[], real Key[] ){
  int h,i,j,k;
  int x;
  for(k=0; (h=ShellIncs[k])>0; k++){
    for(i=h; i<(int)N; i++){
      x=Perm[i];
      for(j=i-h; j>=0 && Key[Perm[j]]<Key[x]; j -= h){ Perm[j+h]=Perm[j]; }
      Perm[j+h]=x;
    }
  }
  assert(SortedRealKey(N,Perm,Key)<=0);
}

bool SelectedRightInt( uint L, uint R, uint K, int A[] ){
  uint i;
  for(i=L; i<K; i++){ if( A[i]>A[K] ) return FALSE; }
  for(i=R; i>K; i--){ if( A[K]>A[i] ) return FALSE; }
  return TRUE;
}

bool SelectedRightReal( uint L, uint R, uint K, real A[] ){
  uint i;
  for(i=L; i<K; i++){ if( A[i]>A[K] ) return FALSE; }
  for(i=R; i>K; i--){ if( A[K]>A[i] ) return FALSE; }
  return TRUE;
}

/* Rearranges A[L..R] so that A[i] <= A[K] <= A[j]  if L<=i<K<j<=R.
 * Then returns A[K].  If R-L+1=N then
 * expected runtime for randomly ordered A[] is 1.5*N compares 
 * asymptotically to find median (and less for other K, finds max in only 1.0*N).
 * In particular if L=0, R=N-1 and K=(N-1)/2 with N odd, then returns median.
 * If N is even, then use K=N/2 to find the upper bimedian, and then
 * the max of A[0..K-1] is the lower bimedian. 
 *   Note: the magic constants 601, 0.5, 0.5, 20, 5, and 5 in the below should
 * be tuned to optimize speed.  They are probably improveable.
 ***/
int FloydRivestSelectInt( uint L, uint R, uint K, int A[] ){
  int I, J, N, S, SD, LL, RR;
  real Z;
  int T,W;
  while( R > L ){
    N = R - L + 1;
    if( N > 601 ){ /* big enough so worth doing FR-type subsampling to find strong splitter */
      I = K - L + 1; 
      Z = log((real)(N));
      S = (int)(0.5 * exp(Z * 2.0/3));
      SD = (int)(0.5 * sqrt(Z * S * (N - S)/N) * SignInt(2*I - N));
      LL = MaxInt(L, K - ((I * S) / N) + SD);
      RR = MinInt(R, K + (((N - I) * S) / N) + SD);
      /* Recursively select inside small subsample to find an element A[K] usually very 
       * near the desired one: */
      FloydRivestSelectInt( LL, RR, K, A );  
    }else if( N > 20 && (int)(5*(K-L)) > N && (int)(5*(R-K)) > N){ 
        /* not big enough to be worth evaluating those expensive logs etc; 
         * but big enough so random splitter is poor; so use median-of-3 */
      I = K-1; if(A[K] < A[I]){ W = A[I]; A[I] = A[K]; A[K] = W; }
      I = K+1; if(A[I] < A[K]){ W = A[I]; A[I] = A[K]; A[K] = W; }
      I = K-1; if(A[K] < A[I]){ W = A[I]; A[I] = A[K]; A[K] = W; }
    } /* otherwise using random splitter (i.e. current value of A[K]) */
    /* now use A[K] to split A[L..R] into two parts... */
    T = A[K]; I = L; J = R;
    W = A[L]; A[L] = A[K]; A[K] = W;
    if( A[R] > T ){ W = A[R]; A[R] = A[L]; A[L] = W; }
    while( I < J ){
      W = A[I]; A[I] = A[J]; A[J] = W; I++; J--;
      while(A[I] < T){ I++; }
      while(A[J] > T){ J--; }
    }
    if( A[L] == T ){ W = A[L]; A[L] = A[J]; A[J] = W; }
    else{ J++; W = A[J]; A[J] = A[R]; A[R] = W; }
    if( J <= (int)K ){ L = J + 1; }
    if( (int)K <= J ){ R = J - 1; }
    /* Now continue on using contracted [L..R] interval... */
  }
  return( A[K] );
}

/* returns twice the median of A[0..N-1] or sum of the bimedians if N even */
int TwiceMedianInt( uint N, int A[] ){ 
  int M,T;
  int i;
  assert(N>0);
  M = FloydRivestSelectInt( 0, N-1, N/2, A );
  assert( SelectedRightInt( 0, N-1, N/2, A ) );
  if((N&1)==0){ /*N is even*/
    T = A[N/2 - 1];
    for(i=N/2 - 2; i>=0; i--){
      if(A[i] > T) T=A[i];
    }
  }else{ T=M; }
  return T+M;
}

/* Rearranges A[L..R] so that A[i] <= A[K] <= A[j]  if L<=i<K<j<=R.
 * Then returns A[K].  If R-L+1=N then
 * expected runtime for randomly ordered A[] is 1.5*N compares 
 * asymptotically to find median (and less for other K, finds max in only 1.0*N).
 * In particular if L=0, R=N-1 and K=(N-1)/2 with N odd, then returns median.
 * If N is even, then use K=N/2 to find the upper bimedian, and then
 * the max of A[0..K-1] is the lower bimedian. 
 ***/
real FloydRivestSelectReal( uint L, uint R, uint K, real A[] ){
  int I, J, N, S, SD, LL, RR;
  real Z;
  real T,W;
  while( R > L ){
    N = R - L + 1;
    if( N > 601 ){ /* big enough so worth doing FR-type subsampling to find strong splitter */
      I = K - L + 1; 
      Z = log((real)(N));
      S = (int)(0.5 * exp(Z * 2.0/3));
      SD = (int)(0.5 * sqrt(Z * S * (N - S)/N) * SignInt(2*I - N));
      LL = MaxInt(L, K - ((I * S) / N) + SD);
      RR = MinInt(R, K + (((N - I) * S) / N) + SD);
      /* Recursively select inside small subsample to find an element A[K] usually very 
       * near the desired one: */
      FloydRivestSelectReal( LL, RR, K, A );  
    }else if( N > 20 && (int)(5*(K-L)) > N && (int)(5*(R-K)) > N){ 
        /* not big enough to be worth evaluating those expensive logs etc; 
         * but big enough so random splitter is poor; so use median-of-3 */
      I = K-1; if(A[K] < A[I]){ W = A[I]; A[I] = A[K]; A[K] = W; }
      I = K+1; if(A[I] < A[K]){ W = A[I]; A[I] = A[K]; A[K] = W; }
      I = K-1; if(A[K] < A[I]){ W = A[I]; A[I] = A[K]; A[K] = W; }
    } /* otherwise using random splitter (i.e. current value of A[K]) */
    /* now use A[K] to split A[L..R] into two parts... */
    T = A[K]; I = L; J = R;
    W = A[L]; A[L] = A[K]; A[K] = W;
    if( A[R] > T ){ W = A[R]; A[R] = A[L]; A[L] = W; }
    while( I < J ){
      W = A[I]; A[I] = A[J]; A[J] = W; I++; J--;
      while(A[I] < T){ I++; }
      while(A[J] > T){ J--; }
    }
    if( A[L] == T ){ W = A[L]; A[L] = A[J]; A[J] = W; }
    else{ J++; W = A[J]; A[J] = A[R]; A[R] = W; }
    if( J <= (int)K ){ L = J + 1; }
    if( (int)K <= J ){ R = J - 1; }
    /* Now continue on using contracted [L..R] interval... */
  }
  return( A[K] );
}

/* returns twice the median of A[0..N-1] or sum of the bimedians if N even */
real TwiceMedianReal( uint N, real A[] ){ 
  real M,T;
  int i;
  assert(N>0);
  M = FloydRivestSelectReal( 0, N-1, N/2, A );
  assert( SelectedRightReal( 0, N-1, N/2, A ) );
  if((N&1)==0){ /*N is even*/
    T = A[N/2 - 1];
    for(i=N/2 - 2; i>=0; i--){
      if(A[i] > T) T=A[i];
    }
  }else{ T=M; }
  return T+M;
}

/*************************** VOTING METHODS:  ********
all are subroutines with a common format - here is the format (which is all subsumed in
the convenient data structure "edata"):
input:
uint NumVoters;
uint NumCands;

uint TopDownPrefs[NumVoters*NumCands];   
Entry x*NumCands+y says the candidate who is the (y-1)th choice of voter x, x=0..NumVoters-1.
Candidates are numbered 0..NumCands-1.

uint CandRankings[NumVoters*NumCands];
Entry x*NumCands+y says the ranking (0=top) of the yth candidate (y=0..NumCands-1)
according to voter x, x=0..NumVoters-1.

real Score[NumVoters*NumCands];
Entry x*NumCands+y says the score (a floating point real from 0 to 1)
of the yth candidate (y=0..NumCands-1) according to voter x, x=0..NumVoters-1.

bool Approve[NumVoters*NumCands];
Entry x*NumCands+y says the Approval (a boolean)
of the yth candidate (y=0..NumCands-1) according to voter x, x=0..NumVoters-1.

bool Approve2[NumVoters*NumCands];
Entry x*NumCands+y says the second-level Approval (a boolean)
of the yth candidate (y=0..NumCands-1) according to voter x, x=0..NumVoters-1.
Used in MCA.  A higher level of approval.

real PerceivedUtility[NumVoters*NumCands];
Entry x*NumCands+y says the utility (a floating point real)
of the yth candidate (y=0..NumCands-1) according to voter x, x=0..NumVoters-1.
Note, this is "honest", i.e. undistorted by strategy, 
although it IS distorted by ignorance.

real Utility[MaxNumCands*MaxNumVoters];
Entry x*NumCands+y says the utility (a floating point real)
of the yth candidate (y=0..NumCands-1) according to voter x, x=0..NumVoters-1.
Note, this is "honest", i.e. undistorted by strategy, AND undistorted by ignorance.

semi-input:  (used as input by some election methods, but computed from the preceding
genuine input by BuildDefeatsMatrix)
int DefeatsMatrix[MaxNumCands*MaxNumCands]; 
DefeatsMatrix[i*NumCands+j] is #voters who ranked i above j  (always nonnegative)

real Armytage[MaxNumCands*MaxNumCands]; 
Armytage[i*NumCands+j] is sum of rating[i]-rating[j] for voters who rated i above j  (always nonnegative)

int MarginsMatrix[MaxNumCands*MaxNumCands]; 
MarginsMatrix[i*NumCands+j] is margin of i over j  (either sign)

real MargArmy[MaxNumCands*MaxNumCands]; 
MargArmy[i*NumCands+j] is Armytage[i*NumCands+j]-Armytage[j*NumCands+i]


output (the returned value):
int Winner;       Number in {0,1,...,NumCands-1}  saying who won. 
                  (Negative number indicates error condition.)
                  Note: tie breaking is random.  A global array RandPerm[0..NumCands-1] is
                  assumed available to help with that task.

CORE METHODS: Plurality, Borda, IRV, SociallyBest, SociallyWorst, RandomWinner, Approval, Range, 
  SmithSet, SchwartzSet, Antiplurality
are "core" voting methods which must always be run (and first), otherwise
some other methods may not work.

side-effects:  can alter these global variables (list follows):
***************************/

 int PlurWinner;
 int AntiPlurWinner;
 int PSecond;
 int RSecond;
 int ASecond;
 int ApprovalWinner;
 int IRVwinner;
 int SmithWinner;
 int RandWinner;
 int SchwartzWinner;
 int RangeWinner;
 int BordaWinner;
 int WorstWinner;
 int BestWinner;
 int RandomUncoveredMemb;
uint RangeGranul;
 int IRVTopLim = BIGINT;
 int CondorcetWinner;   /*negative if does not exist*/
 int TrueCW;   /*cond winner based on undistorted true utilities; negative if does not exist*/
 int CopelandWinner;
 int CopeWinOnlyWinner;
 int SmithIRVwinner;
uint PlurVoteCount[MaxNumCands];
uint AntiPlurVoteCount[MaxNumCands];
uint DabaghVoteCount[MaxNumCands];
 int VFAVoteCount[MaxNumCands];
 int RdVoteCount[MaxNumCands];
 int FavListNext[MaxNumVoters];
 int HeadFav[MaxNumCands];
 int WinCount[MaxNumCands];
 int DrawCt[MaxNumCands];
 int CopeScore[MaxNumCands];
 int LossCount[MaxNumCands];
 int SimmVotesAgainst[MaxNumCands];
 int BeatPathStrength[MaxNumCands*MaxNumCands];
real ArmyBPS[MaxNumCands*MaxNumCands];
uint ApprovalVoteCount[MaxNumCands];
 int UncAAOF[MaxNumCands];
uint MCAVoteCount[MaxNumCands];
real RangeVoteCount[MaxNumCands];
real SumNormedRating[MaxNumCands];
uint RangeNVoteCount[MaxNumCands];
real CCumVoteCount[MaxNumCands];
real MedianRating[MaxNumCands];
real CScoreVec[MaxNumVoters];
 int MedianRank[MaxNumCands];
 int CRankVec[MaxNumVoters];
uint BordaVoteCount[MaxNumCands];
 int NansonVoteCount[MaxNumCands];
uint NauruVoteCount[MaxNumCands];
uint HeismanVoteCount[MaxNumCands];
uint BaseballVoteCount[MaxNumCands];
uint SumOfDefeatMargins[MaxNumCands];
 int WorstDefeatMargin[MaxNumCands];
 int RayDefeatMargin[MaxNumCands];
 int RayBeater[MaxNumCands];
 int ARVictMargin[MaxNumCands];
 int ARchump[MaxNumCands];
real UtilitySum[MaxNumCands];
real UtilityRootSum[MaxNumCands];
uint RandCandPerm[MaxNumCands]; /* should initially contain 0..NumCands-1 */
bool MajApproved[MaxNumCands];
bool Eliminated[MaxNumCands];
bool MDdisquald[MaxNumCands];
bool BSSmithMembs[MaxNumCands];
bool SmithMembs[MaxNumCands];
bool UncoveredSt[MaxNumCands];
bool SchwartzMembs[MaxNumCands];
 int IFav[MaxNumVoters];
uint NauruWt[MaxNumCands];
uint BaseballWt[MaxNumCands];
 int PairApproval[MaxNumCands*MaxNumCands]; 
real SinkRat[MaxNumCands], SinkRow[MaxNumCands], SinkCol[MaxNumCands];
real SinkMat[MaxNumCands*MaxNumCands];
bool CoverMatrix[MaxNumCands*MaxNumCands];
real EigVec[MaxNumCands], EigVec2[MaxNumCands];
bool Rmark[MaxNumCands];
bool rRmark[MaxNumCands];
 int summ[MaxNumCands];
 int Tpath[MaxNumCands*MaxNumCands];
 int Hpotpar[MaxNumCands];
 int Hpar[MaxNumCands];
 int Hroot[MaxNumCands];
uint WoodHashCount[3*MaxNumCands*MaxNumVoters], WoodHashSet[3*MaxNumCands*MaxNumVoters]; 
uint WoodSetPerm[3*MaxNumCands*MaxNumVoters];
uint IBeat[MaxNumCands];

void InitCoreElState(){ /*can use these flags to tell if Plurality() etc have been run*/
  PlurWinner = -1;
  BordaWinner = -1;
  ApprovalWinner = -1;
  RangeWinner = -1;
  IRVwinner = -1;
  BestWinner = -1;
  WorstWinner = -1;
  AntiPlurWinner = -1;
  PSecond = -1;
  RandWinner = -1;
  SmithWinner = -1;
  SchwartzWinner = -1;
  CopeWinOnlyWinner = -1;
  SmithIRVwinner = -1;
  RandomUncoveredMemb = -1;
  IRVTopLim = BIGINT;
}

typedef struct dum1 {
  uint NumVoters;
  uint NumCands;
  uint TopDownPrefs[MaxNumCands*MaxNumVoters];
  uint CandRankings[MaxNumCands*MaxNumVoters];
  real Score[MaxNumCands*MaxNumVoters];
  bool Approve[MaxNumCands*MaxNumVoters];
  bool Approve2[MaxNumCands*MaxNumVoters];
  real PerceivedUtility[MaxNumCands*MaxNumVoters];
  real Utility[MaxNumCands*MaxNumVoters];
  int TrueDefMatrix[MaxNumCands*MaxNumCands];
  int DefeatsMatrix[MaxNumCands*MaxNumCands];
  int MarginsMatrix[MaxNumCands*MaxNumCands];
  real Armytage[MaxNumCands*MaxNumCands];
  int ArmyDef[MaxNumCands*MaxNumCands];
  real MargArmy[MaxNumCands*MaxNumCands];
} edata;

void PrintEdata(FILE *F, edata *E){ /* prints out the edata */
  int v,j;
  fprintf(F, "NumVoters=%d  NumCands=%d\n", E->NumVoters, E->NumCands);
  for(v=0; v<(int)E->NumVoters; v++){
    fprintf(F, "Voter %2d:\n",v);
    fprintf(F, "Utility: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%6.3f", E->Utility[v*E->NumCands + j]);  }
    fprintf(F, "\n");
    fprintf(F, "PercUti: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%6.3f", E->PerceivedUtility[v*E->NumCands + j]);  }
    fprintf(F, "\n");
    fprintf(F, "RangeScore: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%6.3f", E->Score[v*E->NumCands + j]);  }
    fprintf(F, "\n");
    fprintf(F, "CandRank: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%2d", E->CandRankings[v*E->NumCands + j]);  }
    fprintf(F, "\n");
    fprintf(F, "TopDown: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%2d", E->TopDownPrefs[v*E->NumCands + j]);  }
    fprintf(F, "\n");
    fprintf(F, "Approve: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%2d", E->Approve[v*E->NumCands + j] ? 1:0);  }
    fprintf(F, "\n");
    fprintf(F, "Approve2: ");
    for(j=0; j < E->NumCands; j++){  fprintf(F, "%2d", E->Approve2[v*E->NumCands + j] ? 1:0);  }
    fprintf(F, "\n");
  }
  /*???more?*/
  fflush(F);
}

#define EMETH int  /* allows fgrep EMETH IEVS.c to find out what Election methods now available */

void BuildDefeatsMatrix(edata *E
){ /* initializes  E->DefeatsMatrix[], E->MarginsMatrix[], RandCandPerm[], NauruWt[], WinCount[], DrawCt[], CondorcetWinner, CopeWinOnlyWinner, TrueCW */
  int k,i,j,y;
  uint x;
  real t;
  bool CondWin, TrueCondWin;
  MakeIdentityPerm( E->NumCands, RandCandPerm );
  RandomlyPermute( E->NumCands, RandCandPerm );

  assert(E->NumCands <= MaxNumCands);
  x = LCMfact[E->NumCands];
  for(j=E->NumCands -1; j>=0; j--){  NauruWt[j] = x / (j+1); }
  for(j=MinInt(E->NumCands -1,9); j>=0; j--){  BaseballWt[j] = 10-j; }
  BaseballWt[0] = 14;
  ZeroIntArray(  SquareInt(E->NumCands),  E->DefeatsMatrix );
  ZeroRealArray( SquareInt(E->NumCands),  E->Armytage );
  ZeroIntArray(  SquareInt(E->NumCands),  E->ArmyDef );
  ZeroIntArray(  SquareInt(E->NumCands),  E->TrueDefMatrix );
  for(k=0; k<(int)E->NumVoters; k++){
    x = k*E->NumCands;
    for(i=E->NumCands -1; i>=0; i--){
      for(j=i-1; j>=0; j--){
	y = E->CandRankings[x+j]; y -= E->CandRankings[x+i];
	if( y>0 ){
	  E->DefeatsMatrix[i*E->NumCands + j]++; /*i preferred above j*/
	}else{
	  E->DefeatsMatrix[j*E->NumCands + i]++;
	}
	t = E->Score[x+i] - E->Score[x+j];
	if(t > 0.0){
	  E->Armytage[i*E->NumCands + j] += t; /*i preferred above j*/
	  E->ArmyDef[i*E->NumCands + j] ++;
	}else{
	  E->Armytage[j*E->NumCands + i] -= t;
	  E->ArmyDef[j*E->NumCands + i] ++;
	}
	t = E->Utility[x+i] - E->Utility[x+j];
	if(t > 0.0){
	  E->TrueDefMatrix[i*E->NumCands + j] ++;
	}else{
	  E->TrueDefMatrix[j*E->NumCands + i] ++;
	}
      }
    }
  }
  CondorcetWinner = -1;
  TrueCW = -1;
  for(i=E->NumCands -1; i>=0; i--){
    WinCount[i] = DrawCt[i] = 0;
    CondWin = TRUE;
    TrueCondWin = TRUE;
    for(j=(int)E->NumCands -1; j>=0; j--){
      assert( E->DefeatsMatrix[i*E->NumCands+j] <= (int)E->NumVoters );
      assert( E->DefeatsMatrix[i*E->NumCands+j] >= 0 );
      assert( E->DefeatsMatrix[i*E->NumCands+j] + E->DefeatsMatrix[j*E->NumCands+i] <= (int)E->NumVoters );
      y = E->ArmyDef[i*E->NumCands+j]; y -= E->ArmyDef[j*E->NumCands+i];
      if(y>0) E->MargArmy[i*E->NumCands + j] = E->Armytage[i*E->NumCands+j];
      else /*y<=0*/ E->MargArmy[i*E->NumCands + j] = 0;
      y = E->DefeatsMatrix[i*E->NumCands+j]; y -= E->DefeatsMatrix[j*E->NumCands+i];
      E->MarginsMatrix[i*E->NumCands + j] = y;
      assert(i!=j || E->MarginsMatrix[i*E->NumCands + j] == 0);
      if(y > 0) WinCount[i]++;
      if(y==0) DrawCt[i]++;
      if(y<=0 && j!=i){ CondWin = FALSE; } /* if beaten or tied, not a CondorcetWinner by this defn */
      y = E->TrueDefMatrix[i*E->NumCands+j]; y -= E->TrueDefMatrix[j*E->NumCands+i];      
      if(y<=0 && j!=i){ TrueCondWin = FALSE; }
    }
    if( CondWin ) CondorcetWinner = i;  /* i beats all opponents (unique Condorcet winner) */
    if( TrueCondWin ) TrueCW = i;  /* i beats all opponents (unique Condorcet winner) */
  }
  /* find who-you-beat sets: */
  for(i=E->NumCands -1; i>=0; i--){
    x = 0;
    for(j=(int)E->NumCands -1; j>=0; j--){
      if( E->MarginsMatrix[i*E->NumCands + j] > 0 ){
	x |= (1<<j);
      }
    }
    IBeat[i] = x;
  }
  CopeWinOnlyWinner = ArgMaxIntArr( E->NumCands, WinCount, (int*)RandCandPerm );
  ZeroIntArray( SquareInt(E->NumCands), PairApproval );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
       for(k=E->NumCands -1; k>j; k--){
	 y = ((E->Approve[x+j] && E->Approve[x+k]) ? 1:0);
	 PairApproval[j*E->NumCands +k] += y; /* count of voters who approve of both j and k */
	 PairApproval[k*E->NumCands +j] += y;
       }
     }
  }		
}

EMETH SociallyBest(edata *E  /* greatest utility-sum winner */
){ /* side effects: UtilitySum[], BestWinner */
  uint x;
  int i,j;
  ZeroRealArray( E->NumCands, UtilitySum );
  for(i=0; i<(int)E->NumVoters; i++){
    x = i*E->NumCands;
    for(j=E->NumCands -1; j>=0; j--){
      UtilitySum[j] += E->Utility[x+j];
    }
  }		
  BestWinner = ArgMaxRealArr( E->NumCands, UtilitySum, (int*)RandCandPerm );
  for(j=E->NumCands -1; j>=0; j--){ assert( UtilitySum[BestWinner] >= UtilitySum[j] ); }
  return BestWinner;
}

EMETH SociallyWorst(edata *E   /* least utility-sum winner */
){ /* side effects: UtilitySum[], WorstWinner */
  if(BestWinner<0) SociallyBest(E);
  WorstWinner = ArgMinRealArr( E->NumCands, UtilitySum, (int*)RandCandPerm );
  return WorstWinner;
}

EMETH RandomWinner(edata *E){ return (int)RandInt( E->NumCands ); }

EMETH RandomBallot(edata *E){ /*honest top choice of a random voter is elected. Strategyproof.*/
  int winner;
  winner = ArgMaxRealArr( E->NumCands, E->PerceivedUtility + RandInt(E->NumVoters)*E->NumCands, (int*)RandCandPerm );
  return winner;
}

EMETH Hay(edata *E /*Strategyproof. Prob of election proportional to sum of sqrts of [normalized] utilities*/
){ /* side effects: UtilityRootSum[] */
  uint x;
  int i,j,winner;
  real minu, sumrts, t;
  ZeroRealArray( E->NumCands, UtilityRootSum );
  for(i=0; i<(int)E->NumVoters; i++){
    x = i*E->NumCands;
    minu = HUGE;
    for(j=E->NumCands -1; j>=0; j--){
      if(minu > E->Utility[x+j]) minu = E->Utility[x+j];
    }
    sumrts = 0.0;
    for(j=E->NumCands -1; j>=0; j--){
      sumrts += sqrt(E->Utility[x+j] - minu);
    }
    if(sumrts>0.0){
      for(j=E->NumCands -1; j>=0; j--){
	UtilityRootSum[j] += sqrt(E->Utility[x+j] - minu)/sumrts;
      }
    }
  }	
  sumrts = 0.0;
  for(j=E->NumCands -1; j>=0; j--){
    sumrts += UtilityRootSum[j];
  }	
  t = Rand01() * sumrts;
  sumrts = 0.0;
  winner = -1;
  for(j=E->NumCands -1; j>=0; j--){
    sumrts += UtilityRootSum[j];
    if( t<sumrts ){ winner=j; break; }
  }
  assert(winner >= 0);
  return winner;
}

EMETH Plurality(edata *E   /* canddt with most top-rank votes wins */
){ /* side effects: PlurVoteCount[], PlurWinner */
  int i;
  ZeroIntArray( E->NumCands, (int*)PlurVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
    PlurVoteCount[ E->TopDownPrefs[i*E->NumCands+0] ]++;
  }		
  PlurWinner = ArgMaxIntArr( E->NumCands, (int*)PlurVoteCount, (int*)RandCandPerm );
  return PlurWinner;
}

EMETH AntiPlurality(edata *E   /* canddt with fewest bottom-rank votes wins */
		    ){ /* side effects: AntiPlurVoteCount[], AntiPlurWinner */
  int i;
  ZeroIntArray( E->NumCands, (int*)AntiPlurVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
    AntiPlurVoteCount[ E->TopDownPrefs[i*E->NumCands + E->NumCands - 1] ]++;
  }		
  AntiPlurWinner = ArgMinIntArr( E->NumCands, (int*)AntiPlurVoteCount, (int*)RandCandPerm );
  return AntiPlurWinner;
}

/* Plurality needs to have already been run before running Dabagh.  */
EMETH Dabagh(edata *E   /* canddt with greatest Dabagh = 2*#top-rank-votes + 1*#second-rank-votes, wins */
	     ){ /* side effects: DabaghVoteCount[] */
  int i, winner;
  CopyIntArray( E->NumCands, (int*)PlurVoteCount, (int*)DabaghVoteCount );
  ScaleIntVec( E->NumCands, (int*)DabaghVoteCount, 2 );
  for(i=0; i<(int)E->NumVoters; i++){ /*add 2nd-pref votes with weight=1:*/
    DabaghVoteCount[ E->TopDownPrefs[i*E->NumCands +1] ]++;
  }		
  winner = ArgMaxIntArr( E->NumCands, (int*)DabaghVoteCount, (int*)RandCandPerm );
  return winner;
}

/* Plurality needs to have already been run before running VtForAgainst.  */
EMETH VtForAgainst(edata *E   /* canddt with greatest score = #votesFor - #votesAgainst,  wins */
	     ){ /* side effects: VFAVoteCount[] */
  int i, winner, last;
  last = E->NumCands - 1;
  CopyIntArray( E->NumCands, (int*)PlurVoteCount, VFAVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
    VFAVoteCount[ (int)E->TopDownPrefs[i*E->NumCands + last] ]--;
  }		
  winner = ArgMaxIntArr( E->NumCands, VFAVoteCount, (int*)RandCandPerm );
  return winner;
}

EMETH Top2Runoff(edata *E    /* Top2Runoff=top-2-runoff, 2nd round has fully-honest voting */
){ /* side effects: PSecond */
  int i;
  uint offset, pwct=0, wct=0;
  PSecond = -1;
  if(PlurWinner<0) Plurality(E);
  PSecond = Arg2MaxUIntArr( E->NumCands, PlurVoteCount, (int*)RandCandPerm, PlurWinner );
  assert(PSecond>=0);
  for(i=0; i<(int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+PlurWinner] > E->PerceivedUtility[offset+PSecond] ){
       pwct++;
    }else if( E->PerceivedUtility[offset+PlurWinner] < E->PerceivedUtility[offset+PSecond] ){
      wct++;
    }
  }		
  if( pwct > wct || (pwct == wct && RandBool()) ) return(PlurWinner);
  return(PSecond);
}

EMETH VenzkeDisqPlur(edata *E   /* Plurality winner wins unless over 50% vote him bottom; then plurality 2nd-winner wins. */
){ /* side effects: VFAVoteCount[] */
  if(PlurWinner<0) Plurality(E);
  if(AntiPlurWinner<0) AntiPlurality(E);
  if( 2*AntiPlurVoteCount[ PlurWinner ] > E->NumVoters ){
    if(PSecond<0) Top2Runoff(E);
    return PSecond;
  }
  return PlurWinner;
}

EMETH PlurIR(edata *E    /* PlurIR=plur+immediate runoff (give ranking as vote) */
){
  int i;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(PSecond<0) Top2Runoff(E);
  if(PlurWinner<0) Plurality(E);
  if(PSecond<0) Top2Runoff(E);
  assert(PSecond>=0);
  i = E->MarginsMatrix[ PlurWinner*E->NumCands + PSecond ];
  if(i>0) return(PlurWinner);
  if(i<0 || RandBool()) return(PSecond);
  return(PlurWinner);
}

EMETH Borda(edata *E  /* Borda: weighted positional with weights N-1, N-2, ..., 0  if N-canddts */
){ /* side effects: BordaVoteCount[], BordaWinner */
  int i,j,t;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  for(i=E->NumCands -1; i>=0; i--){
    t=0;
    for(j=E->NumCands -1; j>=0; j--){  t += E->MarginsMatrix[i*E->NumCands + j];  }
    BordaVoteCount[i] = t;
  }
  BordaWinner = ArgMaxIntArr( E->NumCands, (int*)BordaVoteCount, (int*)RandCandPerm );  
  return BordaWinner;
}

EMETH Black(edata *E  /* Condorcet winner if exists, else use Borda */
){ 
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CondorcetWinner >= 0) return CondorcetWinner;
  if(BordaWinner<0) Borda(E);
  return BordaWinner;
}

EMETH RandomPair(edata *E){ /*pairwise honest-util winner among 2 random candidates is elected*/
  int x,y;
  x = (int)RandInt(E->NumCands);
  y = (int)RandInt(E->NumCands);
  if( UtilitySum[x] > UtilitySum[y] ) return x;
  if( UtilitySum[x] < UtilitySum[y] ) return y;
  if(RandBool()) return(y);
  return(x);
}

EMETH NansonBaldwin(edata *E  /* repeatedly eliminate Borda loser */
){ /* side effects: NansonVoteCount[], Eliminated[] */
  int i, BordaLoser, rnd, minc, r;
  if(CWSPEEDUP && CondorcetWinner >= 0) return CondorcetWinner;
  if(BordaWinner<0) Borda(E);
  FillBoolArray(E->NumCands, Eliminated, FALSE);
  CopyIntArray(E->NumCands, (int*)BordaVoteCount, NansonVoteCount);
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(rnd=1; rnd < (int)E->NumCands; rnd++){
    BordaLoser = -1;
    minc = BIGINT;
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && NansonVoteCount[r]<minc){
	minc=NansonVoteCount[r]; 
	BordaLoser=r; 
      }
    }
    assert(BordaLoser>=0);
    Eliminated[BordaLoser] = TRUE;
    for(i=E->NumCands -1; i>=0; i--) if(!Eliminated[i]){
      NansonVoteCount[i] -= E->MarginsMatrix[i*E->NumCands + BordaLoser];
    }    
  } /* end of for(rnd) */
  for(i=E->NumCands -1; i>=0; i--){ /* find non-eliminated candidate... */
    if(!Eliminated[i]){
      return i; /*NansonBaldwin winner*/
    }
  }
  return(-1); /*error*/
}

/***
 * Rouse is like Nanson-Baldwin but with an extra level of recursion: it successively
 * pseudo-eliminates the candidate with the highest Borda score until one is left, then it
 * genuinely-eliminates that one from the
 * original list; this step is repeated until a single candidate is left.
 * Slow: O(N^4) steps to find winner for N candidates from pairwise matrix.
 ******/
EMETH Rouse(edata *E  /*like Nanson-Baldwin but with an extra level of recursion BUGGY*/
	    ){ /* side effects: Rmark[], rRmark[] */
  int i,j,k,m,r,highestb,bordsum, maxb, winner;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  FillBoolArray(E->NumCands, rRmark, TRUE); /* nobody eliminated initially */
  for(m= E->NumCands; m>1; m--){ /* NumCands-1 elimination-rounds */
    for(k=1; k<m; k++){  /* m-1 pseudo-elimination rounds */
      maxb = -BIGINT; highestb = -1;
      FillBoolArray(E->NumCands, Rmark, TRUE); /* nobody pseudo-eliminated initially */
      RandomlyPermute( E->NumCands, RandCandPerm );
      for(i=E->NumCands-1; i>=0; i--){
        r = RandCandPerm[i];
        if(rRmark[r] && Rmark[r]){
          bordsum = 0;
	  for(j=E->NumCands -1; j>=0; j--){
	    if(Rmark[j] && rRmark[j]) bordsum += E->MarginsMatrix[r*E->NumCands +j];
	  }
          if(maxb < bordsum){ maxb = bordsum;  highestb = r; }
        }
      }
      assert(highestb >= 0); 
      assert(rRmark[highestb]);
      assert(Rmark[highestb]);
      Rmark[highestb] = FALSE; /* pseudo-eliminate borda-winner */
    }
    /* Find the non-psu-eliminated canddt i: */
    for(i=E->NumCands -1; i>=0; i--){ if(Rmark[i] && rRmark[i]){ break; }} 
    assert(i>=0);
    assert(i<(int)E->NumCands);
    rRmark[i] = FALSE; /* (genuinely) eliminate it */
  }
  winner = -1;
  for(k=0; k< (int)E->NumCands; k++){ if(rRmark[k]){ winner = k;  break; }}
  assert(winner==0 || !rRmark[0]);
  /***???
  if(winner != CondorcetWinner && CondorcetWinner>=0){ *???*
    printf("Rouse aint Condorcet (RW=%d CW=%d):\n", winner, CondorcetWinner);
    PrintEdata(stdout, E);
  }***/
  return winner;
}

EMETH IterCopeland( edata *E  /*iterate Copeland on tied-winner set from previous iter until fixpt*/
 ){ /*Idea of this is due to Alex Small. side effects: summ[] */
  int i, r, j, z, maxc, winner, tiect, oldtiect, x, mxs;
  assert(E->NumCands >= 2);
  if(CWSPEEDUP && CondorcetWinner>=0) return(CondorcetWinner);
  winner = -1; oldtiect = -1;
  ZeroIntArray(E->NumCands, summ);
  RandomlyPermute( E->NumCands, RandCandPerm );
  mxs = 0;
  for(z=E->NumCands; z>0; z--){
    tiect = 0;
    for(i=E->NumCands -1; i>=0; i--) if(summ[i] < mxs){ summ[i] = -BIGINT; }
    maxc = -BIGINT;
    for(i=E->NumCands-1; i>=0; i--){
      r = RandCandPerm[i];
      if(summ[r] >= mxs){
	x = r*E->NumCands;
	for(j=E->NumCands -1; j>=0; j--){
	  if(j!=r && summ[j]>=mxs){
	    summ[r] += 1 + SignInt(E->MarginsMatrix[x+j]);
	  }
	}
	if(summ[r] >= maxc){
	  if(summ[r]==maxc){ tiect++; }
	  else{ maxc=summ[r]; winner=r;  tiect=1; }
	}
      }
    }
    assert(tiect>=1);
    if(tiect<=1 || tiect==oldtiect) return(winner);
    oldtiect = tiect;
  }
  return(-1);
}

EMETH Nauru(edata *E  /* weighted positional with weights 1, 1/2, 1/3, 1/4,... */
){	    /* side effects: NauruVoteCount[] */
  int i,j,x,winner;
  ZeroIntArray( E->NumCands, (int*)NauruVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
       NauruVoteCount[j] += NauruWt[E->CandRankings[x+j]];
     }
  }		
  winner = ArgMaxUIntArr( E->NumCands, NauruVoteCount, (int*)RandCandPerm );  
  return(winner);
}

EMETH HeismanTrophy(edata *E  /* Heisman: weighted positional with weights 3, 2, 1, 0,...,0 */
){  /* side effects: HeismanVoteCount[] */
  int i,j,x,winner;
  ZeroIntArray( E->NumCands, (int*)HeismanVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=MinInt(E->NumCands -1, 2); j>=0; j--){
       HeismanVoteCount[ E->TopDownPrefs[x+j] ] += 3-j;
     }
  }		
  winner = ArgMaxUIntArr( E->NumCands, HeismanVoteCount, (int*)RandCandPerm );  
  return(winner);
}

EMETH BaseballMVP(edata *E  /* weighted positional with weights 14,9,8,7,6,5,4,3,2,1,0,...,0 */
){	    /* side effects: BaseballVoteCount[] */
  int i,j,x,winner;
  ZeroIntArray( E->NumCands, (int*)BaseballVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=MinInt( E->NumCands -1, 9); j>=0; j--){
       BaseballVoteCount[ E->TopDownPrefs[x+j] ] += BaseballWt[j];
     }
  }		
  winner = ArgMaxUIntArr( E->NumCands, BaseballVoteCount, (int*)RandCandPerm );  
  return(winner);
}

EMETH CondorcetLR(edata *E   /* candidate with least sum-of-pairwise-defeat-margins wins */
){ /* side effects:, SumOfDefeatMargins[]  */
   int i,j,t,winner;
   if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
   if(CWSPEEDUP && CondorcetWinner>=0) return(CondorcetWinner);
   for(i=E->NumCands -1; i>=0; i--){
     t = 0;
     for(j=E->NumCands -1; j>=0; j--){  t += PosInt( E->MarginsMatrix[j*E->NumCands+i] );  }
     SumOfDefeatMargins[i] = t;
   }
   winner = ArgMinIntArr( E->NumCands, (int*)SumOfDefeatMargins, (int*)RandCandPerm );
   return winner;
}

EMETH Sinkhorn(edata *E  /* candidate with max Sinkhorn rating (from all-positive DefeatsMatrix+1) */
){ /* side effects SinkRat[], SinkRow[], SinkCol[], SinkMat[] */
  int j,k,winner;
  real t,maxsum,minsum,sum,maxminRatio;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  FillRealArray( E->NumCands, SinkRow, 1.0 );
  FillRealArray( E->NumCands, SinkCol, 1.0 );
  do{
    for(k=0; k < (int)E->NumCands; k++){
      for(j=E->NumCands -1; j>=0; j--){
	SinkMat[k*E->NumCands + j] =
	  SinkRow[k]*SinkCol[j] * (E->DefeatsMatrix[k*E->NumCands + j] + 1.0);
      }
    }
    maxsum = -HUGE; minsum = HUGE;
    for(k=0; k < (int)E->NumCands; k++){
      sum = 0.0; 
      for(j=E->NumCands -1; j>=0; j--){
	sum += SinkMat[k*E->NumCands + j];
      }
      if(minsum > sum) minsum = sum;
      if(maxsum < sum) maxsum = sum;
      SinkRow[k] /= sum;
    }
    maxminRatio = maxsum/minsum;
    maxsum = -HUGE; minsum = HUGE;
    for(k=0; k < (int)E->NumCands; k++){
      sum = 0.0; 
      for(j=E->NumCands -1; j>=0; j--){
	sum += SinkMat[j*E->NumCands + k];
      }
      if(minsum > sum) minsum = sum;
      if(maxsum < sum) maxsum = sum;
      SinkCol[k] /= sum;
    }
    t = maxsum/minsum;
    if( maxminRatio < t ) maxminRatio = t;
  }until(maxminRatio < 1.000003);
  for(j=E->NumCands -1; j>=0; j--){
    SinkRat[j] = SinkCol[j]/SinkRow[j];
  }
  winner = ArgMaxRealArr( E->NumCands, SinkRat, (int*)RandCandPerm );
  return winner;
}

EMETH KeenerEig(edata *E  /* winning canddt has max Frobenius eigenvector entry (of all-positive DefeatsMatrix+1) */
){ /* Side effects EigVec[], EigVec2[] */
  int j,k,winner;
  real t,sum,dist;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  FillRealArray( E->NumCands, EigVec, 1.0 );
  FillRealArray( E->NumCands, EigVec2, 1.0 );
  do{
    for(k=0; k < (int)E->NumCands; k++){
      t = 0;
      for(j=E->NumCands -1; j>=0; j--){
	t += (E->DefeatsMatrix[k*E->NumCands + j] + 1.0) * EigVec[j];
      }
      EigVec2[k] = t;
    }
    sum = SumRealArray( E->NumCands, EigVec2 );
    ScaleRealVec( E->NumCands, EigVec2, 1.0/sum );
    dist = L1Distance(E->NumCands, EigVec, EigVec2);
    CopyRealArray( E->NumCands, EigVec2, EigVec );
  }until( dist < 0.00001 );
  winner = ArgMaxRealArr( E->NumCands, EigVec, (int*)RandCandPerm );
  return winner;
}

EMETH SimpsonKramer(edata *E  /* candidate with mildest worst-defeat wins */
){ /* Side effects: WorstDefeatMargin[] */
  int i,r,x,t,j,winner;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CWSPEEDUP && CondorcetWinner>=0) return(CondorcetWinner);
  for(i=E->NumCands -1; i>=0; i--){
    t = 0;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      x = E->MarginsMatrix[r*E->NumCands+i];
      if(x>t) t=x;
    }
    WorstDefeatMargin[i] = t;
  }
  winner = ArgMinIntArr( E->NumCands, WorstDefeatMargin, (int*)RandCandPerm );
  return winner;
}

EMETH RaynaudElim(edata *E  /* repeatedly eliminate canddt who suffered the worst-margin-defeat */
 ){ /* side effects: Eliminated[], RayDefeatMargin[], RayBeater[] */
  int i, j, x, t, RayLoser, rnd, maxc, r, beater;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CWSPEEDUP && CondorcetWinner >= 0) return CondorcetWinner;
  for(i=E->NumCands -1; i>=0; i--){
    t = -BIGINT; beater = -1;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      x = E->MarginsMatrix[r*E->NumCands+i];
      if(x>t){ t=x; beater=r; }
    }
    assert(beater >= 0);
    RayDefeatMargin[i] = t; /*worst margin of defeat of i, nonpositive if undefeated */
    RayBeater[i] = beater; /*who administered that beating*/
  }
  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(rnd=1; rnd < (int)E->NumCands; rnd++){
    RayLoser = -1;
    maxc = -BIGINT;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && RayDefeatMargin[r]>maxc){
	maxc=RayDefeatMargin[r]; 
	RayLoser=r; 
      }
    }
    assert(RayLoser >= 0);
    if( maxc <= 0 ){ return RayLoser; } /*"loser" is undefeated*/
    Eliminated[RayLoser] = TRUE;
    for(i=E->NumCands -1; i>=0; i--) if(!Eliminated[i] && RayBeater[i]==RayLoser){
      t = -BIGINT; beater = -1;
      RandomlyPermute( E->NumCands, RandCandPerm );
      for(j=E->NumCands -1; j>=0; j--){
	r = RandCandPerm[j];
	if(!Eliminated[r]){
	  x = E->MarginsMatrix[r*E->NumCands+i];
	  if(x>t){ t=x; beater=r; }
	}
      }
      assert(beater >= 0);
      RayDefeatMargin[i] = t;
      RayBeater[i] = beater;
    }    
  } /* end of for(rnd) */
  for(i=E->NumCands -1; i>=0; i--){ /* find non-eliminated candidate... */
    if(!Eliminated[i]){
      return i; /*Raynaud winner*/
    }
  }
  return(-1); /*error*/
}

EMETH ArrowRaynaud(edata *E  /* repeatedly eliminate canddt with smallest {largest margin of victory, which is <=0 if never won} */
){ /* side effects: Eliminated[], ARchump[], ARVictMargin[].  ArrowRaynaud can eliminate a Condorcet Winner in round #1. */
  int i, j, x, t, ARLoser, rnd, minc, r, chump;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  for(i=E->NumCands -1; i>=0; i--){
    t = -BIGINT; chump = -1;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      x = E->MarginsMatrix[i*E->NumCands + r];
      if(x>t){ t=x; chump=r; }
    }
    assert(chump >= 0);
    ARVictMargin[i] = t; /*largest margin of victory of i, nonpositive if never won*/
    ARchump[i] = chump; /*who suffered that beating*/
  }
  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(rnd=1; rnd < (int)E->NumCands; rnd++){
    ARLoser = -1;
    minc = BIGINT;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && ARVictMargin[r]<minc){
	minc=ARVictMargin[r]; 
	ARLoser=r; 
      }
    }
    assert(ARLoser >= 0);
    Eliminated[ARLoser] = TRUE;
    for(i=E->NumCands -1; i>=0; i--) if(!Eliminated[i] && ARchump[i]==ARLoser){
      t = -BIGINT; chump = -1;
      RandomlyPermute( E->NumCands, RandCandPerm );
      for(j=E->NumCands -1; j>=0; j--){
	r = RandCandPerm[j];
	if(!Eliminated[r]){
	  x = E->MarginsMatrix[i*E->NumCands+r];
	  if(x>t){ t=x; chump=r; }
	}
      }
      assert(chump >= 0);
      ARVictMargin[i] = t;
      ARchump[i] = chump;
    }    
  } /* end of for(rnd) */
  for(i=E->NumCands -1; i>=0; i--){ /* find non-eliminated candidate... */
    if(!Eliminated[i]){
      return i; /*ArrowRaynaud winner*/
    }
  }
  return(-1); /*error*/
}

/* O(N^3) algorithm, but it is known how to speed it up to O(N^2) */
EMETH SchulzeBeatpaths(edata *E  /* winner = X so BeatPathStrength over rivals Y exceeds strength from Y */
){ /* Side effects: BeatPathStrength[] */
  int i,j,k,minc,winner;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  for(i=E->NumCands -1; i>=0; i--){
    for(j=E->NumCands -1; j>=0; j--) if(i != j){
      BeatPathStrength[i*E->NumCands +j] = E->MarginsMatrix[i*E->NumCands +j];
    }
  }
  for(i=E->NumCands -1; i>=0; i--){
    for(j=E->NumCands -1; j>=0; j--) if(i != j){
      for(k=0; k<(int)E->NumCands; k++) if(k != j && k != i){
	minc = BeatPathStrength[j*E->NumCands+i];
	if( BeatPathStrength[i*E->NumCands +k] < minc ) minc = BeatPathStrength[i*E->NumCands +k];
	if( BeatPathStrength[j*E->NumCands +k] < minc ) BeatPathStrength[j*E->NumCands +k] = minc;
      }
    }
  }
  for(i=E->NumCands -1; i>=0; i--){
    k = RandCandPerm[i];
    for(j=E->NumCands -1; j>=0; j--) if(k != j){
      if( BeatPathStrength[j*E->NumCands +k] > BeatPathStrength[k*E->NumCands +j] ){
	goto KNOTSCHULZEWINNER;
      }
    }
    winner = k;   return winner;
    KNOTSCHULZEWINNER: ;
  }
  return(-1);
}

/* Note about Smith and Schwartz sets:
BeatPathStrength[k*E->NumCands+j] > 0   for all k in the "Smith Set" and j outside it.
BeatPathStrength[k*E->NumCands+j] >= 0  for all k in the "Schwartz Set" and j outside it.
*******/

void beatDFS( int x, int diff, bool Set[], int Mat[], int N ){
  int i;
  for(i=N-1; i>=0; i--) if(i!=x){
    if( Mat[i*N+x]>=diff ){
      if( !Set[i] ){
	Set[i] = TRUE;
	beatDFS( i, diff, Set, Mat, N );
      }
    }
  }
}

EMETH SmithSet(edata *E  /* Smith set = smallest nonempty set of canddts that pairwise-beat all nonmembers */
){ /* side effects: SmithMembs[] */
  int i,r;
  FillBoolArray(E->NumCands, SmithMembs, FALSE);
  assert(CopeWinOnlyWinner>=0);
  assert(CopeWinOnlyWinner < (int)E->NumCands);
  SmithMembs[CopeWinOnlyWinner] = TRUE;
  beatDFS(CopeWinOnlyWinner, 1, SmithMembs, E->MarginsMatrix, E->NumCands);
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){
    r = RandCandPerm[i];
    if(SmithMembs[r]) return r; /*return random set member*/
  }
  return(-1);
}

EMETH SchwartzSet(edata *E  /* Schwartz set = smallest nonempty set of canddts undefeated by nonmembers */
){ /* side effects: SchwartzMembs[] */
  int i,r;
  FillBoolArray(E->NumCands, SchwartzMembs, FALSE);
  assert(CopeWinOnlyWinner>=0);
  assert(CopeWinOnlyWinner < (int)E->NumCands);
  SchwartzMembs[CopeWinOnlyWinner] = TRUE;
  beatDFS(CopeWinOnlyWinner, 0, SchwartzMembs, E->MarginsMatrix, E->NumCands);
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){
    r = RandCandPerm[i];
    if(SchwartzMembs[r]) return r; /*return random set member*/
  }
  return(-1);
}

/* Uncovered set.  A "covers" B if the candidates A beats pairwise 
are a superset of those B beats pairwise.  "Landau's theorem":
All Copeland winners are members of the uncovered set.
Can do this an order of magnitude faster if MaxNumCands<sizeof(uint)*4 because set inclusion can
be tested in 1 step using wordwide operations??
****/
EMETH SlowUncoveredSet(edata *E /*A "covers" B if A beats a strict?? superset of those B beats.*/
){ /* side effects: UncoveredSt[], CoverMatrix[] */
  int A,B,C,i,r;
  bool cov;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(SchwartzWinner<0) SchwartzSet(E);
  /*find cover relation:*/
  for(A=0; A < (int)E->NumCands; A++){
    for(B=0; B < (int)E->NumCands; B++) if(B!=A){
      cov = TRUE;
      for(C=0; C < (int)E->NumCands; C++){
	if( E->MarginsMatrix[A*E->NumCands + C]<=0 &&
	    E->MarginsMatrix[B*E->NumCands + C]>0  ){
	  cov = FALSE; break;
	}
      }
      CoverMatrix[A*E->NumCands + B] = cov;
    }
    UncoveredSt[A] = TRUE; /*initialization*/
  }
  for(A=0; A < (int)E->NumCands; A++){
    for(B=0; B < A; B++) if(B!=A){
      if( CoverMatrix[A*E->NumCands + B] &&  CoverMatrix[B*E->NumCands + A] ){
	/*enforce strict superset; otherwise coverage both ways would be possible:*/
	CoverMatrix[A*E->NumCands + B] = FALSE;
	CoverMatrix[B*E->NumCands + A] = FALSE;
      }
    }
  }
  /*find UncoveredSt:*/
  for(A=0; A < (int)E->NumCands; A++){
    for(B=0; B < (int)E->NumCands; B++) if(B!=A){
      if( CoverMatrix[B*E->NumCands+A] ){
	UncoveredSt[A] = FALSE; break;
      }
    }
  }
  /*select random uncovered winner:*/
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){
    if( !(UncoveredSt[i]?SchwartzMembs[i]:TRUE) ){
      printf("bozo!\n");
      printf("%d %d %d %d %d %d %d %d %d\n", 
	     E->MarginsMatrix[0*E->NumCands + 0],
	     E->MarginsMatrix[0*E->NumCands + 1],
	     E->MarginsMatrix[0*E->NumCands + 2],
	     E->MarginsMatrix[1*E->NumCands + 0],
	     E->MarginsMatrix[1*E->NumCands + 1],
	     E->MarginsMatrix[1*E->NumCands + 2],
	     E->MarginsMatrix[2*E->NumCands + 0],
	     E->MarginsMatrix[2*E->NumCands + 1],
	     E->MarginsMatrix[2*E->NumCands + 2]  );
      printf("CopeWinOnlyWinner=%d\n",CopeWinOnlyWinner);
      printf("Sc=%d%d%d\n", SchwartzMembs[0], SchwartzMembs[1], SchwartzMembs[2]);
      printf("Un=%d%d%d\n", UncoveredSt[0], UncoveredSt[1], UncoveredSt[2]);
    }
    assert( UncoveredSt[i]?SchwartzMembs[i]:TRUE );
    r = RandCandPerm[i];
    if(UncoveredSt[r]){
      RandomUncoveredMemb = r;
      return r; 
    }
  }
  printf("yikes!\n");
  printf("%d %d %d %d\n", 
	 E->MarginsMatrix[0*E->NumCands + 0],
	 E->MarginsMatrix[1*E->NumCands + 0],
	 E->MarginsMatrix[0*E->NumCands + 1],
	 E->MarginsMatrix[1*E->NumCands + 1] );
  return(-1);
}


/* Uncovered set.  A "covers" B if the candidates A beats pairwise 
are a superset of those B beats pairwise.  "Landau's theorem":
All Copeland winners are members of the uncovered set.
Now an order of magnitude faster assuming MaxNumCands<sizeof(uint)*4 because set inclusion can
be tested in 1 step using wordwide operations.
****/
EMETH UncoveredSet(edata *E /*A "covers" B if A beats a strict superset of those B beats.*/
){ /* side effects: UncoveredSt[], CoverMatrix[] */
  int A,B,i,r;
  if( E->NumCands > 4*sizeof(uint) ){
    printf("UncoveredSet: too many candidates %d to use machine words(%d) to represent sets\n", 
	   E->NumCands, 
	   (int)(4*sizeof(uint)) );
    printf("You could rewrite the code to use uint64s to try allow up to 64 canddts\n");
    exit(EXIT_FAILURE);
  }
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(SchwartzWinner<0) SchwartzSet(E);
  /*find cover relation:*/
  for(A=0; A < (int)E->NumCands; A++){
    for(B=0; B < (int)E->NumCands; B++) if(B!=A){
	CoverMatrix[A*E->NumCands + B] = StrictSuperset(IBeat[A], IBeat[B]);
    }
    UncoveredSt[A] = TRUE; /*initialization*/
  }
  /*find UncoveredSt:*/
  for(A=0; A < (int)E->NumCands; A++){
    for(B=0; B < (int)E->NumCands; B++) if(B!=A){
      if( CoverMatrix[B*E->NumCands+A] ){
	UncoveredSt[A] = FALSE; break;
      }
    }
  }
  /*select random uncovered winner:*/
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){
    if( !(UncoveredSt[i]?SchwartzMembs[i]:TRUE) ){
      printf("bozo! i=%d NumCands=%d\n", i, E->NumCands);
      printf("%d %d %d; %d %d %d; %d %d %d\n", 
	     E->MarginsMatrix[0*E->NumCands + 0],
	     E->MarginsMatrix[0*E->NumCands + 1],
	     E->MarginsMatrix[0*E->NumCands + 2],
	     E->MarginsMatrix[1*E->NumCands + 0],
	     E->MarginsMatrix[1*E->NumCands + 1],
	     E->MarginsMatrix[1*E->NumCands + 2],
	     E->MarginsMatrix[2*E->NumCands + 0],
	     E->MarginsMatrix[2*E->NumCands + 1],
	     E->MarginsMatrix[2*E->NumCands + 2]  );
      printf("CopeWinOnlyWinner=%d\n",CopeWinOnlyWinner);
      printf("Sc=%d%d%d\n", SchwartzMembs[0], SchwartzMembs[1], SchwartzMembs[2]);
      printf("Un=%d%d%d\n", UncoveredSt[0], UncoveredSt[1], UncoveredSt[2]);
    }
    assert( UncoveredSt[i]?SchwartzMembs[i]:TRUE );
    r = RandCandPerm[i];
    if(UncoveredSt[r]){
      RandomUncoveredMemb = r;
      return r; 
    }
  }
  printf("yikes!\n");
  printf("%d %d %d %d\n", 
	 E->MarginsMatrix[0*E->NumCands + 0],
	 E->MarginsMatrix[1*E->NumCands + 0],
	 E->MarginsMatrix[0*E->NumCands + 1],
	 E->MarginsMatrix[1*E->NumCands + 1] );
  return(-1);
}

EMETH Bucklin(edata *E   /* canddt with over half the voters ranking him in top-k wins (for k minimal) */
){ /* side effects: RdVoteCount[] */
  int rnd,i,best,winner;
  winner = -1;
  ZeroIntArray( E->NumCands,  RdVoteCount );
  for(rnd=0; rnd<(int)E->NumCands; rnd++){
    for(i=0; i<(int)E->NumVoters; i++){
       RdVoteCount[ E->TopDownPrefs[i*E->NumCands + rnd] ]++;
    }
    winner = ArgMaxIntArr( E->NumCands, RdVoteCount, (int*)RandCandPerm );
    best = RdVoteCount[winner];
    if(best*2 > (int)E->NumVoters)  break;
  }
  return winner;
}

/******
James Green-Armytage:  Definition of the cardinal-weighted pairwise comparison method
I. Voters rate the candidates, e.g. on a scale from 0 to 100. 
Equal ratings are allowed.   There then is an implied ranking.
II. Tally
1. Determine the direction of the pairwise defeats by using the rankings for a standard pairwise 
comparison tally.
2. Determine the strength of the pairwise defeats by finding the weighted magnitude as follows. 
Suppose candidate A pairwise beats B, and we want to know the strength of the defeat. 
For each voter who ranks A over B, and only for voters who rank A over B, 
subtract their rating of B from their rating of A, to get the rating differential. 
The sum of these individual winning rating differentials is the total weighted magnitude of the defeat. 
(Note that, because voters who rank B over A do not contribute to the weighted magnitude of the A>B defeat, 
it cannot be negative.)
3. Now that the direction of the pairwise defeats have been determined (in step 1) and 
the strength of the defeats have been determined (in step 2), you can choose from a 
variety of Condorcet completion methods to determine the winner. 
I recommend the ranked pairs, beatpath, and river methods.
WDS: IEVS will use beatpaths.  What is good voter strategy in this method?
*********/

EMETH ArmytagePCSchulze(edata *E  /*Armytage pairwise comparison based on Schulze*/
){ /* Side effects: ArmyBPS[] */
  int i,j,k,winner;
  real minc;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  for(i=E->NumCands -1; i>=0; i--){
    for(j=E->NumCands -1; j>=0; j--) if(i != j){
      ArmyBPS[i*E->NumCands +j] = E->MargArmy[i*E->NumCands +j];
    }
  }
  for(i=E->NumCands -1; i>=0; i--){
    for(j=E->NumCands -1; j>=0; j--) if(i != j){
      for(k=0; k<(int)E->NumCands; k++) if(k != j && k != i){
	minc = ArmyBPS[j*E->NumCands+i];
	if( ArmyBPS[i*E->NumCands +k] < minc ) minc = ArmyBPS[i*E->NumCands +k];
	if( ArmyBPS[j*E->NumCands +k] < minc ) ArmyBPS[j*E->NumCands +k] = minc;
      }
    }
  }
  for(i=E->NumCands -1; i>=0; i--){
    k = RandCandPerm[i];
    for(j=E->NumCands -1; j>=0; j--) if(k != j){
      if( ArmyBPS[j*E->NumCands +k] > ArmyBPS[k*E->NumCands +j] ){
	goto KNOTSW;
      }
    }
    winner = k;   return winner;
    KNOTSW: ;
  }
  return(-1);
}

EMETH Copeland(edata *E   /* canddt with largest number of pairwise-wins elected (tie counts as win/2) BUGGY */
){ /* side effects: CopelandWinner, CopeScore[] */
  int i;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CWSPEEDUP && CondorcetWinner>=0){
     CopelandWinner = CondorcetWinner;
     return CopelandWinner;
  }
  for(i=E->NumCands -1; i>=0; i--){ CopeScore[i] = 2*WinCount[i]+DrawCt[i]; }
  CopelandWinner = ArgMaxIntArr( E->NumCands, CopeScore, (int*)RandCandPerm );  
  /* Currently just break ties randomly, return random highest-scorer */
  return CopelandWinner;
}

EMETH SimmonsCond(edata *E  /* winner = X with least sum of top-rank-votes for rivals pairwise-beating X */
){ /* side effects: SimmVotesAgainst[] */
  int i,j,t,winner;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(PlurWinner<0) Plurality(E);
  if(SmithWinner<0) SmithSet(E);
  if(SchwartzWinner<0) SchwartzSet(E);
   if(CWSPEEDUP && CondorcetWinner>=0) return(CondorcetWinner);
   for(i=E->NumCands -1; i>=0; i--){
     t=0;
     for(j=E->NumCands -1; j>=0; j--) if(E->MarginsMatrix[j*E->NumCands+i]>0){ /* j pairwise-beats i */
       t += 3*PlurVoteCount[j] + (SchwartzMembs[j]?1:0) + (SmithMembs[j]?1:0);
       /*Here I am adding 1/3 of a vote if in SmithSet, ditto SchwartzSet, to break Simmons ties*/
     }
     SimmVotesAgainst[i] = t;
   }
   winner = ArgMinIntArr( E->NumCands, SimmVotesAgainst, (int*)RandCandPerm );
   return winner;
}

/*******
 * The following IRV (instant runoff voting) algorithm has somewhat long code, but
 * by using linked lists achieves fast (very sublinear) runtime.  
 * If SmithIRVwinner<0 then It also computes SmithIRVwinner as a side effect.
 * This code can also be used to compute the winner in "Top3" (bastardized) IRV where
 * voters only allowed to indicate their top 3 choices in order..
 * For normal IRV (or SmithIRV) set IRVTopLim=BIGINT before run.
 * For Top3-IRV set IRVTopLim=3 (or any other integer N for TopN-IRV). ***/
EMETH IRV(edata *E   /* instant runoff; repeatedly eliminate plurality loser */
){ /* side effects: Eliminated[], IFav[], RdVoteCount[], FavListNext[], HeadFav[], LossCount[], SmithIRVwinner, IRVwinner  */
  int Iround,i,RdLoser,NextI,j,t,x,minc,r,stillthere,winner;
  assert(E->NumCands <= MaxNumCands);
  if(SmithIRVwinner<0 && IRVTopLim==BIGINT && CopeWinOnlyWinner<0) BuildDefeatsMatrix(E); 
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){ 
    Eliminated[i] = FALSE; 
    HeadFav[i] = -1; /*HeadFav[i] will be the first voter whose current favorite is i*/
    if(SmithIRVwinner<0 && IRVTopLim==BIGINT){ 
      t=0; 
      for(j=E->NumCands -1; j>=0; j--) if(E->MarginsMatrix[j*E->NumCands + i] > 0){ t++; } 
      LossCount[i] = t;  
    }
  } /*end for(i)*/
  ZeroIntArray(E->NumCands, RdVoteCount); 
  ZeroIntArray(E->NumVoters, IFav); 
  /* IFav[i] is the rank of the 1st noneliminated canddt in voter i's topdownpref list (initially 0) */
  FillIntArray(E->NumVoters, FavListNext, -1);
  /* FavListNext is "next" indices in linked list of voters with common current favorite; -1 terminated. */
  if(SmithIRVwinner<0 && IRVTopLim==BIGINT){ SmithIRVwinner = CondorcetWinner; } 
  /* compute vote totals for 1st round and set up forward-linked lists (-1 terminates each list): */
  for(i=E->NumVoters -1; i>=0; i--){
    x = E->TopDownPrefs[i*E->NumCands+IFav[i]]; /* the initial favorite of voter i */
    assert(x >= 0);
    assert(x < (int)E->NumCands);
    RdVoteCount[ x ]++;
    FavListNext[i] = HeadFav[x];
    HeadFav[x] = i;
  }		
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(Iround=1; Iround < (int)E->NumCands; Iround++){ /*perform IRV rounds*/
    RdLoser = -1;
    minc = BIGINT;
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && RdVoteCount[r]<minc){
	minc=RdVoteCount[r]; 
	RdLoser=r; 
      }
    }
    assert(RdLoser>=0);
    assert(RdLoser < (int)E->NumCands);
    Eliminated[RdLoser] = TRUE; /* eliminate RdLoser */
    if(IRVTopLim==BIGINT && SmithIRVwinner < 0){  
      for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){ /* update LossCount[j] */
	t = E->MarginsMatrix[RdLoser*E->NumCands + j];   
	if(t>0){ LossCount[j] --; }   
	if( LossCount[j] <= 0 ){ SmithIRVwinner = j; break; }   
      }   
    }   
    for(i=HeadFav[RdLoser]; i>=0; i=NextI){/*Go thru linked list of voters with favorite=RdLoser, adjust:*/
      j = i*E->NumCands;
      NextI =  FavListNext[i];
      assert(IFav[i] >= 0);
      assert(IFav[i] < (int)E->NumCands);
      assert( E->TopDownPrefs[ j+IFav[i] ] == RdLoser );
      do{ 
	IFav[i]++; x = E->TopDownPrefs[j + IFav[i]]; 
      }while( Eliminated[x] );
      /* x is new favorite of voter i (or ran out of favorites) */
      assert( IFav[i] < (int)E->NumCands );
      assert(x >= 0);
      assert(x < (int)E->NumCands);
      /* update favorite-list: */
      FavListNext[i] = HeadFav[x]; 
      HeadFav[x] = i;
      /* update vote count totals: */
      if(IFav[i] < IRVTopLim){	RdVoteCount[ x ]++;   }
    } /*end for(i)*/
  }  /* end of for(Iround) */
  stillthere = 0;
  if(IRVTopLim >= (int)E->NumCands) IRVwinner = -1;
  winner = -1;
  for(i=E->NumCands -1; i>=0; i--){ /* find non-eliminated candidate... */
    if(!Eliminated[i]){
      winner=i;
      stillthere++;
    }
  }
  if(IRVTopLim >= (int)E->NumCands) IRVwinner=winner;
  assert(stillthere==1);
  return(winner); 
}

EMETH SmithIRV(edata *E  /*  Eliminate plurality loser until unbeaten candidate exists. */
){ /* must be run after IRV. */
  if(IRVwinner<0){
    SmithIRVwinner = -1;
    IRV(E);
  }
  return SmithIRVwinner;
}

EMETH Top3IRV(edata *E  /* Top-3-choices-only IRV */
){ 
  int w;
  IRVTopLim = 3;
  w = IRV(E);
  IRVTopLim = BIGINT;
  return w;
}

EMETH BTRIRV(edata *E  /* Repeatedly eliminate either plur loser or 2nd-loser (whoever loses pairwise) */
){ /* side effects: Eliminated[], IFav[], RdVoteCount[], FavListNext[], HeadFav[], */
  int Iround,x,i,RdLoser,RdLoser2,NextI,j,minc,r;
  assert(E->NumCands <= MaxNumCands);
  if(CWSPEEDUP && CondorcetWinner>=0) return(CondorcetWinner);
  for(i=E->NumCands -1; i>=0; i--){ 
    Eliminated[i] = FALSE; 
    HeadFav[i] = -1; 
  }
  ZeroIntArray(E->NumCands, RdVoteCount); 
  ZeroIntArray(E->NumVoters, IFav);
  /* compute vote totals for 1st round and set up forward-linked lists (-1 terminates each list): */
  for(i=E->NumVoters -1; i>=0; i--){
    assert(IFav[i] >= 0);
    assert(IFav[i] < (int)E->NumCands);
    x = E->TopDownPrefs[i*E->NumCands + IFav[i]]; /* the favorite of voter i */
    assert(x >= 0);
    assert(x < (int)E->NumCands);
    RdVoteCount[ x ]++;
    FavListNext[i] = HeadFav[x];
    HeadFav[x] = i;
  }		
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(Iround=1; Iround<(int)E->NumCands; Iround++){
    RdLoser = -1;
    minc = BIGINT;
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && RdVoteCount[r]<minc){
	minc=RdVoteCount[r]; 
	RdLoser=r; 
      }
    }
    assert(RdLoser>=0);
    RdLoser2 = -1;
    minc = BIGINT;
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      if(!Eliminated[r] && RdVoteCount[r]<minc && r!=RdLoser){
	minc=RdVoteCount[r]; 
	RdLoser2=r; 
      }
    }
    assert(RdLoser2>=0);
    if( E->MarginsMatrix[ RdLoser*E->NumCands + RdLoser2 ] > 0 ) RdLoser = RdLoser2;
    Eliminated[RdLoser] = TRUE; /* eliminate RdLoser */
    for(i=HeadFav[RdLoser]; i>=0; i=NextI){ /* Go thru list of voters with favorite=RdLoser, adjust: */
      j = i*E->NumCands;
      assert( E->TopDownPrefs[j+IFav[i]] == RdLoser );
      assert(IFav[i] >= 0);
      assert(IFav[i] < (int)E->NumCands);
      do{ IFav[i]++; x = E->TopDownPrefs[j + IFav[i]]; }while( Eliminated[x] );
      /* x is new favorite of voter i */
      assert( IFav[i] < (int)E->NumCands );
      NextI =  FavListNext[i];
      /* update favorite-list: */
      FavListNext[i] = HeadFav[x]; 
      HeadFav[x] = i;
      /* update vote count totals: */
      assert(x >= 0);
      assert(x < (int)E->NumCands);
      RdVoteCount[ x ]++;
    }
  }
  for(i=E->NumCands -1; i>=0; i--){ /* find the non-eliminated candidate... */
    if(!Eliminated[i]){
      return i; /*IRV winner*/
    }
  }
  return(-1); /*error*/
}

EMETH Coombs(edata *E /*repeatedly eliminate antiplurality loser (with most bottom-rank votes)*/
){ /*side effects: Eliminated[], IFav[], RdVoteCount[], FavListNext[], HeadFav[] */
  int Iround,i,RdLoser,NextI,j,x,maxc,r,stillthere;
  assert(E->NumCands <= MaxNumCands);
  for(i=E->NumCands -1; i>=0; i--){ 
    Eliminated[i] = FALSE; 
    HeadFav[i] = -1; /*HeadFav[i] will be the first voter whose current most-hated is i*/
  }
  ZeroIntArray(E->NumCands, RdVoteCount); 
  assert(E->NumCands >= 2);
  assert(E->NumCands <= MaxNumVoters);
  FillIntArray(E->NumVoters, IFav, E->NumCands -1); 
  /* IFav[i] is the rank of the last noneliminated canddt in voter i's topdownpref list 
   * (initially NumCands-1) */
  FillIntArray(E->NumVoters, FavListNext, -1);
  /* FavListNext is "next" indices in linked list of voters with common current favorite; -1 terminated. */
  /* compute vote totals for 1st round and set up forward-linked lists (-1 terminates each list): */
  for(i=E->NumVoters -1; i>=0; i--){
    assert(IFav[i] >= 0);
    assert(IFav[i] < (int)E->NumCands);
    x = E->TopDownPrefs[i*E->NumCands + IFav[i]]; /* the initial most-hated of voter i */
    assert(x >= 0);
    assert(x < (int)E->NumCands);
    RdVoteCount[ x ]++; /*antiplurality voting*/
    FavListNext[i] = HeadFav[x];
    HeadFav[x] = i;
  }		
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(Iround=1; Iround < (int)E->NumCands; Iround++){   
    RdLoser = -1;
    maxc = -BIGINT;
    stillthere=0;
    for(i=E->NumCands -1; i>=0; i--){
      r = RandCandPerm[i];
      assert(r >= 0);
      assert(r < (int)E->NumCands);
      if(!Eliminated[r]){
	stillthere++;
	if(RdVoteCount[r]>maxc){
	  maxc=RdVoteCount[r]; 
	  RdLoser=r; 
	}
      }
    }
    assert(RdLoser>=0);
    Eliminated[RdLoser] = TRUE; /* eliminate RdLoser */
    for(i=HeadFav[RdLoser]; i>=0; i=NextI){/*Go thru linked list of voters with favorite=RdLoser, adjust:*/
      j = i*E->NumCands;
      assert( IFav[i]>=0 );
      assert( IFav[i]<= (int)E->NumCands );
      assert( E->TopDownPrefs[ j+IFav[i] ] == RdLoser );
      do{ IFav[i]--; x = E->TopDownPrefs[j+IFav[i]]; }while( Eliminated[x] );
      /* x is new favorite of voter i */
      assert( IFav[i] >= 0 );
      NextI =  FavListNext[i];
      /* update favorite-list: */
      FavListNext[i] = HeadFav[x]; 
      HeadFav[x] = i;
      /* update vote count totals: */
      assert(x>=0);
      assert(x < (int)E->NumCands);
      RdVoteCount[ x ]++;
    }
  }  /* end of for(Iround) */
  for(i=E->NumCands -1; i>=0; i--){ /* find the non-eliminated candidate... */
    if(!Eliminated[i]){
      return i; /*Coombs winner*/
    }
  }
  return(-1); /*error*/
}

EMETH Approval(edata *E   /* canddt with most-approvals wins */
){ /* side effects: ApprovalVoteCount[], ApprovalWinner */
  int i,j,x;
  ZeroIntArray( E->NumCands, (int*)ApprovalVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
        ApprovalVoteCount[j] += E->Approve[x+j];
     }
  }		
  RandomlyPermute( E->NumCands, RandCandPerm );
  ApprovalWinner = ArgMaxUIntArr( E->NumCands, (uint*)ApprovalVoteCount, (int*)RandCandPerm );  
  return(ApprovalWinner);
}

EMETH App2Runoff(edata *E    /*top-2-runoff, 1stRd=approval, 2nd round has fully-honest voting*/
){ /* side effects: ASecond */
  int i;
  uint offset, pwct=0, wct=0;
  ASecond = -1;
  if(ApprovalWinner<0) Approval(E);
  ASecond = Arg2MaxUIntArr( E->NumCands, ApprovalVoteCount, (int*)RandCandPerm, ApprovalWinner );
  assert(ASecond>=0);
  for(i=0; i<(int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+ApprovalWinner] > E->PerceivedUtility[offset+ASecond] ){
       pwct++;
    }else if( E->PerceivedUtility[offset+ApprovalWinner] < E->PerceivedUtility[offset+ASecond] ){
      wct++;
    }
  }		
  if( pwct > wct || (pwct == wct && RandBool()) ) return(ApprovalWinner);
  return(ASecond);
}

EMETH HeitzigDFC(edata *E){ /*winner of honest runoff between ApprovalWinner and RandomBallot winner*/
  int i, Rwnr;
  uint offset, pwct=0, wct=0;
  Rwnr = ArgMaxRealArr( E->NumCands, 
               E->PerceivedUtility + RandInt(E->NumVoters)*E->NumCands, (int*)RandCandPerm );
  if(ApprovalWinner<0) Approval(E);
  for(i=0; i<(int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+ApprovalWinner] > E->PerceivedUtility[offset+Rwnr] ){
       pwct++;
    }else if( E->PerceivedUtility[offset+ApprovalWinner] < E->PerceivedUtility[offset+Rwnr] ){
      wct++;
    }
  }		
  if( pwct > wct || (pwct == wct && RandBool()) ) return(ApprovalWinner);
  return(Rwnr);
}

EMETH HeitzigLFC(edata *E){ /*random canddt who is not "strongly beat" wins; y strongly beats x if Approval(y)>Approval(x) and less than Approval(x)/2 voters prefer x>y.*/
  /*Side effects: Eliminated[] */
  int i,j;
  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(j=E->NumCands -1; j>=0; j--){
    for(i=E->NumCands -1; i>=0; i--){
      if(ApprovalVoteCount[j] > ApprovalVoteCount[i]){
	if(2*E->DefeatsMatrix[i*E->NumCands +j] < (int)ApprovalVoteCount[i]){
          /*Candidate i is "strongly beat"*/
	  Eliminated[i] = TRUE;
	}
      }
    }
  }
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){ /* find random non-eliminated candidate... */
    j = RandCandPerm[i];
    if(!Eliminated[j]){
      return j; /*winner*/
    }
  }
  return(-1); /*error*/
}

/***
Brian Olson's IRNR system:
1. get from each voter a range-style vote-vector, entries in range [-10, +10].
2. normalize so sum of squares equals 1.   (Or in variant, use other power than 2.)
3. sum the normalized vote vectors.
4. disqualify choice with lowest sum.
5. re-normalize (effectively redistributing voters' votes to their remaining choices).
6. back to step 3 until only 1 candidate remains.
7. It wins.
***/
real IRNRPOWER=2.0;
EMETH IRNR(edata *E  /*Brian Olson's voting method described above*/
){ /* side effects: Eliminated[], SumNormedRatings[]*/
  int rd,i,j,x,r,loser;
  real s,t,minc;

  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(rd=E->NumCands; rd>1; rd--){
    ZeroRealArray( E->NumCands, SumNormedRating );
    for(i=0; i<(int)E->NumVoters; i++){
      x = i*E->NumCands;
      s = 0.0;
      for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	t = E->Score[x+j] - 0.5;
	if(t < 0.0) t = -t;
	s += pow(t, IRNRPOWER);
      }
      if(s>0.0){ 
	s = pow(s, -1.0/IRNRPOWER);
	for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	  SumNormedRating[j] += s * E->Score[x+j];
	}
      }
    }
    RandomlyPermute( E->NumCands, RandCandPerm );
    loser = -1;
    minc = HUGE;
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      if(!Eliminated[r]){
	if( SumNormedRating[r] < minc ){
	  minc = SumNormedRating[r];
	  loser = r;
	}
      }
    }
    assert(loser>=0);
    Eliminated[loser] = TRUE;
  }
  for(i=E->NumCands -1; i>=0; i--){ /* find random non-eliminated candidate... */
    j = RandCandPerm[i];
    if(!Eliminated[j]){
      return j; /*winner*/
    }
  }
  return(-1); /*error*/
}

/* Like IRNR based on squares, but when "renomalizing" it now does a TWO-parameter
 * linear transformation so mean=0 and variance=1. */
EMETH IRNRv(edata *E  /*Brian Olson's voting method but with 2-param renorm*/
){ /* side effects: Eliminated[], SumNormedRatings[]*/
  int rd,i,j,x,r,loser,ct;
  real s,t,minc,avg;

  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(rd=E->NumCands; rd>1; rd--){
    ZeroRealArray( E->NumCands, SumNormedRating );
    for(i=0; i<(int)E->NumVoters; i++){
      x = i*E->NumCands;
      s = 0.0; ct=0;
      for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	s += E->Score[x+j]; ct++;
      }
      assert(ct>0);
      avg = s/ct; /*mean*/
      s = 0.0;
      for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	t = E->Score[x+j] - avg;
	s += t*t;
      }
      if(s>0.0){ 
	s = 1.0/sqrt(s);
	for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	  SumNormedRating[j] += s * (E->Score[x+j] - avg);
	}
      }
    }
    RandomlyPermute( E->NumCands, RandCandPerm );
    loser = -1;
    minc = HUGE;
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      if(!Eliminated[r]){
	if( SumNormedRating[r] < minc ){
	  minc = SumNormedRating[r];
	  loser = r;
	}
      }
    }
    assert(loser>=0);
    Eliminated[loser] = TRUE;
  }
  for(i=E->NumCands -1; i>=0; i--){ /* find random non-eliminated candidate... */
    j = RandCandPerm[i];
    if(!Eliminated[j]){
      return j; /*winner*/
    }
  }
  return(-1); /*error*/
}

/* Like IRNR based on squares, but when "renomalizing" it now does a TWO-parameter
 * linear transformation so max=1 and min=0. */
EMETH IRNRm(edata *E  /*Brian Olson's voting method but with 2-param renorm*/
){ /* side effects: Eliminated[], SumNormedRatings[]*/
  int rd,i,j,x,r,loser;
  real s,t,minc,mx,mn;

  FillBoolArray(E->NumCands, Eliminated, FALSE);
  for(rd=E->NumCands; rd>1; rd--){
    ZeroRealArray( E->NumCands, SumNormedRating );
    for(i=0; i<(int)E->NumVoters; i++){
      x = i*E->NumCands;
      mx = -HUGE; mn = HUGE;
      for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	t = E->Score[x+j]; 
	if(t<mn) mn=t;
	if(t>mx) mx=t;
      }
      s = mx-mn;
      if(s>0.0){ 
	s = 1.0/s;
	for(j=E->NumCands -1; j>=0; j--) if(!Eliminated[j]){
	  SumNormedRating[j] += s * (E->Score[x+j] - mn);
	}
      }
    }
    RandomlyPermute( E->NumCands, RandCandPerm );
    loser = -1;
    minc = HUGE;
    for(j=E->NumCands -1; j>=0; j--){
      r = RandCandPerm[j];
      if(!Eliminated[r]){
	if( SumNormedRating[r] < minc ){
	  minc = SumNormedRating[r];
	  loser = r;
	}
      }
    }
    assert(loser>=0);
    Eliminated[loser] = TRUE;
  }
  for(i=E->NumCands -1; i>=0; i--){ /* find random non-eliminated candidate... */
    j = RandCandPerm[i];
    if(!Eliminated[j]){
      return j; /*winner*/
    }
  }
  return(-1); /*error*/
}

EMETH MCA(edata *E  /*canddt with most-2approvals wins if gets >50%, else regular approval-winner wins*/
){ /* side effects: MCAVoteCount[] */
  int i,j,x,winner;
  if(ApprovalWinner<0) Approval(E);
  ZeroIntArray( E->NumCands, (int*)MCAVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
        MCAVoteCount[j] += E->Approve2[x+j];
     }
  }		
  winner = ArgMaxUIntArr( E->NumCands, MCAVoteCount, (int*)RandCandPerm );  
  if(2*MCAVoteCount[winner] > E->NumVoters) return(winner);
  return(ApprovalWinner);
}

/***Chris Benham:
This is my suggested rule for picking the two second-round qualifiers
(for use in a top2 runoff) using Approval for the first round:
  "The two finalists are the Approval winner A; and of those candidates
X who would have more approval than A if ballots that make no
approval distinction between A and X were altered to exclusively
approve X, the second qualifier is the X who is most approved on
ballots that don't approve A."
[There might not be a "second qualifier"; then we just make the approval winner win.]
   This is the "instant" version of a (delayed) 2-round system I
suggested a while back.
  This is a big improvement over simply promoting the two most approved
candidates because (a)it kills the strategy of richer factions trying
to capture both runoff spots by running a pair of clones, and (b) it
makes the method much less vulnerable to "turkey raising" strategy
(because voters in the first round can't have their votes do both of
promote one of their sincerely approved candidates and also
a "turkey").
***/
EMETH Benham2AppRunoff(edata *E  /*description above*/
){ 
  int i,j,y,r,maxc;
  uint jct, awct, offset;
  if(ApprovalWinner<0) Approval(E);
  RandomlyPermute( E->NumCands, RandCandPerm );
  maxc = -BIGINT;
  j = -1;
  for(i=0; i<(int)E->NumCands; i++){
    r = RandCandPerm[i];
    y = PairApproval[ApprovalWinner*E->NumCands +r];
    if( ApprovalVoteCount[r] + y > ApprovalVoteCount[ApprovalWinner] ){
      if( (int)ApprovalVoteCount[r] - y > maxc ){
	maxc = ApprovalVoteCount[r] - y;
	j = r;
      }
    }
  }
  if(j<0) return(ApprovalWinner);
  /* now for honest runoff between ApprovalWinner and j */
  awct=0;
  jct=0;
  for(i=0; i<(int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+ApprovalWinner] > E->PerceivedUtility[offset+j] ){
      awct++;
    }else if( E->PerceivedUtility[offset+ApprovalWinner] < E->PerceivedUtility[offset+j] ){
      jct++;
    }
  }		
  if( awct > jct || (awct == jct && RandBool()) ) return(ApprovalWinner);
  return(j);
}

EMETH Benham2AppRunB(edata *E  /*same as above, except if no 2nd qualifier then just use top 2 approval finishers; always do runoff*/
){ /* side effects: PairApproval[], ASecond */
  uint offset, pwct=0, wct=0;
  int i,j,y,r,maxc;
  uint jct, awct;
  if(ApprovalWinner<0) Approval(E);
  RandomlyPermute( E->NumCands, RandCandPerm );
  maxc = -BIGINT;
  j = -1;
  for(i=0; i<(int)E->NumCands; i++){
    r = RandCandPerm[i];
    y = PairApproval[ApprovalWinner*E->NumCands +r];
    if( ApprovalVoteCount[r] + y > ApprovalVoteCount[ApprovalWinner] ){
      if( (int)ApprovalVoteCount[r] - y > maxc ){
	maxc = ApprovalVoteCount[r] - y;
	j = r;
      }
    }
  }
  if(j<0){ /* no 2nd qualifier*/
    ASecond = -1;
    ASecond = Arg2MaxUIntArr( E->NumCands, ApprovalVoteCount, (int*)RandCandPerm, ApprovalWinner );
    assert(ASecond>=0);
    for(i=0; i<(int)E->NumVoters; i++){
      offset = i*E->NumCands;
      if( E->PerceivedUtility[offset+ApprovalWinner] > E->PerceivedUtility[offset+ASecond] ){
	pwct++;
      }else if( E->PerceivedUtility[offset+ApprovalWinner] < E->PerceivedUtility[offset+ASecond] ){
	wct++;
      }
    }		
    if( pwct > wct || (pwct == wct && RandBool()) ) return(ApprovalWinner);
    return(ASecond);
  }
  /* now for honest runoff between ApprovalWinner and j */
  awct=0;
  jct=0;
  for(i=0; i<(int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+ApprovalWinner] > E->PerceivedUtility[offset+j] ){
      awct++;
    }else if( E->PerceivedUtility[offset+ApprovalWinner] < E->PerceivedUtility[offset+j] ){
      jct++;
    }
  }		
  if( awct > jct || (awct == jct && RandBool()) ) return(ApprovalWinner);
  return(j);
}

EMETH CondorcetApproval(edata *E  /*Condorcet winner if exists, else use Approval*/
){ 
  if(ApprovalWinner<0) Approval(E);
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CondorcetWinner >= 0) return CondorcetWinner;
  return ApprovalWinner;
}

EMETH Range(edata *E    /* canddt with highest average Score wins */
){ /* side effects:   RangeVoteCount[], RangeWinner  */
  int i,j,x;
  ZeroRealArray( E->NumCands, RangeVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
        RangeVoteCount[j] += E->Score[x+j];
     }
  }		
  RangeWinner = ArgMaxRealArr( E->NumCands, RangeVoteCount, (int*)RandCandPerm );  
  return(RangeWinner);
}

EMETH RangeN(edata *E /*highest average rounded Score [rded to integer in 0..RangeGranul-1] wins*/
){ /* side effects:   RangeNVoteCount[], uses global integer RangeGranul  */
  int i,j,x,winner;
  assert(RangeGranul>=2);
  assert(RangeGranul<=10000000);
  ZeroIntArray( E->NumCands, (int*)RangeNVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     for(j=E->NumCands -1; j>=0; j--){
       RangeNVoteCount[j] += (uint)( (E->Score[x+j])*(RangeGranul-0.0000000001) );
     }
  }		
  winner = ArgMaxUIntArr( E->NumCands, RangeNVoteCount, (int*)RandCandPerm );  
  return(winner);
}

EMETH Range2Runoff(edata *E    /*top-2-runoff, 1stRd=range, 2nd round has fully-honest voting*/
){ /* side effects: RSecond */
  int i;
  uint offset, pwct=0, wct=0;
  RSecond = -1;
  if(RangeWinner<0) Range(E);
  RandomlyPermute( E->NumCands, RandCandPerm );
  RSecond = Arg2MaxRealArr( E->NumCands, RangeVoteCount, (int*)RandCandPerm, RangeWinner );
  assert(RSecond>=0);
  for(i=0; i < (int)E->NumVoters; i++){
    offset = i*E->NumCands;
    if( E->PerceivedUtility[offset+RangeWinner] > E->PerceivedUtility[offset+RSecond] ){
       pwct++;
    }else if( E->PerceivedUtility[offset+RangeWinner] < E->PerceivedUtility[offset+RSecond] ){
      wct++;
    }
  }		
  if( pwct > wct || (pwct == wct && RandBool()) ) return(RangeWinner);
  return(RSecond);
}

EMETH ContinCumul(edata *E    /* Renormalize scores so sum(over canddts)=1; then canddt with highest average Score wins */
){ /* side effects:   CCumVoteCount[] */
  int i,j,x,winner;
  real sum;
  ZeroRealArray( E->NumCands, CCumVoteCount );
  for(i=0; i<(int)E->NumVoters; i++){
     x = i*E->NumCands;
     sum = 0.0;
     for(j=E->NumCands -1; j>=0; j--){
        sum += E->Score[x+j];
     }
     if(sum > 0.0) sum = 1.0/sum; else sum=0.0;
     for(j=E->NumCands -1; j>=0; j--){
        CCumVoteCount[j] += sum * E->Score[x+j];
     }
  }		
  winner = ArgMaxRealArr( E->NumCands, CCumVoteCount, (int*)RandCandPerm );  
  return(winner);
}

EMETH TopMedianRating(edata *E    /* canddt with highest median Score wins */
){ /* side effects:   MedianRating[], CScoreVec[]  */
  int i,j,x,winner;
  for(j=E->NumCands -1; j>=0; j--){
    for(i=E->NumVoters -1; i>=0; i--){
      x = i*E->NumCands + j;
      CScoreVec[i] = E->Score[x];
    }
    MedianRating[j] = TwiceMedianReal( E->NumVoters, CScoreVec );
  }
  winner = ArgMaxRealArr( E->NumCands, MedianRating, (int*)RandCandPerm );  
  return(winner);
}

EMETH LoMedianRank(edata *E    /* canddt with best median ranking wins */
){ /* side effects:   MedianRank[], CRankVec[]  */
  int i,j,x,winner;
  for(j=E->NumCands -1; j>=0; j--){
    for(i=E->NumVoters -1; i>=0; i--){
      x = i*E->NumCands + j;
      CRankVec[i] = E->CandRankings[x];
    }
    MedianRank[j] = TwiceMedianInt( E->NumVoters, CRankVec );
    assert( MedianRank[j] >= 0 );
    assert( MedianRank[j] <= 2*((int)E->NumCands - 1) );
  }
  winner = ArgMinIntArr( E->NumCands, MedianRank, (int*)RandCandPerm );  
  return(winner);
}

/* Tideman ranked pairs with Honest voters.
 * Tideman = Condorcet variant in which you
 * pick the A>B comparison with the largest margin and "lock it in".
 * Then you  pick the next largest one available
 * ("available" means: not already used and not creating a cycle),
 * and continue on.  This creates an ordering of the candidates. The
 * topmost in the ordering wins.  It is a bit tricky to spot
 * the cycles as we go...
 * (This code based on email from Blake Cretney bcretney@postmark.net
 * and runs in order N^4 time with N candidates, i.e slow.  It is possible to speed it
 * up to Otilde(N^2) steps with the use of fast data structures...):
 ******************************************/
EMETH TidemanRankedPairs(edata *E  /*lock in comparisons with largest margins not yielding cycle*/
){  /*side effects: Tpath[] is used as a changeable copy of MarginsMatrix.*/
  int i,r,j,winner;
  if(CWSPEEDUP && CondorcetWinner >= 0) return CondorcetWinner;
  CopyIntArray(SquareInt(E->NumCands), E->MarginsMatrix, Tpath);
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=0; i < (int)E->NumCands; i++){ Tpath[i*E->NumCands +i]=BIGINT; }
  /* Whenever a victory
   * is locked in, the appropriate Tpath[] cell is set to BIGINT.
   * pi,pj are used with randperm to give precedence to victories higher
   * in the random permutation (when tie-breaking).
   * This loop finds the next pair (i,j) to lock in: ********/
  for(;;){
    int maxp,oi,oj,pi,pj;
    maxp = -BIGINT; j = -1;
    for(pi=0; pi < (int)E->NumCands; pi++){
      oi=RandCandPerm[pi]; 
      for(pj=pi+1; pj < (int)E->NumCands; pj++){
	oj=RandCandPerm[pj];
	if(   Tpath[oi*E->NumCands +oj]!=BIGINT
	      && Tpath[oj*E->NumCands +oi]!=BIGINT){/*not locked-out*/
	  if(Tpath[oi*E->NumCands +oj]>maxp){ maxp=Tpath[oi*E->NumCands +oj]; i=oi; j=oj; }
	  if(Tpath[oj*E->NumCands +oi]>maxp){ maxp=Tpath[oj*E->NumCands +oi]; i=oj; j=oi; }
	}
      }
    }
    if(maxp == -BIGINT) break;
    assert(j>=0);
    /********* print the pair and its margin:
	       printf("(%d %d) %d\n",i,j,maxp);
    ***********************/
    /*lock in the pair and clobber future no-good pairs:*/
    for(oi=0; oi < (int)E->NumCands; oi++) for(oj=0; oj < (int)E->NumCands; oj++){
      if(Tpath[oi*E->NumCands +i]==BIGINT && Tpath[j*E->NumCands +oj]==BIGINT){
	Tpath[oi*E->NumCands +oj]=BIGINT;
      }
    }
  }
  /* The above code assumes that pairels has been set properly.
   * Tpath[*E->NumCands +] ends up with the winning row having all cells
   * set to BIGINT.  In fact, a complete ranking is given,
   * where Tpath[i*E->NumCands +j]==BIGINT means that i is
   * ranked over j (where i!=j).   So to find the winner: ****/
  winner = -99;
  for(i=0; i < (int)E->NumCands; i++){
    r = RandCandPerm[i];
    for(j=0; j < (int)E->NumCands; j++){
      if(Tpath[r*E->NumCands +j] != BIGINT) break;
    }
    if(j >= (int)E->NumCands){ winner=r; break; }
  }
  assert(winner >= 0);
  return(winner);
}

EMETH HeitzigRiver(edata *E /*http://lists.electorama.com/pipermail/election-methods-electorama.com/2004-October/013971.html*/
){ /* side effects: Hpotpar[], Hpar[], Hrot[] */
  int r,z,i,j,k,pp,oldroot,newroot,maxc;
  if(CWSPEEDUP && CondorcetWinner >= 0) return CondorcetWinner;
  /* find potential (and actual) parents: */
  for(i=0; i < (int)E->NumCands; i++){
    maxc = -BIGINT; pp = -1;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(j=0; j < (int)E->NumCands; j++){
      r = RandCandPerm[j];
      if(r!=i && E->MarginsMatrix[r*E->NumCands +i]>maxc){
        maxc = E->MarginsMatrix[r*E->NumCands +i];
        pp = r;
      }
    }
    assert(pp>=0);
    Hpotpar[i] = pp;
    Hpar[i] = i;
    Hroot[i] = i;
  }
  for(z = E->NumCands; ; ){ /* loop that adds arcs: */
    /* scan tree roots to find best new arc to add: */
    maxc = -BIGINT; k = -1;
    RandomlyPermute( E->NumCands, RandCandPerm );
    for(i=0; i < (int)E->NumCands; i++){
      r = RandCandPerm[i];
      if(Hpar[r]==r){ /* tree root */
	pp = Hpotpar[r];
	if(maxc < E->MarginsMatrix[pp*E->NumCands +r]){
	  maxc = E->MarginsMatrix[pp*E->NumCands +r];
	  k = r;
	}
      }
    }/*end for(i)*/
    assert(k>=0);
    /* add it: */
    pp = Hpotpar[k];
    Hpar[k] = pp;
    /* update roots: */
    newroot = Hroot[pp];
    z--;
    if(z<=1) break;
    oldroot = Hroot[k];
    for(i=0; i < (int)E->NumCands; i++){
      if(Hroot[i] == oldroot) Hroot[i] = newroot;
    }
    /* update potential parent of newroot: */
    maxc = -BIGINT; k = -1;
    for(j=0; j < (int)E->NumCands; j++){
      r = RandCandPerm[j];
      if(Hroot[r]!=newroot && E->MarginsMatrix[r*E->NumCands +newroot]>maxc){
	maxc = E->MarginsMatrix[r*E->NumCands +newroot];
	k = r;
      }
    }
    assert(k>=0);
    Hpotpar[newroot] = k;
  }
  return newroot;
}

EMETH DMC(edata *E  /* eliminate least-approved candidate until unbeaten winner exists */
){ /* side effects: LossCount[] */
  int i,j,t;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(CWSPEEDUP && CondorcetWinner >=0 ) return(CondorcetWinner);
  for(i=E->NumCands -1; i>=0; i--){ 
    t=0;
    for(j=E->NumCands -1; j>=0; j--) if(E->MarginsMatrix[j*E->NumCands + i]>0){ t++; }
    LossCount[i] = t;
  }
  RandomlyPermute( E->NumCands, RandCandPerm );
  IntPermShellSortDown(  E->NumCands, (int*)RandCandPerm,  (int*)ApprovalVoteCount );
  for(i=E->NumCands-1; i>0; i--){
    if( LossCount[i] <= 0 ){  return(i); /*winner*/ }
    for(j=0; j<i; j++){ /* eliminate i and update Losscount[] */
      t = E->MarginsMatrix[i*E->NumCands + j];
      if(t>0){ LossCount[j] --; }
    }
  }
  return(i);
}

void BSbeatDFS( int x, int diff, bool Set[], bool OK[], int Mat[], int N ){
  int i;
  for(i=0; i<N; i++) if(OK[i] && i!=x){
    if( Mat[i*N+x]>=diff ){
      if( !Set[i] ){
	Set[i] = TRUE;
	BSbeatDFS( i, diff, Set, OK, Mat, N );
      }
    }
  }
}

EMETH BramsSanverPrAV(edata *E  /*SJ Brams & MR Sanver: Voting Systems That Combine Approval and Preference,2006*/
){ /* side effects: MajApproved[] */
  int i,j,winner,ctm,t,maxt,CopeWinr,r;
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(ApprovalWinner<0) Approval(E);
  ctm=0;
  for(i=E->NumCands -1; i>=0; i--){
    if(2*ApprovalVoteCount[i] > E->NumVoters){
      MajApproved[i] = TRUE; ctm++;
    }else{  MajApproved[i] = FALSE;  }
  }
  if(ctm<=1){
    /*if exactly 0 or 1 canddt majority-approved, ApprovalWinner wins*/
    return ApprovalWinner;
  }
  assert(ctm>=2);  /*Two or more majority-approved:*/
  for(i=E->NumCands -1; i>=0; i--) if(MajApproved[i]){
    for(j=E->NumCands -1; j>=0; j--) if(j!=i && MajApproved[j]){
      if( E->MarginsMatrix[i*E->NumCands + j] <= 0 ){ goto BADi; }
    }
    winner = i; /*beats-all winner among >=2 majority-approved canddts wins*/
    return winner;
  BADi: ;
  }
  /*now Brams&Sanver want the most-approved member of the 
   *top-cycle among the majority-approved canddts, to win: */
  maxt = 0;
  CopeWinr = -1;
  for(i=E->NumCands -1; i>=0; i--) if(MajApproved[i]){
    t = 0;
    for(j=E->NumCands -1; j>=0; j--) if(j!=i && MajApproved[j]){
      if( E->MarginsMatrix[i*E->NumCands + j] > 0 ){ t++; }
    }
    if(t >= maxt){ maxt=t; CopeWinr=i; }
  }
  assert(CopeWinr >= 0);
  assert(MajApproved[CopeWinr]);
  FillBoolArray(E->NumCands, BSSmithMembs, FALSE);
  BSSmithMembs[CopeWinr] = TRUE;
  BSbeatDFS(CopeWinr, 1, BSSmithMembs, MajApproved, E->MarginsMatrix, E->NumCands);
  assert(BSSmithMembs[CopeWinr]);
  RandomlyPermute( E->NumCands, RandCandPerm );
  winner = -1; 
  maxt = 0;
  for(i=E->NumCands -1; i>=0; i--){
    r = RandCandPerm[i];
    if(BSSmithMembs[r] && ApprovalVoteCount[r]>maxt){ maxt=ApprovalVoteCount[r]; winner=r; }
  }
  return winner;
}

EMETH MDDA(edata *E  /* approval-count winner among canddts not majority-defeated (or all, if all maj-defeated) */
){ /* side effects: MDdisquald[] */
  int i,j,r,dqct,thresh,maxc,winner;
  /*if(CWSPEEDUP && CondorcetWinner >=0 ) return(CondorcetWinner); valid??*/
  if(CopeWinOnlyWinner<0) BuildDefeatsMatrix(E);
  if(ApprovalWinner<0) Approval(E);
  dqct=0;
  thresh = (E->NumVoters)/2;
  for(i=E->NumCands -1; i>=0; i--){ 
    MDdisquald[i] = FALSE;
    for(j=E->NumCands -1; j>=0; j--) if(E->DefeatsMatrix[j*E->NumCands + i] > thresh){
      MDdisquald[i] = TRUE; 
      dqct++;
      break;
    }
  }
  if( dqct >= (int)E->NumCands ){ /*If all disqualified, none are. */
    return(ApprovalWinner);
  }
  winner = -1;
  maxc = 0;
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=E->NumCands -1; i>=0; i--){
    r = RandCandPerm[i];
    if(!MDdisquald[r] && (int)ApprovalVoteCount[r] >= maxc){
      maxc=ApprovalVoteCount[r];
      winner=r; 
    }
  }
  return(winner);
}

/***Forest Simmons 27 Feb 2007:
UncAAO stands for Uncovered, Approval, Approval Opposition.  
Here's how it works:
For each candidate X, if X is uncovered, then let f(X)=X,
else let f(X) be the candidate against which X has the least approval 
opposition, among those candidates that cover X.
["Approval opposition" of X against Y is the number of ballots on which 
X but not Y is approved.]
Start with the approval winner A and apply the function f repeatedly 
until the output equals the input.  This "fixed point" of f is the 
method winner.
  This method requires a tally of both pairwise approval and pairwise 
ordinal information, but both are efficiently summable in N by N 
matrices, where N is the number of candidates.
  This method (UncAAO) is monotone, clone free, and always picks from the 
uncovered set, which is a subset of the Smith set.
  Zero info strategy is sincere.
Even perfect info incentives for burial and betrayal are practically  nil.
   As near as I can tell, the only bad thing about the method is the 
"tyranny of the majority" problem shared by most, if not all, deterministic methods.
*****/
EMETH UncAAO(edata *E  
){ 
  int i,j,ff,AppOpp,MnAO,r,winner;
  if(ApprovalWinner<0) Approval(E);
  if(RandomUncoveredMemb<0) UncoveredSet(E);
  RandomlyPermute( E->NumCands, RandCandPerm );
  for(i=(int)E->NumCands -1; i>=0; i--){
    if( UncoveredSt[i] ){
      UncAAOF[i] = i;
    }else{
      MnAO = BIGINT;
      ff = -1;
      for(j=(int)E->NumCands -1; j>=0; j--){
	r = RandCandPerm[j];
	if( CoverMatrix[r*E->NumCands+i]  ){
	  AppOpp = ApprovalVoteCount[r] - PairApproval[r*E->NumCands+i];
	  if(AppOpp < MnAO){ MnAO = AppOpp; ff = r; }
	}
      }
      assert(ff >= 0);
      UncAAOF[i] = ff;
    }
  }
  winner = ApprovalWinner;
  do{ winner =  UncAAOF[winner]; }until( winner == UncAAOF[winner] );
  return(winner);
}


/*  Woodall-DAC:
1. Rank-order ballots.
2. A voter "acquiesces"' to a set of candidates if he does not rank any
candidate outside the set higher than any inside the set.
(Every voter acquiesces to the full candidate-set.)
3. Sort all possible sets from most acquiescing voters to fewest.
Going down the list, disqualify every candidate not found in each set (i.e. take
set intersections) unless that
would disqualify all remaining candidates (i.e. would result in the empty set).
Continue until only one candidate
is not disqualified; he is the winner.
*****************************************/
EMETH WoodallDAC(edata *E  /*Woodall: Monotonocity of single-seat preferential election rules, Discrete Applied Maths 77 (1997) 81-98.*/
){ /*side effects: WoodHashCount[],  WoodHashSet[], WoodSetPerm[] */
  /* Hash Tab entries contain counter and set-code which is a single machine word. */
   int v,c,k,r;
   uint offset, s, x, h, numsets;
   if( E->NumCands > 4*sizeof(uint) ){
     printf("WoodallDAC: too many candidates %d to use machine words(%d) to represent sets\n", 
	    E->NumCands, 
	    (int)(4*sizeof(uint)) );
     printf("You could rewrite the code to use uint64s to try allow up to 64 canddts\n");
     exit(EXIT_FAILURE);
   }
   for(v=ARTINPRIME-1; v>=0; v--){ WoodHashCount[v] = 0; WoodHashSet[v] = 0; }
   for(v=E->NumVoters-1; v>=0; v--){
      s = 0;
      offset = v*E->NumCands;
      for(c=0; c < (int)E->NumCands; c++){
         s |= (1<<(E->TopDownPrefs[offset+c]));
         h = s%ARTINPRIME;
	 assert( !EmptySet(s) );
	 for(;;){
	   x = WoodHashSet[h];
	   if( EmptySet(x) ){ /* insert new set s into hash table */
	     WoodHashSet[h] = s;  WoodHashCount[h] = 1;  break;
	   }else if( x==s ){ /* already there so increment count */
	     WoodHashCount[h]++;  break;
	   }
	   h++; /* hash table collision so walk ("linear probing") */
	 }
      }
   }
   numsets=0;
   for(v=ARTINPRIME-1; v>=0; v--){
     if( !EmptySet(WoodHashSet[v]) ){
       WoodSetPerm[numsets] = v;
       numsets++;
     }
   }
   assert(numsets <= E->NumCands * E->NumVoters);
   IntPermShellSortDown( numsets, (int*)WoodSetPerm, (int*)WoodHashCount );
   s = WoodHashSet[WoodSetPerm[0]];
   assert( !EmptySet(s) );
   if(SingletonSet(s)){  goto DONE;  }
   for(k=1; k<(int)numsets; k++){ /*decreasing order!*/
      h = WoodSetPerm[k];
      x = s & WoodHashSet[h];
      if(!EmptySet(x)){
	s = x;
	if(SingletonSet(s)){  goto DONE;  }
      }
   }
 DONE: ;
   /*printf("C%d/%d\n", CardinalitySet(s), E->NumCands);*/
   /*It is extremely rare that s is not a singleton set.  In fact may never happen.*/
   RandomlyPermute( E->NumCands, RandCandPerm );
   for(k=E->NumCands -1; k>=0; k--){
     r = RandCandPerm[k];
     if( (s>>r)&1 ){
       return r; /* return random set-element */
     }
   }
   return(-1); /*failure*/
}

/******** uber-routines which package many voting methods into one: **********/

void PrintMethName( int WhichMeth, bool Padding ){
  switch(WhichMeth){
  case(0) : printf("SociallyBest"); if(Padding) PrintNSpaces(3); break;
  case(1) : printf("SociallyWorst"); if(Padding) PrintNSpaces(2); break;
  case(2) : printf("RandomWinner"); if(Padding) PrintNSpaces(3); break;
  case(3) : printf("Plurality"); if(Padding) PrintNSpaces(5); break;
  case(4) : printf("Borda"); if(Padding) PrintNSpaces(8); break;
  case(5) : printf("IRV"); if(Padding) PrintNSpaces(10); break;
  case(6) : printf("Approval"); if(Padding) PrintNSpaces(9); break;
  case(7) : printf("Range"); if(Padding) PrintNSpaces(8); break;
  case(8) : printf("SmithSet"); if(Padding) PrintNSpaces(8); break;
  case(9) : printf("SchwartzSet"); if(Padding) PrintNSpaces(6); break;
  case(10) : printf("AntiPlurality"); if(Padding) PrintNSpaces(3); break;  
  case(11) : printf("Top2Runoff"); if(Padding) PrintNSpaces(3); break;
    /****** above methods were "Core"; below are optional *****/
  case(12) : printf("CondorcetLR"); if(Padding) PrintNSpaces(7); break;
  case(13) : printf("SimpsonKramer"); if(Padding) PrintNSpaces(3); break;
  case(14) : printf("Bucklin"); if(Padding) PrintNSpaces(7); break;  
  case(15) : printf("Copeland"); if(Padding) PrintNSpaces(6); break; 
  case(16) : printf("SimmonsCond"); if(Padding) PrintNSpaces(3); break;
  case(17) : printf("SmithIRV"); if(Padding) PrintNSpaces(8); break;  
  case(18) : printf("BTRIRV"); if(Padding) PrintNSpaces(7); break;  
  case(19) : printf("DMC"); if(Padding) PrintNSpaces(10); break;  
  case(20) : printf("Dabagh"); if(Padding) PrintNSpaces(7); break;  
  case(21) : printf("VtForAgainst"); if(Padding) PrintNSpaces(3); break;  
  case(22) : printf("SchulzeBeatpaths"); if(Padding) PrintNSpaces(0); break;
  case(23) : printf("PlurIR"); if(Padding) PrintNSpaces(9); break;
  case(24) : printf("Black"); if(Padding) PrintNSpaces(8); break;
  case(25) : printf("RandomBallot"); if(Padding) PrintNSpaces(1); break;
  case(26) : printf("RandomPair"); if(Padding) PrintNSpaces(3); break;
  case(27) : printf("NansonBaldwin"); if(Padding) PrintNSpaces(0); break;
  case(28) : printf("Nauru"); if(Padding) PrintNSpaces(8); break;  
  case(29) : printf("TopMedianRating"); if(Padding) PrintNSpaces(0); break;  
  case(30) : printf("LoMedianRank"); if(Padding) PrintNSpaces(3); break;  
  case(31) : printf("RaynaudElim"); if(Padding) PrintNSpaces(4); break;  
  case(32) : printf("ArrowRaynaud"); if(Padding) PrintNSpaces(3); break;  
  case(33) : printf("Sinkhorn"); if(Padding) PrintNSpaces(7); break;  
  case(34) : printf("KeenerEig"); if(Padding) PrintNSpaces(7); break;  
  case(35) : printf("MDDA"); if(Padding) PrintNSpaces(12); break;  
  case(36) : printf("VenzkeDisqPlur"); if(Padding) PrintNSpaces(2); break;  
  case(37) : printf("CondorcetApproval"); if(Padding) PrintNSpaces(0); break;  
  case(38) : printf("UncoveredSet"); if(Padding) PrintNSpaces(4); break;  
  case(39) : printf("BramsSanverPrAV"); if(Padding) PrintNSpaces(3); break;  
  case(40) : printf("Coombs"); if(Padding) PrintNSpaces(12); break;  
  case(41) : printf("Top3IRV"); if(Padding) PrintNSpaces(11); break;  
  case(42) : printf("ContinCumul"); if(Padding) PrintNSpaces(7); break;  
  case(43) : printf("IterCopeland"); if(Padding) PrintNSpaces(6); break;  
  case(44) : printf("HeitzigRiver"); if(Padding) PrintNSpaces(6); break;  
  case(45) : printf("MCA"); if(Padding) PrintNSpaces(15); break;  
  case(46) : printf("Range3"); if(Padding) PrintNSpaces(12); break;  
  case(47) : printf("Range10"); if(Padding) PrintNSpaces(11); break;  
  case(48) : printf("HeismanTrophy"); if(Padding) PrintNSpaces(2); break;  
  case(49) : printf("BaseballMVP"); if(Padding) PrintNSpaces(4); break;  
  case(50) : printf("App2Runoff"); if(Padding) PrintNSpaces(5); break;  
  case(51) : printf("Range2Runoff"); if(Padding) PrintNSpaces(3); break;  
  case(52) : printf("HeitzigDFC"); if(Padding) PrintNSpaces(5); break;  
  case(53) : printf("ArmytagePCSchulze"); if(Padding) PrintNSpaces(0); break;  
  case(54) : printf("Hay"); if(Padding) PrintNSpaces(12); break;  
  case(55) : printf("HeitzigLFC"); if(Padding) PrintNSpaces(5); break;  
  case(56) : printf("Benham2AppRunoff"); if(Padding) PrintNSpaces(0); break;  
  case(57) : printf("Benham2AppRunB"); if(Padding) PrintNSpaces(0); break;  
  case(58) : printf("WoodallDAC"); if(Padding) PrintNSpaces(5); break;  
  case(59) : printf("UncAAO"); if(Padding) PrintNSpaces(9); break;  
    /****** below methods are "Slow": *****/
  case(NumFastMethods+0) : printf("TidemanRankedPairs"); if(Padding) PrintNSpaces(0); break;  
  case(NumFastMethods+1) : printf("IRNR2"); if(Padding) PrintN