/******************************************************************************/
/* */
/* COMPUTE THE MINIMUM AND MAXIMUM ODD ELEMENTS IN A CYCLE */
/* 10/29/10 (dkc) */
/* */
/* This C program computes the minimum element in a cycle using Halbeisen */
/* and Hungerbuhler's formula and computes the maximum odd element using an */
/* analogous formula. "l" is the number of elements in the cycle (where the */
/* element following an odd element i is (3*i+1)/2) and "n" is the number of */
/* odd elements in the cycle. */
/* */
/******************************************************************************/
#include <math.h>
/*****************************************************************************/
/* */
/* CARRY */
/* 01/12/07 (dkc) */
/* */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum) {
unsigned int c;
if ((a&b)>>31!=0)
c=1;
else {
if ((a|b)>>31==0)
c=0;
else {
if (sum>=0x80000000)
c=0;
else
c=1;
}
}
return c;
}
/*****************************************************************************/
/* */
/* LEFT-MOST BIT DETECTION */
/* 01/12/07 (dkc) */
/* */
/*****************************************************************************/
unsigned int lmbd(unsigned int mode, unsigned int a) {
unsigned int i,mask,count;
if (mode==0)
a=~a;
if (a==0)
return(32);
count=0;
mask=0x80000000;
for (i=0; i<32; i++) {
if ((a&mask)!=0)
break;
count+=1;
mask>>=1;
}
return(count);
}
/*****************************************************************************/
/* */
/* 128/64 BIT DIVIDE (UNSIGNED) */
/* 01/12/07 (dkc) */
/* */
/*****************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
unsigned int lmbd(unsigned int mode, unsigned int a);
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int *quotient, unsigned int d2,
unsigned int d3) {
unsigned int i,d0,d1,dshift,ashift,count,flag;
unsigned int shift,c,c0,c1,c2,temp,temp0,temp1,temp2,temp3;
if (d2==0) {
if ((a0==0)&&(a1==0)&&(a2==0)&&(a3<d3)) {
*quotient=0;
*(quotient+1)=0;
*(quotient+2)=0;
*(quotient+3)=0;
return;
}
}
else {
if ((a0==0)&&(a1==0)&&(a2<d2)) {
*quotient=0;
*(quotient+1)=0;
*(quotient+2)=0;
*(quotient+3)=0;
return;
}
}
dshift=lmbd(1,d2);
if (d2==0)
dshift+=lmbd(1,d3);
dshift+=64;
ashift=lmbd(1,a0);
if (a0==0)
ashift+=lmbd(1,a1);
if ((a0|a1)==0)
ashift+=lmbd(1,a2);
if ((a0|a1|a2)==0)
ashift+=lmbd(1,a3);
shift=dshift-ashift;
count=shift+1;
d0=0;
d1=0;
if (shift<32) {
if (shift!=0) {
d1=d2>>(32-shift);
d2=(d2<<shift)|(d3>>(32-shift));
}
else
d1=0; // added to get MSVC to work
d3=d3<<shift;
flag=3;
shift=32-shift;
}
else {
shift=shift-32;
d1=d2;
d2=d3;
d3=0;
if (shift<32) {
if (shift!=0) {
d0=d1>>(32-shift);
d1=(d1<<shift)|(d2>>(32-shift));
}
else
d0=0; // added to get MSVC to work
d2=d2<<shift;
flag=2;
shift=32-shift;
}
else {
shift=shift-32;
d0=d1;
d1=d2;
d2=0;
if (shift<32) {
if (shift!=0) // added to get MSVC to work
d0=(d0<<shift)|(d1>>(32-shift));
d1=d1<<shift;
flag=1;
shift=32-shift;
}
else {
shift=shift-32;
d0=d1;
d1=0;
d0=d0<<shift;
flag=0;
shift=32-shift;
}
}
}
d0=~d0;
d1=~d1;
d2=~d2;
d3=~d3;
temp=d3+1;
c=carry(d3,1,temp);
d3=temp;
temp=d2+c;
c=carry(d2,c,temp);
d2=temp;
temp=d1+c;
c=carry(d1,c,temp);
d1=temp;
d0=d0+c;
for (i=0; i<count; i++) {
temp3=a3+d3;
c2=carry(a3,d3,temp3);
temp=a2+c2;
c1=carry(a2,c2,temp);
temp2=temp+d2;
c1+=carry(temp,d2,temp2);
temp=a1+c1;
c0=carry(a1,c1,temp);
temp1=temp+d1;
c0+=carry(temp,d1,temp1);
temp0=a0+d0+c0;
if ((temp0>>31)==0) {
a0=temp0<<1;
if ((temp1>>31)!=0)
c=1;
else
c=0;
a0=a0+c;
a1=temp1<<1;
if ((temp2>>31)!=0)
c=1;
else
c=0;
a1=a1+c;
a2=temp2<<1;
if ((temp3>>31)!=0)
c=1;
else
c=0;
a2=a2+c;
a3=temp3<<1;
a3=a3+1;
}
else {
a0=a0<<1;
if ((a1>>31)!=0)
c=1;
else
c=0;
a0=a0+c;
a1=a1<<1;
if ((a2>>31)!=0)
c=1;
else
c=0;
a1=a1+c;
a2=a2<<1;
if ((a3>>31)!=0)
c=1;
else
c=0;
a2=a2+c;
a3=a3<<1;
}
}
shift=shift-1;
if (flag==3) {
a0=0;
a1=0;
a2=0;
a3=a3<<shift;
a3=a3>>shift;
}
else {
if (flag==2) {
a0=0;
a1=0;
a2=a2<<shift;
a2=a2>>shift;
}
else {
if (flag==1) {
a0=0;
a1=a1<<shift;
a1=a1>>shift;
}
else {
a0=a0<<shift;
a0=a0>>shift;
}
}
}
*quotient=a0;
*(quotient+1)=a1;
*(quotient+2)=a2;
*(quotient+3)=a3;
return;
}
/*****************************************************************************/
/* */
/* LEFT-SHIFT 32*N BITS */
/* 10/16/09 (dkc) */
/* */
/*****************************************************************************/
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned
int n) {
unsigned int i;
for (i=0; i<n; i++)
*(b+i)=*(a+i);
while (shift>31) {
for (i=0; i<n-1; i++)
*(b+i)=*(b+i+1);
*(b+n-1)=0;
shift=shift-32;
}
if (shift==0)
return;
for (i=0; i<n-1; i++)
*(b+i)=(*(b+i)<<shift)|(*(b+i+1)>>(32-shift));
*(b+n-1)=*(b+n-1)<<shift;
return;
}
/*****************************************************************************/
/* */
/* RIGHT-SHIFT 32*N BITS */
/* 04/02/10 (dkc) */
/* */
/*****************************************************************************/
void shiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned
int n) {
unsigned int i;
for (i=0; i<n; i++)
*(b+i)=*(a+i);
while (shift>31) {
for (i=n-1; i>0; i--)
*(b+i)=*(b+i-1);
*b=0;
shift=shift-32;
}
if (shift==0)
return;
for (i=n-1; i>0; i--)
*(b+i)=(*(b+i)>>shift)|(*(b+i-1)<<(32-shift));
*b=*b>>shift;
return;
}
/******************************************************************************
* *
* NORMALIZE N WORDS *
* 04/03/10 (dkc) *
* *
******************************************************************************/
unsigned int normn(unsigned int *a, unsigned int n) {
unsigned int i,shift,mask;
i=0;
while ((*(a+i)==0)&&(i<n))
i=i+1;
shift=32*i;
if (i==n)
return shift;
mask=0x80000000;
while ((*(a+i)&mask)==0) {
shift=shift+1;
mask=mask>>1;
}
return shift;
}
/******************************************************************************
* *
* N-WORD ADD *
* 04/19/10 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void addn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i,s;
unsigned int c[8192];
if (n>8192)
return;
for (i=n-1; i>0; i--) {
s=*(a+i)+*(b+i);
c[i]=carry(*(a+i),*(b+i),s);
*(b+i)=s;
}
*b=*a+*b;
for (i=n-2; i>0; i--) {
s=*(b+i)+c[i+1];
c[i]+=carry(*(b+i),c[i+1],s);
*(b+i)=s;
}
*b=*b+c[1];
return;
}
/******************************************************************************
* *
* N-WORD SUBTRACT *
* 04/19/10 (dkc) *
* *
******************************************************************************/
unsigned int carry(unsigned int a, unsigned int b, unsigned int sum);
void subn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i,s,d;
unsigned int c[8192];
if (n>8192)
return;
for (i=0; i<n; i++) {
s=*(b+i);
*(b+i)=~s;
}
d=1;
for (i=n-1; i>0; i--) {
s=*(b+i)+d;
d=carry(*(b+i),d,s);
*(b+i)=s;
}
*b=*b+d;
for (i=n-1; i>0; i--) {
s=*(a+i)+*(b+i);
c[i]=carry(*(a+i),*(b+i),s);
*(b+i)=s;
}
*b=*a+*b;
for (i=n-2; i>0; i--) {
s=*(b+i)+c[i+1];
c[i]+=carry(*(b+i),c[i+1],s);
*(b+i)=s;
}
*b=*b+c[1];
return;
}
/******************************************************************************
* *
* 32*N-BIT COPY *
* 04/02/10 (dkc) *
* *
******************************************************************************/
void copyn(unsigned int *a, unsigned int *b, unsigned int n) {
unsigned int i;
for (i=0; i<n; i++)
*(b+i)=*(a+i);
return;
}
/******************************************************************************
* *
* SET THE LAST OF N WORDS *
* 04/02/10 (dkc) *
* *
******************************************************************************/
void setn(unsigned int *a, unsigned int b, unsigned int n) {
unsigned int i;
for (i=0; i<n-1; i++)
*(a+i)=0;
*(a+n-1)=b;
return;
}
/******************************************************************************
* *
* "OR" N WORDS TOGETHER *
* 04/02/10 (dkc) *
* *
******************************************************************************/
unsigned int orn(unsigned int *a, unsigned int n) {
unsigned int i,b;
b=*a;
for (i=1; i<n; i++)
b=b|*(a+i);
return b;
}
unsigned int halbhung(unsigned int l, unsigned int n, unsigned int *M,
unsigned int *N, unsigned int *sv, unsigned int *A,
unsigned int *B, unsigned int *C, unsigned int *D,
unsigned int *L, unsigned int *S, unsigned int m) {
unsigned int h,i,j,o,a,b,offset,count,temp,maxoff;
unsigned int flag1,sign;
int g;
count=0; // odd element count
setn(D, 0, m); // clear maximum odd element
//
// compute rotated parity vector using ceiling function
//
maxoff=0;
for (offset=0; offset<l; offset++) {
for (j=1; j<=l; j++) {
a=j*n;
if (a==((a/l)*l))
a=a/l;
else
a=(a/l)+1;
b=(j-1)*n;
if (b==((b/l)*l))
b=b/l;
else
b=(b/l)+1;
temp=j+offset;
if (temp>l)
temp=temp-l;
sv[temp-1]=a-b;
}
if (sv[0]==0) // continue if even element
continue;
//
// compute odd element
//
setn(S, 0, m);
temp=0;
for (j=1; j<=l; j++) {
if (sv[j-1]!=0) {
setn(A, 1, m);
for (h=0; h<(n-temp-1); h++) {
copyn(A, B, m);
addn(A, A, m);
addn(B, A, m);
}
lshiftn(A, B, j-1, m);
addn(B, S, m);
if ((S[0]&0xc0000000)!=0) {
return(1);
}
}
temp=temp+sv[j];
}
//
// find maximum odd element
//
copyn(S, A, m);
subn(D, A, m);
if ((A[0]&0x80000000)!=0) {
copyn(S, D, m);
maxoff=offset;
}
//
// save odd element
//
if (offset==0)
copyn(S, L, m);
count=count+1; // increment odd element count
}
if (count!=n) {
return(2);
}
//
// compute 2^(K+L)-3^K
//
setn(A, 1, m);
for (i=0; i<n; i++) { // 3**K
copyn(A, B, m);
addn(A, A, m);
addn(B, A, m);
}
setn(B, 1, m);
lshiftn(B, C, l, m); // 2**(K+L)
subn(C, A, m); // 2**(K+L)-3**K
sign=0;
if ((A[0]&0x7fffffff)!=0) { // |2**(K+L)-3**K|
setn(B, 0, m);
subn(B, A, m);
sign=1;
}
//
copyn(L, S, m);
if ((S[0]&0xffc00000)!=0) {
return(3);
}
lshiftn(S, C, 2, 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(C, S, g, m);
shiftn(A, B, g, m);
flag1=g;
}
else {
copyn(C, S, m);
copyn(A, B, m);
}
o=orn(S, m-4);
if ((o|(S[m-4]&0xc0000000))!=0) {
return(4);
}
div128_64(S[m-4],S[m-3],S[m-2],S[m-1],M,B[m-2],B[m-1]); // minimum/(2**l-3**n)
//
copyn(D, S, m);
if ((S[0]&0xffc00000)!=0) {
return(5);
}
lshiftn(S, C, 2, m);
if (g>0)
shiftn(C, S, g, m);
else
copyn(C, S, m);
o=orn(S, m-4);
if ((o|(S[m-4]&0xc0000000))!=0) {
return(6);
}
div128_64(S[m-4],S[m-3],S[m-2],S[m-1],N,B[m-2],B[m-1]); // maximum/(2**l-3**n)
return(0);
}