]> www.ginac.de Git - cln.git/commitdiff
* src/integer/conv/cl_I_from_digits.cc: Made input of all numbers in
authorRichard Kreckel <kreckel@ginac.de>
Wed, 2 Nov 2005 23:13:39 +0000 (23:13 +0000)
committerRichard Kreckel <kreckel@ginac.de>
Wed, 2 Nov 2005 23:13:39 +0000 (23:13 +0000)
        non-power-of-two base much faster.
        * tests/test_I_io.cc: New file...
        * tests/Makefile.in, tests/test_I.cc: ...used here.

ChangeLog
src/integer/conv/cl_I_from_digits.cc
tests/Makefile.in
tests/test_I.cc
tests/test_I_io.cc [new file with mode: 0644]

index 662167322a40dd5da2e4f82fc1f1d1820b42be8c..18e8e19a31523c4d63f306360abb786b7544862b 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,4 +1,11 @@
-2004-10-22  Richard B. Kreckel  <kreckel@ginac.de>
+2005-11-02  Richard B. Kreckel  <kreckel@ginac.de>
+
+       * src/integer/conv/cl_I_from_digits.cc: Made input of all numbers in
+       non-power-of-two base much faster.
+       * tests/test_I_io.cc: New file...
+       * tests/Makefile.in, tests/test_I.cc: ...used here.
+
+2005-10-22  Richard B. Kreckel  <kreckel@ginac.de>
 
        * Version 1.1.10 released.
 
index 1f49bbdf65a36e6b1edea0abafad3f0c5b6cb38d..7e101ac39175e4bd9013955a5a354138675d6150 100644 (file)
 \r
 namespace cln {\r
 \r
+static const cl_I digits_to_I_base2 (const char * MSBptr, uintL len, uintD base)\r
+{\r
+       // base is a power of two: write the digits from least significant\r
+       // to most significant into the result NUDS. Result needs\r
+       // 1+ceiling(len*log(base)/(intDsize*log(2))) or some more digits\r
+       CL_ALLOCA_STACK;\r
+       var uintD* erg_MSDptr;\r
+       var uintC erg_len;\r
+       var uintD* erg_LSDptr;\r
+       var int b = (base==2 ? 1 : base==4 ? 2 : base==8 ? 3 : base==16 ? 4 : /*base==32*/ 5);\r
+       num_stack_alloc(1+(len*b)/intDsize,,erg_LSDptr=);\r
+       erg_MSDptr = erg_LSDptr; erg_len = 0;\r
+       var uintD d = 0;  // resulting digit\r
+       var int ch_where = 0;  // position of ch inside d\r
+       while (len > 0) {\r
+               var uintB ch = *(const uintB *)(MSBptr+len-1); // next character\r
+               if (ch!='.') { // skip decimal point\r
+                       // Compute value of ch ('0'-'9','A'-'Z','a'-'z'):\r
+                       ch = ch - '0';\r
+                       if (ch > '9'-'0') { // not a digit?\r
+                               ch = ch+'0'-'A'+10;\r
+                               if (ch > 'Z'-'A'+10) {// not an uppercase letter?\r
+                                       ch = ch+'A'-'a'; // must be lowercase!\r
+                               }\r
+                       }\r
+                       d = d | (uintD)ch<<ch_where;\r
+                       ch_where = ch_where+b;\r
+                       if (ch_where >= intDsize) {\r
+                           // d is ready to be written into the NUDS:\r
+                           lsprefnext(erg_MSDptr) = d;\r
+                           ch_where = ch_where-intDsize;\r
+                           d = (uintD)ch >> b-ch_where;  // carry\r
+                           erg_len++;\r
+                       }\r
+               }\r
+               len--;\r
+       }\r
+       if (d != 0) {  // is there anything left over?\r
+               lsprefnext(erg_MSDptr) = d;\r
+               ++erg_len;\r
+       }\r
+       return NUDS_to_I(erg_MSDptr,erg_len);\r
+}\r
+\r
+// For each base b in [2..36], power_table[b-2] contains the largest exponent e\r
+// such that b^e<2^intDsize, i.e. floor(log(2^intDsize-1,b)).\r
+static const uintC power_table [36-2+1] = {\r
+#if (intDsize==8)\r
+       /* base  2..7  */           7,  5,  3,  3,  3,  2,\r
+       /* base  8..15 */   2,  2,  2,  2,  2,  2,  2,  2,\r
+       /* base 16..23 */   1,  1,  1,  1,  1,  1,  1,  1,\r
+       /* base 24..31 */   1,  1,  1,  1,  1,  1,  1,  1,\r
+       /* base 32..36 */   1,  1,  1,  1,  1\r
+#endif\r
+#if (intDsize==16)\r
+       /* base  2..7  */          15, 10,  7,  6,  6,  5,\r
+       /* base  8..15 */   5,  5,  4,  4,  4,  4,  4,  4,\r
+       /* base 16..23 */   3,  3,  3,  3,  3,  3,  3,  3,\r
+       /* base 24..31 */   3,  3,  3,  3,  3,  3,  3,  3,\r
+       /* base 32..36 */   3,  3,  3,  3,  3\r
+#endif\r
+#if (intDsize==32)\r
+       /* base  2..7  */          31, 20, 15, 13, 12, 11,\r
+       /* base  8..15 */  10, 10,  9,  9,  8,  8,  8,  8,\r
+       /* base 16..23 */   7,  7,  7,  7,  7,  7,  7,  7,\r
+       /* base 24..31 */   6,  6,  6,  6,  6,  6,  6,  6,\r
+       /* base 32..36 */   6,  6,  6,  6,  6\r
+#endif\r
+#if (intDsize==64)\r
+       /* base  2..7  */          63, 40, 31, 27, 24, 22,\r
+       /* base  8..15 */  21, 20, 19, 18, 17, 17, 16, 16,\r
+       /* base 16..23 */  15, 15, 15, 15, 14, 14, 14, 14,\r
+       /* base 24..31 */  13, 13, 13, 13, 13, 13, 13, 12,\r
+       /* base 32..36 */  12, 12, 12, 12, 12\r
+#endif\r
+};\r
+\r
+static const cl_I digits_to_I_baseN (const char * MSBptr, uintL len, uintD base)\r
+{\r
+       // base is not a power of two: Add digits one by one. Result nees\r
+       // 1+ceiling(len*log(base)/(intDsize*log(2))) or some more digits.\r
+       CL_ALLOCA_STACK;\r
+       var uintD* erg_MSDptr;\r
+       var uintC erg_len;\r
+       var uintD* erg_LSDptr;\r
+       var uintL need = 1+floor(len,intDsize*256); // > len/(intDsize*256) >=0\r
+       switch (base) { // multiply need with ceiling(256*log(base)/log(2)):\r
+               case 2: need = 256*need; break;\r
+               case 3: need = 406*need; break;\r
+               case 4: need = 512*need; break;\r
+               case 5: need = 595*need; break;\r
+               case 6: need = 662*need; break;\r
+               case 7: need = 719*need; break;\r
+               case 8: need = 768*need; break;\r
+               case 9: need = 812*need; break;\r
+               case 10: need = 851*need; break;\r
+               case 11: need = 886*need; break;\r
+               case 12: need = 918*need; break;\r
+               case 13: need = 948*need; break;\r
+               case 14: need = 975*need; break;\r
+               case 15: need = 1001*need; break;\r
+               case 16: need = 1024*need; break;\r
+               case 17: need = 1047*need; break;\r
+               case 18: need = 1068*need; break;\r
+               case 19: need = 1088*need; break;\r
+               case 20: need = 1107*need; break;\r
+               case 21: need = 1125*need; break;\r
+               case 22: need = 1142*need; break;\r
+               case 23: need = 1159*need; break;\r
+               case 24: need = 1174*need; break;\r
+               case 25: need = 1189*need; break;\r
+               case 26: need = 1204*need; break;\r
+               case 27: need = 1218*need; break;\r
+               case 28: need = 1231*need; break;\r
+               case 29: need = 1244*need; break;\r
+               case 30: need = 1257*need; break;\r
+               case 31: need = 1269*need; break;\r
+               case 32: need = 1280*need; break;\r
+               case 33: need = 1292*need; break;\r
+               case 34: need = 1303*need; break;\r
+               case 35: need = 1314*need; break;\r
+               case 36: need = 1324*need; break;\r
+               default: NOTREACHED\r
+       }\r
+       // Now we have need >= len*log(base)/(intDsize*log(2)).\r
+       need += 1;\r
+       // Add digits one by one:\r
+       num_stack_alloc(need,,erg_LSDptr=);\r
+       // erg_MSDptr/erg_len/erg_LSDptr is a NUDS, erg_len < need.\r
+       erg_MSDptr = erg_LSDptr; erg_len = 0;\r
+       while (len > 0) {\r
+               var uintD newdigit = 0;\r
+               var uintC chx = 0;\r
+               var uintD factor = 1;\r
+               while (chx < power_table[base-2] && len > 0) {\r
+                       var uintB ch = *(const uintB *)MSBptr; MSBptr++; // next character\r
+                       if (ch!='.') { // skip decimal point\r
+                               // Compute value of ('0'-'9','A'-'Z','a'-'z'):\r
+                               ch = ch-'0';\r
+                               if (ch > '9'-'0') { // not a digit?\r
+                                       ch = ch+'0'-'A'+10;\r
+                                       if (ch > 'Z'-'A'+10) {// not an uppercase letter?\r
+                                               ch = ch+'A'-'a';  // must be lowercase!\r
+                                       }\r
+                               }\r
+                               factor = factor*base;\r
+                               newdigit = base*newdigit+ch;\r
+                               chx++;\r
+                       }\r
+                       len--;\r
+               }\r
+               var uintD carry = mulusmall_loop_lsp(factor,erg_LSDptr,erg_len,newdigit);\r
+               if (carry!=0) {\r
+                       // need to extend NUDS:\r
+                       lsprefnext(erg_MSDptr) = carry;\r
+                       erg_len++;\r
+               }\r
+       }\r
+       return NUDS_to_I(erg_MSDptr,erg_len);\r
+}\r
+\r
 const cl_I digits_to_I (const char * MSBptr, uintL len, uintD base)\r
 {\r
-      CL_ALLOCA_STACK;\r
-      var uintD* erg_MSDptr;\r
-      var uintC erg_len;\r
-      var uintD* erg_LSDptr;\r
-      if ((base & (base-1)) == 0) {\r
-        // Fast path for powers of two: write the digits from least\r
-        // significant to most significant into the result NUDS.\r
-        var int b = (base==2 ? 1 : base==4 ? 2 : base==8 ? 3 : base==16 ? 4 : /*base==32*/ 5);\r
-        num_stack_alloc(1+(len*b)/intDsize,,erg_LSDptr=);\r
-        erg_MSDptr = erg_LSDptr; erg_len = 0;\r
-        var uintD d = 0;  // resulting digit\r
-        var int ch_where = 0;  // position of ch inside d\r
-        while (len > 0)\r
-          { var uintB ch = *(const uintB *)(MSBptr+len-1); // next character\r
-            if (!(ch=='.')) // skip decimal point\r
-              { // Compute value of ch ('0'-'9','A'-'Z','a'-'z'):\r
-                ch = ch - '0';\r
-                if (ch > '9'-'0') // not a digit?\r
-                  { ch = ch+'0'-'A'+10;\r
-                    if (ch > 'Z'-'A'+10) // not an uppercase letter?\r
-                      { ch = ch+'A'-'a'; } // must be lowercase!\r
-                  }\r
-                d = d | (uintD)ch<<ch_where;\r
-                ch_where = ch_where+b;\r
-                if (ch_where >= intDsize) {\r
-                  // d is ready to be written into the NUDS:\r
-                  lsprefnext(erg_MSDptr) = d;\r
-                  ch_where = ch_where-intDsize;\r
-                  d = (uintD)ch >> b-ch_where;  // carry\r
-                  erg_len++;\r
-                }\r
-              }\r
-            len--;\r
-          }\r
-        if (d != 0) {  // is there anything left over?\r
-          lsprefnext(erg_MSDptr) = d;\r
-          ++erg_len;\r
-        }\r
-      } else {\r
-        // Slow path: Add digits one by one for non-powers of two.\r
-        // Platz fürs Ergebnis:\r
-        // 1+ceiling(len*log(base)/(intDsize*log(2))) or some more digits\r
-        var uintL need = 1+floor(len,intDsize*256); // > len/(intDsize*256) >=0\r
-        switch (base) // multiply need with ceiling(256*log(base)/log(2)):\r
-          {\r
-          //case 2: need = 256*need; break;\r
-            case 3: need = 406*need; break;\r
-          //case 4: need = 512*need; break;\r
-            case 5: need = 595*need; break;\r
-            case 6: need = 662*need; break;\r
-            case 7: need = 719*need; break;\r
-          //case 8: need = 768*need; break;\r
-            case 9: need = 812*need; break;\r
-            case 10: need = 851*need; break;\r
-            case 11: need = 886*need; break;\r
-            case 12: need = 918*need; break;\r
-            case 13: need = 948*need; break;\r
-            case 14: need = 975*need; break;\r
-            case 15: need = 1001*need; break;\r
-          //case 16: need = 1024*need; break;\r
-            case 17: need = 1047*need; break;\r
-            case 18: need = 1068*need; break;\r
-            case 19: need = 1088*need; break;\r
-            case 20: need = 1107*need; break;\r
-            case 21: need = 1125*need; break;\r
-            case 22: need = 1142*need; break;\r
-            case 23: need = 1159*need; break;\r
-            case 24: need = 1174*need; break;\r
-            case 25: need = 1189*need; break;\r
-            case 26: need = 1204*need; break;\r
-            case 27: need = 1218*need; break;\r
-            case 28: need = 1231*need; break;\r
-            case 29: need = 1244*need; break;\r
-            case 30: need = 1257*need; break;\r
-            case 31: need = 1269*need; break;\r
-          //case 32: need = 1280*need; break;\r
-            case 33: need = 1292*need; break;\r
-            case 34: need = 1303*need; break;\r
-            case 35: need = 1314*need; break;\r
-            case 36: need = 1324*need; break;\r
-            default: NOTREACHED\r
-          }\r
-        // Now we have need >= len*log(base)/(intDsize*log(2)).\r
-        need += 1;\r
-        // Add digits one by one:\r
-        num_stack_alloc(need,,erg_LSDptr=);\r
-        erg_MSDptr = erg_LSDptr; erg_len = 0;\r
-        while (len > 0)\r
-          { // erg_MSDptr/erg_len/erg_LSDptr is a NUDS, erg_len < need.\r
-            var uintB ch = *(const uintB *)MSBptr; MSBptr++; // next character\r
-            if (!(ch=='.')) // skip decimal point\r
-              { // Compute value of ('0'-'9','A'-'Z','a'-'z'):\r
-                ch = ch - '0';\r
-                if (ch > '9'-'0') // not a digit?\r
-                  { ch = ch+'0'-'A'+10;\r
-                    if (ch > 'Z'-'A'+10) // not an uppercase letter?\r
-                      { ch = ch+'A'-'a'; } // must be lowercase!\r
-                  }\r
-                // multiply erg with base and add ch:\r
-               {var uintD carry = mulusmall_loop_lsp(base,erg_LSDptr,erg_len,ch);\r
-                if (!(carry==0))\r
-                  // need to extend NUDS:\r
-                  { lsprefnext(erg_MSDptr) = carry; erg_len++; }\r
-              }}\r
-            len--;\r
-          }\r
-      }\r
-      return NUDS_to_I(erg_MSDptr,erg_len);\r
+       if ((base & (base-1)) == 0) {\r
+               return digits_to_I_base2(MSBptr, len, base);\r
+       } else {\r
+               // This is quite insensitive to the breakeven point.\r
+               // On a 1GHz Athlon I get approximately:\r
+               //   base  3: breakeven == 15000\r
+               //   base 10: breakeven ==  5000\r
+               //   base 36: breakeven ==  2000\r
+               if (len>50000/base)\r
+                       // Divide-and-conquer:\r
+                       return digits_to_I(MSBptr,len/2,base)*expt_pos(base,len-len/2)\r
+                             +digits_to_I(MSBptr+len/2,len-len/2,base);\r
+               else\r
+                       return digits_to_I_baseN(MSBptr, len, base);\r
+       }\r
 }\r
 \r
 }  // namespace cln\r
index 18802f097aa8d4829d415c21351043d0355cee74..dc6bf808668693eac8cdb8819ac8845afbee5651 100644 (file)
@@ -56,7 +56,7 @@ MODULES_tests = tests \
                 test_I_gcd test_I_xgcd \
                 test_I_ash test_I_evenp test_I_oddp test_I_lognot test_I_logand test_I_logandc1 test_I_logandc2 test_I_logior test_I_logorc1 test_I_logorc2 test_I_logxor test_I_lognand test_I_lognor test_I_logeqv test_I_boole test_I_logbitp test_I_logtest test_I_ldb test_I_ldbtest test_I_mkf test_I_dpb test_I_dpf test_I_logcount test_I_ilength test_I_ord2 test_I_power2p \
                 test_I_isqrt test_I_sqrtp \
-                test_I_GV \
+                test_I_io test_I_GV \
                 test_MI \
                 test_MI_canonhom test_MI_plus test_MI_minus test_MI_mul test_MI_recip test_MI_div test_MI_expt \
                 test_nt \
index 6073620813d2204d8f2e4c68440db9514ae2456a..8e5ea63f494dd4604a71ebe70f4ef469a005b42c 100644 (file)
@@ -43,6 +43,7 @@ extern int test_I_power2p (int iterations);
 extern int test_I_isqrt (int iterations);
 extern int test_I_sqrtp (int iterations);
 // Miscellaneous.
+extern int test_I_io (int iterations);
 extern int test_I_GV (int iterations);
 
 #define RUN(tester,iterations)  \
@@ -96,6 +97,7 @@ int test_I (int iterations)
        RUN(test_I_isqrt,iterations);
        RUN(test_I_sqrtp,iterations);
        // Miscellaneous.
+       RUN(test_I_io,iterations);
        RUN(test_I_GV,iterations);
        return error;
 }
diff --git a/tests/test_I_io.cc b/tests/test_I_io.cc
new file mode 100644 (file)
index 0000000..cbc4b1e
--- /dev/null
@@ -0,0 +1,17 @@
+#include "test_I.h"
+#include <cln/input.h>
+#include <sstream>
+
+int test_I_io (int iterations)
+{
+       int error = 0;
+       for (int i = iterations; i > 0; i--) {
+               cl_I a = testrandom_I();
+               int base = iterations % (36-1) + 2;
+               cl_read_flags rflags = {syntax_integer, lsyntax_standard, base};
+               stringstream buf;
+               print_integer(buf, base, a);
+               ASSERT1(a == read_integer(buf, rflags), a);
+       }
+       return error;
+}