[ prog / sol / mona ]

prog


/prog/ challenge #5 - Code Archaeology

2 2021-02-14 21:32

The Swing-Factorial
https://archive.tinychan.net/read/prog/1301799586/13-15

#lang racket
(require srfi/1)
(require math/number-theory rnrs/arithmetic/bitwise-6)
(define small-odd-swing
  (vector 1 1 1 3 3 15 5 35 35 315 63 693 231 3003 429 6435 6435 109395
          12155 230945 46189 969969 88179 2028117 676039 16900975 1300075
          35102025 5014575 145422675 9694845 300540195 300540195))
(define (sieve n)
  (let*
      ((limit (round (/ (+ 1 (inexact->exact (round (sqrt n))))
                        2)))
       (m (round (/ (+ 1 n) 2)))
       (ar (make-vector m 1)))
    (for ([i (in-range 1 limit)])
         (when (= (vector-ref ar i) 1)
           (let
               ((p (+ (* 2 i) 1)))
             (for ([j (in-range (* (+ p 1) i) m p)])
                  (vector-set! ar j 0)))))
    (let ((result
           (for/vector
            ([(x i) (in-indexed (in-vector ar))]
             #:when (> x 0))
            (+ 1 (* 2 i)))))
      (vector-set! result 0 2)
      (vector->list result))))

(define (Π f n m #:threads (c (processor-count)))
  (let* ((t (lambda (a)
              (lambda ()
                (let/ec k
                  (do ((x (+ n a) (+ x c))
                       (r 1 (* r (f x))))
                    ((> x m) (k r)))))))
         (xs (do ((x (sub1 c) (sub1 x))
                  (xs '() (cons (future (t x)) xs)))
               ((zero? x) xs)))
         (x ((t 0))))
    (foldr * x (map touch xs))))

(define (swing n)
  (if (< n 33) (vector-ref small-odd-swing n)
      (let* ((primes (cdr (sieve n)))
             (rootn (integer-sqrt n))
             (primesa (take-while (curryr <= rootn) primes))
             (primesb (filter (λ (x) (odd? (quotient n x)))
                                 (take-while (curryr <= (quotient n 3))
                                             (drop-while (curryr <= rootn) primes)))))
        (apply
         *
         (append
          (map
           (λ (prime)
             (do ((q (quotient n prime) (quotient q prime))
                  (p 1 (if (odd? q) (* p prime) p)))
               ((zero? q) p)))
           primesa)
          (take-while (curryr <= n)
                      (drop-while (curryr <= (quotient n 2)) primes))
          primesb)))))

(define (rec-fact n)
  (if (< n 2) 1
      (* (expt (rec-fact (quotient n 2)) 2) (swing n))))

(define (fact n)
  (if (< n 20) (Π values 1 n)
      (arithmetic-shift (rec-fact n) (- n (bitwise-bit-count n)))))
> (time (fact 1000))
cpu time: 1 real time: 0 gc time: 0
402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
> (collect-garbage)
> (time (void (fact 100000)))
cpu time: 143 real time: 143 gc time: 80
> (collect-garbage)
> (time (void (fact 1000000)))
cpu time: 7306 real time: 7333 gc time: 3850
8


VIP:

do not edit these