;;; This R7RS program tests the efficiency of Marc Nieper-Wißkirchen's
;;; semantics for the R6RS x|p notation.
;;;
;;; It has been tested in Larceny and in Sagittarius.
;;; (To run with Sagittarius v0.8.4, one assertion must be disabled.)
;;;
;;; This is an R7RS program because the R7RS (small) language
;;; includes procedures that can be used for timing benchmarks.
;;; It should not be difficult to translate this into an R6RS
;;; program that uses some auxiliary library for timing.

;;; Marc's preferred semantics says x|p should evaluate to
;;; the result of (inexact (approximate #ex p)),
;;; where the approximate procedure is defined as in
;;;
;;; https://srfi-email.schemers.org/srfi-77/msg/25583008/

;;; Marc's definition of that approximate procedure is written
;;; in R6RS Scheme, so the following R7RS program imports R6RS
;;; libraries.  (Both Larceny and Sagittarius allow that.)
;;; If your implementation of the R7RS doesn't allow R6RS
;;; libraries to be imported, then you have to choose between
;;; translating the following library into R7RS or translating
;;; the main program that follows it into R6RS.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(library (local approximate)
  (export approximate)
  (import (rnrs base)
          (rnrs arithmetic bitwise))

;;; On 04 Sep 2024, Marc Nieper-Wißkirchen posted this code
;;; to the SRFI 77 discussion list.
;;;
;;; https://srfi-email.schemers.org/srfi-77/msg/25583008/

(define approximate
  (lambda (x p)
    (assert (and (exact? x) (rational? x)))
    (assert (and (exact? p) (integer? p) (positive? p)))
    (let ([a (numerator x)]
          [b (denominator x)])
      (let ([sa (bitwise-length a)]
            [sb (bitwise-length b)])
        (let-values
            ([(aa bb)
              (if (>= sa sb)
                  (values a (bitwise-arithmetic-shift-left b (- sa sb)))
                  (values (bitwise-arithmetic-shift-left a (- sb sa)) b))])
          (let-values
              ([(bb sa)
                (if (>= aa bb)
                    (values (bitwise-arithmetic-shift-left bb 1) (+ sa 1))
                    (values bb sa))])
            (let ([aa (bitwise-arithmetic-shift-left aa (+ p 1))]
                  [sb (+ sb p 1)])
              (let-values ([(aa bb) (div-and-mod aa bb)])
                (let
                    ([aa
                      (cond
                       [(not (bitwise-bit-set? aa 0))
                        aa]
                       [(not (zero? bb))
                        (if (positive? aa)
                            (+ aa 1)
                            (- aa 1))]
                       [(not (bitwise-bit-set? aa 1))
                        (if (positive? aa)
                            (- aa 1)
                            (+ aa 1))]
                       [(positive? aa)
                        (+ aa 1)]
                       [else
                        (- aa 1)])])
                  (* aa (expt 2 (- sa sb))))))))))))

) ; end of (local approximate)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(import (scheme base)
        (scheme time)
        (scheme write)
        (scheme process-context)
        (local approximate))

;;; Returns the value of x|p as given by Marc's preferred semantics.

(define (x!p x p)
  (inexact (approximate x p)))

;;; Will Clinger's preferred semantics says x|p should evaluate to
;;; the result of (inexact x) unless p is less than the precision
;;; of (inexact x), in which case x|p might evaluate to an inexact
;;; real approximation to x whose precision is less than p.

;;; Returns the value of x|p as given by Will's preferred semantics,
;;; assuming the precision of (inexact x) is the only precision
;;; available for inexact reals.

(define (x~p x p)
  (inexact x))

(define (assert boolean)
  (if (not boolean)
      (begin (display "Assertion failed.\n")
             (exit))))

;;; In the above, I think Marc and I are both assuming that for
;;; large values of p, x|p evaluates to a floating point value
;;; with the same precision as (inexact x).  Let's check that.

;;; The following assertion fails in Sagittarius 0.8.4, apparently
;;; because x|128 is read with less precision than plain old x.
;;; I assume that's just a bug in that old version of Sagittarius.
;;; Rather than install a newer version, I simply remove this
;;; assertion when testing with Sagittarius.

(assert (= (inexact #e1234567890987654321.0123456789)
           1234567890987654321.0123456789|128))            ; x|p

;;; In the test cases that follow, I will assume (inexact x)
;;; returns a floating point value with 53 bits of precision,
;;; as in IEEE double precision.  Let's check that.

(assert (= (inexact #e1.0e23)
           99999999999999991611392))

;;; Let's compare the efficiency of Marc's preferred semantics
;;; against Will's preferred semantics.

(define (compare-efficiency make-test-cases n msg semantics1 semantics2)
  (let ((test-cases (make-test-cases n)))
    (define (timing-loop semantics test-cases y)
      (if (null? test-cases)
          y
          (timing-loop semantics
                       (cdr test-cases)
                       (semantics (car test-cases)))))
    (let* ((t1start (current-jiffy))
           (y1 (timing-loop semantics1 test-cases 0.0))
           (t1end (current-jiffy))
           (t2start (current-jiffy))
           (y2 (timing-loop semantics2 test-cases 0.0))
           (t2end (current-jiffy))
           (time1 (- t1end t1start))
           (time2 (- t2end t2start))
           (time1/time2 (/ (round (* 10 (/ time1 time2))) 10.0)))
      ;; To discourage compilers from optimizing away the call to timing-loop:
      (assert (positive? (+ y1 y2)))
      (display "In this particular implementation of Scheme, with ")
      (display msg)
      (display ",\nMarc's preferred semantics is ")
      (display time1/time2)
      (display " times as slow as Will's preferred semantics.\n\n"))))

;;; Returns n randomly selected exact numbers between 1 and 1 million.

(define (make-test-cases n)
  (do ((n n (- n 1))
       (test-cases '() (cons (make-test-case) test-cases)))
      ((= n 0)
       test-cases)))

;;; Returns n copies of the exact integer 1.

(define (make-test-cases1 n)
  (make-list n 1))

;;; For this purpose, the random number generator need not be of high quality.

(define make-test-case
  (let ((seed1 5000000)
        (seed2 1205936125)
        (seed3 6528514201)
        (modulus1 3000000)
        (modulus2 4125015)
        (modulus3 7821551)
        (m1 10251013)
        (m2 55013547)
        (m3 98215133)
        (a1 1025505029)
        (a2 53512053)
        (a3 87012924))
    (lambda ()
      (set! seed1
            (modulo (+ a1 (* m1 seed1))
                    modulus1))
      (set! seed2
            (modulo (+ a2 (* m2 seed2))
                    modulus2))
      (set! seed3
            (modulo (+ a3 (* m3 seed3))
                    modulus3))
      (let ((x (+ seed1
                  (* 13 (/ seed2 10000000))
                  (/ seed3 1000000000))))
        (if (< 1 x 1000000)
            x
            (make-test-case))))))

(compare-efficiency make-test-cases1
                    10000
                    "x = 1/10 and p = 54"
                    (lambda (x) (x!p x 54))
                    (lambda (x) (x~p x 54)))

(compare-efficiency make-test-cases1
                    10000
                    "x = 1/10 and p = 20"
                    (lambda (x) (x!p x 20))
                    (lambda (x) (x~p x 20)))

(compare-efficiency make-test-cases1
                    10000
                    "x = 1/10 and p = 3"
                    (lambda (x) (x!p x 3))
                    (lambda (x) (x~p x 3)))

(compare-efficiency make-test-cases
                    10000
                    "random x and p = 54"
                    (lambda (x) (x!p x 54))
                    (lambda (x) (x~p x 54)))

(compare-efficiency make-test-cases
                    10000
                    "random x and p = 20"
                    (lambda (x) (x!p x 20))
                    (lambda (x) (x~p x 20)))

(compare-efficiency make-test-cases
                    10000
                    "random x and p = 3"
                    (lambda (x) (x!p x 3))
                    (lambda (x) (x~p x 3)))