/******************************************************************************/
/* */
/* COMPUTE UPPER BOUNDS OF MINIMA */
/* 04/22/10 (dkc) */
/* */
/* This C program computes upper bounds of minima where n is prime. */
/* */
/******************************************************************************/
#include <math.h>
void mul3shft(unsigned int *c, unsigned int b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
far unsigned int hipop[200]; // n values in second population
far unsigned int hipop3[2]; // n values in third population
far double max0[200],min0[200]; // maximum and minimum x values
far double max1[200],min1[200];
far double max2[200],min2[200];
far double max3[200],min3[200];
far double max4[200],min4[200];
far unsigned int max4i[200];
far unsigned int output[200];
far unsigned int hist0[200],hist1[200],hist2[200],hist3[200],hist4[200];
far unsigned int error[4]; // error array
far double outx[18]; // maximum and minimum x values
far unsigned int table[667]={ // prime look-up table
7, 11, 13, 17, 19, 23, 29, 31,
37, 41, 43, 47, 53, 59, 61, 67, 71, 73,
79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
131, 137, 139, 149, 151, 157, 163, 167, 173, 179,
181, 191, 193, 197, 199, 211, 223, 227, 229, 233,
239, 241, 251, 257, 263, 269, 271, 277, 281, 283,
293, 307, 311, 313, 317, 331, 337, 347, 349, 353,
359, 367, 373, 379, 383, 389, 397, 401, 409, 419,
421, 431, 433, 439, 443, 449, 457, 461, 463, 467,
479, 487, 491, 499, 503, 509, 521, 523, 541, 547,
557, 563, 569, 571, 577, 587, 593, 599, 601, 607,
613, 617, 619, 631, 641, 643, 647, 653, 659, 661,
673, 677, 683, 691, 701, 709, 719, 727, 733, 739,
743, 751, 757, 761, 769, 773, 787, 797, 809, 811,
821, 823, 827, 829, 839, 853, 857, 859, 863, 877,
881, 883, 887, 907, 911, 919, 929, 937, 941, 947,
953, 967, 971, 977, 983, 991, 997, 1009, 1013, 1019,
1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087,
1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153,
1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229,
1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291, 1297,
1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381,
1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453,
1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523,
1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583, 1597,
1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663,
1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823,
1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901,
1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993,
1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063,
2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131,
2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221,
2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287, 2293,
2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371,
2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437,
2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539,
2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617, 2621,
2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689,
2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749,
2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833,
2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903, 2909,
2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001,
3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079, 3083,
3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187,
3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259,
3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343,
3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433,
3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517,
3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581,
3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659,
3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733,
3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823,
3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911,
3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989, 4001,
4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073,
4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139, 4153,
4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241,
4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327,
4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421,
4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493, 4507,
4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591,
4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663,
4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759,
4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831, 4861,
4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943,
4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003};
int main () {
unsigned int maxn=4000000; // maximum n value
unsigned int i,j,k,a,b,p;
double x,ln2,ln3,maxx,minx;
double x0,x1,x2,x3;
double maxx0,maxx1,maxx2,maxx3,maxx4,maxx5,maxx6,maxx7;
double minx0,minx1,minx2,minx3,minx4,minx5,minx6,minx7;
unsigned int iter0,iter1,iter2,iter3;
unsigned int save0,save1,save2,save3;
unsigned int last0,last1,last2,last3;
unsigned int delta,count,limit,popcnt,temp,n;
unsigned int X[198140];
if (maxn<=165392) // set number of words
n=8192;
else {
if (maxn<=330784)
n=16384;
else {
if (maxn<=661576)
n=32768;
else {
if (maxn<=1323152)
n=65536;
else {
if (maxn<=2646304)
n=131072;
else {
if (maxn<=4000000)
n=198140;
else{
error[0]=1; // not enough internal memory
n=198140;
}
}
}
}
}
}
for (i=0; i<18; i++) // clear maximum/minimum array
outx[i]=0.0;
maxx=-655360.0; // set maximum
minx=65536.0; // set minimum
maxx0=-65536.0;
minx0=65536.0;
maxx1=-65536.0;
minx1=65536.0;
maxx2=-65536.0;
minx2=65536.0;
maxx3=-65536.0;
minx3=65536.0;
maxx4=-65536.0;
minx4=65536.0;
maxx5=-65536.0;
minx5=65536.0;
maxx6=-65536.0;
minx6=65536.0;
maxx7=-65536.0;
minx7=65536.0;
for (i=0; i<200; i++) { // clear histogram
hist0[i]=0;
hist1[i]=0;
hist2[i]=0;
hist3[i]=0;
hist4[i]=0;
output[i]=0;
hipop[i]=0;
max0[i]=0.0;
min0[i]=65536.0*65536.0*65536.0*65536.0;
max1[i]=0.0;
min1[i]=65536.0*65536.0*65536.0*65536.0;
max2[i]=0.0;
min2[i]=65536.0*65536.0*65536.0*65536.0;
max3[i]=0.0;
min3[i]=65536.0*65536.0*65536.0*65536.0;
}
count=0; // clear counter
popcnt=0; // clear counter
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
iter0=0; // period count
iter1=0;
iter2=0;
iter3=0;
save0=0; // prime from last period
save1=0;
save2=0;
save3=0;
last0=0; // last prime
last1=0;
last2=0;
last3=0;
x0=0.0; // last minimum
x1=0.0;
x2=0.0;
x3=0.0;
ln2=log(2); // log(2)
ln3=log(3); // log(3)
a=0; // clear a
b=0; // clear b
setn(X, 0, n);
X[0]=0x10000000; // x=1/2
for (i=0; i<maxn; i++) { // maximum n value
if (((X[0]+0xe0000000)&0x80000000)==0) { // x>1
mul3shft(X, 2, n); // x=(3/4)x
a=a+1;
}
else { // x<1
mul3shft(X, 1, n); // x=(3/2)x
b=b+1;
}
if (i>=(maxn-200)) { // check for overflow
output[count]=X[n-10];
count=count+1;
}
error[2]=a;
error[3]=b;
j=i+1; // n
if (j==1) // check if n=1
goto bskip;
if ((j&1)==0) // check if 2 divides n
goto bskip;
if (j==(j/3)*3) // check if 3 divides n
goto bskip;
if (j==(j/5)*5) // check if 5 divides n
goto bskip;
if (j>table[666]) // check if n is greater than last prime
goto yskip;
for (k=0; k<667; k++) { // check if n is in the prime LUT
if (j==table[k])
goto askip;
if (table[k]>j)
goto bskip;
}
goto bskip;
yskip:
limit=(unsigned int)(sqrt((double)j)+20.0);
for (k=0; k<667; k++) {
p=table[k]; // load prime
if (p>limit) // check limit
goto askip;
if (j==(j/p)*p) // check if n is divisible by prime
goto bskip;
}
askip:
x=(double)(a+1)/((double)(2*a+b+1)*ln2-(double)(a+b)*ln3); // upper bound of minimum
if (x>maxx) // compare to maximum
maxx=x;
if (x<minx) // compare to minimum
minx=x;
if ((j+1)==((j+1)/12)*12) { // check if 12 divides n+1
if ((x0>0.0)&&(x<0.0)) { // check for new period
delta=(last0-save0)/12; // difference in n values
temp=save0;
save0=last0; // save old n value
if (iter0!=0) {
if (delta<200) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x0)
max0[delta]=x0;
if (min0[delta]>x0)
min0[delta]=x0;
if (delta<75) { // maximum/minimum in lower population
if (x0>maxx0)
maxx0=x0;
if (x<minx0)
minx0=x;
}
else { // maximum/minimum in upper population
if (x0>maxx4)
maxx4=x0;
if (x<minx4)
minx4=x;
if (popcnt<100) {
hipop[popcnt*2]=temp;
hipop[popcnt*2+1]=save0;
popcnt=popcnt+1;
}
}
if (delta>120) {
hipop3[0]=temp;
hipop3[1]=save0;
}
}
iter0=iter0+1; // increment period count
}
last0=j; // save n value
x0=x; // last x value
}
if ((j+5)==((j+5)/12)*12) { // check if 12 divides n+5
if ((x1>0.0)&&(x<0.0)) { // check for new period
delta=(last1-save1)/12; // difference in n values
temp=save1;
save1=last1; // save old n value
if (iter1!=0) {
if (delta<200) // histogram differences
hist1[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max1[delta]<x1)
max1[delta]=x1;
if (min1[delta]>x1)
min1[delta]=x1;
if (delta<75) { // maximum/minimum in lower population
if (x1>maxx1)
maxx1=x1;
if (x<minx1)
minx1=x;
}
else { // maximum/minimum in upper population
if (x1>maxx5)
maxx5=x1;
if (x<minx5)
minx5=x;
if (popcnt<100) {
hipop[popcnt*2]=temp;
hipop[popcnt*2+1]=save1;
popcnt=popcnt+1;
}
}
if (delta>120) {
hipop3[0]=temp;
hipop3[1]=save1;
}
}
iter1=iter1+1; // increment period count
}
last1=j; // save n value
x1=x; // last x value
}
if ((j+7)==((j+7)/12)*12) { // check if 12 divides n+7
if ((x2>0.0)&&(x<0.0)) { // check for new period
delta=(last2-save2)/12; // difference in n values
temp=save2;
save2=last2; // save old n value
if (iter2!=0) {
if (delta<200) // histogram differences
hist2[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max2[delta]<x2)
max2[delta]=x2;
if (min2[delta]>x2)
min2[delta]=x2;
if (delta<75) { // maximum/minimum in lower population
if (x2>maxx2)
maxx2=x2;
if (x<minx2)
minx2=x;
}
else { // maximum/minimum in upper population
if (x2>maxx6)
maxx6=x2;
if (x<minx6)
minx6=x;
if (popcnt<100) {
hipop[popcnt*2]=temp;
hipop[popcnt*2+1]=save2;
popcnt=popcnt+1;
}
}
if (delta>120) {
hipop3[0]=temp;
hipop3[1]=save2;
}
}
iter2=iter2+1; // increment period count
}
last2=j; // save n value
x2=x; // last x value
}
if ((j+11)==((j+11)/12)*12) { // check if 12 divides n+11
if ((x3>0.0)&&(x<0.0)) { // check for new period
delta=(last3-save3)/12; // difference in n values
temp=save3;
save3=last3; // save old n value
if (iter3!=0) {
if (delta<200) // histogram differences
hist3[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max3[delta]<x3)
max3[delta]=x3;
if (min3[delta]>x3)
min3[delta]=x3;
if (delta<75) { // maximum/minimum in lower population
if (x3>maxx3)
maxx3=x3;
if (x<minx3)
minx3=x;
}
else { // maximum/minimum in upper population
if (x3>maxx7)
maxx7=x3;
if (x<minx7)
minx7=x;
if (popcnt<100) {
hipop[popcnt*2]=temp;
hipop[popcnt*2+1]=save3;
popcnt=popcnt+1;
}
}
if (delta>120) {
hipop3[0]=temp;
hipop3[1]=save3;
}
}
iter3=iter3+1; // increment period count
}
last3=j; // save n value
x3=x; // last x value
}
bskip:
j=0;
}
for (i=0; i<200; i++) { // sum up histograms
hist4[i]=hist0[i]+hist1[i]+hist2[i]+hist3[i];
max4[i]=max0[i];
min4[i]=min0[i];
}
for (i=0; i<200; i++) {
if (max4[i]<max1[i])
max4[i]=max1[i];
if (min4[i]>min1[i])
min4[i]=min1[i];
if (max4[i]<max2[i])
max4[i]=max2[i];
if (min4[i]>min2[i])
min4[i]=min2[i];
if (max4[i]<max3[i])
max4[i]=max3[i];
if (min4[i]>min3[i])
min4[i]=min3[i];
}
for (i=0; i<200; i++)
max4i[i]=(unsigned int)(max4[i]/65536.0);
outx[0]=maxx;
outx[1]=maxx0;
outx[2]=maxx1;
outx[3]=maxx2;
outx[4]=maxx3;
outx[5]=maxx4;
outx[6]=maxx5;
outx[7]=maxx6;
outx[8]=maxx7;
outx[9]=minx;
outx[10]=minx0;
outx[11]=minx1;
outx[12]=minx2;
outx[13]=minx3;
outx[14]=minx4;
outx[15]=minx5;
outx[16]=minx6;
outx[17]=minx7;
zskip:
return(0);
}