[plt-scheme] Some numeric functions in Scheme

Hi,

Following a previous thread, this is just a couple of function to
compute some number theory and some other stuff in scheme, this is just
some loose code that I developed through the years. Some functions were
written for fun, others were written so that I could learn a particular
concept. I cannot assure you they have no bugs since I've not tested
them (except some tries to some some exercises or to compute simple
values) but they're really fun. 

Next is code for number-theory functions I developed some years ago:
;; Number theoretic functions
(module numbertheory
  mzscheme
  (require (lib "list.ss"))
   
  (define sieve
    (lambda (n)
      (define filter-vec
        (lambda (vec)
          (let recur ((l null) (i 0))
            (cond ((= i (- (vector-length vec) 1))
                   (reverse l))
                  ((zero? (vector-ref vec i))
                   (recur (cons i l) (+ i 1)))
                  (else (recur l (+ i 1)))))))
      (cddr (filter-vec (let cut ((vec (make-vector n)) (cutn 2))
                          (if (<= cutn
                                  (sqrt n))
                              (begin0
                                (do ((j (* cutn 2) (+ cutn j)))
                                  ((>= j (- n 1)) vec)
                                  (vector-set! vec j 1))
                                (cut vec (+ cutn 1)))))))))
  (define sq (lambda (x) (* x x)))
   
  (define expnt
    (lambda (x n)
      (define expnt-int
        (lambda (x n)
          (cond ((even? n) (sq (expnt-int x (/ n 2))))
                ((= n 1) x)
                (else
                 (* x (expnt-int x (- n 1)))))))
      (if (and (integer? x)
               (integer? n))
          (expnt-int x n)
          (expt x n))))
   
  (define expnt-mod
    (lambda (x n m)
      (cond ((even? n) (modulo (sq (expnt-mod x (/ n 2) m)) m))
            ((= n 1) x)
            (else
             (modulo (* x (expnt-mod x (- n 1) m)) m)))))
   
  (define map-expnt-mod
    (lambda (l n m)
      (map (lambda (x)
             (expnt-mod x n m))
           l)))
   
   
  ;Returns a list with three elements:
  ;the gcd of a and b,
  ;the multiplicative coeficient of (max a b)
  ;and the multiplicative coeficient of (min a b)
  (define ext-gcd
    (lambda (a b)
      (define ext-gcd-aux
        (lambda (a b q x1 x2 y1 y2)
          (let* ((id (quotient a b))
                 (r (- a (* b id))))
            (if (zero? r)
                (list b x2 y2)
                (ext-gcd-aux b
                             r
                             id
                             x2
                             (- x1 (* id x2))
                             y2
                             (- y1 (* id y2)))))))
      (let ((res (ext-gcd-aux (max a b) (min a b) 0 1 0 0 1)))
        (if (< a b)
            (begin
              (list (first res) (third res) (second res)))
            res))))
   
  (define field-inverse
    (lambda (x n)
      (if (> x n)
          (error 'field-inverse "there isn't any ~a in Z/~aZ." x n)
          (let ((val (ext-gcd n x)))
            (if (not (= (first val) 1))
                (error 'field-inverse "~a is not prime with ~a, so it
has no inverse in Z/~aZ." x n n)
                (third val))))))
   
  (define mod
    (lambda (a b)
      (if (and (positive? a)
               (positive? b))
          (modulo a b)
          (error 'mod "both arguments have to be positive."))))
   
  (define field-generator?
    (lambda (g q)
      (define fg-aux
        (lambda (i vec)
          (let ((val (expnt-mod g i q)))
            (cond ((= i q) #t)
                  ((not (zero? (vector-ref vec val)))
                   #f)
                  (else
                   (vector-set! vec val 1)
                   (fg-aux (+ i 1) vec))))))
      (fg-aux 1 (make-vector q))))
   
  (define field-generator
    (lambda (q)
      (let recur ((n 1))
        (cond ((field-generator? n q)
               n)
              ((= n (- q 1))
               (error "No generators found for Z/~aZ." q))
              (else (recur (+ n 1)))))))
   
  ; Solves systems of congruences using the chinese remainder theorem
  ; Receives n the number of eqs and l the list of lists ((a1 m1) (a2
m2) ... (a3 m3))
  ; with system x = a mod m
  (define crt
    (lambda (l)
      (call-with-values
       (var-vectors (length l) l)
       (lambda (a m) (crt-vector a m)))))
  (define var-vectors
    (lambda (n l)
      (make-recur n n l (make-vector n) (make-vector n))))
  (define make-recur
    (lambda (n nn
               nl
               avec
               mvec)
      (if (null? nl)
          (values avec mvec)
          (begin
            (vector-set! avec (- nn n) (first (car nl)))
            (vector-set! mvec (- nn n) (second (car nl)))
            (make-recur n (+ 1 nn) (cdr nl) avec mvec)))))
   
  (define crt-vector
    (lambda (a m)
      (let ((bigm (multiply-all m (vector-length m)))
            (solution 0))
        (do ((i 0 (+ 1 i)))
          ((= i (vector-length m)) solution)
          (set! solution (+ solution
                            (* (vector-ref a i)
                               (/ bigm (vector-ref m i))
                               (second (ext-gcd (/ bigm (vector-ref m
i))
                                                (vector-ref m
i))))))))))
   
  (define multiply-all
    (lambda (v l)
      (let loop ((prod 1)
                 (i 0))
        (if (= i l) prod
            (loop (* (vector-ref v i) prod) (+ 1 i))))))
   
  (provide
   (all-defined-except crt-vector multiply-all)))


Next, some functions to calculate derivatives.
;;
;; Symbolic derivatives
;;
 
(require (lib "list.ss"))
(require (lib "trace.ss"))
 
(load "simplify.ss")
 
(define sder
  (lambda (f v)
    (if (not (list? f))
        (cond ((eq? f v) 1)
              (else 0))
        (let ((op (first f))
              (args (rest f)))
          (cond ((eq? op '+)
                 (cons '+
                       (map (lambda (x) (sder x v)) args)))
                ((eq? op '-)
                 (cons '-
                       (map (lambda (x) (sder x v)) args)))
                ((eq? op '/)
                 `(/ (- (* ,(sder (first args) v) ,(second args)) (*
,(first args) ,(sder (second args) v)))
                     (pow ,(second args) 2)))
                ((eq? op '*)
                 `(+ (* ,(sder (first args) v) ,(second args))
                     (* ,(first args) ,(sder (second args) v))))
                (else
                 `,((namespace-variable-value (string->symbol
(string-append "der-" (symbol->string op)))) args v)))))))
 
(define der-pow
  (lambda (f v)
    (if (> (length f) 2)
        (error "pow: should receive 2 arguments of power, got f(~a) =
(pow ~a)" v f)
        (let ((base (first f))
              (exponent (second f)))
          `(* ,exponent (pow ,base (- ,exponent 1)) ,(sder base v))))))
 
(define der-sin
  (lambda (f v)
    (if (> (length f) 1)
        (error "sin: should receive 1 argument of sine, got f(~a) = (sin
~a)" v f)
        `(* ,(sder (first f) v)
            (cos ,(first f))))))
 
(define der-cos
  (lambda (f v)
    (if (> (length f) 1)
        (error "cos: should receive 1 argument of cosine, got f(~a) =
(cos ~a)"  v f)
        `(* ,(sder (first f) v)
            (- (sin ,(first f)))))))
 
(define der-tan
  (lambda (f v)
    (if (> (length f) 1)
        (error "tan: should receive 1 argument of tangent, got f(~a) =
(tan ~a)"  v f)
        `(* ,(sder (first f) v)
            (pow (sec ,(first f)) 2)))))
 
(define der-cotan
  (lambda (f v)
    (if (> (length f) 1)
        (error "cotan: should receive 1 argument of cotangent, got f(~a)
= (cotan ~a)"  v f)
        `(* ,(sder (first f) v)
            (- (pow (cosec ,(first f)) 2))))))
 
(define der-sec
  (lambda (f v)
    (if (> (length f) 1)
        (error "sec: should receive 1 argument of secant, got f(~a) =
(sec ~a)"  v f)
        `(* ,(sder (first f) v)
            (sec ,(first f))
            (tan ,(first f))))))
 
(define der-cosec
  (lambda (f v)
    (if (> (length f) 1)
        (error "cosec: should receive 1 argument of cosecant, got f(~a)
= (cosec ~a)"  v f)
        `(* ,(sder (first f) v)
            (- (cosec ,(first f)))
            (cotan ,(first f))))))
 
(define der-arcsin
  (lambda (f v)
    (if (> (length f) 1)
        (error "arcsin: should receive 1 argument of arc-sine, got f(~a)
= (arcsin ~a)" v f)
        `(/ ,(sder (first f) v)
            (sqrt (- 1 (pow ,(first f) 2)))))))
 
(define der-arccos
  (lambda (f v)
    (if (> (length f) 1)
        (error "arccos: should receive 1 argument of arc-cosine, got
f(~a) = (arccos ~a)" v f)
        `(- (/ ,(sder (first f) v)
               (sqrt (- 1 (pow ,(first f) 2))))))))
 
(define der-arctan
  (lambda (f v)
    (if (> (length f) 1)
        (error "arctan: should receive 1 argument of arc-tangent, got
f(~a) = (arctan ~a)" v f)
        `(/ ,(sder (first f) v)
            (+ 1 (pow ,(first f) 2))))))
 
(define der-arccotan
  (lambda (f v)
    (if (> (length f) 1)
        (error "arccotan: should receive 1 argument of arc-cotangent,
got f(~a) = (arccotan ~a)" v f)
        `(- (/ ,(sder (first f) v)
               (+ 1 (pow ,(first f) 2)))))))
 
(define der-arcsec
  (lambda (f v)
    (if (> (length f) 1)
        (error "arcsec: should receive 1 argument of arc-secant, got
f(~a) = (arcsec ~a)" v f)
        `(/ ,(sder (first f) v)
            (* ,(first f) (sqrt (- (pow ,(first f) 2) 1)))))))
 
(define der-arccosec
  (lambda (f v)
    (if (> (length f) 1)
        (error "arccosec: should receive 1 argument of arc-cosecant, got
f(~a) = (arccosec ~a)" v f)
        `(- (/ ,(sder (first f) v)
            (* ,(first f) (sqrt (- (pow ,(first f) 2) 1))))))))
 
(define der-ln
  (lambda (f v)
    (if (> (length f) 1)
        (error "ln: should receive 1 argument of natural logarithm, got
f(~a) = (ln ~a)" v f)
        `(/ ,(sder (first f) v)
            (first f)))))
         
         
(trace sder)


 I think the concept to use a namespace search is very nice.
Extensibility is awesome. Just add a new function with the name
der-<functionname> and then when trying to compute the derivative of
that function, scheme handles the search of the respective function.
It's very, very nice, however, it does not simplify its output and
that's a shame! :D

Well, I hope to use some of those functions and ideas in my new Computer
Algebra System which is being developed in scheme :D. Very, very nice.

I'll post some code just to try it out asap.

Cheers,
-- 

Paulo J. Matos : pocm [_at_] mega . ist . utl . pt
Instituto Superior Tecnico - Lisbon
Computer and Software Eng. - A.I.
 - > http://mega.ist.utl.pt/~pocm
---
        -> God had a deadline...
                So, he wrote it all in Lisp!