!r6rs

(import (rnrs base)
        (rnrs arithmetic bitwise)
        (rnrs control)
        (rnrs io simple))

;;; This is a test of Marc Nieper-Wißkirchen's claims 1 and 2.
;;;
;;; It has been tested in MzScheme, Vicare, Larceny, and Sagittarius.
;;; (To run with Sagittarius v0.8.4, one assertion must be disabled.)

;;; Marc's Claim 1 amounts to saying the R6RS requires x|p
;;; to 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/
;;;
;;; I refer to that Claim 1 as Marc's preferred semantics.

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

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

;;; Marc's Claim 2 says that when p > q, x|p should evaluate
;;; to the same inexact number as x|q.

(define (claim2-holds-for-this-example? x p q)
  (or (not (> p q))
      (= (x!p x p)
         (x!p x q))))

;;; 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))

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

;;; 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))))))))))))

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

;;; Let's see how frequently Marc's Claim 2 holds for this system.

(define (test-claim2 n p q)
  (let ((test-cases (make-test-cases n)))
    (do ((test-cases test-cases (cdr test-cases))
         (failed 0 (if (claim2-holds-for-this-example? (car test-cases) p q)
                       failed
                       (+ failed 1))))
        ((null? test-cases)
         (if (= failed 0)
             (begin
               (display "Marc's Claim 2 was found to be true for each of ")
               (display n)
               (display " random test cases.\n"))
             (begin
               (display "Marc's Claim 2 was found to be false for ")
               (display (/ (round (* 1000 (/ failed n))) 10.0))
               (display "% of ")
               (display n)
               (display " random test cases.\n"))))))
  (newline))

;;; 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)))

;;; 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
            (mod (+ a1 (* m1 seed1))
                 modulus1))
      (set! seed2
            (mod (+ a2 (* m2 seed2))
                 modulus2))
      (set! seed3
            (mod (+ a3 (* m3 seed3))
                 modulus3))
      (let ((x (+ seed1
                  (* 13 (/ seed2 10000000))
                  (/ seed3 1000000000))))
        (if (< 1 x 1000000)
            x
            (make-test-case))))))

(test-claim2 5000 54 53)