A C-code for permutation polynomials mod 2^n

In three recent threads [1] initiated by me, I have directly or
indirectly referred to the use of permutation polynomials mod 2^32.
In order to ease the work of persons needing an implementation of such
polynomials, I am taking the liberty to post my own C-code below. The
function invpermpoly is not required for issues of the aforementioned
threads, but could be of some interest on other occasions. (It is
used, though, in the demonstration part of the code below.)

M. K. Shen


[1] A poorman's block encryption algorithm, 11.03.2010
Update of my old idea on random number generation, 20.03.2010
A poorman's stream encryption algorithm, 21.03.2010


#include <stdio.h>
#include <stdlib.h>
#include <assert.h>

// permpoly.c Permutation polynomials mod 2^n
// Version 1.0 Release date: 21.3.2010.
// Please send comments and critiques to mok-kong.shen@xxxxxxxxxxxx

typedef unsigned long int ULI;

// Assumed number of bits of ULI (please check on 64-bit OS and
// eventually modify the following statement)
const int bnULI=32;

// Array permpoly contains m coefficients of a polynomial P(x) of
// degree m-1 in ascending order, m>=2. modulus is 2^n, n>=2.
// We check creteria given in:
// V. Anashin, A. Khrennikov, Applied Algebraic Dynamics, Berlin,
// 2009, p283:
// c3 + c5 + c7 + ... = 2*c2 mod 4
// c4 + c6 + c8 + ... = c1 + c2 - 1 mod 4
// c1 = 1 mod 2
// c0 = 1 mod 2
// When these criteria are satisfied, sequence x_(i+1)=P(x_i) is full
// period. If this full period property is not required, then the
// criteria may be weakened to those given in R. L. Rivest, Finite
// Fields and Appl. 7 (2001), pp.287-292:
// c1 = 1 mod 2
// c2 + c4 + c6 + ... = 0 mod 2
// c3 + c6 + c8 + ... = 0 mod 2
// In this case, the procedures checkpermpoly and correctpermpoly below
// have to be modified correspondingly, while the procedures
// evalpolymod and invpermpoly remain valid.
// To use permpoly as PRNG, it may be desirable to avoid coefficients
// lying within 5% of both ends of the modulus for practical reasons
// (a literature attributes this advice to Knuth). In this case, input
// to checkpermpoly should have already satisfied that condition.

void checkpermpoly(ULI *permpoly, int m, int n, int *err)
{ int i;
ULI s1,s2,s3,s4,ff;
assert(m>=2 && n>=2 && n<=bnULI);
if (n<bnULI)
{ ff=(~0)<<n;
for (i=0; i<m; i++) if (permpoly[i]&ff)
{ printf("Warning: Coefficients out of ranges\n"); break;
s1= m==2 ? 0 : -2*permpoly[2];
for (i=3; i<m; i+=2) s1+=permpoly[i]; s1&=3;
s2= m==2 ? -permpoly[1]+1 : -permpoly[1]-permpoly[2]+1;
for (i=4; i<m; i+=2) s2+=permpoly[i]; s2&=3;
s3=(permpoly[1]-1)&1; s4=(permpoly[0]-1)&1;
if (s1!=0 || s2!=0 || s3!=0 || s4!=0)
{ printf("checkpermpoly failed: %lu %lu %lu %lu\n",s1,s2,s3,s4);
else *err=0;

// correctpermpoly may be called after checkpermpoly fails, i.e. when
// the user-given polynomial doesn't satisfy the criteria mentioned
// above. If a given permpoly fails checkpermpoly, there can in general
// be many different ways of correction. The following is only one
// particular way of correction, in that we impose permpoly[1]=1 mod 4
// and permpoly[2]=0 mod 4. Further we correct permpoly[3] and
// permpoly[4] and leave the others unchanged.
// After call, permpoly is modified and satisfies the criteria. User
// should preferably examine whether this particular type of correction
// is acceptable for purposes of his application.

void correctpermpoly(ULI *permpoly, int m, int n)
{ int i,j;
ULI s,ff;
ff= n<bnULI ? ~((~0)<<n) : ~0;
s=permpoly[1]&3; if (s!=1) permpoly[1]=(permpoly[1]+5-s)&ff;
if (m>2)
{ s=permpoly[2]&3; if (s!=0) permpoly[2]=(permpoly[2]+4-s)&ff;
for (i=3; i<=4; i++) if (i<m)
{ s=permpoly[i]; for (j=i+2; j<m; j+=2) s+=permpoly[j]; s&=3;
if (s!=0) permpoly[i]=(permpoly[i]+4-s)&ff;

// evalpolymod evaluates a general polynomial poly of degree m-1 mod 2^n
// with argument x.

ULI evalpolymod(ULI *poly, int m, int n, ULI x)
{ int i;
ULI t=poly[m-1];
for (i=m-1; i>0; i--) t=t*x+poly[i-1];
if (n<bnULI) return t&(~((~0)<<n)); else return t;

// Given y, invpermpoly solves permpoly(x)=y mod 2^n for x.
// permpoly is assumed to have passed checkpermpoly.

ULI invpermpoly(ULI *permpoly, int m, int n, ULI y)
{ int i;
ULI x;
// First derivative of permpoly(x)=1 mod 2. (See also reference above
// p.306.) Hensel's lifting lemma leads to the iteration
// x_(i+1)=x_i-g(x_i) mod 2^(i+1), with g(x)=permpoly(x)-y. Since the
// solution is unique, we may compute in all steps mod 2^bnLUI and
// simply do a final mod 2^n.
x = ((permpoly[0]-y)&1)==0 ? 0 : 1;
for (i=2; i<=n; i++) x-=evalpolymod(permpoly,m,bnULI,x)-y;
if (n<bnULI) return x&(~((~0)<<n)); else return x;

// The following code lines are for demonstration only and show the
// correctness of the implementation.

void showvec(int n, ULI *vec)
{ int i;
int n1=4;
for (i=0; i<n; i++)
{ if (i!=0 && i%n1==0) printf("\n");
printf("%lu ",vec[i]);

int main()
{ int m,n,err;
ULI x,y,xn,xa,xb,xmax;

//Test example, not satisfying criteria
ULI permpoly[4]={1,3,2,5}; m=4; n=8;


if (err!=0)
{ printf("Coefficients error\n");
if (err!=0)
{ printf("error in correctpermpoly\n"); return(1);
printf("Corrected coefficients:\n");

xmax= n<bnULI ? ~((~0)<<n) : ~0;

printf("Input range of checking: xa xb in [0,%lu]\n",xmax);
if (xb>xmax)
{ printf("xb out of range\n"); goto uu;

x=xa; err=0;

{ y=evalpolymod(permpoly,m,n,x); xn=invpermpoly(permpoly,m,n,y);
// xn should equal x
if (xn!=x)
{ err=1; printf("Deviations: *****");
printf("x=%lu, y=%lu, xn=%lu\n",x,y,xn);
//This is needed if xb==xmax and n=bnULI
if (x==xmax) break;

if (err)
{ printf("permpoly and invpermpoly don't check\n"); return(2);
printf("Check for the given range o.k.\n");

Relevant Pages