// Fac(p) and Binom(p,q) are just made for comparison of the results with the method in poly.cpp, // p! = p.(p-1).(p-2) .. .3.2.1 #macro Fac(p) ( #if (p<2) 1 #local Result = 1; #else #local Result = p; p #local Count=p-1; #while (Count>1) #local Result = Result * Count; *Count #local Count=Count-1; #end #end #debug "\n Macro Fac(p):" #declare StrD = concat("\np = ",str(p,0,0)," Fac(p) = ",str(Result,0,0) ) #debug StrD ) #end //#macro Fac(p) //---------------------------------------------------------------------------------------------------------------------------------- /* ( p ) p! ( ) = -------- ( q ) q!(p-q)! */ #macro Binom(p,q) ( Fac(p) /Fac(q) /Fac(p-q) ) #end //#macro Binom(p,q) //---------------------------------------------------------------------------------------------------------------------------------- /* **************** Translation of 'binomial(int n, int r)' in poly.cpp to POV-Ray SDL *************************************** */ #declare Stack1=array[400]; // Contains the factors (r+1) -> n. (The remaining factors after n!/r!) #declare Stack2=array[100]; // Contains the prime factors of the numbers 2 -> (n-r). // Principle: The factors in Stack1 are divided by the prime factors of // 2 to (n-r) that are repeatedly put in Stack2. // The subroutines 'factor1' and 'factor_out' produce the prime factors // in Stack2. // Stack1 is then searched through for a divisible factor and divided. // This is repeated until all primes in Stack2 are done. // This is repeated from 2 to (n-r). // The remaining factors in Stack1 are multiplied together. /* ******************************************************************************************************************************** static int factor_out(int n, int i, int *c, int *s) { while (!(n % i)) // (n % i) is remainder of integer division, if divisible then not(zero) = one = true. { n /= i; // n = n/i ,factor found and divide by factor. s[(*c)++] = i; // factors (i) are saved in array *s[], with *c = 'count' = number of factors. } return(n); } ********************************************************************************************************************************** */ #macro Factor_Out(n, i, c, Stack2) ( #debug "\n Macro Factor_Out(n, i, c, Stack2):" #declare StrC = concat("\n n = ",str(n,0,0)," i = ",str(i,0,0)," c = ",str(c,0,0) ) #debug StrC #while (mod(n,i)=0) #declare n = div(n,i); #declare Stack2[c] = i; #declare StrC = concat("\n n = ",str(n,0,0)," Stack2[",str(c,0,0), "] = ",str(Stack2[c],0,0) ); #debug StrC #declare c = c+1; #end //while (mod(n,i)=0) n #declare StrC = concat("\n c = ",str(c,0,0) ); #debug StrC #debug "\n End Macro Factor_Out." ) #end //macro Factor_Out(n, i, c, Stack2) /* ******************************************************************************************************************************** static void factor1(int n, int *c, int *s) { int i,k; n = factor_out(n, 2, c, s); // Remove factors 2 from n, save in *s[], *c = 'count' = number of factors. k = (int)sqrt((DBL)n) + 1; // k = (Biggest is sqrt(n)+1 ). for (i = 3; (n > 1) && (i <= k); i += 2) // Now any odd factors. { if (!(n % i)) // (n % i) is remainder of integer division, if divisible then not(zero) = one = true. { n = factor_out(n, i, c, s); // Remove factors i from n, save in *s[], *c ='count' = number of factors. k = (int)sqrt((DBL)n)+1; // Make new k. } } if (n > 1) { s[(*c)++] = n; } // Add last prime to *s[] if bigger then 1. } ********************************************************************************************************************************** */ #macro Factor1(n, c, Stack2) #debug "\n Macro Factor1(n, c, Stack2):" #declare StrB = concat("\n n = ",str(n,0,0)," c = ",str(c,0,0) ) #debug StrB #declare n = Factor_Out(n, 2, c, Stack2); #local k = int(sqrt(n)+1); #local i = 3; #declare StrB = concat("\n n = ",str(n,0,0)," sqrt(n) = ",str(k,0,0), " i = ",str(i,0,0) ); #debug StrB #while ((n>1)&(i<=k)) #if (mod(n,i)=0) #declare n = Factor_Out(n, i, c, Stack2); #local k = int(sqrt(n)+1); #declare StrB = concat("\n n = ",str(n,0,0)," sqrt(n) = ",str(k,0,0) ); #debug StrB #end //if (mod(n,i)=0) #local i = i+2; #end //while ((n>1)&(i<=k)) #if (n>1) #declare Stack2[c] = n; #declare StrB = concat("\n n = ",str(n,0,0), " Stack2[",str(c,0,0),"] =",str(Stack2[c],0,0) ); #debug StrB #declare c = c+1; #declare StrB = concat("\n c = ",str(c,0,0) ); #debug StrB #end //if (n>1) #debug "\n End Macro Factor1." #end //macro Factor1(n, c, Stack2) /* ******************************************************************************************************************************** static long binomial(int n, int r) { int h,i,j,k,l; unsigned long result; static int stack1[BINOMSIZE], stack2[BINOMSIZE]; // BINOMSIZE = 40 if ((n < 0) || (r < 0) || (r > n)) { result = 0L; } // Prevent rare mistakes. else { if (r == n) { result = 1L; } // Prevent unnecessary calculations. else { if ((r < 15) && (n < 15)) // Use precalculated values. { result = (long)binomials[n][r]; } else // For values that are not precomputed. { j = 0; // 'j' gives last index in Stack1. for (i = r + 1; i <= n; i++) // Fill Stack1 with: r+1, r+2, .. , n-1, n. { stack1[j++] = i; } // r! can be scratched away against n!. Rest = (r+1).(r+2)...(n-1).n. for (i = 2; i <= (n-r); i++) // i = 2, 3, 4, .., (n-r). { h = 0; // 'h' gives last index in Stack2. factor1(i, &h, stack2); // Find all prime factors of i and put them in Stack2. for (k = 0; k < h; k++) // Walk through Stack2 until h. { for (l = 0; l < j; l++) // Walk through Stack1 until j. { if (!(stack1[l] % stack2[k])) // If Stack1[l] is divisible by Stack2[k] then divide Stack1[l]. { stack1[l] /= stack2[k]; // Divide the number in Stack1 by factor in Stack2. goto l1; // Jump over error message. } } } // Error if we get here // Debug_Info("Failed to factor\n"); // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! l1:; // If we jump to this position only the first input of Stack2 is done // , more entrees in Stack2 are skipped!!!!!!!!!!!!!!!!!! // l1: must be placed one bracket back, so Stack2 is run completely!!! // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! } result = 1; // Multiply the remaining factors together. for (i = 0; i < j; i++) { result *= stack1[i]; } } } } return(result); } ********************************************************************************************************************************** */ #macro Binomial(n, r) #debug "\nMacro Binomial(n,r):" #declare StrA = concat("\nn = ",str(n,0,0)," r = ",str(r,0,0) ); #debug StrA #if ((n<0) | (r<0) | (r>n)) 0 #else #if (r = n) 1 #else // **************************************************************************************************************** #if (r < div(n,2)) // This is not in the original file but a great speed improvement !!! #local r = n-r; // Can be done because Binomials are symmetric. #end // **************************************************************************************************************** #declare StrA = concat("\nn = ",str(n,0,0)," r = ",str(r,0,0) ); #debug StrA #local j = 0; #local i = r+1; #debug "\n Fill Stack1 with (r+1) => n:" #while (i<=n) #declare Stack1[j] = i; #declare StrA = concat("\ni = ",str(i,0,0), " Stack1[",str(j,0,0),"] = ",str(Stack1[j],0,0) ); #debug StrA #local j = j+1; #local i = i+1; #end //while (i<=n) #local i = 2; #declare StrA = concat("\n\nWalking i=2 => i = (n-r) = ",str(n-r,0,0) ); #debug StrA #while (i <= (n-r)) #local h = 0; #debug "\n Factor i:" #declare StrA = concat(" i = ",str(i,0,0) ); #debug StrA #local ii = i; // Prefents reset of i by Factor1(). That's why I needed all that debugging !!! // This problem does not take place in the cpp-code I presume. Factor1(ii, h, Stack2) #declare StrA = concat("\n i = ",str(i,0,0)," h = ",str(h,0,0) ); #debug StrA #local k = 0; #debug "\n Start walking through Stack2: (n-r)!" #declare StrA = concat("\n k = ",str(k,0,0)," h = ",str(h,0,0) ); #debug StrA // ********************************************************************************************************************* // #local BBB = 1; // BBB is a substitute for the 'goto l1:' in the cpp-code. // #while ((k < h) & BBB) // According to the poly.cpp-file here BBB should be used, but that's wrong !!! // If put here you skip the rest of Stack2 when the first prime factor is done. // ********************************************************************************************************************* #while (k < h) // (k< h) and BBB in the next stage is correct. Then all prime factors of // Stack2 are done. #local l = 0; #local BBB = 1; // So here initiation of BBB (BBB = 1;) #declare StrA = concat("\n k = ",str(k,0,0)," l = ",str(l,0,0)," BBB = 1" ); #debug StrA #debug "\n Start walking through Stack1: n!/r!" #while ((l < j) & BBB ) // Here the escape. #declare StrA = concat("\n l = ",str(l,0,0), " Stack1[",str(l,0,0),"] = ",str(Stack1[l],0,0), " k = ",str(k,0,0)," Stack2[",str(k,0,0),"] = ",str(Stack2[k],0,0) ); #debug StrA #if (mod(Stack1[l],Stack2[k])=0) #declare Stack1[l] = div(Stack1[l],Stack2[k]); #declare StrA = concat("\n Stack1[l] = ",str(Stack1[l],0,0), " Stack2[k] = ",str(Stack2[k],0,0)," BBB = 0" ); #debug StrA #local BBB = 0; // Raise the escape flag. (= goto l1: in cpp) #end #local l = l+1; #if (l>=j) #debug "\n End walking through Stack1" #end #end //while ((l=h) #debug "\n End walking through Stack2" #end #end //while (k(n-r)) #debug "\nEnd walking i=2 => (n-r)" #end #end //while (i <= (n-r)) #local i = 0; 1 //---------------------------------------------------------------------------------------------------------------- #declare Result = 1; // This is for the debugging only! //---------------------------------------------------------------------------------------------------------------- #while (i < j) // Multiply the remaining factors together. *Stack1[i] //---------------------------------------------------------------------------------------------------------------- #declare Result = Result*Stack1[i]; // For debugging #declare StrA = concat("\ni = ",str(i,0,0), " Stack1[",str(i,0,0),"] = ",str(Stack1[i],0,0), " Result = ",str(Result,0,0) ); #debug StrA //---------------------------------------------------------------------------------------------------------------- #local i = i+1; #end #end //if (r=n) else #end //if ((n<0) | (r<0) | (r>n)) else #debug "\nEnd Macro Binomial(n,r).\n" #end //macro Binomial(n, r) #declare n = 18; //n #declare r = 13; //r #declare Str2 = concat( "\n Binomial(",str(n,0,0),",",str(r,0,0),") = ",str(Binomial(n,r),0,0),"\n\n" ); #declare Str1 = concat( "\n Binom(",str(n,0,0),",",str(r,0,0),") = ",str(Binom(n,r),0,0) ); #debug Str1 #debug Str2