# Estimate the number of discrepancy 2 sequences

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;

namespace EDP
{
class Program
{
static double[] answers = { 1, 2, 4, 6, 12, 18, 28, 44, 88, 100, 152, 240, 370, 556, 882, 750, 1500, 2250, 2784, 4284, 6438, 6062, 9526, 14856, 22944, 26164, 39528, 35122, 54800, 80940, 81326, 122422, 244844, 234934, 356154, 309068, 388042, 589796, 900000, 813466, 1212450, 1837030, 1882194, 2921946, 4544342, 2274560, 3542738, 5495686, 8436986, 9597362, 11352364, 10876536, 16040144, 23626898, 24060696, 22332908, 34665316, 32813078, 48774494, 77333978, 79086932, 118815026, 181723488, 101003862, 202007724, 171112060, 175332604, 266556200, 388397590, 379477372, 355489350, 544231036, 685850834, 1035438526, 1600139990, 967443184, 1435672238, 1335563798, 1305891128, 1985108178, 2984972718, 2245417744, 3449035210, 5401519132, 5272055840, 4470731784, 7013199284, 6544305958, 10138992314, 15284191798, 9839916650, 8579138578, 13353765596, 12170187396, 18839898968, 17144676512, 26322363104, 39457407576 };
static void Run()
{
double logLower = 0;
double logUpper = 0;
double log = 0;
double epsilon = 0.000001;
for (int idx = 1; idx <= answers.Length - 1; idx++)
{
int HAPs = 0;
for (int i = 1; i <= idx; i++)
{
if (i != idx && idx % i == 0 && (((idx / i) - 1) % 2) == 0) HAPs++;
}

double goodTransition = 0.5;

// generalize again: foreach HAP f that constrains you now, go back a step
// in that HAP and foreach HAP g that constrains THAT step, calculate its
// bias towards or away from 0 by saying that f had wlog bias +1 just before there
// and doing some independence-assumptions prob on where its +1 went
// in (or not in) g

// (note: same loop as above, combine)
// check each f for factor of idx and whether its HAP would be on an odd step at idx
for (int f = 1; f <= idx; f++)
{
// does HAPf constrain us?
if (f != idx && idx % f == 0 && (((idx / f) - 1) % 2) == 0)
{
// it does, here's its last step
int idxm1 = idx - f;
// check each g for factor of idxm1 and whether its HAP would be on
// an odd step at idxm1
for (int g = 1; g <= idxm1; g++)
{
// does HAPg constrain us?
if (g != idxm1 && idxm1 % g == 0 && (((idxm1 / g) - 1) % 2) == 0)
{
// it does! set up some numbers
int k = ((idx / f) - 1) / 2;
int m = ((idxm1 / g) - 1) / 2;
int mult = lcm(f, g);
int s = ((2 * k - 1) * f) / mult;
int ns = 2 * m - s;
double sharedPlus = ((double)k) / (2 * k - 1);

// that's wrong, fix it!
double nonsharedPlus = 0.5;
if (ns > 0) nonsharedPlus = ((idxm1 - 2 * k - 1) / 2.0) / (idxm1 - 2 * k);

// this is wrong too, fix it!
double plusses = (s * sharedPlus + ns * nonsharedPlus) / (s + ns);

// this is wrong too, not sure how to fix it yet though
double plus2 = bin(2*m, m + 1) * Math.Pow(plusses, m + 1) * Math.Pow(1 - plusses, 2*m - (m + 1));
double plus0 = bin(2*m, m) * Math.Pow(plusses, m) * Math.Pow(1 - plusses, 2*m - (m));
double plusn2 = bin(2*m, m - 1) * Math.Pow(plusses, m - 1) * Math.Pow(1 - plusses, 2*m - (m - 1));
// normalize
double norm = plus2 + plus0 + plusn2;
plus2 /= norm;
plus0 /= norm;
plusn2 /= norm;

// apply to find the probability of a bad transition (away from 0)
// (wlog f had bias 1 so we use the -2 and a 50/50 of the 0 to get a +,
//  which is the direction away from 0 with bias 1)
double badProb = plusn2 + 0.5 * plus0;
}
}
}
}

double middleProb = goodTransition;
double edgeProb = (1 - goodTransition) / 2;

double p1 = Math.Pow(middleProb + edgeProb, HAPs) - Math.Pow(middleProb, HAPs);
double p2 = Math.Pow(middleProb, HAPs);
double multiplier = p1 + p1 + 2 * p2;
double lowerEV = Math.Log(2 * Math.Floor(Math.Exp(logLower + Math.Log(multiplier)) / 2 + epsilon));
double upperEV = Math.Log(2 * Math.Ceiling(Math.Exp(logUpper + Math.Log(multiplier)) / 2 - epsilon));
double EV = log + Math.Log(multiplier);
logLower = lowerEV;
logUpper = upperEV;
log = EV;
if (idx < answers.Length) Console.WriteLine(idx + ":\t" + HAPs + "\t" + Math.Exp(lowerEV) + "\t" + Math.Exp(upperEV) + "\t" + answers[idx] / Math.Exp(upperEV));
else Console.WriteLine(idx + ":\t" + HAPs + "\t" + Math.Exp(lowerEV) + "\t" + Math.Exp(upperEV) + "\t" + Math.Exp(EV));
}
return;
}
static int lcm(int A, int B)
{
int AB = A * B;
while (A > 0)
{
if (A < B)
{
int t = A;
A = B;
B = t;
}
A = A % B;
}
return AB / B;
}
static double hyp(int N, int m, int n, int k)
{
return bin(m, k) * bin(N - m, n - k) / bin(N, n);
}
static double bin(int N, int k)
{
double r = 1;
for (int i = 0; i < k; i++)
{
r *= N - i;
r /= i + 1;
}
return r;
}
static void Main(string[] args)
{
Run();
return;
}
}
}