/* $Id: hexcuberoots.c $ */

/* Compute the cube root of a number in hexadecimal format. */
/*
Use this to compute the first 64 bits of the fractional parts of the cube roots 
of the first eighty prime numbers as used for the SHA-384 and SHA-512 constants 
in section 4.2.3 of FIPS PUB 180-4 "Secure Hash Standard (SHS)", March 2012.
(and for the sixty-four 32-bit SHA-224 and SHA-256 ones, too, in section 4.2.2)
*/

/*
 * $Date: 2013-06-05 12:53Z $
 * $Revision: 1.0.0 $
 * $Author: dai $
 */

/**************************** COPYRIGHT NOTICE ********************************
  This code was originally written by David Ireland and is copyright (C)
  2013 DI Management Services Pty Ltd <www.di-mgt.com.au>. It is provided
  "as is" with no warranties. Use at your own risk. You must make your own
  assessment of its accuracy and suitability for your own purposes. It is
  not to be altered or distributed, except as part of an application. You
  are free to use it in any application, provided this copyright notice is
  left unchanged.  
************************* END OF COPYRIGHT NOTICE ****************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include "bigd.h"

/* Compile with bigd.c, bigdigits.c, bigd.h, bigdigits.h and bigdtypes.h available
   for free from http://di-mgt.com.au/bigdigits.html
*/

char *cbrthex(const char *numstr)
{	/* Returns pointer to static string, or the empty string if error */
	static char buf[1024];
	BIGD root;
	BIGD R, S, r0;
	BIGD A, B;
	unsigned int d, triple;
	const char *cp;

	/* Easier to deal with strings of digits rather than a large number */
	if (strlen(numstr) % 3)
	{
		printf("ERROR: number of digits must be a multiple of 3\n");
		return "";
	}

	root = bdNew();
	R = bdNew();
	S = bdNew();
	r0 = bdNew();
	A = bdNew();
	B = bdNew();

	// root = R = 0;
	bdSetZero(root);
	bdSetZero(R);

	for (cp = numstr; *cp; cp += 3)
	{
		/* Get next triple of digits */
		char block[4];
		block[0] = *cp;
		block[1] = *(cp+1);
		block[2] = *(cp+2);
		block[3] = '\0';
		triple = strtol(block, NULL, 16);
		//printf("Triple = %d\n", triple);
		//R = R * 16 * 16 * 16 + triple;
		bdShortMult(R, R, 16*16*16);
		bdShortAdd(R, R, triple);
		/* Compute constants A and B */
		// A = 3 * 16 * 16 * root * root
		bdSquare(A, root);
		bdShortMult(A, A, 3 * 16 * 16);
		// B = 3 * 16 * root
		bdShortMult(B, root, 3 * 16);
		bdSetZero(r0);
		/* Find largest value d in [0,15] that satisfies the inequality
		(3 * 16 * 16 * root * root + 3 * 16 * root * d + d * d) * d <= R
		*/
		for (d = 1; d < 16; d++)
		{
			//if ((A + B*d + d*d)*d > R)
			bdShortMult(S, B, d);	// S = B * d
			bdAdd(S, S, A);			// S += A
			bdShortAdd(S, S, d*d);	// S += d*d
			bdShortMult(S, S, d);	// S *= d
			if (bdCompare(S, R) > 0)
				break;
			bdSetEqual(r0, S);	/* Keep previous value for later */
		}
		d--;	/* We went one too far */
		//R -= (3 * 16 * 16 * root * root + 3 * 16 * root * d + d * d) * d;
		assert(bdCompare(R, r0) >= 0);
		bdSubtract(R, R, r0);

		//root = 16 * root + d;
		bdShortMult(root, root, 16);
		bdShortAdd(root, root, d);

	}
	/* Encode root in hex form */
	bdConvToHex(root, buf, sizeof(buf));

	bdFree(&root);
	bdFree(&R);
	bdFree(&S);
	bdFree(&r0);
	bdFree(&A);
	bdFree(&B);

	return buf;
}


int main(void)
{
	char *ptrbuf;
	unsigned int primes[] = {	/* The first eighty primes */
		2, 3, 5, 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, 
	};
	int nprimes = sizeof(primes)/sizeof(primes[0]);
	unsigned int prime;
	int i;
	char strbuf[64];
	char *padding = "000000000000000000000000000000000000000000000000";

	/* Do a check on known cubes */
	ptrbuf = cbrthex("001000000000000000000000000000000000000000000000000");
	printf("CHECK: cbrt(1)=%s\n", &ptrbuf[0]);
	ptrbuf = cbrthex("008000000000000000000000000000000000000000000000000");
	printf("CHECK: cbrt(8)=%s\n", &ptrbuf[0]);
	ptrbuf = cbrthex("07D000000000000000000000000000000000000000000000000");
	printf("CHECK: cbrt(125)=%s\n", &ptrbuf[0]);

	/* 002 = prime in hex + 48 x '0' digits => 16 hex digits in cube root => 16 x 4 = 64 bits */
	ptrbuf = cbrthex("002000000000000000000000000000000000000000000000000");
	/* Skip the first digit to get fractional part only */
	printf("%2d: frac_cbrt(2)=\t%s\n", 1, &ptrbuf[1]);

	/* Compose similar strings for each subsequent prime */
	for (i = 1; i < nprimes; i++)
	{
		prime = primes[i];
		sprintf(strbuf, "%03X%s", prime, padding);
		//printf("[%s]\n", strbuf);
		ptrbuf = cbrthex(strbuf);
		printf("%2d: frac_cbrt(%d)=\t%s\n", i+1, prime, &ptrbuf[1]);
	}

}