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, 3449035210, 2245417744, 5401519132, 5272055840, 4470731784, 7013199284, 6544305958, 10138992314, 15284191798, 9839916650, 8579138578, 13353765596, 12170187396, 18839898968 };
        static void Run()
        {
            double logLower = 0;
            double logUpper = 0;
            double log = 0;
            double epsilon = 0.000001;
            for (int idx = 1; idx <= 1125; idx++)
            {
                int HAPs = 0;
                for (int i = 1; i <= idx; i++)
                {
                    if (i != idx && idx % i == 0 && (((idx / i) - 1) % 2) == 0) HAPs++;
                }

                double badTransition = 0.5;
                double goodTransition = 0.5;
                double oddsBad = badTransition / goodTransition;

                // 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 extraEvidence = (badProb / (1 - badProb)) * (goodTransition / badTransition);
                                oddsBad *= extraEvidence;
                            }
                        }
                    }
                }

                badTransition = oddsBad / (1 + oddsBad);
                goodTransition = 1 - badTransition;

                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));
            }
            Console.ReadLine();
            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;
        }
    }
}