/* pack.c: Illustrate an efficient "information packing" function.
*
* By Terry R. McConnell
*
* Let W = {0,1,2,...} be the set of whole numbers. A recursive function
* C: WxW --> W is an information packing function if it is 1-1 and there
* are recursive functions T ("tail") and H ("head") such that
*
* H(C(x,y)) = x and T(C(x,y)) = y.
*
* Such sets of functions are of interest in recursion theory because they
* allow functions of many variables to be treated as functions of one
* variable.
*
* One example of an information packing function is C(x,y) = [(x+y)^2 +
* (x+y) + 2x]/2. (Note that the expression in brackets is always even.) This
* example is particularly efficient because it maps WxW onto W. This
* program illustrates this example by implementing the head and tail functions
* as recursive, or primitive recursive functions. The head and tail are
* documented below with their implementations.
*
* To see that C(x,y) is a bijection, it suffices to show 2C(x,y) is a
* bijection from WxW onto the even integers. Introduce the auxiliary
* function v(n) = n^2 + n, which is clearly strictly increasing onto
* a certain subset of the even whole numbers. It is easy to see that
* with n = x+y,
*
* 2C(x,y) = v(n) + 2x < v(n+1), 0 <= x <= n, n = 0,1,...
*
* and over the stated range, (n,x) covers the set of all even
* whole numbers exactly once. The desired result follows, since (x,y) ->
* (x+y,x) is 1-1.
*
* The function C, and variants of it, appear in a number of texts on the theory
* of computation. See, e.g., problem 10.3-4 in Marvin Minsky, Computation:
* finite and infinite machines, Prentice Hall, Englewood Cliffs, N.J., 1967.
*
* Run with -h option to see usage information.
*
*/
/* compile: cc -o pack pack.c
Use -D_SHORT_STRINGS if your compiler does not support multiline
string constants.
These are very slow and are included only to prove a point:
Use -DPRIMITIVE_RECURSIVE to implement all functions as primitive
recursive functions.
Use -DRECURSIVE to implement all functions as recursive functions.
*/
/* Set these appropriately for your own platform: */
#define _GNU_SOURCE
#define USE_LONG_LONG /* or USE_LONG or USE_INT */
#include
#include
#include
#ifdef USE_LONG_LONG
#define CONVERT atoll
#define TYPEMAX ULONG_LONG_MAX
#define PRINT_AS "llu" /* conversion after % in printf */
typedef unsigned long long whole;
#endif
#ifdef USE_LONG
#define CONVERT atol
#define TYPEMAX ULONG_MAX
#define PRINT_AS "lu"
typedef unsigned long whole;
#endif
#ifdef USE_INT
#define CONVERT atoi
#define TYPEMAX UINT_MAX
#define PRINT_AS "u"
typedef unsigned int whole;
#endif
#define VERSION "1.0"
#define USAGE "pack [ -1 -h -v] [x y]"
#ifndef _SHORT_STRINGS
#define HELP "\n\npack [ -1 -h -v ] [x y]\n\n\
Pack information in whole numbers x and y into a single whole number z and \n\
print the result. \n\n\
-1: Unpack the information in the whole number z and print result as x y.\n\
-v: Print version number and exit. \n\
-h: Print this helpful information. \n\n"
#else
#define HELP USAGE
#endif
/* The implementation of the packing function, C(x,y). */
whole pack( whole x, whole y)
{
return ((x+y)*(x+y) + (x+y) + 2*x)/2;
/* Well, that was pretty easy, wasn't it? */
}
/* The auxiliary function V */
whole V(whole z){
return z*z + z;
}
/* U, an integer valued "inverse" of V. I.e., U(x) is the unique whole
* number n that satisfies
*
* (*) V(n) <= x < V(n+1).
*
* */
#ifdef PRIMITIVE_RECURSIVE
/* Implement U as a primitive recursive function. I claim that the
* following is a recursion for U:
*
* U(x+1) = U(x) + I(U^2(x) + 3U(x) + 1 - x)
* U(0) = 0
*
* where I is the indicator of the set {0} and the subtraction has
* a - b = 0 when a < b. It suffices to prove that the function defined by
* this recursion satisfies (*). We do this by induction on x. For x = 0 this
* is obvious, so suppose (*) holds for n = U(x). We must show that
*
* V(U(x+1)) <= x+1 < V(U(x+1)+1) (**)
*
* Case 1: U^2(x) + 3U(x) + 1 = x. Then U(x+1) = U(x) + 1, so
* V(U(x+1)) = U^2(x) + 3U(x) + 1 = x < x+1 so the left side of (**) holds
* with <. OTOH V(U(x+1)+1) = V(U(x)+2) = U^2(x) + 6U(x) + 4 which is clearly
* > x+1.
*
* Case 2: U^2(x) + 3U(x) + 1 > x (This is the only alternative allowed
* by the inductive hypothesis.) Then U(x+1) = U(x) so V(U(x+1)) =
* V(U(x)) <= x by the inductive hypothesis. OTOH, V(U(x+1)+1) = V(U(x)+1)
* = U^2(x) + 3U(x) + 2 > x + 1.
*
* For large z this is very slow and consumes lots of stack. */
whole U(whole z){
whole w;
if(z == 0)return 0;
w = U(z-1);
if(w*w + 3*w + 1 == z)return w + 1;
return w;
}
#elif defined( RECURSIVE )
/* We use the "mu" operator, i.e., linear search for the largest n such
* V(n) <= z. For large z this is slow. */
whole U(whole z){
whole i = 0;
while(V(i+1) <= z)i++;
return i;
}
#else
/* Implement U `efficiently'. Since this may be passed the largest possible
* number representable, we must be careful not to generate anything larger
* than z during the calculation. We use a bisection method and are
* careful about the possibililty of overflow.
*/
whole U(whole z){
whole top = z, bot = 0,m,n;
while(1){
/* Compute the integer average of top and bot */
if((top%2)&&(bot%2))
m = 1 + top/2 + bot/2;
else m = top/2 + bot/2;
/* Are we done ? */
if((top==m)||(bot==m))break;
/* Head off the possibility that V(m) will overflow. */
/* First consider the possibility that m*m overflows. */
n = m*m;
if((n % m) || (n == m) || (n == 0)) {
top = m;
continue;
}
/* OK, m*m does not overflow, but V(m) still could. */
if(m > TYPEMAX - n){
top = m;
continue;
}
/* OK, V(m) = n + m really is what it says it is. */
if(n + m > z)
top = m;
else bot = m;
}
return m;
}
#endif
whole head(whole z){
return (2*z - V(U(2*z)))/2;
}
whole tail(whole z){
return U(V(U(2*z))) - head(z);
}
int
main(int argc, char **argv)
{
int j=0;
whole x,y,z,w,t;
/* Process command line */
while(++j < argc){
if(argv[j][0] == '-')
switch(argv[j][1]){
case '1':
z = (whole)CONVERT(argv[j+1]);
j++;
if(z > TYPEMAX/2)fprintf(stderr,"pack: value is too big. May not unpack correctly.\n");
printf("%"PRINT_AS" %"PRINT_AS"\n",head(z),tail(z));
return 0;
case 'v':
case 'V':
printf("%s\n",VERSION);
exit(0);
case '?':
case 'h':
case 'H':
printf("%s\n",HELP);
#if !defined(RECURSIVE) && !defined(PRIMITIVE_RECURSIVE)
t = U(TYPEMAX);
w = TYPEMAX - V(t);
if(2*t > w)t = t-1;
printf("Maximum x + y = %"PRINT_AS".\n",t);
#endif
printf("Maximum value handled by -1 option: z = %"PRINT_AS".\n\n",TYPEMAX/2);
return 0;
default:
fprintf(stderr,"pack: unkown option %s\n",
argv[j]);
return 1;
}
else break;
}
if(j < argc)
x = (whole)CONVERT(argv[j++]);
else {
fprintf(stderr,"%s\n",USAGE);
exit(1);
}
if(j < argc)
y = (whole)CONVERT(argv[j++]);
else {
fprintf(stderr,"%s\n",USAGE);
exit(1);
}
/* Check to make sure the arguments won't produce an encoding that's
* too big. */
t = pack(x,y);
#if !defined(RECURSIVE) && !defined(PRIMITIVE_RECURSIVE)
if((x != head(t))||(y != tail(t)))
fprintf(stderr,"pack: warning: x or y is too big! Will not unpack correctly.\n");
#endif
printf("%"PRINT_AS"\n",t);
return 0;
}