/* (c) 2009 Julian D. A. Wiseman Released under GNU General Public License. Problem: http://www.jdawiseman.com/papers/easymath/open_questions.html#stars Please inform author of progress: contact details at http://www.jdawiseman.com/author.html Also see http://developer.apple.com/documentation/Darwin/Reference/ManPages/man3/math.3.html#//apple_ref/doc/man/3/math which, alas, does not refer to an extended-precision cos function. */ #include #include #include #include void hsort( void *base, size_t n, size_t size, int (*cmp)(const void *, const void *) ); int InSeries(unsigned long n1, unsigned long m1, unsigned long n2, unsigned long m2) ; inline double InnerRadius(unsigned long n, unsigned long m); void OutputMathematicaCode(FILE *fp_Mathematica, unsigned long n1, unsigned long m1, unsigned long n2, unsigned long m2); static double diff; /* saves time of repeated creation */ const double pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798 ; typedef struct { unsigned long n; unsigned long m; double IR; } star ; inline double InnerRadius(unsigned long n, unsigned long m) {return( cos((pi*m)/n) / cos((pi*(m-1))/n) );}; int CompareInnerRadii(const star * const s1, const star * const s2) { diff = s1->IR - s2->IR ; return( diff>0 ? 1 : (diff<0 ? -1 : 0) ); }; /* CompareInnerRadii */ inline int InSeries(unsigned long n1, unsigned long m1, unsigned long n2, unsigned long m2) { /* (6i-2)/i and (18i-6)/(6i-2) or (6i-4)/i and (18i-12)/(6i-3) */ if( m1 < m2 ) { return( (m2==n1 && n2==m2*3 && n1==6*m1-2) || (m2==n1+1 && n1==m1*6-4 && n2==m2*3-3) ); } else { return( (m1==n2 && n1==m1*3 && n2==6*m2-2) || (m1==n2+1 && n2==m2*6-4 && n1==m1*3-3) ); } }; /* InSeries */ void OutputMathematicaCode(FILE *fp_Mathematica, unsigned long n1, unsigned long m1, unsigned long n2, unsigned long m2) { fprintf( fp_Mathematica, "Print[{ %lu, %lu, %lu, %lu, N[ (Cos[Pi %lu/%lu]/Cos[Pi (%lu-1)/%lu]) - (Cos[Pi %lu/%lu]/Cos[Pi (%lu-1)/%lu]) , 100] }] \n", n1, m1, n2, m2, m1, n1, m1, n1, m2, n2, m2, n2 ); } /* OutputMathematicaCode */ int main() { star *s; unsigned long ArrayLength, k, k_prime, ArrayNumUsed, m, m_start, MaxNumPossibleMatches, NumPossibleMatches; signed long ArrayIndex; double IR, IR_upper, IR_lower, Epsilon; char filename_Mathematica[255]; FILE *fp_Mathematica; for(;;) { printf("How many stars in array? "); scanf ("%lu",&ArrayLength); printf("\nWorking with ArrayLength=%lu\n", ArrayLength); s = (star*) malloc( sizeof(star) * ArrayLength ) ; if( s == NULL ) {printf("malloc failed. Try a new length.\n");} else {printf("malloc succeeded.\n"); break;} } printf("Upper limit on number of possible matches? "); scanf("%lu",&MaxNumPossibleMatches); Epsilon = 0; printf("Epsilon enlarged at k ="); for( k=2; k<244609L; k++ ) /* k not having its usual meaning here; doing the role of i */ { /* (6i-2)/i and (18i-6)/(6i-2) or (6i-4)/i and (18i-12)/(6i-3) */ IR = fabs( InnerRadius(6*k-2,k) - InnerRadius(18*k-6,6*k-2) ); if( Epsilon < IR ) {Epsilon = IR; printf(" %lu",k);} IR = fabs( InnerRadius(6*k-4,k) - InnerRadius(18*k-12,6*k-3) ); if( Epsilon < IR ) {Epsilon = IR; printf(" %lu",k);} } /* k */ Epsilon *= 1.1; /* Arbitrary small increase */ printf(", so Epsilon=%le\n", Epsilon); sprintf(filename_Mathematica,"matching_stars_Mathematica_%lu.txt", time(0)%100000000L ); if(NULL==(fp_Mathematica=fopen(filename_Mathematica, "w"))) {fprintf(stderr,"fopen problem with '%s'.",filename_Mathematica); exit(EXIT_FAILURE);} printf("Output file = %s\n",filename_Mathematica); NumPossibleMatches = 0; for( k_prime=5 ; ; k_prime++ ) { /* Searching for matching inner radii between (k_prime-1)/(k_prime+1) and k_prime/(k_prime+2) */ IR_lower = ((double)(k_prime-1)) / (k_prime+1) ; IR_upper = ((double)k_prime) / (k_prime+2) ; ArrayNumUsed = 0; printf("k_prime=%lu, IR_lower=%0.17lf, IR_upper=%0.17lf\n", k_prime, IR_lower, IR_upper); m_start=2; for( k=1 ; k IR_lower ) { if( ArrayNumUsed >= ArrayLength ) goto ArrayFull; s[ArrayNumUsed].n = m_start*2+k ; s[ArrayNumUsed].m = m_start ; s[ArrayNumUsed].IR = IR; ArrayNumUsed++ ; for( m=m_start+1 ; ; m++ ) { IR = InnerRadius(m*2+k, m) ; if( IR <= IR_lower ) break; if( ArrayNumUsed >= ArrayLength ) goto ArrayFull; s[ArrayNumUsed].n = m*2+k ; s[ArrayNumUsed].m = m ; s[ArrayNumUsed].IR = IR; ArrayNumUsed++ ; } /* m */ } /* first one below ceiling is above floor */ } /* k */ /* printf("\n"); printf("k_prime=%lu\n", k_prime); printf("IR_lower=%0.17lf\n", IR_lower); printf("IR_upper=%0.17lf\n", IR_upper); printf("ArrayNumUsed=%lu\n", ArrayNumUsed); printf("n\tk\tm\tIR\n"); for( ArrayIndex=0; ArrayIndex= MaxNumPossibleMatches ) goto Finished; } } /* closer than Epsilon */ } /* ArrayIndex */ k = k_prime-1; /* Generally ArrayIndex points to an item smaller than the current IR */ ArrayIndex = ArrayNumUsed - 1 ; for( m=m_start+1 ; ArrayIndex >= 0 ; m++ ) { IR = InnerRadius(m*2+k, m) ; if( fabs(IR - s[ArrayIndex].IR) < Epsilon ) { if( ! InSeries(s[ArrayIndex].n, s[ArrayIndex].m, m*2+k, m ) ) { printf( "Match? k_prime=%lu: %lu/%lu and %lu/%lu, inner radii %0.17lf and %0.17lf (natural).\n", k_prime, s[ArrayIndex].n, s[ArrayIndex].m, m*2+k, m, s[ArrayIndex].IR, IR ); OutputMathematicaCode(fp_Mathematica, s[ArrayIndex].n, s[ArrayIndex].m, m*2+k, m); if( ++NumPossibleMatches >= MaxNumPossibleMatches ) goto Finished; } } /* closer than Epsilon */ if( ArrayIndex < ArrayNumUsed-1 ) { if( fabs(s[ArrayIndex+1].IR - IR) < Epsilon ) { if( ! InSeries(s[ArrayIndex+1].n, s[ArrayIndex+1].m, m*2+k, m ) ) { printf( "Match? k_prime=%lu: %lu/%lu and %lu/%lu, inner radii %0.17lf and %0.17lf (reverse).\n", k_prime, s[ArrayIndex+1].n, s[ArrayIndex+1].m, m*2+k, m, s[ArrayIndex+1].IR, IR ); OutputMathematicaCode(fp_Mathematica, s[ArrayIndex+1].n, s[ArrayIndex+1].m, m*2+k, m); if( ++NumPossibleMatches >= MaxNumPossibleMatches ) goto Finished; } } /* closer than Epsilon */ } /* not last element of array */ for( ; ArrayIndex >= 0 ; ) if( s[ArrayIndex].IR > IR ) ArrayIndex-- ; else break; } /* loop m, test ArrayIndex */ } /* k_prime */ Finished: ; ArrayFull: ; free(s); return(EXIT_SUCCESS); }