/******************************************************************************/
/* */
/* COMPUTE M(l,m)/(2^l-3^m)FOR AN M-CYCLE */
/* 04/03/10 (dkc) */
/* */
/* This C program computes M(l,m)/(2^l-3^m) for an M-cycle. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int *quotient, unsigned int d2,
unsigned int d3);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
void subn(unsigned int *c, unsigned int *d, unsigned int n);
void copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
unsigned int normn(unsigned int *a, unsigned int n);
unsigned int orn(unsigned int *a, unsigned int n);
unsigned int ab[476*2]={ // l and n values
0x0003, 0x002,
0x0005, 0x003,
0x0008, 0x005,
0x000B, 0x007,
0x000D, 0x008,
0x0010, 0x00A,
0x0013, 0x00C,
0x0015, 0x00D,
0x0018, 0x00F,
0x001B, 0x011,
0x001E, 0x013,
0x0020, 0x014,
0x0023, 0x016,
0x0026, 0x018,
0x0028, 0x019,
0x002B, 0x01B,
0x002E, 0x01D,
0x0031, 0x01F,
0x0033, 0x020,
0x0036, 0x022,
0x0039, 0x024,
0x003B, 0x025,
0x003E, 0x027,
0x0041, 0x029,
0x0044, 0x02B,
0x0046, 0x02C,
0x0049, 0x02E,
0x004C, 0x030,
0x004E, 0x031,
0x0051, 0x033,
0x0054, 0x035,
0x0056, 0x036,
0x0059, 0x038,
0x005C, 0x03A,
0x005F, 0x03C,
0x0061, 0x03D,
0x0064, 0x03F,
0x0067, 0x041,
0x0069, 0x042,
0x006C, 0x044,
0x006F, 0x046,
0x0072, 0x048,
0x0074, 0x049,
0x0077, 0x04B,
0x007A, 0x04D,
0x007C, 0x04E,
0x007F, 0x050,
0x0082, 0x052,
0x0085, 0x054,
0x0087, 0x055,
0x008A, 0x057,
0x008D, 0x059,
0x008F, 0x05A,
0x0092, 0x05C,
0x0095, 0x05E,
0x0098, 0x060,
0x009A, 0x061,
0x009D, 0x063,
0x00A0, 0x065,
0x00A2, 0x066,
0x00A5, 0x068,
0x00A8, 0x06A,
0x00AA, 0x06B,
0x00AD, 0x06D,
0x00B0, 0x06F,
0x00B3, 0x071,
0x00B5, 0x072,
0x00B8, 0x074,
0x00BB, 0x076,
0x00BD, 0x077,
0x00C0, 0x079,
0x00C3, 0x07B,
0x00C6, 0x07D,
0x00C8, 0x07E,
0x00CB, 0x080,
0x00CE, 0x082,
0x00D0, 0x083,
0x00D3, 0x085,
0x00D6, 0x087,
0x00D9, 0x089,
0x00DB, 0x08A,
0x00DE, 0x08C,
0x00E1, 0x08E,
0x00E3, 0x08F,
0x00E6, 0x091,
0x00E9, 0x093,
0x00EC, 0x095,
0x00EE, 0x096,
0x00F1, 0x098,
0x00F4, 0x09A,
0x00F6, 0x09B,
0x00F9, 0x09D,
0x00FC, 0x09F,
0x00FE, 0x0A0,
0x0101, 0x0A2,
0x0104, 0x0A4,
0x0107, 0x0A6,
0x0109, 0x0A7,
0x010C, 0x0A9,
0x010F, 0x0AB,
0x0111, 0x0AC,
0x0114, 0x0AE,
0x0117, 0x0B0,
0x011A, 0x0B2,
0x011C, 0x0B3,
0x011F, 0x0B5,
0x0122, 0x0B7,
0x0124, 0x0B8,
0x0127, 0x0BA,
0x012A, 0x0BC,
0x012D, 0x0BE,
0x012F, 0x0BF,
0x0132, 0x0C1,
0x0135, 0x0C3,
0x0137, 0x0C4,
0x013A, 0x0C6,
0x013D, 0x0C8,
0x0140, 0x0CA,
0x0142, 0x0CB,
0x0145, 0x0CD,
0x0148, 0x0CF,
0x014A, 0x0D0,
0x014D, 0x0D2,
0x0150, 0x0D4,
0x0152, 0x0D5,
0x0155, 0x0D7,
0x0158, 0x0D9,
0x015B, 0x0DB,
0x015D, 0x0DC,
0x0160, 0x0DE,
0x0163, 0x0E0,
0x0165, 0x0E1,
0x0168, 0x0E3,
0x016B, 0x0E5,
0x016E, 0x0E7,
0x0170, 0x0E8,
0x0173, 0x0EA,
0x0176, 0x0EC,
0x0178, 0x0ED,
0x017B, 0x0EF,
0x017E, 0x0F1,
0x0181, 0x0F3,
0x0183, 0x0F4,
0x0186, 0x0F6,
0x0189, 0x0F8,
0x018B, 0x0F9,
0x018E, 0x0FB,
0x0191, 0x0FD,
0x0194, 0x0FF,
0x0196, 0x100,
0x0199, 0x102,
0x019C, 0x104,
0x019E, 0x105,
0x01A1, 0x107,
0x01A4, 0x109,
0x01A6, 0x10A,
0x01A9, 0x10C,
0x01AC, 0x10E,
0x01AF, 0x110,
0x01B1, 0x111,
0x01B4, 0x113,
0x01B7, 0x115,
0x01B9, 0x116,
0x01BC, 0x118,
0x01BF, 0x11A,
0x01C2, 0x11C,
0x01C4, 0x11D,
0x01C7, 0x11F,
0x01CA, 0x121,
0x01CC, 0x122,
0x01CF, 0x124,
0x01D2, 0x126,
0x01D5, 0x128,
0x01D7, 0x129,
0x01DA, 0x12B,
0x01DD, 0x12D,
0x01DF, 0x12E,
0x01E2, 0x130,
0x01E5, 0x132,
0x01E8, 0x134,
0x01EA, 0x135,
0x01ED, 0x137,
0x01F0, 0x139,
0x01F2, 0x13A,
0x01F5, 0x13C,
0x01F8, 0x13E,
0x01FA, 0x13F,
0x01FD, 0x141,
0x0200, 0x143,
0x0203, 0x145,
0x0205, 0x146,
0x0208, 0x148,
0x020B, 0x14A,
0x020D, 0x14B,
0x0210, 0x14D,
0x0213, 0x14F,
0x0216, 0x151,
0x0218, 0x152,
0x021B, 0x154,
0x021E, 0x156,
0x0220, 0x157,
0x0223, 0x159,
0x0226, 0x15B,
0x0229, 0x15D,
0x022B, 0x15E,
0x022E, 0x160,
0x0231, 0x162,
0x0233, 0x163,
0x0236, 0x165,
0x0239, 0x167,
0x023B, 0x168,
0x023E, 0x16A,
0x0241, 0x16C,
0x0244, 0x16E,
0x0246, 0x16F,
0x0249, 0x171,
0x024C, 0x173,
0x024E, 0x174,
0x0251, 0x176,
0x0254, 0x178,
0x0257, 0x17A,
0x0259, 0x17B,
0x025C, 0x17D,
0x025F, 0x17F,
0x0261, 0x180,
0x0264, 0x182,
0x0267, 0x184,
0x026A, 0x186,
0x026C, 0x187,
0x026F, 0x189,
0x0272, 0x18B,
0x0274, 0x18C,
0x0277, 0x18E,
0x027A, 0x190,
0x027D, 0x192,
0x027F, 0x193,
0x0282, 0x195,
0x0285, 0x197,
0x0287, 0x198,
0x028A, 0x19A,
0x028D, 0x19C,
0x028F, 0x19D,
0x0292, 0x19F,
0x0295, 0x1A1,
0x0298, 0x1A3,
0x029A, 0x1A4,
0x029D, 0x1A6,
0x02A0, 0x1A8,
0x02A2, 0x1A9,
0x02A5, 0x1AB,
0x02A8, 0x1AD,
0x02AB, 0x1AF,
0x02AD, 0x1B0,
0x02B0, 0x1B2,
0x02B3, 0x1B4,
0x02B5, 0x1B5,
0x02B8, 0x1B7,
0x02BB, 0x1B9,
0x02BE, 0x1BB,
0x02C0, 0x1BC,
0x02C3, 0x1BE,
0x02C6, 0x1C0,
0x02C8, 0x1C1,
0x02CB, 0x1C3,
0x02CE, 0x1C5,
0x02D1, 0x1C7,
0x02D3, 0x1C8,
0x02D6, 0x1CA,
0x02D9, 0x1CC,
0x02DB, 0x1CD,
0x02DE, 0x1CF,
0x02E1, 0x1D1,
0x02E3, 0x1D2,
0x02E6, 0x1D4,
0x02E9, 0x1D6,
0x02EC, 0x1D8,
0x02EE, 0x1D9,
0x02F1, 0x1DB,
0x02F4, 0x1DD,
0x02F6, 0x1DE,
0x02F9, 0x1E0,
0x02FC, 0x1E2,
0x02FF, 0x1E4,
0x0301, 0x1E5,
0x0304, 0x1E7,
0x0307, 0x1E9,
0x0309, 0x1EA,
0x030C, 0x1EC,
0x030F, 0x1EE,
0x0312, 0x1F0,
0x0314, 0x1F1,
0x0317, 0x1F3,
0x031A, 0x1F5,
0x031C, 0x01F6,
0x031F, 0x01F8,
0x0322, 0x01FA,
0x0325, 0x01FC,
0x0327, 0x01FD,
0x032A, 0x01FF,
0x032D, 0x0201,
0x032F, 0x0202,
0x0332, 0x0204,
0x0335, 0x0206,
0x0337, 0x0207,
0x033A, 0x0209,
0x033D, 0x020B,
0x0340, 0x020D,
0x0342, 0x020E,
0x0345, 0x0210,
0x0348, 0x0212,
0x034A, 0x0213,
0x034D, 0x0215,
0x0350, 0x0217,
0x0353, 0x0219,
0x0355, 0x021A,
0x0358, 0x021C,
0x035B, 0x021E,
0x035D, 0x021F,
0x0360, 0x0221,
0x0363, 0x0223,
0x0366, 0x0225,
0x0368, 0x0226,
0x036B, 0x0228,
0x036E, 0x022A,
0x0370, 0x022B,
0x0373, 0x022D,
0x0376, 0x022F,
0x0379, 0x0231,
0x037B, 0x0232,
0x037E, 0x0234,
0x0381, 0x0236,
0x0383, 0x0237,
0x0386, 0x0239,
0x0389, 0x023B,
0x038B, 0x023C,
0x038E, 0x023E,
0x0391, 0x0240,
0x0394, 0x0242,
0x0396, 0x0243,
0x0399, 0x0245,
0x039C, 0x0247,
0x039E, 0x0248,
0x03A1, 0x024A,
0x03A4, 0x024C,
0x03A7, 0x024E,
0x03A9, 0x024F,
0x03AC, 0x0251,
0x03AF, 0x0253,
0x03B1, 0x0254,
0x03B4, 0x0256,
0x03B7, 0x0258,
0x03BA, 0x025A,
0x03BC, 0x025B,
0x03BF, 0x025D,
0x03C2, 0x025F,
0x03C4, 0x0260,
0x03C7, 0x0262,
0x03CA, 0x0264,
0x03CD, 0x0266,
0x03CF, 0x0267,
0x03D2, 0x0269,
0x03D5, 0x026B,
0x03D7, 0x026C,
0x03DA, 0x026E,
0x03DD, 0x0270,
0x03DF, 0x0271,
0x03E2, 0x0273,
0x03E5, 0x0275,
0x03E8, 0x0277,
0x03EA, 0x0278,
0x03ED, 0x027A,
0x03F0, 0x027C,
0x03F2, 0x027D,
0x03F5, 0x027F,
0x03F8, 0x0281,
0x03FB, 0x0283,
0x03FD, 0x0284,
0x0400, 0x0286,
0x0403, 0x0288,
0x0405, 0x0289,
0x0408, 0x028B,
0x040B, 0x028D,
0x040E, 0x028F,
0x0410, 0x0290,
0x0413, 0x0292,
0x0416, 0x0294,
0x0418, 0x0295,
0x041B, 0x0297,
0x041E, 0x0299,
0x0420, 0x029A,
0x0423, 0x029C,
0x0426, 0x029E,
0x0429, 0x02A0,
0x042B, 0x02A1,
0x042E, 0x02A3,
0x0431, 0x02A5,
0x0433, 0x02A6,
0x0436, 0x02A8,
0x0439, 0x02AA,
0x043C, 0x02AC,
0x043E, 0x02AD,
0x0441, 0x02AF,
0x0444, 0x02B1,
0x0446, 0x02B2,
0x0449, 0x02B4,
0x044C, 0x02B6,
0x044F, 0x02B8,
0x0451, 0x02B9,
0x0454, 0x02BB,
0x0457, 0x02BD,
0x0459, 0x02BE,
0x045C, 0x02C0,
0x045F, 0x02C2,
0x0462, 0x02C4,
0x0464, 0x02C5,
0x0467, 0x02C7,
0x046A, 0x02C9,
0x046C, 0x02CA,
0x046F, 0x02CC,
0x0472, 0x02CE,
0x0474, 0x02CF,
0x0477, 0x02D1,
0x047A, 0x02D3,
0x047D, 0x02D5,
0x047F, 0x02D6,
0x0482, 0x02D8,
0x0485, 0x02DA,
0x0487, 0x02DB,
0x048A, 0x02DD,
0x048D, 0x02DF,
0x0490, 0x02E1,
0x0492, 0x02E2,
0x0495, 0x02E4,
0x0498, 0x02E6,
0x049A, 0x02E7,
0x049D, 0x02E9,
0x04A0, 0x02EB,
0x04A3, 0x02ED,
0x04A5, 0x02EE,
0x04A8, 0x02F0,
0x04AB, 0x02F2,
0x04AD, 0x02F3,
0x04B0, 0x02F5,
0x04B3, 0x02F7,
0x04B6, 0x02F9,
0x04B8, 0x02FA,
0x04BB, 0x02FC,
0x04BE, 0x02FE,
0x04C0, 0x02FF,
0x04C3, 0x0301,
0x04C6, 0x0303,
0x04C8, 0x0304,
0x04CB, 0x0306,
0x04CE, 0x0308,
0x04D1, 0x030A,
0x04D3, 0x030B,
0x04D6, 0x030D,
0x04D9, 0x030F,
0x04DB, 0x0310,
0x04DE, 0x0312,
0x04E1, 0x0314,
0x04E4, 0x0316,
0x04E6, 0x0317,
0x04E9, 0x0319,
0x04EC, 0x031B,
0x04EE, 0x031C,
0x04F1, 0x031E,
0x04F4, 0x0320,
0x04F7, 0x0322,
0x04F9, 0x0323,
0x04FC, 0x0325,
0x04FF, 0x0327,
0x0501, 0x0328};
int main () {
unsigned int m=512; // number of words
int g,h,i,j,l,n,a,b;
unsigned int A[1024],B[1024],C[4],S[1024],T[1024],flag,flag1,o;
FILE *Outfp;
Outfp = fopen("out18a.dat","w");
setn(T, 0, m);
for (i=0; i<476; i++) {
l=ab[2*i]; // length of loop
n=ab[2*i+1]; // number of odd elements in loop
setn(S, 0, m); // clear sum
for (j=1; j<=l; j++) {
a=j*n;
if (a==((a/l)*l)) // a=(j*n)/l
a=a/l;
else
a=(a/l)+1;
b=(j-1)*n;
if (b==((b/l)*l)) // b=((j-1)*n)/l
b=b/l;
else
b=(b/l)+1;
if ((a-b)>0x20000000) {
printf("error: a-b is too large \n");
goto zskip;
}
setn(A, a-b, m);
for (h=0; h<(n-a); h++) { // (a-b)*3**(n-a)
copyn(A, B, m);
addn(A, A, m);
addn(B, A, m);
}
lshiftn(A, B, j-1, m); // (a-b)*2**(j-1)*3**(n-a)
addn(B, S, m); // sum+=(a-b)*2**(j-1)*3**(n-a)
if ((S[0]&0xc0000000)!=0) {
printf("error A: sum too large \n");
goto zskip;
}
}
setn(A, 1, m);
for (h=0; h<n; h++) { // 3**n
copyn(A, B, m);
addn(A, A, m);
addn(B, A, m);
if ((A[0]&0xc0000000)!=0) {
printf("error: n=%d \n",n);
goto zskip;
}
}
if ((unsigned int)l>=(32*m)) {
printf("error: l=%d \n",l);
goto zskip;
}
setn(B, 1, m);
lshiftn(B, B, l, m); // 2**l
subn(B, A, m); // 2**l-3**n
flag=0;
if ((A[0]&0x80000000)!=0) { // negate 2**l-3**n
setn(B, 0, m);
subn(B, A, m);
flag=1;
}
copyn(T, B, m);
subn(A, B, m);
if ((B[0]&0x80000000)!=0) {
printf("not monotonic \n");
fprintf(Outfp,"not monotonic \n");
// goto zskip;
}
copyn(A, T, m);
if ((S[0]&0xff000000)!=0) {
printf("error: sum too big \n");
goto zskip;
}
lshiftn(S, S, 8, m);
flag1=0;
o=normn(A, m); // count to normalize divisor
g=(int)(32*m)-(int)o-62; // right shift amount
if (g>0) {
shiftn(S, S, g, m);
shiftn(A, A, g, m);
flag1=g;
}
o=orn(S, m-4);
if ((o|(S[m-4]&0xc0000000))!=0) {
printf("error B: sum too large \n");
goto zskip;
}
div128_64(S[m-4],S[m-3],S[m-2],S[m-1],C,A[m-2],A[m-1]); // sum/(2**l-3**n)
shiftn(C, C, 8, 4);
if ((C[0]|C[1]|C[2])!=0) {
printf("warning: x too large \n");
goto zskip;
}
printf("l=%d, n=%d, sign=%d, shift=%d, x=%d %d %d %d \n",l,n,flag,flag1,C[0],C[1],C[2],C[3]);
fprintf(Outfp,"l=%d, n=%d, sign=%d, shift=%d, x=%d %d %d %d \n",l,n,flag,flag1,C[0],C[1],C[2],C[3]);
}
zskip:
fclose(Outfp);
return(0);
}