#include <stdio.h>

#define max(a,b)(((a)>(b))?(a):(b))

/*
This script is used to find product bases whose orthogonality graphs have a given number of parties each with three sets of two states that are equal on that party.
This script finds that it is impossible to have 5 or more such parties, which plays a crucial role in the proof of case 3(b) of Lemma 6 of the paper "The Minimum Size of Qubit Unextendible Product Bases". This script can also show that it is impossible to have 3 or more such parties if we also require that those three sets of two states have common endpoints, which is important in cases 4(b) and 4(c) of the same paper.

The script works by recursively testing all possible configurations of 3 sets of 2 equal states on the requested number of parties. The fillParty function does the majority of the difficult work. i[0] and i[1] are vertices in one set of states that are equal to each other on the current party, i[2] and i[3] form another such set of vertices, and i[4] and i[5] form the final such set of vertices.

Author: Nathaniel Johnston (nathaniel@njohnston.ca)
Version: 1.00
URL: http://www.njohnston.ca/publications/qubit-upbs/code/
Last Updated: Feb. 5, 2013
*/

unsigned short int maxV[100],two[100][6];
unsigned int ct,case4,s;
FILE *fp;

void fillParty(unsigned short int curP);
unsigned short int checkUnextend(unsigned short int curP);

int main()
{
    printf("This tool will try to help you find (or prove non-existence of) a UPB of size 4k+3 in (C^2)^{otimes 4k}.\nSpecifically, this script looks for arrangements of three sets of two vectors being equal to each other on s parties that could potentially lead to an unextendible product basis.\n\nPlease enter s (s >= 3): ");
    scanf("%d",&s);
    printf("\nIf you wish to perform the default search (required in case 3(b) of Lemma 6), enter \"0\". If you wish to perform the search that has all sets of 2 equal states having common endpoints (required in cases 4(b) and 4(c) of Lemma 6), enter \"1\".\n\nType 0 or 1: ");
    scanf("%d",&case4);

    /* WLOG, we can assume that on the first party, the three groupings of two identical qubits happen on vertices 0-1, 2-3, and 4-5. */
    two[1][0] = 0;
    two[1][1] = 1;
    two[1][2] = 2;
    two[1][3] = 3;
    two[1][4] = 4;
    two[1][5] = 5;
    maxV[1] = 5;
    ct = 0;

    if(case4 == 1){
        fp = fopen("upb_case4_output.txt", "w");
    }else{
        fp = fopen("upb_case3b_output.txt", "w");
    }
    fillParty(2);
    fclose(fp);

    if(ct > 0){
        printf("\nDone searching! Found %d valid arrangements (note that many of these will be equivalent under relabelling parties and vertices). Output written to file \"upb_case3b_output.txt\".\n",ct);
    }else{
        printf("\nDone searching! There is no arrangement of sets of 2 equal states satisfying the specified requirements.\n");
    }
    getchar();
    return 0;
}

void fillParty(unsigned short int curP)
{
    unsigned short int i[6];
    unsigned short int lLim2,lLim4,uLim0,uLim2,uLim3,uLim4,uLim5;
    unsigned int j;
    if(curP == 2 || case4 == 1){
        uLim0 = 0; /* Upper limit of 0 because of WLOG arguments in this case. */
    }else{
        uLim0 = 3; /* Upper limit of 3 instead of 5 because there must be room for vertices 2 and 4 after vertex 0. */
    }

    /* Start by placing vertices 0, 2, and 4. Then place vertices 1, 3, and 5. */
    for(i[0]=0; i[0]<=uLim0; i[0]++){
        if(case4 == 1){
            lLim2 = 2; uLim2 = 2;
        }else{
            lLim2 = i[0]+1; uLim2 = 4;
        }
        for(i[2]=lLim2; i[2]<=uLim2; i[2]++){ /* Upper limit of 4 instead of 5 because there must be room for vertex 4 after vertex 2. */
            if(case4 == 1){
                lLim4 = 4; uLim4 = 4;
            }else{
                lLim4 = i[2]+1; uLim4 = 5;
            }
            for(i[4]=lLim4; i[4]<=uLim4; i[4]++){
                for(i[1]=i[0]+1; i[1]<=maxV[curP-1]+1; i[1]++){
                    if(i[1] != i[2] && i[1] != i[4]){ /* Make sure vertex 1 isn't the same as vertex 2 or 4 (the for loop bounds make sure it's not the same as vertex 0). */
                        uLim3 = maxV[curP-1]+1;
                        uLim5 = maxV[curP-1]+1;
                        if(i[1] == maxV[curP-1]+1){uLim3++;uLim5++;}
                        for(i[3]=i[2]+1; i[3]<=uLim3; i[3]++){
                            if(i[3] != i[1] && i[3] != i[4]){ /* Make sure vertex 3 isn't the same as vertex 1 or 4 (the for loop bounds make sure it's not the same as vertex 0 or 2). */
                                if(i[3] == uLim3){uLim5++;}

                                /* If on the second party, ignore cases where no group of 2 contains vertex 4 or 5 (these cases are isomorphic to cases where no group of two contains vertex 2 or 3). */
                                if(curP == 2 && (i[1]<=3 || i[1]>=6) && (i[3]<=3 || i[3]>=6) && i[2]<=3 && i[4]<=3){uLim5 = 5;}

                                for(i[5]=i[4]+1; i[5]<=uLim5; i[5]++){
                                    if(i[5] != i[1] && i[5] != i[3]){ /* Make sure vertex 5 isn't the same as vertex 1 or 3 (the for loop bounds make sure it's not the same as vertex 0, 2, or 4). */
                                        two[curP][0] = i[0];
                                        two[curP][1] = i[1];
                                        two[curP][2] = i[2];
                                        two[curP][3] = i[3];
                                        two[curP][4] = i[4];
                                        two[curP][5] = i[5];
                                        maxV[curP] = max(maxV[curP-1],max(i[0],max(i[1],max(i[2],max(i[3],max(i[4],i[5]))))));

                                        if(curP <= 2 || (curP < s && checkUnextend(curP)==1)){
                                            fillParty(curP+1);
                                        }else if(curP == s && checkUnextend(curP)==1){
                                            fprintf(fp,"convert(convert([");
                                            if(ct == 0){
                                                printf("\nFound a valid configuration of two equal states:");
                                            }
                                            for(j=1; j<=s; j++){
                                                fprintf(fp,"{{%d,%d},{%d,%d},{%d,%d}}",two[j][0]+1,two[j][1]+1,two[j][2]+1,two[j][3]+1,two[j][4]+1,two[j][5]+1);
                                                if(j<s){
                                                    fprintf(fp,",");
                                                }
                                                if(ct == 0){
                                                    printf("\n%d,%d   %d,%d   %d,%d",two[j][0],two[j][1],two[j][2],two[j][3],two[j][4],two[j][5]);
                                                }
                                            }
                                            if(ct == 0){
                                                printf("\n\nAny further configurations that are found will be written to the output file, but not displayed here.\n");
                                            }
                                            fprintf(fp,"],multiset),set),\n");
                                            ct++;
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
    }
}

unsigned short int checkUnextend(unsigned short int curP)
{
    unsigned short int i[3],p[3];

    for(p[0]=1; p[0]<=curP-2; p[0]++){
        for(p[1]=p[0]+1; p[1]<=curP-1; p[1]++){
            for(p[2]=p[1]+1; p[2]<=curP; p[2]++){
                for(i[0]=0; i[0]<=4; i[0]+=2){
                    for(i[1]=0; i[1]<=4; i[1]+=2){
                        /* Check if there's some overlap between the first and second party. */
                        if(two[p[0]][i[0]]!=two[p[1]][i[1]] && two[p[0]][i[0]]!=two[p[1]][i[1]+1] && two[p[0]][i[0]+1]!=two[p[1]][i[1]] && two[p[0]][i[0]+1]!=two[p[1]][i[1]+1]){
                            for(i[2]=0; i[2]<=4; i[2]+=2){
                                /* Check if there's some overlap between the third party and one of the other parties. */
                                if(two[p[2]][i[2]]!=two[p[0]][i[0]] && two[p[2]][i[2]]!=two[p[0]][i[0]+1] && two[p[2]][i[2]]!=two[p[1]][i[1]] && two[p[2]][i[2]]!=two[p[1]][i[1]+1] && two[p[2]][i[2]+1]!=two[p[0]][i[0]] && two[p[2]][i[2]+1]!=two[p[0]][i[0]+1] && two[p[2]][i[2]+1]!=two[p[1]][i[1]] && two[p[2]][i[2]+1]!=two[p[1]][i[1]+1]){
                                    return 0;
                                }
                            }
                        }
                    }
                }
            }
        }
    }

    return 1;
}