Xenharmonic Wiki

User:VIxen/Code: Difference between revisions

VisualWikitext
VIxen (talk | contribs)
autopatrolled, emailconfirmed
155 edits
DR triad finder code
 
VIxen (talk | contribs)
autopatrolled, emailconfirmed
155 edits
m no sense in if...break when I can add the condition to the for-clause, smh
 
(One intermediate revision by the same user not shown)
Line 71: Line 71:
int printIsodiffTriads(){
int printIsodiffTriads(){
/* Find triads in ETs that satisfy the frequency relation bottom + top == 2 * middle.
/* Find triads in ETs that satisfy the frequency relation bottom + top == 2 * middle.
The root of x^{top} - 2 * x^{mid} + 1 is found by the Newton-Raphson method
The root of x^{top} - 2 * x^{mid} + 1 is found by the Newton-Raphson method. */
after dividing the polynomial by x-1. */


     const int MINDEG = 12;
     const int MINDEG = 12;
Line 93: Line 92:
           to omit boring 0-(2k+1)-4k triads for large k */  
           to omit boring 0-(2k+1)-4k triads for large k */  


         for (mid = top - (top/6) - 3; mid >= minmid; mid--){
         for (mid = top - (top/6) - 3; (mid >= minmid) && (x > MINX); mid--){
         /* mid starts a bit lower than top to omit boring triads like 0-4k-4(k+1) */
         /* mid starts a bit lower than top to omit boring triads like 0-4k-4(k+1) */


             do{ /* Step of Newton-Raphson. TODO: avoid summing lots of powers,
             do{ /* Step of Newton-Raphson */
                  differentiate (x^k-1)/(x-1) like a fraction instead. */
                 pwr = pow(x, mid);
                 pwr = x;
                 val = 1. - 2. * pwr;
                 val = -x - 1.;
                 der = -2. * mid * pwr / x;
                 der = -1.;
                pwr *= pow(x, top - mid);
                for (p = 2; p < mid; p++){
                val += pwr;
                    pwr *= x;
                der += top * pwr / x;
                    val -= pwr;
                    der -= pwr * p;
                }
                for (; p < top; p++){
                    pwr *= x;
                    val += pwr;
                    der += pwr * p;
                }
                 x -= val / der;
                 x -= val / der;
             } while (val > 1e-5);
             } while (val > 1e-5);
Line 116: Line 107:
             val = log2(x);
             val = log2(x);
             fprintf(res, "%d,%d,%f,%.3f,%.3f\n", mid, top, 1200.*val, 1./val, val * top);
             fprintf(res, "%d,%d,%f,%.3f,%.3f\n", mid, top, 1200.*val, 1./val, val * top);
            if (x < MINX) break;
         }
         }
     }
     }
Retrieved from "https://en.xen.wiki/w/User:VIxen/Code"