// Program to calculate Solinas's Joint Sparse Form

#include <iostream>

using namespace std;


// This is OK for Windows..

#define u64type unsigned __int64  

// in Linux
// #define u64type unsigned long long

// To output in binary in C++ I need this function

void output_in_binary(u64type x)
{
	int j,len,i=0;
	char s[80],w;
	for (j=0;j<64;j++) s[j]='0';
	
	while (x>0)
	{
		if (x%2==0) w='0';
		else w='1';
		s[i++]=w;
		x/=2;
	}
	for (i=0,j=63;i<=j;i++,j--)
	{
		w=s[i];
		s[i]=s[j];
		s[j]=w;
	}
	s[64]=0;
	cout << s << endl;
}

// Solinas's JSF algorithm

void jsf(u64type &k0,u64type &k1,u64type &u0p,u64type &u0m,u64type &u1p,u64type &u1m)
{
	int j,u0,u1,d0,d1,i0,i1;
	u64type w2,w3,b;
	d0=d1=0;

	w2=k0;
	w3=k1;
	u0p=u0m=u1p=u1m=0;

	for(j=0;;j++)
	{
		if (w2==0 && d0==0 && w3==0 && d1==0) break;  // finished
		i0=(w2+d0)&0x7;
		i1=(w3+d1)&0x7;
		if (i0%2==0) u0=0;
		else
		{
			u0=2-(i0%4);
			if ((i0==3 || i0==5) && i1%4==2) u0=-u0;
		}
		if (i1%2==0) u1=0;
		else
		{
			u1=2-(i1%4);
			if ((i1==3 || i1==5) && i0%4==2) u1=-u1;
		}
		b=((u64type)1<<j);
		if (u0>0) u0p+=b;
		if (u0<0) u0m+=b;
		if (u1>0) u1p+=b;
		if (u1<0) u1m+=b;

		if (2*d0==1+u0) d0=1-d0;
		if (2*d1==1+u1) d1=1-d1;
		
		w2/=2;
		w3/=2;
	}

}

int main()
{

	u64type k0,k1,u0p,u0m,u1p,u1m;
	
	cout << "Input k0 and k1 in Hex" << endl;
	cout << "k0= "; cin >> hex >> k0;
	cout << "k1= "; cin >> hex >> k1;

	cout << "k0,k1=" << endl;


	output_in_binary(k0);
	output_in_binary(k1);

// Call Solinas's JSF algorithm

	jsf(k0,k1,u0p,u0m,u1p,u1m);

	cout << "u0p, u0m, u1p, u1m" << endl;
	output_in_binary(u0p);
	output_in_binary(u0m);
	output_in_binary(u1p);
	output_in_binary(u1m);

	return 0;
}