From e7eb602118e3dec6ad84c9e8fcf1a64a8d399a3f Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Fri, 19 May 2000 12:55:08 +0000 Subject: [PATCH] - Rearranged break-even points to better match present-day CPUs whenever GMP3 is used. --- src/base/digitseq/cl_DS_div.cc | 98 ++++++++++---- src/base/digitseq/cl_DS_mul.cc | 234 ++++++++++++++++++++------------- 2 files changed, 214 insertions(+), 118 deletions(-) diff --git a/src/base/digitseq/cl_DS_div.cc b/src/base/digitseq/cl_DS_div.cc index de88316..d9e8d5d 100644 --- a/src/base/digitseq/cl_DS_div.cc +++ b/src/base/digitseq/cl_DS_div.cc @@ -12,33 +12,75 @@ #include "cl_N.h" #include "cl_abort.h" -// We observe the following timings: -// Time for dividing 2*N digits by N digits, on a i486 33 MHz running Linux: -// N standard Newton -// 10 0.0003 0.0012 -// 25 0.0013 0.0044 -// 50 0.0047 0.0125 -// 100 0.017 0.037 -// 250 0.108 0.146 -// 500 0.43 0.44 -// 1000 1.72 1.32 -// 2500 11.2 4.1 -// 5000 44.3 9.5 -// 10000 187 20.6 -// -----> Newton faster for N >= 550. -// Time for dividing 3*N digits by N digits, on a i486 33 MHz running Linux: -// N standard Newton -// 10 0.0006 0.0025 -// 25 0.0026 0.0103 -// 50 0.0092 0.030 -// 100 0.035 0.089 -// 250 0.215 0.362 -// 500 0.85 1.10 -// 1000 3.44 3.21 -// 2500 23.3 7.9 -// 5000 89.0 15.6 -// 10000 362 33.1 -// -----> Newton faster for N >= 740. +// We observe the following timings in seconds: +// Time for dividing a 2*n word number by a n word number: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n standard Newton standard Newton +// 10 0.000010 0.000024 0.000036 0.000058 +// 30 0.000026 0.000080 0.00012 0.00027 +// 100 0.00018 0.00048 0.00084 0.0016 +// 300 0.0013 0.0028 0.0062 0.0090 +// 1000 0.014 0.019 0.064 0.066 <-(~2200) +// 2000 0.058 0.058 <-(~2000) 0.26 0.20 +// 3000 0.20 0.11 0.57 0.24 +// 10000 2.3 0.50 6.7 1.2 +// 30000 24.4 1.2 62.0 2.8 +// Time for dividing a 3*n word number by a n word number: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n standard Newton standard Newton +// 10 0.000013 0.000040 0.000063 0.00011 +// 30 0.000046 0.00018 0.00024 0.00062 +// 100 0.00035 0.0012 0.0016 0.0040 +// 300 0.0027 0.0071 0.012 0.021 +// 1000 0.029 0.047 0.13 0.16 +// 2000 0.12 0.14 <-(~2200) 0.51 0.45 <-(~1600) +// 3000 0.40 0.22 1.1 0.52 +// 10000 4.5 0.76 13.2 2.0 +// 30000 42.0 2.8 123.0 6.0 +// Time for dividing m digits by n digits: +// OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, +// Machine: P-III/450MHz Machine: EV5/300MHz: +// n Newton faster for: Newton faster for: +// 2-400 never never +// 600 never 670 n + { if (n < 900) + return cl_false; + else + if (n < 2200) + return (cl_boolean)((m >= n+50) && (m < 2*n-600)); + else + return (cl_boolean)(m >= n+30); + } +#else +// Use the old default values from CLN version <= 1.0.3 as a crude estimate. +// They came from the timings for dividing m digits by n digits on an i486/33: +// Dividing 2*N digits by N digits: Dividing 3*N digits by N digits: +// N standard Newton standard Newton +// 10 0.0003 0.0012 0.0006 0.0025 +// 25 0.0013 0.0044 0.0026 0.0103 +// 50 0.0047 0.0125 0.0092 0.030 +// 100 0.017 0.037 0.035 0.089 +// 250 0.108 0.146 0.215 0.362 +// 500 0.43 0.44 <-(~550) 0.85 1.10 +// 1000 1.72 1.32 3.44 3.21 <-(~740) +// 2500 11.2 4.1 23.3 7.9 +// 5000 44.3 9.5 89.0 15.6 +// 10000 187 20.6 362 33.1 // Time for dividing m digits by n digits: // n = 2,3,5,10,25,50,100,250: Newton never faster. // n = 400: Newton faster for m >= 440, m < 600 @@ -52,7 +94,6 @@ // n = 2000: Newton faster for m >= 2020 // n = 2500: Newton faster for m >= 2520 // n = 5000: Newton faster for m >= 5020 -// Break-even-point. When in doubt, prefer to choose the standard algorithm. static inline cl_boolean cl_recip_suitable (uintL m, uintL n) // m > n { if (n < 500) return cl_false; @@ -62,6 +103,7 @@ else return (cl_boolean)(m >= n+20); } +#endif // Dividiert zwei Unsigned Digit sequences durcheinander. // UDS_divide(a_MSDptr,a_len,a_LSDptr, b_MSDptr,b_len,b_LSDptr, &q,&r); diff --git a/src/base/digitseq/cl_DS_mul.cc b/src/base/digitseq/cl_DS_mul.cc index 285c642..be332ef 100644 --- a/src/base/digitseq/cl_DS_mul.cc +++ b/src/base/digitseq/cl_DS_mul.cc @@ -119,55 +119,55 @@ } while (len > 0); } -// Karatsuba-Multiplikation: O(n^(log 3 / log 2)) +// Karatsuba-multiplication: O(n^(log 3 / log 2)) static void mulu_karatsuba (const uintD* sourceptr1, uintC len1, const uintD* sourceptr2, uintC len2, uintD* destptr); static void mulu_karatsuba_square (const uintD* sourceptr, uintC len, uintD* destptr); #include "cl_DS_mul_kara.h" - // karatsuba_threshold = Länge, ab der die Karatsuba-Multiplikation bevorzugt - // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen. - // Als Test dient (progn (time (! 5000)) nil), das viele kleine und einige - // ganz große Multiplikationen durchführt. Miß die Runtime. - // Unter Linux mit einem 80486: Auf einer Sparc 2: - // threshold time in 0.01 sec. - // 5 125 127 - // 6 116 117 - // 7 107 110 - // 8 101 103 - // 9 99 102 - // 10 98 100 - // 11 97 100 - // 12 96 99 - // 13 97 99 - // 14 97 100 - // 15 97 99 - // 16 98 100 - // 17 98 100 - // 18 98 100 - // 19 98 101 - // 20 99 102 - // 25 103 105 - // 30 109 111 - // 40 115 118 - // 50 122 125 - // 70 132 134 - // 100 151 152 - // 150 164 167 - // 250 183 187 - // 500 203 205 - // 1000 203 205 - // (clisp)(cln) - // Das Optimum scheint bei karatsuba_threshold = 12 zu liegen. - // Da das Optimum aber vom Verhältnis - // Zeit für uintD-Multiplikation / Zeit für uintD-Addition - // abhängt und die gemessenen Zeiten auf eine Unterschreitung des Optimums - // empfindlicher reagieren als auf eine Überschreitung des Optimums, - // sind wir vorsichtig und wählen einen Wert etwas über dem Optimum: + // karatsuba_threshold = length, from which on Karatsuba-multiplication is a + // gain and will be preferred. The break-even point is determined from + // timings. The test is (progn (time (! 5000)) nil), which does many small + // and some very large multiplications. The measured runtimes are: + // OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, + // Machine: P-III/450MHz Machine: EV5/300MHz: + // threshold time in 0.01 sec. time in 0.01 sec. + // 5 3.55 2.29 + // 10 2.01 1.71 + // 15 1.61 1.61 + // 20 1.51 1.60 <- + // 25 1.45 1.63 + // 30 1.39 1.66 + // 35 1.39 <- 1.67 + // 40 1.39 1.71 + // 45 1.40 1.75 + // 50 1.41 1.78 + // 55 1.41 1.79 + // 60 1.44 1.84 + // 65 1.44 1.85 + // 70 1.43 1.85 + // 75 1.45 1.89 + // 80 1.47 1.91 + // 90 1.51 1.96 + // 100 1.53 1.97 + // 150 1.62 2.13 + // 250 1.75 2.19 + // 500 1.87 2.17 + // 1000 1.87 2.18 + // 2000 1.88 2.17 + // The optimum appears to be between 20 and 40. But since that optimum + // depends on the ratio time(uintD-mul)/time(uintD-add) and the measured + // times are more sensitive to a shift towards lower thresholds we are + // careful and choose a value at the upper end: +#if CL_USE_GMP + const unsigned int cl_karatsuba_threshold = 35; +#else const unsigned int cl_karatsuba_threshold = 16; + // (In CLN version <= 1.0.3 cl_karatsuba_threshold was always 16) +#endif -#if 0 // Lohnt sich nicht +#if 0 // Doesn't seem to be worth the effort // FFT-Multiplikation nach Nussbaumer: O(n log n log log n) #include "cl_DS_mul_nuss.h" @@ -246,58 +246,112 @@ // FFT-Multiplikation in Z/pZ: O(n^1.29) #include "cl_DS_mul_fftm.h" - // fftm_threshold = Länge, ab der die FFT-Multiplikation mod m bevorzugt - // wird. Der Break-Even-Point bestimmt sich aus Zeitmessungen. - // Multiplikation zweier N-Wort-Zahlen unter - // Linux mit einem 80486: Solaris, Sparc 10/20: + // fftm_threshold = length, from which on FFT multiplication mod m is a gain + // and will be preferred. The break-even point is determined from timings. + // The times to multiply two N-limb numbers are: + // OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, + // Machine: P-III/450MHz Machine: EV5/300MHz: // N kara fftm (time in sec.) kara fftm - // 1000 0.36 0.54 0.08 0.10 - // 5000 4.66 2.48 1.01 0.51 - // 25000 61.1 13.22 13.23 2.73 - // 32500 91.0 27.5 20.0 5.8 - // 35000 102.1 27.5 21.5 5.6 - // 50000 183 27.6 40.7 5.6 - // Multiplikation zweier N-Wort-Zahlen unter - // Linux mit einem 80486: Solaris, Sparc 10/20: - // N kara fftm (time in sec.) kara fftm - // 1000 0.36 0.54 0.08 0.10 - // 1260 0.52 0.50 0.11 0.10 - // 1590 0.79 0.51 0.16 0.10 - // 2000 1.09 1.07 0.23 0.21 - // 2520 1.57 1.08 0.33 0.21 - // 3180 2.32 1.08 0.50 0.21 - // 4000 3.29 2.22 0.70 0.41 - // 5040 4.74 2.44 0.99 0.50 - // N1 N2 kara fftm (time in sec.) kara fftm - // 1250 1250 0.51 0.50 0.11 0.10 - // 1250 1580 0.70 0.50 0.15 0.10 - // 1250 2000 0.89 0.51 0.18 0.10 - // 1250 2250 0.99 0.51 0.21 0.10 - // 1250 2500 1.08 1.03 <--- 0.22 0.21 - // 1250 2800 1.20 1.07 0.26 0.21 - // 1250 3100 1.35 1.07 0.28 0.21 - // Es gibt also noch Werte von (len1,len2) mit 1250 <= len1 <= len2, bei - // denen "kara" schneller ist als "fftm", aber nicht viele und dort auch - // nur um 5%. Darum wählen wir ab hier die FFT-Multiplikation. - const unsigned int cl_fftm_threshold = 1250; // muß stets >= 6 sein (sonst Endlosrekursion!) + // 1000 0.005 0.016 0.018 0.028 + // 1500 0.009 0.012 0.032 0.028 + // 2000 0.015 0.025 0.053 0.052 <- + // 2500 0.022 0.026 0.067 0.052 + // 3000 0.029 0.027 <- 0.093 0.053 + // 3500 0.035 0.037 0.12 0.031 + // 4000 0.045 0.028 0.16 0.12 + // 5000 0.064 0.050 0.20 0.11 + // 7000 0.110 0.051 0.37 0.20 + // 10000 0.19 0.11 0.61 0.26 + // 20000 0.59 0.23 1.85 0.55 + // 30000 1.10 0.25 3.79 0.56 + // 50000 2.52 1.76 8.15 1.37 + // 70000 4.41 2.30 14.09 2.94 + // 100000 7.55 1.53 24.48 2.96 + // More playing around with timings reveals that there are some values where + // FFT multiplication is somewhat slower than Karatsuba, both for len1==len2 + // and also if len1= 6 (else infinite recursion) +#else + // Use the old default value from CLN version <= 1.0.3 as a crude estimate. + const unsigned int cl_fftm_threshold = 1250; +#endif // This is the threshold for multiplication of equally sized factors. // When the lengths differ much, the threshold varies: - // len2 = 3000 len1 >= 800 - // len2 = 3500 len1 >= 700 - // len2 = 4000 len1 >= 580 - // len2 = 4500 len1 >= 430 - // len2 = 5000 len1 >= 370 - // len2 = 5500 len1 >= 320 - // len2 = 6000 len1 >= 500 - // len2 = 7000 len1 >= 370 - // len2 = 8000 len1 >= 330 - // len2 = 9000 len1 >= 420 - // len2 =10000 len1 >= 370 - // len2 =11000 len1 >= 330 - // len2 =12000 len1 >= 330 - // len2 =13000 len1 >= 350 + // OS: Linux 2.2, intDsize==32, OS: TRU64/4.0, intDsize==64, + // Machine: P-III/450MHz Machine: EV5/300MHz: + // len2 = 3000 len1 >= 2600 len1 >= 800 + // len2 = 4000 len1 >= 1500 len1 >= 700 + // len2 = 5000 len1 >= 1100 len1 >= 600 + // len2 = 6000 len1 >= 1300 len1 >= 700 + // len2 = 7000 len1 >= 1100 len1 >= 600 + // len2 = 8000 len1 >= 900 len1 >= 500 + // len2 = 9000 len1 >= 1300 len1 >= 600 + // len2 = 10000 len1 >= 1100 len1 >= 500 + // len2 = 11000 len1 >= 1000 len1 >= 500 + // len2 = 12000 len1 >= 900 len1 >= 700 + // len2 = 13000 len1 >= 900 len1 >= 500 + // len2 = 14000 len1 >= 900 len1 >= 600 + // Here are the timigs from CLN version <= 1.0.3: + // // len2 = 3000 len1 >= 800 + // // len2 = 3500 len1 >= 700 + // // len2 = 4000 len1 >= 580 + // // len2 = 4500 len1 >= 430 + // // len2 = 5000 len1 >= 370 + // // len2 = 5500 len1 >= 320 + // // len2 = 6000 len1 >= 500 + // // len2 = 7000 len1 >= 370 + // // len2 = 8000 len1 >= 330 + // // len2 = 9000 len1 >= 420 + // // len2 =10000 len1 >= 370 + // // len2 =11000 len1 >= 330 + // // len2 =12000 len1 >= 330 + // // len2 =13000 len1 >= 350 // Let's choose the following condition: +#if CL_USE_GMP + const unsigned int cl_fftm_threshold1 = 600; +#else + // Use the old default values from CLN version <= 1.0.3 as a crude estimate. const unsigned int cl_fftm_threshold1 = 330; +#endif const unsigned int cl_fftm_threshold2 = 2*cl_fftm_threshold; // len1 > cl_fftm_threshold1 && len2 > cl_fftm_threshold2 // && len1 >= cl_fftm_threshold1 + cl_fftm_threshold/(len2-cl_fftm_threshold1)*(cl_fftm_threshold-cl_fftm_threshold1). @@ -314,8 +368,8 @@ } return cl_false; } - -#if 0 // Lohnt sich nicht + +#if 0 // Doesn't seem to be worth the effort // FFT-Multiplikation über den komplexen Zahlen. #include "cl_DS_mul_fftc.h"