Guile-hoot/module/hoot/numbers.scm

;;; Bytevectors
;;; Copyright (C) 2024 Igalia, S.L.
;;;
;;; Licensed under the Apache License, Version 2.0 (the "License");
;;; you may not use this file except in compliance with the License.
;;; You may obtain a copy of the License at
;;;
;;;    http://www.apache.org/licenses/LICENSE-2.0
;;;
;;; Unless required by applicable law or agreed to in writing, software
;;; distributed under the License is distributed on an "AS IS" BASIS,
;;; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
;;; See the License for the specific language governing permissions and
;;; limitations under the License.

;;; Commentary:
;;;
;;; Bytevectors.
;;;
;;; Code:

(library (hoot numbers)
  (export + * - /
          < <= = >= >

          1+ 1-

          abs
          floor
          ceiling
          round
          truncate
          number?
          complex?
          real?
          rational?
          integer?
          exact-integer?
          exact?
          inexact?

          finite?
          infinite?
          nan?
          inexact
          exact
          quotient
          remainder
          modulo
          even?
          odd?
          numerator
          denominator
          exact-integer-sqrt
          floor/

          floor-quotient
          floor-remainder
          truncate/
          truncate-quotient
          truncate-remainder
          gcd lcm

          max min

          negative?
          positive?
          zero?
          sin
          cos
          tan
          asin
          acos

          atan
          sqrt
          log
          exp

          rationalize
          square
          expt

          make-rectangular
          make-polar
          magnitude
          angle
          real-part
          imag-part)
  (import (hoot primitives)
    (hoot values)
    (hoot eq)
    (hoot not)
    (hoot errors)
    (hoot match)
    (hoot bitwise))

  (define (1+ x) (%+ x 1))
  (define (1- x) (%- x 1))

  (define-syntax-rule (define-associative-eta-expansion f %f)
    (define f
      (case-lambda
       (() (%f))
       ((x) (%f x))
       ((x y) (%f x y))
       ((x y . z) (apply f (%f x y) z)))))

  (define-associative-eta-expansion * %*)
  (define-associative-eta-expansion + %+)

  (define-syntax-rule (define-sub/div-eta-expansion f %f zero)
    (begin
      (define %generic
        (case-lambda
         ((y) (%f zero y))
         ((x y) (%f x y))
         ((x y . z) (apply %generic (%f x y) z))))
      (define-syntax f
        (lambda (stx)
          (syntax-case stx ()
            ((_ . x) #'(%f . x))
            (f (identifier? #'f) #'%generic))))))

  (define-sub/div-eta-expansion - %- 0)
  (define-sub/div-eta-expansion / %/ 1)

  (define-syntax-rule (define-comparison-expansion f %f)
    (begin
      (define %generic
        (case-lambda
         ((a b) (%f a b))
         ((a b . c)
          (let lp ((res (%f a b)) (a b) (c c))
            (match c
              (() res)
              ((b . c)
               (lp (and (%f a b) res) b c)))))))
      (define-syntax f
        (lambda (stx)
          (syntax-case stx ()
            ((_ x y . z) #'(%f x y . z))
            (f (identifier? #'f) #'%generic))))))

  (define-comparison-expansion < %<)
  (define-comparison-expansion <= %<=)
  (define-comparison-expansion = %=)
  (define-comparison-expansion >= %>=)
  (define-comparison-expansion > %>)

  (define (number? x) (%number? x))
  (define (complex? x) (%complex? x))
  (define (real? x) (%real? x))
  (define (rational? x) (%rational? x))
  (define (integer? x) (%integer? x))
  (define (exact-integer? x) (%exact-integer? x))
  (define (exact? x) (%exact? x))
  (define (inexact? x) (%inexact? x))

  (define (abs x) (%abs x))
  (define (floor x) (%floor x))
  (define (ceiling x) (%ceiling x))
  (define (round x) (%floor (+ x 0.5)))
  (define (truncate x)
    (check-type x real? 'truncate)
    (if (exact? x)
        (if (integer? x)
            x
            (truncate-quotient (numerator x) (denominator x)))
        (%inline-wasm
         '(func (param $x f64) (result f64)
                (f64.trunc (local.get $x)))
         x)))

  ;; Unlike R7RS, these only operate on real numbers.
  (define (infinite? x)
    (or (= x +inf.0) (= x -inf.0)))
  (define (nan? x)
    (not (= x x)))
  (define (finite? x)
    (and (not (infinite? x)) (not (nan? x))))

  (define (inexact x) (%inexact x))
  (define (exact x)
    (cond
     ((exact? x) x)
     (else
      (check-type x finite? 'exact)
      (%inline-wasm
       '(func (param $x f64)
              (result (ref eq))
              (call $f64->exact (local.get $x)))
       x))))

  (define (quotient x y) (%quotient x y))
  (define (remainder x y) (%remainder x y))
  (define (modulo x y) (%modulo x y))

  (define (even? x) (zero? (logand x 1)))
  (define (odd? x) (not (even? x)))

  (define (numerator x)
    (cond
     ((exact-integer? x) x)
     ((exact? x)
      (%inline-wasm
       '(func (param $x (ref $fraction))
              (result (ref eq))
              (struct.get $fraction $num (local.get $x)))
       x))
     (else (inexact (numerator (exact x))))))
  (define (denominator x)
    (cond
     ((exact-integer? x) 1)
     ((exact? x)
      (%inline-wasm
       '(func (param $x (ref $fraction))
              (result (ref eq))
              (struct.get $fraction $denom (local.get $x)))
       x))
     (else (inexact (denominator (exact x))))))
  (define (exact-integer-sqrt n)
    ;; FIXME: There's a compiler bug that makes this procedure return
    ;; junk when this exact-integer? check is enabled.
    ;;
    (check-type n exact-integer? 'exact-integer-sqrt)
    (assert (>= n 0) 'exact-integer-sqrt)
    (let loop ((x n) (y (quotient (+ n 1) 2)))
      (if (< y x)
          (loop y (quotient (+ y (quotient n y)) 2))
          (values x (- n (* x x))))))

  (define (floor/ x y)
    (values (floor-quotient x y) (floor-remainder x y)))
  ;; Adapted from the SRFI-141 reference implementation
  (define (floor-quotient x y)
    (check-type x integer? 'floor-quotient)
    (check-type y integer? 'floor-quotient)
    (cond
     ((and (negative? x) (negative? y))
      (quotient (- x) (- y)))
     ((negative? x)
      (let ((x (- x)))
        (call-with-values (lambda () (truncate/ x y))
          (lambda (q r)
            (if (zero? r)
                (- q)
                (1- (- q)))))))
     ((negative? y)
      (let ((y (- y)))
        (call-with-values (lambda () (truncate/ x y))
          (lambda (q r)
            (if (zero? r)
                (- q)
                (1- (- q)))))))
     (else (quotient x y))))
  (define (floor-remainder x y) (modulo x y))

  (define (truncate/ x y)
    (values (truncate-quotient x y)
            (truncate-remainder x y)))
  (define (truncate-quotient x y) (quotient x y))
  (define (truncate-remainder x y) (remainder x y))

  (define (%binary-gcd x y)
    (check-type x integer? 'gcd)
    (check-type y integer? 'gcd)
    (let ((result
           (%inline-wasm
            '(func (param $x (ref eq)) (param $y (ref eq))
                   (result (ref eq))
                   (call $gcd (local.get $x) (local.get $y)))
            (exact x)
            (exact y))))
      (if (or (inexact? x) (inexact? y))
          (inexact result)
          result)))
  (define-syntax %gcd
    (syntax-rules ()
      ((_) 0)
      ((_ x) x)
      ((_ x y) (%binary-gcd x y))))

  (define (%binary-lcm x y)
    (check-type x integer? 'lcm)
    (check-type y integer? 'lcm)
    (let* ((exact-x (exact x))
           (exact-y (exact y))
           (result (if (and (eqv? exact-x 0) (eqv? exact-y 0))
                       0
                       (quotient (abs (* exact-x exact-y))
                                 (gcd exact-x exact-y)))))
      (if (or (inexact? x) (inexact? y))
          (inexact result)
          result)))
  (define-syntax %lcm
    (syntax-rules ()
      ((_) 1)
      ((_ x) x)
      ((_ x y) (%binary-lcm x y))))

  (define-syntax-rule (define-associative-eta-expansion f %f)
    (define f
      (case-lambda
       (() (%f))
       ((x) (%f x))
       ((x y) (%f x y))
       ((x y . z) (apply f (%f x y) z)))))
  (define-associative-eta-expansion gcd %gcd)
  (define-associative-eta-expansion lcm %lcm)

  (define max
    (case-lambda
     ((x) x)
     ((x y) (if (> x y) x y))
     ((x y . z) (apply max (max x y) z))))
  (define min
    (case-lambda
     ((x) x)
     ((x y) (if (< x y) x y))
     ((x y . z) (apply min (min x y) z))))

  (define (negative? x) (< x 0))
  (define (positive? x) (> x 0))
  (define (zero? x) (= x 0))

  (define (sin x) (%sin x))
  (define (cos x) (%cos x))
  (define (tan x) (%tan x))
  (define (asin x) (%asin x))
  (define (acos x) (%acos x))
  (define atan
    (case-lambda
     ((x) (%atan x))
     ((x y) (%atan x y))))
  (define (sqrt x) (%sqrt x))

  (define* (log x #:optional y)
    (define (%log x)
      (%inline-wasm
       '(func (param $x (ref eq)) (result (ref eq))
              (call $log (local.get $x)))
       x))
    (if y
        (/ (%log x)
           (%log y))
        (%log x)))

  (define (exp x)
    (define (%exp x)
      (%inline-wasm
       '(func (param $x (ref eq)) (result (ref eq))
              (call $exp (local.get $x)))
       x))
    (%exp x))

  ;; Adapted from the comments for scm_rationalize in libguile's numbers.c
  (define (rationalize x y)
    (check-type x rational? 'rationalize)
    (check-type y rational? 'rationalize)
    (define (exact-rationalize x eps)
      (let ((n1  (if (negative? x) -1 1))
            (x   (abs x))
            (eps (abs eps)))
        (let ((lo (- x eps))
              (hi (+ x eps)))
          (if (<= lo 0)
              0
              (let loop ((nlo (numerator lo)) (dlo (denominator lo))
                         (nhi (numerator hi)) (dhi (denominator hi))
                         (n1 n1) (d1 0) (n2 0) (d2 1))
                (let-values (((qlo rlo) (floor/ nlo dlo))
                             ((qhi rhi) (floor/ nhi dhi)))
                  (let ((n0 (+ n2 (* n1 qlo)))
                        (d0 (+ d2 (* d1 qlo))))
                    (cond ((zero? rlo) (/ n0 d0))
                          ((< qlo qhi) (/ (+ n0 n1) (+ d0 d1)))
                          (else (loop dhi rhi dlo rlo n0 d0 n1 d1))))))))))
    (if (and (exact? x) (exact? y))
        (exact-rationalize x y)
        (inexact (exact-rationalize (exact x) (exact y)))))

  (define (square x) (* x x))

  (define (expt x y)
    (check-type x number? 'expt)
    (check-type y number? 'expt)
    (cond
     ((eqv? x 0)
      (cond ((zero? y) (if (exact? y) 1 1.0))
            ((positive? y) (if (exact? y) 0 0.0))
            (else +nan.0)))
     ((eqv? x 0.0)
      (cond ((zero? y) 1.0)
            ((positive? y) 0.0)
            (else +nan.0)))
     ((exact-integer? y)
      (if (< y 0)
          (/ 1 (expt x (abs y)))
          (let lp ((y y)
                   (result 1))
            (if (= y 0)
                result
                (lp (1- y) (* x result))))))
     (else (exp (* y (log x))))))

  ;; (scheme complex)
  ;; Adapted from Guile's numbers.c
  (define (make-rectangular real imag)
    (check-type real real? 'make-rectangular)
    (check-type imag real? 'make-rectangular)
    (if (eq? imag 0)
        real
        (%inline-wasm
         '(func (param $real f64) (param $imag f64) (result (ref eq))
                (struct.new $complex
                            (i32.const 0)
                            (local.get $real)
                            (local.get $imag)))
         (inexact real) (inexact imag))))
  (define (make-polar mag ang)
    (check-type mag real? 'make-polar)
    (check-type ang real? 'make-polar)
    (cond
     ((eq? mag 0) 0)
     ((eq? ang 0) mag)
     (else
      (%inline-wasm
       '(func (param $mag f64) (param $ang f64) (result (ref eq))
              (local $f0 f64) (local $f1 f64)
              (local.set $f0 (call $fcos (local.get $ang)))
              (local.set $f1 (call $fsin (local.get $ang)))
              ;; If sin/cos are NaN and magnitude is 0, return a complex
              ;; zero.
              (if (ref eq)
                  (i32.and (call $f64-is-nan (local.get $f0))
                           (call $f64-is-nan (local.get $f1))
                           (f64.eq (local.get $mag) (f64.const 0.0)))
                  (then (struct.new $complex
                                    (i32.const 0)
                                    (f64.const 0.0)
                                    (f64.const 0.0)))
                  (else (struct.new $complex
                                    (i32.const 0)
                                    (f64.mul (local.get $mag) (local.get $f0))
                                    (f64.mul (local.get $mag) (local.get $f1))))))
       (inexact mag) (inexact ang)))))
  (define (magnitude z)
    (cond
     ((real? z) (abs z))
     (else
      (check-type z complex? 'magnitude)
      (let ((r (real-part z))
            (i (imag-part z)))
        (sqrt (+ (* r r) (* i i)))))))
  (define (angle z)
    (cond
     ((real? z)
      (if (negative? z)
          (atan 0.0 -1.0)
          0.0))
     (else
      (check-type z complex? 'angle)
      (atan (imag-part z) (real-part z)))))
  (define (real-part z)
    (cond
     ((real? z) z)
     (else
      (check-type z complex? 'real-part)
      (%inline-wasm
       '(func (param $z (ref $complex)) (result f64)
              (struct.get $complex $real (local.get $z)))
       z))))
  (define (imag-part z)
    (cond
     ((real? z) 0.0)
     (else
      (check-type z complex? 'real-part)
      (%inline-wasm
       '(func (param $z (ref $complex)) (result f64)
              (struct.get $complex $imag (local.get $z)))
       z))))

  )