[-*- Mode: emacs-lisp -*-]

[
 Copyright 2012 Ken Takusagawa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
]

["Can all the pierpont primes be listed, like the Mersennes? Looks like not."]
["This program pauses at the beginning to generate a large array"]

[-----grammar
 (gz start (prog))
 (gz prog (module))
 (gz module (id language-pragma-opt exports imports topdecl-star
                ::pr("[[language-pragma-opt]]module [[id]] [[exports]] where{\n[[imports]]\n"
                     "[[topdecl-star('',';\n','\n')]]}\n\n")))
 (gz language-pragma ( :language-pragma f id-non-star j
                                        ::pr
                                        ("{-# LANGUAGE [[id-non-star('',',','')]] #-}\n")))
 (gz exports (f export-star j
                ::pr("[[export-star('(',', ',')')]]")) (:export-everything ::pr("")))
 (gz export (id)
     (f :module-export id j ::pr("module [[id]]"))
     )
 (gz imports (f import-star j
                ::pr("[[import-star('',';\n',';\n')]]")))
 (gz import(id ::pr( "import [[id]]"))
     (f :qualified (original-name ::is id) (new-name ::is id) j
        ::pr ( "import qualified [[original-name]] as [[new-name]]"))
     (f :specific id id-non-star j
        ::pr ("import [[id]][[id-non-star('(',',',')')]]"))
     (f :hiding import id-non-star j
        ::pr("[[import]] hiding[[id-non-star('(',',',')')]]")
        [dunno if this will work for complicated cases]
        )
     
     )

 (gz type-class (f decl-mark (class-name ::is id) :type-class context-opt  id-non-plus f type-class-decl-star j j
                   ::pr("class [[context-opt]][[class-name]] [[id-non-plus('',' ','')]] where{\n"
                        "[[type-class-decl-star('',';\n','\n')]]}")))
 (gz topdecl (decl) (data)(type-synonym)(newtype)(instance)(type-class))
 (gz type-class-decl (type-signature)(decl))
 (gz type-signature (f :tysig name ret-type-and-params j
                       ::pr("[[name]] :: [[ret-type-and-params]]")))
 (gz instance (f :instance (type ::is id) (name ::is simpletype)  decls j
                 ::pr ("instance [[type]] ([[name]]) where [[decls]]")))
 [(gz instance (f :instance context-opt (type ::is id) simpletype-plus :x decls j
                  ::pr ("instance [[context-opt]][[type]] [[simpletype-plus('(',')(',')')]] where [[decls]]")))]
 (gz newtype (f decl-mark (name ::is id) :newtype type-vars-opt  constr-or-wrap deriving-opt j
                ::pr("newtype [[name]] [[type-vars-opt]] = [[constr-or-wrap{my_name}]][[deriving-opt]]")))
 (gz constr-or-wrap ::gets "tr-id*{name}" (constr)(wrap-constr ::pr("[[wrap-constr{name}]]")))
 (gz wrap-constr ::gets "tr-id*{name}"
     ( :wrap id
             ::pr(
                  (::c "name->print();")
                  " {un_"
                  (::c "name->print();")
                  " :: [[id]]}")))

 (gz deriving (:deriving f id-non-plus j
                         ::pr(" deriving [[id-non-plus('(',', ',')')]]")))
 (gz id-non (id))
 (gz type-synonym (f decl-mark id :type-synonym type-vars-opt type j
                     ::pr("type [[id]] [[type-vars-opt]] = [[type]]")))
 (gz data (f decl-mark id :data type-vars-opt constrs deriving-opt j ::pr("data [[id]] [[type-vars-opt]] = [[constrs]][[deriving-opt]]")))
 (gz simpletype (f id-non-plus j ::pr ("[[id-non-plus('',' ','')]]")))
 (gz type-vars (:args f id-non-star j ::pr ("[[id-non-star('',' ','')]]")))
 (gz constrs(constr-star ::pr("[[constr-star('',' | ','')]]"))

     )
 (gz field-type-and-param (f param type j
                             ::pr("[[param]] :: [[type]]")))
 (gz type-and-param ( f param type j ::pr("[[type]]")))
 (gz constr(positional-constructor) (field-label-constructor) )
 (gz field-label-constructor(f type-ctor :field field-type-and-param-star j
                               ::pr("[[type-ctor]][[field-type-and-param-star('{',', ','}')]]")))

 (gz decls (  decl-star  ::pr("{[[decl-star('\n',';\n','\n')]]}\n")))
 (gz context (:context f a-context-plus j ::pr
                       ("[[a-context-plus('(',', ',')')]] => ")))
 (gz a-context [(f (type ::is id) id-non-plus j
                   ::pr("[[type]] [[id-non-plus('',' ','')]]"))]
     (f (class ::is id) type-plus j
        ::pr("[[class]] [[type-plus('(',')(',')')]]"))
     )
 (gz forall (:forall f id-non-plus j
                     ::pr("forall [[id-non-plus('',' ','')]] . ")))
 (gz ret-type-and-params
     (type f type-and-param-star j forall-opt context-opt
           ::pr("[[forall-opt]][[context-opt]][[type-and-param-star('',' \x2d> ','')]]"
                (::c "if(my_type_and_param_star->v.size()>0)out(' \x2d> ');")
                "[[type]]")))

 (gz decl (f decl-mark name :fun haddock-opt ret-type-and-params  expr j
             ::pr("[[haddock-opt]]"
                  "[[name]] :: [[ret-type-and-params]];\n"
                  "[[name]]"
                  (::c "for(many_trees::const_iterator pos = my_ret_type_and_params->my_type_and_param_star->v.begin();pos!= my_ret_type_and_params->my_type_and_param_star->v.end();++pos){"
                       "const tr_type_and_param* t=dynamic_cast<const tr_type_and_param*>(*pos);"
                       "assert(t);"
                       "out(' ');"
                       "t->my_param->print();" "}")
                  " = [[expr]]"))
     (f decl-mark name :fun :no-sig ret-type-and-params  expr j
        ::pr("[[name]]"
             (::c "for(many_trees::const_iterator pos = my_ret_type_and_params->my_type_and_param_star->v.begin();pos!= my_ret_type_and_params->my_type_and_param_star->v.end();++pos){"
                  "const tr_type_and_param* t=dynamic_cast<const tr_type_and_param*>(*pos);"
                  "assert(t);"
                  "out(' ');"
                  "t->my_param->print();" "}")
             " = [[expr]]"))
     (f decl-mark name :simple  expr j
        ::pr("[[name]] = [[expr]]"))
     )
 (gz name (id))
 (gz positional-constructor ["this one is sketchy"]
     (type-ctor ::pr("[[type-ctor]]"))
     (f   type-ctor typepls-opt j
          ::pr("[[type-ctor]][[typepls-opt]]"))
     (f :tuple type-plus j
        ::pr("[[type-plus('(',', ',')')]]"))
     )
 (gz pattern
     (id)
     (f  pattern-ctor pattern-star j
         ::pr ("([[pattern-ctor]] [[pattern-star('',' ','')]])"))

     (f pattern-ctor :fpat f fpat-star j j
        ::pr ("[[pattern-ctor]][[fpat-star('{',', ','}')]]"))

     (f :ptuple pattern-plus j [pattern-plus cuz :nil exists for empty lists]
        ::pr("[[pattern-plus('(',', ',')')]]"))

     (f :plist pattern-plus j [pattern-plus cuz :nil exists for empty lists]
        ::pr("[[pattern-plus('\x5b',', ','\x5d')]]"))

     (f :pchar astring j ::pr("(\x27[[astring]]\x27)"))

     (f :pstring astring j ::pr("\x22[[astring]]\x22"))

     (f :as id pattern j ::pr("[[id]]@[[pattern]]"))

     )
 (gz pattern-ctor (id) (:cons ::pr ("(:)")) (:nil ::pr ("[]"))
     (:paren id  ::pr("([[id]])")) ["workaround for Qualified parenthesized constructors"]
     ( :capture id type ::pr("[[id]] :: [[type]]"))
     )

 (gz fpat (f (variable ::is id) pattern j
             ::pr("[[variable]] = [[pattern]]")))

 (gz type (f :fn ret-type-and-params j ::pr ("([[ret-type-and-params]])"))
     (:inforall f id-non-plus j type  ::pr("(forall [[id-non-plus('',' ','')]] . [[type]])"))
     (:unit ::pr("()"))
     (positional-constructor))
 (gz typepls (paren-type-plus))
 (gz paren-type (type ::pr( "([[type]])"))
     (f :strict type j ::pr("!([[type]])"))
     (f :generic id j ::pr (" [[id]] ")))
 (gz type-ctor(id)(:list ::pr ("[]"))(:nondet ::pr ("[]")))
 (gz param (pattern))
 (gz qastring  (astring ::pr("\x22[[astring]]\x22")))

 (gz expr (id) (:mcons ::pr ("(:)")) [(:nil ::pr ("[]"))]
     (f :chain astring expr-plus j ::pr 
        ([("[[expr-plus('(',' op ',')')]]")]
         "("
         (::c "for(many_trees::const_iterator pos=my_expr_plus->v.begin();pos!=my_expr_plus->v.end();++pos){")
         (::c "if(pos!=my_expr_plus->v.begin()){")
         "[[astring]]"
         (::c "}(*pos)->print();}")
         ")"
        ))
     (f :join expr-plus j ::pr("[[expr-plus('(',' >>= ',')')]]"))
     (f :cc expr-star j ::pr ("[[expr-star('(',' ++ ',')')]]"))

     (f :rpipe expr-plus j ::pr[("[[expr-star('(',' $ ',')')]]")]
        (["http;//gcc.gnu.org/bugzilla/show_bug.cgi?id=11729"]
         (::c "for(many_trees::reverse_iterator pos = my_expr_plus->v.rbegin();"
              "pos!=my_expr_plus->v.rend();++pos){")
         "("
         (::c "(*pos)->print();" "}")
         (::c "for(many_trees::const_iterator pos = my_expr_plus->v.begin();"
              "pos!=my_expr_plus->v.end();++pos){")
         ")"
         (::c  "}")
         )
        )
     (f :rcompose expr-plus j ::pr
        ("("
         (::c "for(many_trees::reverse_iterator pos = my_expr_plus->v.rbegin();"
              "pos!=my_expr_plus->v.rend();++pos){")
         (::c "if(pos!=my_expr_plus->v.rbegin())")
         " . "
         (::c "(*pos)->print();" "}")
         ")"
         )
        )

     (qastring)
     (f :lit astring j ::pr("[[astring]]"))

     (f :ty type expr j ::pr("([[expr]] :: [[type]])"))
     (f (fun-name ::is expr) expr-star j
        ::pr ("([[fun-name]][[expr-star(' ',' ','')]])"))
     (f :do  stmt-star  j ::pr("(do{\n[[stmt-star('  ','\n  ','\n')]]})"))
     (f :case expr  alt-star  j
        ::pr("(case [[expr]] of {\n[[alt-star('  ',';\n  ','\n')]]})"))
     (f :case expr alt-star :else (underbar ::is expr) j
        ["the else is there so the grammar does not have a reduce/reduce conflict"]
        ::pr("(case [[expr]] of {\n[[alt-star('  ',';\n  ',';\n')]]"
             "  _ -> [[underbar]]\n"
             "})"))
     (f :lcase alt-star j
        ::pr("(\x5clambda_case_var ->"
             "case lambda_case_var of {\n"
             "[[alt-star('  ',';\n  ','\n')]]})"))
     (f :lcase alt-star :else (underbar ::is expr) j
        ::pr("(\x5clambda_case_var ->"
             "case lambda_case_var of {\n"
             "[[alt-star('  ',';\n  ',';\n')]]"
             "  _ -> [[underbar]]\n"
             "})"))

     (f :let decl-star expr j
        ::pr("(let {[[decl-star('\n',';\n','\n')]]}\n in [[expr]])"))
     (f :rlet expr decl-star j
        ::pr("(let {[[decl-star('\n',';\n','\n')]]}\n in [[expr]])"))
     (f :cfd expr assignments-star j
        ::pr("([[expr]][[assignments-star('{',', ','}')]])"))
     (f :mlist expr-star j ::pr("[[expr-star('\x5b',', ','\x5d')]]"))
     (f :cons-list expr-star j ::pr("[[expr-star('(',':',')')]]"))
     (f :mtuple expr-star j ::pr("[[expr-star('(',', ',')')]]"))
     (:nothing  ::pr ("()"))

     (f :lambda name ret-type-and-params expr j
        ::pr("(let {[[name]] :: [[ret-type-and-params]];\n"
             "[[name]]"
             (::c "for(many_trees::const_iterator pos = my_ret_type_and_params->my_type_and_param_star->v.begin();pos!= my_ret_type_and_params->my_type_and_param_star->v.end();++pos){"
                  "const tr_type_and_param* t=dynamic_cast<const tr_type_and_param*>(*pos);"
                  "assert(t);"
                  "out(' ');"
                  "t->my_param->print();" "}")
             " = [[expr]]} in [[name]])"))
     (f :lambda-simple id-non expr j
        [recommended only for reordering arguments to functions
                     and other simple expressions]
        [only one variable to keep it simple]
        ::pr ("(\x5c[[id-non]] -> [[expr]])"))
     (f :field-edit expr field-edit-plus j
        [plus is required by Haskell]
        ::pr("([[expr]][[field-edit-plus('{',',','}')]])"))
     )
 (gz field-edit (f id expr j ::pr ("[[id]] = [[expr]]")))

 (gz assignments (f id expr j ::pr("[[id]] = [[expr]]")))
 (gz stmt (expr ::pr("[[expr]];"))
     (f ":=" pattern type expr j ::pr("[[pattern]] :: [[type]] <- [[expr]];"))
     (f :dlet decl-star j ::pr ("let {[[decl-star('\n',';\n','\n')]]};"))
     )
 (gz alt (f  pattern expr-or-gpat j
             ::pr("[[pattern]][[expr-or-gpat]]")))
 (gz expr-or-gpat
     (expr ::pr ("-> [[expr]]"))
     (where-opt :gpats pred-expr-star [silly lookahead limitation]
        ::pr ("\n[[pred-expr-star('','','')]] [[where-opt]]"
              )))
 (gz pred-expr ( f (pred ::is expr) (do ::is expr) j
                   ::pr ("| [[pred]]\n -> [[do]]\n")))
 (gz where (:where decls ::pr ("where [[decls]]")))

 (gz decl-mark (":"))
 (gz haddock ( :doc f docline-star j 
                    ::pr ("[[docline-star('\n\x7b\x2d |','\n',' \x2d\x7d\n')]]")))
 (gz docline (astring))
 ]



Main :language-pragma
(
 ScopedTypeVariables
 [PatternSignatures]
 GeneralizedNewtypeDeriving
 [NoMonomorphismRestriction]
 RankNTypes
 )
(main)
(
 [ (:hiding Prelude mapM)]
 (:specific System.Environment getArgs)
 (:specific Data.List groupBy sortBy)
 Primality
 IO
 Data.Array.IArray
 Data.Array.Unboxed
 (:specific Control.Exception assert)
 (:specific Control.Monad when)
 (:specific Debug.Trace trace)
 )

(: main :fun (IO :unit) ()
   (:join
    getArgs
    (:lcase
     ((:plist(:pstring "nothing"))(return :nothing))
     ((:plist(:pstring "sortPs"))(:rpipe sortPs (mapM- print)))
     ((:plist(:pstring "numbers"))(:rpipe sortPs (map number) (mapM- print)))
     ((:plist(:pstring "large")n)(:rpipe n read getPs (mapM- print)))
     ((:plist(:pstring "era1"))(:rpipe sortPs (groupBy same-era) (mapM- print)))
     ((:plist(:pstring "era-limit")n m) (:rpipe sortPs (groupBy same-era) (drop (read n)) (take (read m))
                                                (mapM- single-print2)))
     ((:plist(:pstring "era")n)(:do
                                linebuffering
                                (:rpipe sortPs (groupBy same-era) (drop (read n))
                                        (mapM- single-print2))))
     ((:plist(:pstring "era10")n part)(:do
                                      linebuffering
                                      (:rpipe sortPs (groupBy same-era) (drop (read n))
                                              (filter (cpu-divide (read part)))
                                              (mapM- single-print2))))
     ((:plist(:pstring "div-io"))(:rpipe sortPs (mapM- div-io)))
     ((:plist(:pstring "test-lucas")d length)(:rpipe sortPs (drop (read d)) (take (read length)) (mapM- test-lucas-pierpont)))
     ((:plist(:pstring "time")d amount)(:rpipe sortPs (drop (read d)) (take (read amount)) (filter time-pierpont) length print))
     ((:plist(:pstring "nth")n )(:rpipe (!! sortPs (read n)) getlog print))
     
     [((:plist n)(:rpipe
                  n
                  read
                  run1))])
    
    ))
(: P :data (P Integer Integer) :deriving (Show Eq))
(: getlog :fun Double (((P (:capture two Integer)(:capture three Integer))P))
   (:rpipe
    three
    fromInteger
    (* log2-3)
    (+ (fromInteger two))
    ))
(: log2-3 :fun Double ()
   (/ (log 3)(log 2)))
(: ordP :fun Ordering ((x P)(y P))
   (compare (getlog x)(getlog y)))
(: getPs1 :fun (:list P)((max Integer))
   (:do
    (:= x Integer (enumFromTo 1 max))
    (:= y Integer (enumFromTo 0 (floor(/ (fromInteger max) log2-3))))
    (return (P x y))))
(: getPs :fun (:list P)((max Integer))
   (:rpipe
    max getPs1
    (sortBy ordP)
    (takeWhile (/= (P max 0)))))
    
(: five-isPrimeP :fun Bool ((x P))
   (miller_rabin_1
    (:rpipe
     x number
     (+ 1)
     ) 5))

(: isPrimeP :fun Bool ((x P))
   (:rpipe
    x
    isPrimeP-miller-rabin
    (&& (all-check2 x))
    ))
(: isPrimeP-miller-rabin :fun Bool ((x P))
   (:rpipe x number (+ 1) miller-rabin-isPrime))
(: linebuffering :fun (IO :unit)()
   (hSetBuffering stdout LineBuffering)
   )
(: run1 :fun (IO :unit)((n Integer))
   (:do
    linebuffering
   [(:rpipe
    (find-lt (getPs n) (P n 0))
    print)]
   [(single-print (getPs n) n)]
   (mapM- (single-print (getPs n)) (enumFromTo 2 n))
   ))
   
(: find-lt :fun P ((ps(:list P))(target P))
   (:rpipe
    ps
    (takeWhile (/= target))
    reverse
    (dropWhile (:rcompose isPrimeP not))
    head))
(: single-print :fun (IO :unit)((ps(:list P))(i Integer))
   (:do
    (putStr (:cc (show i) ": "))
    (:rpipe
     (find-lt ps (P i 0))
     print)
    ))
   
(: inc2 :fun P (((P(:capture two Integer)(:capture three Integer))P))
   (P (+ 1 two) three))
(: inc3 :fun P (((P(:capture two Integer)(:capture three Integer))P))
   (P two (+ 1 three)))

(: sortPs :fun (:list P)()
   (:mcons (P 1 0) ["must have 2*x to make n+1 odd"]
           (merge
            (map inc2 sortPs)
            (map inc3 sortPs)
            )))
(: merge :fun (:list P)((x(:list P))(y(:list P)))
   (:case (ordP (head x)(head y))
          (LT (:mcons (head x) (merge (tail x)y)))
          (GT (:mcons (head y) (merge x (tail y))))
          (EQ (:mcons (head x) (merge (tail x)(tail y))))
          ))
(: same-era :fun Bool ((x P)(y P))
   (:let
    (: era :fun Int ((i P))
       (:rpipe i getlog floor))
    (== (era x)(era y))))
(: maybe-head :fun (Maybe a)((l(:list a)))
   (:case l
          ((:nil)Nothing)
          :else (:rpipe l head Just)))
(: single-print2 :fun (IO :unit) ((g(:list P)))
   (print
    (:mtuple(:rpipe g head)
            (:rpipe g reverse (dropWhile (:rcompose isPrimeP not)) maybe-head))
    ))
(: number :fun Integer (((P(:capture two Integer)(:capture three Integer))P))
   (gen-number two three))
(: gen-number :fun Integer ((two Integer)(three Integer))
   (* (^ 2 two)(^ 3 three)))
(: Aii :type-synonym (Array (:tuple Int Int) Bool))

(: div-n-array :fun Aii ((n Integer))
   (:let (: minus :fun Int () (:rpipe n fromInteger pred pred))
         (array (:mtuple(:mtuple 0 0)
                        (:mtuple minus minus))
                (:do
                 (:= a Int (enumFromTo 0 minus))
                 (:= b Int (enumFromTo 0 minus))
                 (return (:mtuple (:mtuple a b)
                                  (:rpipe
                                   (gen-number (toInteger a)(toInteger b))
                                   (+ 1)
                                   ((flip mod)n)
                                   (/= 0))))))))
(: Div-array :type-synonym (:tuple Aii Int))
(: make-div-array :fun Div-array ((n Integer))
   (:mtuple (div-n-array n) (fromInteger n)))
(: all-div :fun (:list Div-array) ()
   (:rpipe :nothing primes (drop 2)
           (take 40) (map make-div-array)))

(: div-check :fun Bool (((P(:capture two Integer)(:capture three Integer))P)
                        ((:ptuple (:capture darray Aii)(:capture n Int))Div-array))
   (:let
    (: m :fun Int ((x Integer))
       (mod (fromInteger x)(pred n)))
    (! darray (:mtuple (m two)(m three)))
    ))


(: div5 :fun Bool ((p P))
   (:rpipe
    all-div
    head
    (div-check p)
    ))

(: all-check :fun Bool ((p P))
   (:rpipe
    all-div
    (map (div-check p))
    and))

(: div-io :fun (IO :unit) ((n P))
   (:case
    (== (/= 0 (mod (+ 1 (number n))5))
        (div5 n))
    (False (putStrLn (:cc "differs " (show n))))
    :else [(putStrLn "same")]
    (return :nothing)
    ))


(: lucas-primality-test-1factor :fun Bool ((test Integer)(base Integer)(factor Integer))
   (:case
    (divMod (pred test)factor)
    ((:ptuple(:capture q Integer)(:capture remainder Integer))
     ((assert(== 0 remainder))
      (not-mod-1 base q test)
      ))))

(: not-mod-1 :fun Bool ((base Integer)(exponent Integer)(modulus Integer))
   (/= 1 (modular-exponentiation base exponent modulus)))

(: Lucas-result :data (Prime-lucas-result)(Composite-lucas-result)(Unknown-lucas-result)
   :deriving (Show Eq))
(: lucas-primality-test1 :fun Lucas-result ((test Integer)(factors (:list Integer))(base Integer))
   (:case (miller-rabin-1 test base)
          (False Composite-lucas-result)
          :else
          (:case (and (map (lucas-primality-test-1factor test base) factors))
                 (True Prime-lucas-result)
                 (False Unknown-lucas-result)
                 )))
(: lucas-primality-with-witness :fun (:tuple Integer Bool) ((test Integer)(factors(:list Integer)))
   (:let
    (: loop :fun (:tuple Integer Bool) ((base Integer))
       (:case (lucas-primality-test1 test factors base)
              (Prime-lucas-result (:mtuple base True))
              (Composite-lucas-result (:mtuple base False))
              (Unknown-lucas-result (loop (+ 1 base)))
              ))
    (loop 2)))

(: lucas-pierpont :fun (:tuple Integer Bool) ((p P))
   (lucas-primality-with-witness (+ 1 (number p))
                                 (pierpont-get-factors p)))

(: optimized-lucas-pierpont :fun (:tuple Integer Bool) ((p P))
   (optimized-lucas-primality-with-witness (+ 1 (number p))
                                 (pierpont-get-factors p)))

[(: opt2-lucas-pierpont :fun (:tuple Integer Bool) ((p P))
   (lucas-primality-with-witness (+ 1 (number p))
                                 (optimized-pierpont-get-factors p)))]

(: pierpont-get-factors :fun (:list Integer) (((P(:capture two Integer)(:capture three Integer))P))
   (++
    (:case two
           (0 (:mlist))
           :else (:mlist 2))
    (:case three
           (0 (:mlist))
           :else (:mlist 3))))

[
this does not work, because the (n-1)/2 might be 1 due to repeated squarings of
-1, as with Fermat numbers and base 2.
(: optimized-pierpont-get-factors :fun (:list Integer) (((P _(:capture three Integer))P))
   ["two is not necessary because it automatically gets tested in Miller-Rabin"]
   (:case three
          (0 (:mlist))
          :else (:mlist 3)))]

(: test-lucas-pierpont :fun (IO :unit) ((p P))
   (when (/= (isPrimeP-miller-rabin p)
             (:rpipe p optimized-lucas-pierpont snd))
     (putStrLn (:cc "fail " (show p)))))

(: opt-lucas :fun Lucas-result ((test Integer)(factors (:list Integer))(base Integer))
   (:let
    (: all-factors :fun Integer () (product factors))
    (: remaining :fun (:list Integer) ()
       (:do
        (:= f Integer factors)
        (return (div all-factors f))))
    (:case
     (divMod (pred test) all-factors)
     ((:ptuple(:capture pre-multiplied Integer)(:capture r Integer))
      ((assert (== 0 r))
       ([trace
        (:cc "pre-multiplied " (show pre-multiplied))]
        (:case (miller-rabin-1 test base)
              (False Composite-lucas-result)
              :else
              (:let
               (: almost-there :fun Integer ()
                  (modular-exponentiation base pre-multiplied test))
               (:case (and (:do
                            (:= leftover Integer remaining)
                            [(trace (:cc "not-mod-1 " (show almost-there) " " (show leftover) " " (show test))
                                   (return :nothing))]
                            (return (not-mod-1 almost-there leftover test))
                            ))
                      (True Prime-lucas-result)
                      (False Unknown-lucas-result)))
              )))))))
(: base-sequence :fun (:list Integer)()
   (drop 2 (primes :nothing)))
(: optimized-lucas-primality-with-witness :fun (:tuple Integer Bool) ((test Integer)(factors(:list Integer)))
   (:let
    (: loop :fun (:tuple Integer Bool) (((:cons(:capture base Integer)(:capture more-bases (:list Integer)))
                                         (:list Integer)))
       (:case (opt-lucas test factors base)
              (Prime-lucas-result (:mtuple base True))
              (Composite-lucas-result (:mtuple base False))
              (Unknown-lucas-result (loop more-bases))
              ))
    (loop base-sequence)))
(: time-pierpont :fun Bool ((p P))
   (:rpipe p [optimized-lucas-pierpont snd]
           [isPrimeP-miller-rabin]
           lucas-pierpont snd
           [all-check]
           )
   )
(: powers-mod :fun (:list Int) ((base Int) (modulus Int))
   (:let
    (: f :fun Int ((x Int))
       (mod (* x base) modulus)
       ) 
    (:mcons 1 (map f (powers-mod base modulus))))
   )
(: Adiv :type-synonym (UArray Int Int))
(: single-divisibility-array :fun Adiv ((base Int) (modulus Int))
   (:let (: minus :fun Int () (:rpipe modulus pred pred))
         (listArray (:mtuple 0 minus)
                    (powers-mod base modulus)
                    )
         ))
(: pair-divisibility-array :fun (:tuple Adiv Adiv) ((modulus Int))
   (:mtuple (single-divisibility-array 2 modulus)
            (single-divisibility-array 3 modulus)))
(: all-pair-divisibilities :fun (:list (:tuple Adiv Adiv)) ()
   (:rpipe :nothing primes (drop 2)
           (take 300) (map fromInteger) (map pair-divisibility-array))
   )
(: div-check-pair :fun Bool (((P(:capture two Integer)(:capture three Integer))P)
                             ((:ptuple (:capture atwo Adiv) (:capture athree Adiv))
                              (:tuple Adiv Adiv)))
   (:let
    (: modulus :fun Int ()
       (:rpipe atwo bounds snd (+ 2))
       )
    (: m :fun Int ((i Integer))
       (mod (fromInteger i) (pred modulus)))
    (/= 0
        (mod
         (+ 1
            (* (! atwo (m two)) (! athree (m three))))
         modulus
         )))
   )
(: all-check2 :fun Bool ((p P))
   (:rpipe all-pair-divisibilities
           (map (div-check-pair p))
           and))
(: cpu-divide :fun Bool ((part Integer) (g (:list P)))
   (:let
    (: answer :fun Bool ()
       (== part
           (mod (:case (head g)
                       ((P (:capture two Integer)0)two))
                10)))
    (:rpipe 
     answer 
     [ (trace (:cc "cpu-divide " (show g) " " (show answer)))])
))