﻿ /*******************************
```/******************************************************************************/
/*									      */
/*  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) {
if (mode==0)
a=~a;
if (a==0)
return(32);
count=0;
for (i=0; i<32; i++) {
break;
count+=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) {
i=0;
while ((*(a+i)==0)&&(i<n))
i=i+1;
shift=32*i;
if (i==n)
return shift;
shift=shift+1;
}
return shift;
}
/******************************************************************************
*									      *
*  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);
}
lshiftn(A, B, j-1, 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);
}
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);
}
```