/******************************************************************************/
/* */
/* COMPUTE UPPER BOUNDS OF MINIMA */
/* 08/25/10 (dkc) */
/* */
/* This C program computes upper bounds of minima where l is prime. This */
/* program is for use on the C64. */
/* */
/******************************************************************************/
#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 XG[198140]; // global copy of X
far double max0[400],min0[400]; // maximum and minimum upper bounds of minima
far unsigned int maxi[400]; // maximum upper bounds of minima (integer)
far unsigned int output[400]; // check for overflow of X array
far unsigned int hist0[400]; // histogram of x values
far unsigned int error[6]; // error array
far double outx[2]; // maximum and minimum upper bounds of minima
far unsigned int outp[2]; // maximum and minimum l 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,maxp,minp;
double x,ln2,ln3,maxx,minx;
double x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17;
unsigned int iter0,iter1,iter2,iter3,iter4,iter5,iter6,iter7,iter8,iter9,iter10;
unsigned int iter11,iter12,iter13,iter14,iter15,iter16,iter17;
unsigned int save0,save1,save2,save3,save4,save5,save6,save7,save8,save9,save10;
unsigned int save11,save12,save13,save14,save15,save16,save17;
unsigned int last0,last1,last2,last3,last4,last5,last6,last7,last8,last9,last10;
unsigned int last11,last12,last13,last14,last15,last16,last17;
unsigned int delta,count,limit,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<2; i++) { // clear maximum/minimum array
outx[i]=0.0;
outp[i]=0;
}
maxx=-655360.0; // set maximum
minx=65536.0; // set minimum
for (i=0; i<400; i++) { // clear histogram
hist0[i]=0;
output[i]=0;
max0[i]=0.0;
min0[i]=65536.0*65536.0*65536.0*65536.0;
maxi[i]=0;
}
count=0; // clear counter
error[0]=0; // clear error array
error[1]=0;
error[2]=0;
error[3]=0;
error[4]=0;
error[5]=0;
iter0=0; // period count
iter1=0;
iter2=0;
iter3=0;
iter4=0;
iter5=0;
iter6=0;
iter7=0;
iter8=0;
iter9=0;
iter10=0;
iter11=0;
iter12=0;
iter13=0;
iter14=0;
iter15=0;
iter16=0;
iter17=0;
save0=0; // prime from last period
save1=0;
save2=0;
save3=0;
save4=0;
save5=0;
save6=0;
save7=0;
save8=0;
save9=0;
save10=0;
save11=0;
save12=0;
save13=0;
save14=0;
save15=0;
save16=0;
save17=0;
last0=0; // last prime
last1=0;
last2=0;
last3=0;
last4=0;
last5=0;
last6=0;
last7=0;
last8=0;
last9=0;
last10=0;
last11=0;
last12=0;
last13=0;
last14=0;
last15=0;
last16=0;
last17=0;
x0=0.0; // last minimum
x1=0.0;
x2=0.0;
x3=0.0;
x4=0.0;
x5=0.0;
x6=0.0;
x7=0.0;
x8=0.0;
x9=0.0;
x10=0.0;
x11=0.0;
x12=0.0;
x13=0.0;
x14=0.0;
x15=0.0;
x16=0.0;
x17=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-400)) { // check for overflow
output[count]=X[n-10];
count=count+1;
}
error[2]=a;
error[3]=b;
j=2*a+b+1; // l
if (j==1) // check if l=1
goto bskip;
if ((j&1)==0) // check if 2 divides l
goto bskip;
if (j==(j/3)*3) // check if 3 divides l
goto bskip;
if (j==(j/5)*5) // check if 5 divides l
goto bskip;
if (j>table[666]) // check if l is greater than last prime
goto yskip;
for (k=0; k<667; k++) { // check if l 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 l is divisible by prime
goto bskip;
}
askip:
error[4]=j; // save l
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;
maxp=j;
}
if (x<minx) { // compare to minimum
minx=x;
minp=j;
}
if ((j+1)==((j+1)/54)*54) { // check if 54 divides l+1
if ((x0>0.0)&&(x<0.0)) { // check for new period
delta=(last0-save0)/54; // difference in l values
save0=last0; // save old l value
if (iter0!=0) {
if (delta<400) // 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;
}
iter0=iter0+1; // increment period count
}
last0=j; // save l value
x0=x; // last x value
}
if ((j+5)==((j+5)/54)*54) { // check if 54 divides l+5
if ((x1>0.0)&&(x<0.0)) { // check for new period
delta=(last1-save1)/54; // difference in l values
save1=last1; // save old l value
if (iter1!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x1)
max0[delta]=x1;
if (min0[delta]>x1)
min0[delta]=x1;
}
iter1=iter1+1; // increment period count
}
last1=j; // save l value
x1=x; // last x value
}
if ((j+7)==((j+7)/54)*54) { // check if 54 divides l+7
if ((x2>0.0)&&(x<0.0)) { // check for new period
delta=(last2-save2)/54; // difference in l values
save2=last2; // save old l value
if (iter2!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x2)
max0[delta]=x2;
if (min0[delta]>x2)
min0[delta]=x2;
}
iter2=iter2+1; // increment period count
}
last2=j; // save l value
x2=x; // last x value
}
if ((j+11)==((j+11)/54)*54) { // check if 54 divides l+11
if ((x3>0.0)&&(x<0.0)) { // check for new period
delta=(last3-save3)/54; // difference in l values
save3=last3; // save old l value
if (iter3!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x3)
max0[delta]=x3;
if (min0[delta]>x3)
min0[delta]=x3;
}
iter3=iter3+1; // increment period count
}
last3=j; // save l value
x3=x; // last x value
}
if ((j+13)==((j+13)/54)*54) { // check if 54 divides l+13
if ((x4>0.0)&&(x<0.0)) { // check for new period
delta=(last4-save4)/54; // difference in l values
save4=last4; // save old l value
if (iter4!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x4)
max0[delta]=x4;
if (min0[delta]>x4)
min0[delta]=x4;
}
iter4=iter4+1; // increment period count
}
last4=j; // save l value
x4=x; // last x value
}
if ((j+17)==((j+17)/54)*54) { // check if 54 divides l+17
if ((x5>0.0)&&(x<0.0)) { // check for new period
delta=(last5-save5)/54; // difference in l values
save5=last5; // save old l value
if (iter5!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x5)
max0[delta]=x5;
if (min0[delta]>x5)
min0[delta]=x5;
}
iter5=iter5+1; // increment period count
}
last5=j; // save l value
x5=x; // last x value
}
if ((j+19)==((j+19)/54)*54) { // check if 54 divides l+19
if ((x6>0.0)&&(x<0.0)) { // check for new period
delta=(last6-save6)/54; // difference in l values
save6=last6; // save old l value
if (iter6!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x6)
max0[delta]=x6;
if (min0[delta]>x6)
min0[delta]=x6;
}
iter6=iter6+1; // increment period count
}
last6=j; // save l value
x6=x; // last x value
}
if ((j+23)==((j+23)/54)*54) { // check if 54 divides l+23
if ((x7>0.0)&&(x<0.0)) { // check for new period
delta=(last7-save7)/54; // difference in l values
save7=last7; // save old l value
if (iter7!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x7)
max0[delta]=x7;
if (min0[delta]>x7)
min0[delta]=x7;
}
iter7=iter7+1; // increment period count
}
last7=j; // save l value
x7=x; // last x value
}
if ((j+25)==((j+25)/54)*54) { // check if 54 divides l+25
if ((x8>0.0)&&(x<0.0)) { // check for new period
delta=(last8-save8)/54; // difference in l values
save8=last8; // save old l value
if (iter8!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x8)
max0[delta]=x8;
if (min0[delta]>x8)
min0[delta]=x8;
}
iter8=iter8+1; // increment period count
}
last8=j; // save l value
x8=x; // last x value
}
if ((j+29)==((j+29)/54)*54) { // check if 54 divides n+29
if ((x9>0.0)&&(x<0.0)) { // check for new period
delta=(last9-save9)/54; // difference in l values
save9=last9; // save old l value
if (iter9!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x9)
max0[delta]=x9;
if (min0[delta]>x9)
min0[delta]=x9;
}
iter9=iter9+1; // increment period count
}
last9=j; // save l value
x9=x; // last x value
}
if ((j+31)==((j+31)/54)*54) { // check if 54 divides n+31
if ((x10>0.0)&&(x<0.0)) { // check for new period
delta=(last10-save10)/54; // difference in l values
save10=last10; // save old l value
if (iter10!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x10)
max0[delta]=x10;
if (min0[delta]>x10)
min0[delta]=x10;
}
iter10=iter10+1; // increment period count
}
last10=j; // save l value
x10=x; // last x value
}
if ((j+35)==((j+35)/54)*54) { // check if 54 divides l+35
if ((x11>0.0)&&(x<0.0)) { // check for new period
delta=(last11-save11)/54; // difference in l values
save11=last11; // save old l value
if (iter11!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x11)
max0[delta]=x11;
if (min0[delta]>x11)
min0[delta]=x11;
}
iter11=iter11+1; // increment period count
}
last11=j; // save l value
x11=x; // last x value
}
if ((j+37)==((j+37)/54)*54) { // check if 54 divides l+37
if ((x12>0.0)&&(x<0.0)) { // check for new period
delta=(last12-save12)/54; // difference in l values
save12=last12; // save old l value
if (iter12!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x12)
max0[delta]=x12;
if (min0[delta]>x12)
min0[delta]=x12;
}
iter12=iter12+1; // increment period count
}
last12=j; // save l value
x12=x; // last x value
}
if ((j+41)==((j+41)/54)*54) { // check if 54 divides l+41
if ((x13>0.0)&&(x<0.0)) { // check for new period
delta=(last13-save13)/54; // difference in l values
save13=last13; // save old l value
if (iter13!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x13)
max0[delta]=x13;
if (min0[delta]>x13)
min0[delta]=x13;
}
iter13=iter13+1; // increment period count
}
last13=j; // save l value
x13=x; // last x value
}
if ((j+43)==((j+43)/54)*54) { // check if 54 divides l+43
if ((x14>0.0)&&(x<0.0)) { // check for new period
delta=(last14-save14)/54; // difference in l values
save14=last14; // save old l value
if (iter14!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x14)
max0[delta]=x14;
if (min0[delta]>x14)
min0[delta]=x14;
}
iter14=iter14+1; // increment period count
}
last14=j; // save l value
x14=x; // last x value
}
if ((j+47)==((j+47)/54)*54) { // check if 54 divides l+47
if ((x15>0.0)&&(x<0.0)) { // check for new period
delta=(last15-save15)/54; // difference in l values
save15=last15; // save old l value
if (iter15!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x15)
max0[delta]=x15;
if (min0[delta]>x15)
min0[delta]=x15;
}
iter15=iter15+1; // increment period count
}
last15=j; // save l value
x15=x; // last x value
}
if ((j+49)==((j+49)/54)*54) { // check if 54 divides l+49
if ((x16>0.0)&&(x<0.0)) { // check for new period
delta=(last16-save16)/54; // difference in l values
save16=last16; // save old l value
if (iter16!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x16)
max0[delta]=x16;
if (min0[delta]>x16)
min0[delta]=x16;
}
iter16=iter16+1; // increment period count
}
last16=j; // save l value
x16=x; // last x value
}
if ((j+53)==((j+53)/54)*54) { // check if 54 divides l+53
if ((x17>0.0)&&(x<0.0)) { // check for new period
delta=(last17-save17)/54; // difference in l values
save17=last17; // save old l value
if (iter17!=0) {
if (delta<400) // histogram differences
hist0[delta]+=1;
else {
error[0]=3;
error[1]=delta;
goto zskip;
}
if (max0[delta]<x17)
max0[delta]=x17;
if (min0[delta]>x17)
min0[delta]=x17;
}
iter17=iter17+1; // increment period count
}
last17=j; // save l value
x17=x; // last x value
}
bskip:
j=0;
}
outx[0]=maxx;
outx[1]=minx;
outp[0]=maxp;
outp[1]=minp;
for (i=0; i<400; i++)
maxi[i]=(unsigned int)(max0[i]/65536.0);
for (i=0; i<198140; i++)
XG[i]=X[i];
zskip:
return(0);
}