;;;; srfi-27-numbers.scm
;;;; Kon Lovett, Feb '10

(module srfi-27-numbers

  (;export
    ;
    check-integer
    #;check-cardinal-integer
    check-positive-integer
    #;check-real
    #;check-nonzero-real
    #;check-nonnegative-real
    #;check-positive-real
    #;check-real-open-interval
    #;check-real-closed-interval
    check-real-precision
    #;check-real-unit
    ;
    random-large-integer
    random-large-real
    ;
    native-real-precision?)

  (import
    (except scheme
      <= < zero? positive? negative?
      + * - / quotient expt
      integer? real?
      exact->inexact inexact->exact
      floor)
    chicken
    (only numbers
      <= < zero? positive? negative?
      + * - / quotient expt
      integer? real?
      exact->inexact inexact->exact
      floor)
    (only type-checks
      check-real)
    (only type-errors
      error-argument-type error-open-interval error-closed-interval))

  (require-library numbers type-errors)

  (declare
    (not usual-integrations
      <= < zero? positive? negative?
      + * - / quotient expt
      integer? real?
      exact->inexact inexact->exact) )

;;;

;;

(define (check-integer loc obj #!optional argnam)
  (unless (integer? obj)
    (error-argument-type loc obj "integer" argnam)) )

#;
(define (check-cardinal-integer loc obj #!optional argnam)
  (unless (and (integer? obj) (<= 0 obj))
    (error-argument-type loc obj "cardinal-integer" argnam)) )

(define (check-positive-integer loc obj #!optional argnam)
  (unless (and (integer? obj) (positive? obj))
    (error-argument-type loc obj "positive-integer" argnam)) )

;;

#;
(define (check-real loc obj #!optional argnam)
  (unless (real? obj)
    (error-argument-type loc obj "real" argnam)) )

#;
(define (check-nonzero-real loc obj #!optional argnam)
  (unless (and (real? obj) (not (zero? obj)))
    (error-argument-type loc obj "nonzero-real" argnam)) )

#;
(define (check-nonnegative-real loc obj #!optional argnam)
  (unless (and (real? obj) (not (negative? obj)))
    (error-argument-type loc obj "nonnegative-real" argnam)) )

#;
(define (check-positive-real loc obj #!optional argnam)
  (unless (and (real? obj) (positive? obj))
    (error-argument-type loc obj "positive-real" argnam)) )

;;

(define (check-real-open-interval loc obj mn mx #!optional argnam)
  (check-real loc obj argnam)
  (unless (< mn obj mx)
    (error-open-interval loc obj mn mx argnam)) )

#;
(define (check-real-closed-interval loc obj mn mx #!optional argnam)
  (check-real loc obj argnam)
  (unless (<= mn obj mx)
    (error-closed-interval loc obj mn mx argnam)) )

(define (check-real-precision loc obj #!optional argnam)
  (check-real-open-interval loc obj 0 1 argnam) )

#;
(define (check-real-unit loc obj #!optional argnam)
  (check-real-closed-interval loc obj 0 1 argnam) )

;;;

; 'max  - maximum range of "core" random integer generator (maybe inexact)
; 'm    - exact maximum range (passed so conversion performed by caller)

; Chicken 'numbers' egg determines precision on arg1 so
; convert result of the "base" random integer (which may return a
; flonum due to integer representation issues) into an exact integer
; to prevent infinity generation when the value of n * m is very large.
; (the conversion of the large exact integer to a flonum may result in
; +inf)

(define (random-power rndint state max m k) ; n = m^k, k >= 1
  (do ((k k (fx- k 1))
       (n
         (inexact->exact (rndint state max))
         (+ (inexact->exact (rndint state max)) (* n m))) )
      ((fx= 1 k) n) ) )

; Large Integers
; ==============
;
; To produce large integer random deviates, for n > m, we first
; construct large random numbers in the range {0..m^k-1} for some
; k such that m^k >= n and then use the rejection method to choose
; uniformly from the range {0..n-1}.

(define (random-large-integer rndint state max m n) ; n > m
  (do ((k 2 (fx+ k 1))
       (mk (* m m) (* mk m)))
      ((<= n mk)
        (let* ((mk-by-n (quotient mk n))
               (a (* mk-by-n n)) )
          (let loop ()
            (let ((x (random-power rndint state max m k)))
              (if (< x a)
                (quotient x mk-by-n)
                (loop) ) ) ) ) ) ) )

; Multiple Precision Reals
; ========================
;
; To produce multiple precision reals we produce a large integer value
; and convert it into a real value. This value is then normalized.
; The precision goal is prec <= 1/(m^k + 1), or 1/prec - 1 <= m^k.
; If you know more about the floating point number types of the
; Scheme system, this can be improved.

(define (random-large-real rndint state max m prec)
  (do ((k 1 (fx+ k 1))
       (u (- (/ 1 prec) 1) (/ u m)) )
      ((<= u 1)
        (exact->inexact
          (/
            (+ 1 (random-power rndint state max m k))
            (+ 1 (expt m k)))) ) ) )

;;;

(define (native-real-precision? prec max)
  (or
    (not prec)
    (<= (- (floor (/ 1 prec)) 1) max)) )

) ;module srfi-27-numbers