#
# cycpolcheck.py -- by Antti Karttunen (his-firstname.his-surname@gmail.com)
# Written May 03 2006 and placed in Public Domain.
#
# This is a small Python-program to post-filter the solution-candidates to
# a puzzle presented in New Scientist (some issue of August, September 2005 ?)
# or equivalently, to MuddMath "Bernoff's Puzzler"
# which can be found on the page 12/20 of the following PDF-file:
# http://www.math.hmc.edu/muddmath/v4n1.pdf
# (or page 10/18 of the original). 

# See September 11 - 14, 2005 postings on the SeqFan mailing list
# by Brendan McKay, Edwin Clark, T. D. Noe, Hugo Pfoertner, et al.

# This program is based on the idea given by T. D. Noe
# on his May 2 2006 SeqFan posting:
# "The computations can be done without floating-point arithmetic by using
# polynomial GCD.  As mentioned in a previous e-mail, for a permutation
# (w1,w2,..wn) to be a solution, the polynomial P(x)=w1+w2*x+...+wn*x^(n-1)
# must have e^(2 pi I/n) as a root.  This means that P(x) and the nth
# cyclotomic polynomial have a common zero or, equivalently,
# gcd(P(x),Phi(n,x)) = Phi(n,x) because Phi(n,x) is irreducible. This means
# Phi(n,x) divides P(x) or, equivalently, P(x)=0 (mod Phi(n,x)). These are
# all computable via integer operations."

# Currently the program (function checkfile) assumes that
# its input is in the format produced by Hugo Pfortner's
# Fortran-program.

# Usage:
# from cycpolcheck import *
# checkfile('bc15.txt') # (for example...)
# or checkdir('.','txt') # if running checkfile for every *.txt file in this directory.

# Program uses the poly1d class from SciPy-module, which can be downloaded
# from http://www.scipy.org/ and its documentation found at:
# http://www.tau.ac.il/~kineret/amit/scipy_tutorial/

import os
import re
from math import *
from scipy import poly1d



# Positions of the positive and negative coefficients
# in cyclotomic polynomials, shown in binary:
#
# http://www.research.att.com/~njas/sequences/A063697
# http://www.research.att.com/~njas/sequences/A063699
# 
# 	A063697		A063699
# 0 	10	 	0
# 1	10		1
# 2	11		0
# 3	111		0
# 4	101		0
# 5	11111		0
# 6	101		10
# 7	1111111		0
# 8	10001		0
# 9	1001001		0
# 10	10101		1010
# 11	11111111111	0
# 12	10001		100
# 13	1111111111111	0
# 14	1010101		101010
# 15	100101001	10010010
# 16	100000001	0
# 17	11111111111111111	0
# 18	1000001		1000
# 19	1111111111111111111	0
# 20	100010001	1000100
# 21	1001001001001	100100010010
# 22	10101010101	1010101010

# Same positions in decimal (for n=0..30):

A063696 = [0,2,3,7,5,31,5,127,17,73,21,2047,17,8191,85,297,
           257,131071,65,524287,273,4681,1365,8388607,257,
           1082401,5461,262657,4369,536870911,387]

A063698 = [0,1,0,0,0,0,2,0,0,0,10,0,4,0,42,146,0,0,8,0,68,
           2322,682,0,16,0,2730,0,1092,0,56]



def floor_log2(n):
  '''See http://www.research.att.com/~njas/sequences/A000523'''
  if n < 1: return(0)
  else: return(int(log(n)/log(2)))

def make_plusminusone_coeff_polynom(poscoeffs,negcoeffs):
  '''Creates a polynomial of +-1 coefficients from binary masks poscoeffs and negcoeffs.'''
  coeffs = [((poscoeffs>>i)&1)-((negcoeffs>>i)&1)
            for i in range(1+max(floor_log2(poscoeffs),floor_log2(negcoeffs)))]
  coeffs.reverse()
  return(poly1d(coeffs))

CyclotomicPolynomials = map(make_plusminusone_coeff_polynom,A063696,A063698)

ZeroPolynom = poly1d([0])

def nthcycpoldivides(n,pol):
  '''Returns True if the n:th Cyclotomic Polynomial divides Polynom pol'''
  (quot,rema) = pol/CyclotomicPolynomials[n]
  return(rema == ZeroPolynom)

# We search for lines like this:
#  0.130126187626938E+04    1    1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

# The output is written to similarly named file,
# with the extension replaced by a new extension ".sol"

def checkfile(filename):
  '''Opens file "filename" for reading and similar file with extension ".sol" for writing.'''
  infp = open(filename,'r')
  try:
    lastdot = filename.rindex('.')
    outfilename = filename[:lastdot] + ".sol"
  except ValueError: # There were no dot in the filename
    outfilename = filename + ".sol"

  outfp = open(outfilename,'w')
  permlinepat = re.compile(r'^\s*(-?0\.[-+0-9E]*)\s+(.*)')
  sols = 0
  firstmax = 0
  for line in infp.xreadlines(): 
    m = permlinepat.match(line)
    if(m):
      sum = m.group(1) # Or is it norm, i.e. the distance from center?
      coeffs = map(int,(m.group(2)).split())
      maxnow = max(coeffs[1:]) # First term after 'sum', skip it, Hugo knows what it is.
      if(0 == firstmax): firstmax = maxnow
      elif(maxnow != firstmax):
        print "Skipping the following illegal line because its max. coeff. is",maxnow,"not",firstmax
        print line

      p = poly1d(coeffs[1:])
      if(nthcycpoldivides(maxnow,p)):
        sols += 1
        output = "Sol " + str(sols) + ": " + line.lstrip()
        outfp.write(output)

  print filename,sols,"solutions found."
  outfp.close()


def checkdir(path,suff):
  '''Executes checkdir for all the files ending with "suff" in directory "path"'''
  olddir = os.getcwd()
  os.chdir(path)
  allfiles = os.listdir('.')
  sufflen = len(suff)
  files = filter(lambda(s): (suff == s[-sufflen:]),allfiles)
  for file in files: checkfile(file)
  os.chdir(olddir)


# This module ends here.