;;; mathematical utilities Copyright Igor Rivin 2004,2005, all rights reserved.

(define sgn (lambda (x) 
              (cond (( > x 0) 1)
                    ((< x 0) -1)
                    (else 0))))

(define pi (acos -1.0))


;recursive horner's rule  
(define horneval 
  (lambda (coeflist var)
    (if (null? coeflist)
        0
        (+ (car coeflist) (* var (horneval (cdr coeflist) var))))))

; the iterative version of horner's rule
(define horneval2
  (lambda (coeflist var)
    (cond ((null? coeflist) 0)
          ((null? (cdr coeflist)) (car coeflist))
          (else 
           (let ((thefunc 
                  (lambda (x1 x2)
                    (+ (* var x1) x2))))
             (fold thefunc (reverse coeflist) 0))))))



;; returns all permutations of a list l.

(define allperms
  (lambda (l)
    (if (null? l) (list l)
        (apply append (map (lambda (x)
               (map (lambda (y)
                      (cons (car x) y))
                    (allperms (cdr x))))(n-1sets l))))))

;; computes the list of k-subsets of a list l.

(define subsets 
  (lambda (l k)
    (if (> k (length l))
        '()
        (if (= k 0) (list '())
            (append (subsets (cdr l) k) 
                   (map (lambda (x) (cons (car l) x)) (subsets (cdr l) (- k 1))))))))



(define randperm
  (lambda (ar) 
    (let ((swap (lambda (i j)
                  (let ((tmp (vector-ref ar i)))
                    (vector-set! ar i (vector-ref ar j))
                    (vector-set! ar j tmp))))
          (len (vector-length ar)))
      (let araux ((thelen len))
        (if (= thelen 0)
            ar
            (let ((randind (random thelen)))
              (swap (- thelen 1) randind)
              (araux (- thelen 1))))))))



;; computes the product of two matrices m1 and m2 given as lists of lists

(define mat-mult
  (lambda (m1 m2)
    (outer ldot m1 (transpose2 m2))))


;; fixed point of a function. The arguments are func (the function) the starting value and the tolerance for the equality comparison.
(define fixed-point 
  (lambda (func start tol)
    (let composer ((curval start) (newval (func start)))
      (if (< (abs (- curval newval)) tol)
          newval
          (composer newval (func newval))))))



; rotate a list one place to the left.

(define simplerot 
  (lambda (l)
    (append (cdr l) (list (car l)))))

; all cyclic permutations of a list

(define n-1sets 
  (lambda (l)
    (letrec ((allrots (lambda (l)
                        (if (null? l) (list l)
                            (let ((oldcar (car l)))
                              (let rotator ( (lastrot l) (theres (list l)))
                                (let ((rotated (simplerot lastrot)))
                                  (if (eq? (car rotated) oldcar)
                                      theres
                                      (rotator rotated (cons rotated theres))))))))))
      (allrots l))))

;; Schur product of two matrices

(define schurprod
  (lambda (mat1 mat2)
    (map (lambda (x y) (map (lambda (z w) (* z w)) x y)) mat1 mat2)))

;; L2 norm of a list of numbers
              
(define ll2norm
  (lambda (l)
    (sqrt (apply + (map (lambda (x) (* x x)) l)))))

;; makes a random nxm matrix of zeros and ones (entries uniformly distributed).
(define rand01
  (lambda (n m)
    (outer (lambda (x y) (random (+ x y)))
           (constlist n 1)
           (constlist m 1))))

(define zerosearch ;; works for monotonically INCREASING functions
  (lambda (func high low start eps)
    (let theloop ((newguess start) (lowguess low) (highguess high)(oldguess low))
      (if (< (abs (- oldguess newguess)) eps) 
          newguess
          (let ((startval (func newguess)))
            (cond ((> startval 0)
                   (theloop (* 0.5 (+ lowguess newguess)) lowguess newguess newguess))
                   ((< startval 0) 
                    (theloop (* 0.5 (+ highguess newguess)) newguess highguess newguess))
                   (else newguess)))))))