/*
 * Get Victor Shoup's NTL-library (version 5.3.1) from http://shoup.net/ntl/
 *
 * And after installing the library, compile this program as
 * g++ -I/usr/local/include -L/usr/local/lib GF2Xfacs.c -o GF2Xfacs -lntl -lm
 *
 * where g++ --version gives:
 * g++ (GCC) 3.2.2
 * Copyright (C) 2002 Free Software Foundation, Inc.
 *
 * for the compiler what I have used now in January 2004.
 *
 */

#include <NTL/GF2XFactoring.h>

NTL_CLIENT


typedef unsigned long int ULLI;

/*
   n =  1:     x + 1                    11
   n =  2:    x² + 1                   101
   n =  3:    x³ + 1                  1001
   n =  4:   x^4 + 1                 10001
   n =  5:   x^5 + 1                100001
   n =  6:   x^6 + 1               1000001
   n =  7:   x^7 + 1              10000001 (1 byte)
   n =  8:   x^8 + 1             100000001 (2 bytes)
   n =  9:   x^9 + 1            1000000001 (2 bytes)
   n = 10:  x^10 + 1           10000000001 (2 bytes)

   n = 15:  x^15 + 1      1000000000000001 (2 bytes)
   n = 16:  x^16 + 1     10000000000000001 (3 bytes)

   Note: we should have: (list 3 3 137438953471 137438953471)
   instead of (3 3 4294967295 4294967295)
   for (Xn_1factor 74)
   I.e. use bignums, not the integers!
 */


void initXn_1(GF2X *p_z2p,ULLI n) /* n >= 1 */
{
    int i;
    unsigned char bytes[65];

    /* n divided by eight tells the number of intermediate zero-bytes needed. */
    for(i=0; i < (n >> 3);)
     {
       bytes[i++] = 0;
     }

    bytes[i] = 0;

    /* i now (n >> 3) */
    bytes[0] = 1; /* The constant coefficient 1. */

    /* Then the leading coefficient: */
    bytes[i++] |= (1 << (n & 7)); /* I.e. shift 1 to position (n mod 8). */

    /* i now (n >> 3)+1, the total number of bytes. */

    *p_z2p = GF2XFromBytes(bytes,(long)i);

    p_z2p->normalize();  /* Not necessary here?. */
}


void initFromBinexp(GF2X *p_z2p,ULLI n)
{
    int i = 0;
    unsigned char bytes[65];

    while(0 != n)
     {
       bytes[i] = (n&255);
       i++;
       n >>= 8;
     }

    *p_z2p = GF2XFromBytes(bytes,(long)i);

/*  p_z2p->normalize();  Not necessary here. */

}

ULLI GF2XtoBinexp(GF2X p_z2p)
{
    ULLI n = 0;
    int i;

    unsigned char bytes[65];
    long int n_bytes = NumBytes(p_z2p);
 
    BytesFromGF2X(bytes,p_z2p,n_bytes);

    for(i=n_bytes; i > 0; )
     {
       n <<= 8;
       n |= bytes[--i];
     }

    return(n);
}


ZZ GF2XtoZZ(GF2X p_z2p)
{
    ULLI n = 0;
    int i;

    unsigned char bytes[65];
    long int n_bytes = NumBytes(p_z2p);
 
    BytesFromGF2X(bytes,p_z2p,n_bytes);
    return(ZZFromBytes(bytes,n_bytes));
}


int main(int argc, char **argv)
{
   int i,upto_n;
   GF2X z2p;
   vec_pair_GF2X_long factors;
   const char *vec_name = "GF2X_FACVEC";
   int compute_Xn_1_factorizations = 0;

   if(argc <= 1) { cerr << "Usage: GF2Xfacs upto_n\n"; exit(1); }
   else { upto_n = atol(argv[1]); }

   if(upto_n < 0)
    {
      vec_name = "Xn_1_FACVEC";
      upto_n = -upto_n;
      compute_Xn_1_factorizations = 1;
    }

   cout << "(define "; cout << vec_name; cout << " (make-vector ";
   cout << (upto_n+1);  cout << "))\n\n";

   for(i=1; i <= upto_n; i++)
    {
      int j,u;
      if(compute_Xn_1_factorizations) { initXn_1(&z2p,i); }
      else { initFromBinexp(&z2p,i); }

      CanZass(factors, z2p);
      cout << "(vector-set! "; cout << vec_name; cout << " ";
      cout << i; cout << " (sort (list";
      u=factors.length();
      for(j=0; j<u; j++)
       {
         int h = factors[j].b;
         ZZ n = GF2XtoZZ(factors[j].a); /* Was: ULLI n = GF2XtoBinexp(factors[j].a); */
         while(h--) { cout << " "; cout << n; }
       }
      cout << ") <))\n";
    }
}