From: antispam@fricas.org   
      
   albert@spenarnc.xs4all.nl wrote:   
   > In article ,   
   > wrote:   
   >>Using the standard program dc on linux to test for primality.   
   >>This hinges on | that does modular exponentiation.   
   >>   
   >>ciforth idiosyncrasies:   
   >> Uses &c to denote characters   
   >> Uses ^c to denote control characters   
   >> Uses &, in large numbers   
   >> Uses RDROP   
   >> Uses real strings "ORANG UTAN"   
   >> WANT announces a facility that you need   
   >> REGRESS introduces a test, use it as an example.   
   >> The latter two can be regarded as comment   
   >>   
   >>You cannot get around the prefix h that introduce huge numbers.   
   >>You could however replace   
   >> '' h10 '' by '' 10 >S ''   
   >>for single precision numbers.   
   >>   
   >>A working /dev/random is used, but you could replace RAND   
   >>by a suitable random number generation.   
   >>Even `` 117 CONSTANT RAND '' would work here.   
   >>   
   >>Communication to dc is established using the following interface   
   >>With this dummy file it is supposed to compile.   
   >>------------------ 8< --------------------------   
   >>   
   >>\ At this point we have available   
   >>\ command (sc -- ) Get a string plus cr to dc.   
   >>\ cmd1 (char -- ) Get a single char command to dc.   
   >>\ h PREFIX : number preceded by h is put on the dc stack.   
   >>\ get : get one string from dc, possibly containing &\..   
   >>\ Define dummies and kill the regression test.   
   >>: command 2DROP ;   
   >>: cmd1 DROP ;   
   >>: h NAME 2DROP ; IMMEDIATE PREFIX   
   >>: get 0 0 ;   
   >>: REGRESS POSTPONE \ ; IMMEDIATE   
   >>   
   >>------------------ 8< --------------------------   
   >>   
   >>This is the actual program.   
   >>h is the operator in dc   
   >>h is a huge number to dc.   
   >>   
   >>------------------ 8< --------------------------   
   >>\ Copyright (2026): Albert van der Horst {by GNU Public License}   
   >>   
   >>\ An interface with the program dc that can handle h-numbers (huge).   
   >>   
   >>WANT BOUNDS TRUE REGRESS   
   >>   
   >>INCLUDE pipe.frt   
   >>   
   >>: $>S command ;   
   >>: >S 0 <# #S #> command ;   
   >>: S>$ "ps_" command get ;   
   >>: S> S>$ EVALUATE ;   
   >>: S>d 0.0 S>$ >NUMBER 2DROP ;   
   >>   
   >>REGRESS 1234 >S S> S: 1234   
   >>REGRESS "1234" $>S S>$ EVALUATE S: 1234   
   >>   
   >>: +h &+ cmd1 ;   
   >>: -h &- cmd1 ;   
   >>: *h &* cmd1 ;   
   >>: /h &/ cmd1 ;   
   >>: MODh &% cmd1 ;   
   >>: /MODh &~ cmd1 &r cmd1 ;   
   >>: **h &^ cmd1 ;   
   >>: */h "s_ * l_ /" command ;   
   >>: |h &| cmd1 ;   
   >>   
   >>REGRESS h13 h3 +h S> S: 16   
   >>REGRESS h13 h3 -h S> S: 10   
   >>REGRESS h13 h3 *h S> S: 39   
   >>REGRESS h13 h3 /h S> S: 4   
   >>REGRESS h13 h3 MODh S> S: 1   
   >>REGRESS h13 h3 /MODh S> S> S: 4 1   
   >>REGRESS h13 h2 **h S> S: 169   
   >>REGRESS h13 h2 h70 |h S> S: 29   
   >>REGRESS h12 h2 h6 */h S> S: 4   
   >>   
   >>: DROPh "s_" command ;   
   >>: SWAPh &r cmd1 ;   
   >>: DUPh &d cmd1 ;   
   >>: OVERh "sAsBlBlAlB" command ;   
   >>REGRESS h1 h2 h3 h4 SWAPh S> S> S> S> S: 3 4 2 1   
   >>REGRESS h1 h2 h3 h4 DROPh S> S> S> S: 3 2 1   
   >>REGRESS h1 h2 h3 h4 DUPh S> S> S> S> S> S: 4 4 3 2 1   
   >>REGRESS h1 h2 h3 h4 OVERh S> S> S> S> S> S: 3 4 3 2 1   
   >>   
   >>: .h S>$ TYPE ;   
   >>REGRESS h1 h2 +h .h S:   
   >>   
   >>\ On dc the .S (f) command fails if the stack is empty.   
   >>\ So first DEPTH is added (z) and later removed (s_).   
   >>: .Sh "zfs_ %n" command get TYPE ;   
   >>   
   >>\ Testing huge primes.   
   >>"/dev/random" 0 OPEN-FILE THROW CONSTANT devrand   
   >>: RAND _ DSP@ 7 devrand READ-FILE THROW 7 <> 1001 ?ERROR ;   
   >>REGRESS RAND . S:   
   >>   
   >>   
   >>: inithex "16o" command "16i" command ;   
   >>   
   >>1,000,000,000 CONSTANT infinity   
   >>   
   >>\ For X leave odd X' and POWER of 2. (H:X -- H:X' ) ( -- p)   
   >>\ X = X' * 2^p   
   >>: even-norm infinity 0 DO   
   >> DUPh h2 MODh S> 1 = IF I LEAVE THEN h2 /h   
   >> LOOP ;   
   >>REGRESS h26 even-norm S> S: 1 13   
   >>   
   >>\ The number on the DC stack is one. (Z is how many digits).   
   >>: one? DUPh "Z" command S> 1 > IF DROPh FALSE ELSE S> 1 = THEN ;   
   >>REGRESS h1234 one? S: FALSE   
   >>REGRESS h0 one? S: FALSE   
   >>REGRESS h1 one? S: TRUE   
   >>   
   >>\ The number on the DC stack is zero.   
   >>: zero? DUPh "Z" command S> 1 > IF DROPh FALSE ELSE S> 0 = THEN ;   
   >>REGRESS h1234 zero? S: FALSE   
   >>REGRESS h0 zero? S: TRUE   
   >>REGRESS h1 zero? S: FALSE   
   >>   
   >>\ For dc containing a number: "it IS probably prime using one probe."   
   >>\ (H:n -- ) ( -- flag)   
   >>: ?Rabbin DUPh "sp" command \ local p   
   >> h1 -h even-norm >R   
   >> RAND >S SWAPh "lp" command |h   
   >> DUPh one? IF DROPh RDROP TRUE EXIT THEN   
   >>   
   >> BEGIN DUPh h1 +h "lp" command -h zero? IF DROPh RDROP TRUE EXIT THEN   
   >> R@ WHILE R> 1- >R   
   >> DUPh *h "lp" command MODh \ Square   
   >> REPEAT DROPh RDROP FALSE ;   
   >>REGRESS h1001 ?Rabbin S: FALSE   
   >>REGRESS h1013 ?Rabbin S: TRUE   
   >>   
   >>: test 1,000,000,000,000,000,100 1,000,000,000,000,000,001 DO   
   >> I 0 <# #S #> command ?Rabbin IF I . CR THEN 2 +LOOP ;   
   >>\ Expect 03 09 31 79   
   >>   
   >>: post-google   
   >> h10 h100 **h   
   >> BEGIN DUPh ?Rabbin 0= WHILE h1 +h REPEAT   
   >> .h ;   
   >>   
   >>------------------ 8< --------------------------   
   >>   
   >>Examples:   
   >>Finding primes inside a range:   
   >>   
   >>S[ ] OK test   
   >>1000000000000000003   
   >>1000000000000000009   
   >>1000000000000000031   
   >>1000000000000000079   
   >>   
   >>Find the first prime past google (10^100)   
   >>S[ ] OK post-google   
   >>100000000000000000000000000000000000000000000000000000000000000000000\   
   >>00000000000000000000000000000267   
   >>   
   >>P.S. using the rabbin oracle, see Knuth TAOC.   
   >>   
   > Great disappointment. I thought that dc and python use the same   
   > infinite precision libraries.   
   > Apparently not. Using dc I am substantial slower than python.   
      
   If you want fast bignum arithmetic use gmp. The relevant functions   
   are: 'mpn_mul' (multiplication), 'mpn_tdiv_qr' (quotient with reminder),   
   'mpn_gcd' (GCD) and possibly 'mpn_add' and 'mpn_sub'. For the last two   
   it is not hard to write reasonably good implementation in assembly,   
   and for big numbers most execution time goes into multiplicative   
   operations, that is why they are not critical. Other are pretty complex   
   to get good speed. Of course, you can use more functions to get   
   somewhat better speed.   
      
   Note that converting to text or from text has similar cost to   
   multiplication, so if you use text interface for each operation   
   you will _not_ get full speed.   
      
   --   
    Waldek Hebisch   
      
   --- SoupGate-Win32 v1.05   
    * Origin: you cannot sedate... all the things you hate (1:229/2)   
|