xyzzy電卓半端移植+拡張その3 (ひとまず区切り)

xyzzyのナイスな電卓をコマンドラインから使えるように+科学定数を組み込む+単位が(申し訳程度に)扱えるようにする、の第3回。

sigmaが使えるようにしたのと関数が定義できるように直したのとあとこまごました修正。これでhttp://xyzzy.s53.xrea.com/wiki/index.php?Calc-mode%A4%CEManual にある関数はだいたい動くはず。setも使える。これ以上はもうやめておく。大体移植は出来た…のかな*1。こういうのを移植とかじゃなくて書けるのはすごいなあ。

Nos@ubuntu-colunux:~$ calc
$ sigma(j,1,10,5 j J)
275 kg m^2 s^-2

式自体にあんまり意味はないけど \Sigma_{j=1}^{10} 5j[J]相当。似たような感じでproductというのも使えるようにしてある。

数値積分まがいのこともできなくもない。

$ f(x) = exp(- (x^2))/sqrt(pi)
FUNCTION
$ dx = 0.1
0.1
$ sigma(i,-1/dx,1/dx,dx*f(i*dx))
0.8627640914878524179
$ sigma(i,-3/dx,3/dx,dx*f(i*dx))
0.9999841793541379511
$ sigma(i,-5/dx,5/dx,dx*f(i*dx))
0.99999999999911749787

整合が取れていれば単位が入っていてもいけるが面倒かも。

$ f(x) = exp(-(x/m)^2)/sqrt(pi)*m
FUNCTION
$ dx = 0.1m
0.1m
$ sigma(i,-5m/dx,5m/dx,dx*f(i*dx))
0.99999999999911749787 m^2

ビット演算とかもできる。そこら辺はhttp://xyzzy.s53.xrea.com/wiki/index.php?Calc-mode%A4%CEManual 参照。同じ動作をしている、はず。

定数に関する注意点としては

ぐらい? 何が入っているかに関してはunits.lの末尾を参照。mmとか定義してないけどm*2が定義されているのでSI接頭辞を勝手に解釈して動く。ミリリットルのmlとかも同様。

使い方としては第1回と同じく下の2つを適当なディレクトリに放り込んで

 alias calc='(cd ~/path/to/the/directory/; rlwrap clisp calc.l)'

と.bashrcにでも書いておく。

前回まで:

以下コード

SI接頭辞のExaが抜けていてそれ以降がずれていたので修正。
calc.l

;xyzzyに付属するcalc.lを移植してみたもの
;主な変更点(うろ覚え):
;case sensitiveなので sin -> |sin| とかした
;calc-primary-exprのcond節でlparenが捕まっていたので条件の順序を変えた
;
;calc-stringに文字列を渡すと計算する。
;infix-stringに渡すと評価直前のS式の形で返す。
;#$ に続けて計算式を打つとその行を計算式として受け取って計算する。

; 単位とか使えるようにしてみた。

(Defpackage :calc
  (:use :cl)
  (:export #:calc-string #:infix-string #:start-calc))

(defvar *calc-package* (find-package :calc))
(in-package :cl-user)

(defmacro while (p &rest body)
  `(loop while ,p do
	,@body))


(defconstant *calc-token-alist*
  '((#\+ . +) (#\- . -) (#\* . *) (#\/ . /) (#\% . rem) (#\= . =) (#\\ . truncate)
    (#\, . comma) (#\^ . expt) (#\! . !) (#\~ . lognot)
    (#\< (#\< . calc-lshift)) (#\> (#\> . calc-rshift))
    (#\& . logand) (#\| . logior) (#\` . logxor) (#\( . lparen) (#\) . rparen)))
(labels ((set-token-name (c s l)
	   (cond ((listp l)
		  (dolist (x l)
		    (set-token-name (car x) (cons c s) (cdr x))))
		 (s
		  (setf (get l 'calc-operator)
			(format nil "~{~A~}~A" (reverse s) (if c c ""))))
		 (t
		  (setf (get l 'calc-operator) c)))))
  (dolist (x *calc-token-alist*)
    (set-token-name (car x) nil (cdr x))))

(defvar *calc-bits* nil)
(defvar *calc-radix* 10)
(defvar *calc-unsigned* nil)
(defvar *calc-ratio* 'ratio)

(defvar *calc-token* nil)

(load "units.l")

;;;;
(setf (symbol-function 'calc::|gcd|) #'gcd)
(setf (symbol-function 'calc::|lcm|) #'lcm)
(setf (symbol-function 'calc::|exp|) #'exp)
(setf (symbol-function 'calc::|expt|) #'expt)
(setf (symbol-function 'calc::|pow|) #'expt)
(setf (symbol-function 'calc::|log|) #'log)
(setf (symbol-function 'calc::|log10|) #'(lambda (x) (log (float x 0d0) 10d0)))
(setf (symbol-function 'calc::|sqrt|) #'sqrt)
(setf (symbol-function 'calc::|isqrt|) #'isqrt)
(setf (symbol-function 'calc::|abs|) #'abs)
(setf (symbol-function 'calc::|sin|) #'sin)
(setf (symbol-function 'calc::|cos|) #'cos)
(setf (symbol-function 'calc::|tan|) #'tan)
(setf (symbol-function 'calc::|asin|) #'asin)
(setf (symbol-function 'calc::|acos|) #'acos)
(setf (symbol-function 'calc::|atan|) #'atan)
(setf (symbol-function 'calc::|atan2|) #'(lambda (y x) (atan y x)))
(setf (symbol-function 'calc::|sinh|) #'sinh)
(setf (symbol-function 'calc::|cosh|) #'cosh)
(setf (symbol-function 'calc::|tanh|) #'tanh)
(setf (symbol-function 'calc::|asinh|) #'asinh)
(setf (symbol-function 'calc::|acosh|) #'acosh)
(setf (symbol-function 'calc::|atanh|) #'atanh)
(setf (symbol-function 'calc::|floor|) #'floor)
(setf (symbol-function 'calc::|ceil|) #'ceiling)
(setf (symbol-function 'calc::|ceiling|) #'ceiling)
(setf (symbol-function 'calc::|trunc|) #'truncate)
(setf (symbol-function 'calc::|truncate|) #'truncate)
(setf (symbol-function 'calc::|round|) #'round)
(setf (symbol-function 'calc::|rem|) #'rem)
(setf (symbol-function 'calc::|mod|) #'mod)
(setf (symbol-function 'calc::|ffloor|) #'ffloor)
(setf (symbol-function 'calc::|fceiling|) #'fceiling)
(setf (symbol-function 'calc::|fceil|) #'fceiling)
(setf (symbol-function 'calc::|ftruncate|) #'ftruncate)
(setf (symbol-function 'calc::|ftrunc|) #'ftruncate)
(setf (symbol-function 'calc::|fround|) #'fround)
(setf (symbol-function 'calc::|float|) #'(lambda (x) (float x 0d0)))
(setf (symbol-function 'calc::|int|) #'truncate)
(setf (symbol-function 'calc::|ratio|) #'rationalize)
(setf (symbol-function 'calc::|complex|) #'complex)
(setf (symbol-function 'calc::|phase|) #'phase)
(setf (symbol-function 'calc::|cis|) #'cis)
(setf (symbol-function 'calc::|min|) #'min)
(setf (symbol-function 'calc::|max|) #'max)
(setf (symbol-function 'calc::|conjugate|) #'conjugate)
(setf (symbol-function 'calc::|signum|) #'signum)
(setf (symbol-function 'calc::|sign|) #'signum)
(setf (symbol-function 'calc::|realpart|) #'realpart)
(setf (symbol-function 'calc::|real|) #'realpart)
(setf (symbol-function 'calc::|imagpart|) #'imagpart)
(setf (symbol-function 'calc::|imag|) #'imagpart)
(setf (symbol-function 'calc::|numerator|) #'numerator)
(setf (symbol-function 'calc::|num|) #'numerator)
(setf (symbol-function 'calc::|denominator|) #'denominator)
(setf (symbol-function 'calc::|den|) #'denominator)
(setf (symbol-function 'calc::|ash|) #'ash)
(setf (symbol-function 'calc::|shift|) #'ash)
(setf (symbol-function 'calc::|random|) #'random)
(setf (symbol-function 'calc::|rand|) #'random)


;;;;
(defun calc::|fibonacci| (n)
  (check-type n integer)
  (cond ((not (plusp n))
	 (error 'range-error :datum n))
	((= n 1) 1)
	((= n 2) 1)
	(t
	 (do* ((f1 1 f2)
	       (f2 1 f)
	       (f (+ f1 f2) (+ f1 f2))
	       (x 3 (+ x 1)))
	     ((>= x n) f)))))

(setf (symbol-function 'calc::|fib|) #'calc::|fibonacci|)

(defun calc-fact (n)
  (check-type n integer)
  (when (minusp n)
    (error 'range-error :datum n))
  (do* ((i 1 (+ i 1))
	(r 1 (* i r)))
      ((>= i n) r)))

(setf (symbol-function 'calc::|fact|) #'calc-fact)


(defmacro calc::|sigma| (var from to sexp)
	  `(calc::with-unit-calculation
	       (do* ((,var ,from (1+ ,var))
		     (#1=#:sum (calc::zero-of-unit ,sexp))) ;次元が必要
		    ((> ,var ,to) #1#)
		 (setq #1#  (+ #1# ,sexp)))))
(defmacro calc::|sum| (var from to sexp)
	  `(calc::|sigma| ,var ,from ,to ,sexp))

(defmacro calc::|product| (var from to sexp)
	  `(calc::with-unit-calculation
	       (do* ((,var ,from (1+ ,var))
		     (#1=#:product 1))
		    ((> ,var ,to) #1#)
		 (setq #1#  (* #1# ,sexp)))))


(setf (symbol-function 'calc-lshift) #'ash)
(defun calc-rshift (x y) (ash x (- y)))

;;

(defun calc-token-string (x)
  (or (and (symbolp x) (get x 'calc-operator)) x ""))


(defun calc-get-token (c s fn)
  (with-output-to-string (so)
    (when c
      (write-char c so))
    (loop
      (setq c (read-char s nil nil))
      (unless c
	(return))
      (unless (funcall fn c)
	(unread-char c s)
	(return))
      (write-char c so))))


(defun calc-next-token (s)
  (loop
    (let ((c (read-char s nil nil)))
      (cond ((null c)
	     (return nil))
	    ((let ((tok (assoc c *calc-token-alist*)))
	       (when tok
		 (return (loop
			   (let ((al (cdr tok)))
			     (unless (listp al)
			       (return al))
			     (setq c (read-char s nil nil))
			     (setq tok (assoc c al))
			     (unless tok
			       (setq tok (assoc nil al))
			       (cond ((null tok)
				      (if c
					  (error "Syntax error: ~c" c)
					(error "Unexpected EOL")))
				     (c
				      (unread-char c s))))))))))
	    ((or (digit-char-p c)
		 (char= c #\.))
	     (when (char= c #\0)
	       (let ((c (read-char s nil nil)))
		 (unless c
		   (return 0))
		 (let ((base (cond ((char-equal c #\x) 16)
				   ((char-equal c #\b) 2)
				   ((digit-char-p c 8)
				    (unread-char c s)
				    8)
				   (t
				    (unread-char c s)
				    nil))))
		   (when base
		     (return (parse-integer
			      (calc-get-token nil s
					      #'(lambda (c)
						  (digit-char-p c base)))
			      :radix base))))))
	     (let ((state 'leadnum)
		   (invalid-p nil))
	       (flet ((flonum (c)
			(setq state
			      (cond ((digit-char-p c)
				     (case state
				       ((leadnum trailnum expnum) state)
				       (dot 'trailnum)
				       ((expchar expsign) 'expnum)))
				    ((char= c #\.)
				     (case state
				       (leadnum 'dot)
				       (t (setq invalid-p state))))
				    ((or (char-equal c #\e)
					 (char-equal c #\d))
				     (case state
				       ((leadnum dot trailnum) 'expchar)
				       ((expchar expsign)
					(setq invalid-p state))
				       (expnum nil)))
				    ((or (char= c #\+)
					 (char= c #\-))
				     (case state
				       ((leadnum dot trailnum expnum) nil)
				       (expchar 'expsign)
				       (t (setq invalid-p state))))
				    (t
				     (case state
				       ((leadnum dot trailnum expnum) nil)
				       (t (setq invalid-p state))))))))
		 (flonum c)
		 (let* ((tok (read-from-string (calc-get-token c s #'flonum))))
		   (when (or invalid-p (not (numberp tok)))
		     (error "Invalid format: ~a" tok))
		   (return tok)))))
;	    ((find c " \t\n\r\f"))
	    ((find c (list #\space #\tab #\newline #\return #\page)))
	    ((or (alpha-char-p c)
		 (char= c #\_))
	     (return (intern (calc-get-token c s #'(lambda (c)
	     					     (or (alphanumericp c)
	     						 (char= c #\_))))

			      *calc-package*)))

	    (t
	     (error "Syntax error: ~c" c))))))

(defmacro calc-assoc-left (name next operators)
  `(defun ,name ()
     (let ((left (,next)))
       (while (,(if (listp operators) 'member 'eq) (car *calc-token*) ',operators)
	 (setq left `(,(pop *calc-token*) ,left ,(,next))))
       left)))

(defmacro calc-assoc-right (name next operators)
  `(defun ,name ()
     (let ((left (,next)))
       (if (,(if (listp operators) 'member 'eq) (car *calc-token*) ',operators)
	   `(,(pop *calc-token*) ,left ,(,name))
	 left))))


(defun calc-expr-list ()
  (when (eq (car *calc-token*) 'rparen)
    (pop *calc-token*)
    (return-from calc-expr-list nil))
  (let ((list nil))
    (loop
      (push (calc-expr) list)
      (cond ((eq (car *calc-token*) 'comma)
	     (pop *calc-token*))
	    ((eq (car *calc-token*) 'rparen)
	     (pop *calc-token*)
	     (return (nreverse list)))
	    (t
	     (error ") expected."))))))

(defun calc-primary-expr ()
  (let ((tok (car *calc-token*)))
    (cond ((numberp tok)
	   (pop *calc-token*)
	   tok)
	  ((eq (car *calc-token*) 'lparen)
	   (pop *calc-token*)
	   (prog1
	       (calc-expr)
	     (unless (eq (pop *calc-token*) 'rparen)
	       (error ") expected."))))
	  ((and (symbolp tok)
		(eq (symbol-package tok) *calc-package*))
	   (pop *calc-token*)
	   (cond ((eq (car *calc-token*) 'lparen)
		  (pop *calc-token*)
		  `(,tok ,@(calc-expr-list)))
		 (t tok)))
	  (t (error "Primary expected: ~A" (calc-token-string tok))))))

(defun calc-postfix-expr ()
  (let ((expr (calc-primary-expr)))
    (while (eq (car *calc-token*) '!)
      (pop *calc-token*)
      (setq expr `(calc-fact ,expr)))
    expr))

(defun calc-unary-expr ()
  (if (member (car *calc-token*) '(+ - lognot))
      `(,(pop *calc-token*) ,(calc-unary-expr))
    (calc-postfix-expr)))

(calc-assoc-right calc-power-expr calc-unary-expr expt)

(defun calc-multiplicative-expr ()
  (let ((left (calc-power-expr)))
    (loop
      (let ((ope (car *calc-token*)))
	(cond ((null ope)
	       (return left))
	      ((member ope '(* / rem truncate))
	       (setq left `(,(pop *calc-token*) ,left ,(calc-power-expr))))
	      ((or (and (eq ope 'lparen)
			(not (symbolp left)))
		   (not (and (symbolp ope)
			     (get ope 'calc-operator))))
	       (setq left `(* ,left ,(calc-power-expr))))
	      (t
	       (return left)))))))

(calc-assoc-left calc-additive-expr calc-multiplicative-expr (+ -))
(calc-assoc-left calc-shift-expr calc-additive-expr (calc-lshift calc-rshift))
(calc-assoc-left calc-logand-expr calc-shift-expr logand)
(calc-assoc-left calc-logxor-expr calc-logand-expr logxor)
(calc-assoc-left calc-logior-expr calc-logxor-expr logior)

(defun calc-assign-expr ()
  (let ((left (calc-logior-expr)))
    (cond ((eq (car *calc-token*) '=)
	   (pop *calc-token*)
	   (let ((right (calc-assign-expr)))
	     (cond ((symbolp left)
		    `(setq ,left ,right))
		   ((and (listp left)
			 (symbolp (car left))
			 (eq (symbol-package (car left)) *calc-package*))
		    (dolist (x (cdr left))
		      (or (symbolp x)
			  (error "Syntax error: ~a" (calc-token-string x))))
;		    (princ right)
		    (setf (symbol-function (car left)) 
			  (eval `(lambda ,(cdr left) ,(calc::with-unit-calculation% right))))
		    ''function)
		   (t
		    (error "left operand must be l-value")))))
	  (t left))))

(defun calc-expr ()
  (calc-assign-expr))

(defun calc-options (args)
  (cond (args
	 (while args
	   (let ((var (pop args))
		 (val (when (eq (car args) '=)
			(pop args)
			(pop args))))
	     (cond ((eq var 'calc::|bits|)
		    (setq *calc-bits*
			  (cond ((member val '(16 32 64)) val)
				((or (null val) (eq val 'calc::|unlimit|)) nil)
				(t (error "Invalid bits: ~A" (calc-token-string val))))))
		   ((eq var 'calc::|radix|)
		    (setq *calc-radix*
			  (cond ((member val '(2 8 10 16)) val)
				((eq val 'calc::|bin|) 2)
				((eq val 'calc::|oct|) 8)
				((eq val 'calc::|dec|) 10)
				((eq val 'calc::|hex|) 16)
				(t (error "Invalid radix: ~A" (calc-token-string val))))))
		   ((eq var 'calc::|signed|)
		    (setq *calc-unsigned* nil))
		   ((eq var 'calc::|unsigned|)
		    (setq *calc-unsigned* t))
		   ((eq var 'calc::|ratio|)
		    (setq *calc-ratio*
			  (cond ((eq val 'calc::|int|) 'calc::|int|)
				((eq val 'calc::|ratio|) 'calc::|ratio|)
				((eq val 'calc::|float|) 'calc::|float|)
				((null val) *calc-ratio*)
				(t (error "Invalid ratio: ~A" (calc-token-string val))))))
		   (t
		    (error "Unknown option: ~A" (calc-token-string var))))))
	 nil)
	(t
	 (format nil "bits=~A radix=~A ~A ratio=~A"
		 (if *calc-bits* *calc-bits* '|unlimit|)
		 *calc-radix* (if *calc-unsigned* '|unsigned| '|signed|) *calc-ratio*))))

(defun calc::infix-string (string)
  "中置記法の数式を文字列として与えるとS式に組み直して返す"
  (let ((*read-default-float-format* 'double-float)
	(*calc-token* (let ((token nil))
			(with-input-from-string (s string)
			  (do ((tok (calc-next-token s) (calc-next-token s)))
			      ((null tok) (nreverse token))
			    (push tok token))))))
    (if (eq (car *calc-token*) 'calc::|set|)
	(calc-options (cdr *calc-token*))
      (let ((expr (and *calc-token* (calc-expr))))
	(when *calc-token*
	  (error "unexpected ~A." (calc-token-string (car *calc-token*))))
	expr
	))))

(defun calc::calc-string (string)
  "中置記法の数式を文字列として与えるとその結果を返す"
  (eval (calc::with-unit-calculation% (calc::infix-string string))))

;; (defun calc::one-line-reader (stream ch n)
;;   (declare (ignore ch n))
;;   (let ((chars))
;;     (do ((curr (read-char stream)
;; 	       (read-char stream)))
;; 	((char= #\newline curr))
;;       (push curr chars))
;;     (coerce (nreverse chars) 'string)))

;; (defun calc::one-line-calculator (stream ch n)
;;   (declare (ignore ch n))
;;   (let ((chars))
;;     (do ((curr (read-char stream)
;; 	       (read-char stream)))
;; 	((char= #\newline curr))
;;       (push curr chars))
;;     (calc::calc-string (coerce (nreverse chars) 'string))))

;; (set-dispatch-macro-character #\# #\$ #'calc::one-line-calculator)

(defun calc::start-calc-quiet ()
  "電卓replもどき"
  (if (ignore-errors
	(princ "$ ")
	(do ((str (read-line) (read-line)))
	    ((member str '("bye" "quit") :test #'equal) t)
	  (let ((result (calc::calc-string str)))
	    (when result
	      (princ result)
	      (princ #\newline)))
	  (princ "$ ")))
      (quit)
      (calc::start-calc-quiet)))

(defun calc::start-calc-with-eval ()
  "式を入力すると答えを表示するループ。
\"eval\"を受け取ると次に入力されたS式を評価する。
\"bye\"、\"quit\"、EOFを受け取ると終了。"
  (labels ((read-line_ () (read-line *standard-input* nil 'eof)))
    (princ "$ ")
    (do ((str (read-line_) (read-line_)))
	((or (member str '("bye" "quit") :test #'equal)
	     (eq str 'eof))
	 t)
      (if
       (equal str "eval")
       (progn
	 (princ "eval: ")
	 (princ (eval (read-from-string (read-line_))))
	 (princ #\newline))
       (let ((result))
	 (handler-case 
	     (let ((*read-default-float-format* 'double-float))
	       (setf result (calc::calc-string str)))
	   (error (c)
	     (setf result (let ((*package* *calc-package*)) c))))
	 (when result
	   (calc::print-unit result *calc-bits* *calc-radix* *calc-unsigned* *calc-ratio*)
	   (princ #\newline))))
      (princ "$ "))))


(calc::start-calc-with-eval)

units.l

(in-package :calc)

;(defvar *fundamental-units* '(|kg| |m| |s| |C| |K|))
(defvar *fundamental-units* '(|kg| |m| |s| |C| |K| |mol|))


;;
(defun mkstr (&rest args)
  (with-output-to-string (s)
    (dolist (a args)(princ a s))))
(defun symb (&rest args)
  (values (intern (apply #'mkstr args))))

(defvar *si-prefix-alist*
  `((|y| . ,(expt 10 -24))
    (|z| . ,(expt 10 -21))
    (|a| . ,(expt 10 -18))
    (|f| . ,(expt 10 -15))
    (|p| . ,(expt 10 -12))
    (|n| . ,(expt 10 -9))
    (|u| . ,(expt 10 -6))
    (|m| . ,(expt 10 -3))
    (|c| . ,(expt 10 -2))
    (|d| . ,(expt 10 -1))
    (|K| . ,(expt 10 3))
    (|k| . ,(expt 10 3))
    (|h| . ,(expt 10 2))
    (|M| . ,(expt 10 6))
    (|G| . ,(expt 10 9))
    (|T| . ,(expt 10 12))
    (|P| . ,(expt 10 15))
    (|E| . ,(expt 10 18))
    (|Z| . ,(expt 10 21))
    (|Y| . ,(expt 10 24))))		;da とか使わないよね?

(defstruct (unit (:constructor make-unit))
  "単位や定数を格納する。値と次元を表すリストとからなる。"
;  name
  (value 1)
  (dimensions (zero-dimension-list)))

(defun zero-dimension-list ()
  "無次元であることを表すリストを返す"
  (make-list (length *fundamental-units*) :initial-element 0))

(defun dimension (fundamental-unit u)
  "uのfundamental-unitで指定された次元がいくつかを返す。fundamental-unitは*fundamental-unit*の要素である必要がある。"
  (nth (position fundamental-unit *fundamental-units*) (unit-dimensions u)))

(dolist (u *fundamental-units*)
  (setf (symbol-value u)
;	(make-unit :name (symbol-name u))
	(make-unit)
	)
  (setf (nth (position u *fundamental-units*) (unit-dimensions (symbol-value u))) 1))

(defun zero-dimension-p (u)
  "与えられた単位の構造体が無次元量か"
  (equal (unit-dimensions u) (make-list (length *fundamental-units*) :initial-element 0)))

(defun zero-dimensional-unit-to-number (u)
  "単位の構造体でかつ次元が0なら単なる数値にして返す。そうでなければそのまま返す。"
    (if (and (unit-p u) (zero-dimension-p u))
	(unit-value u)
	u))

(defun multiply-unit (&rest args)
  "*を単位構造体が扱えるようにしたもの。結果が無次元になったら数値として返す。"
  (let ((units   (remove-if-not #'unit-p args))
	(numbers (remove-if     #'unit-p args)))
    (if (null units)
	(apply #'* numbers)
	(zero-dimensional-unit-to-number 
	 (make-unit :value
		    (* (apply #'* numbers) (apply  #'* (mapcar #'unit-value units)))
		    :dimensions 
		    (apply #'mapcar #'+ (mapcar #'unit-dimensions units)))))))

(defun invert-unit (arg)
  "逆数"
  (if (unit-p arg)
      (make-unit :value (/ (unit-value arg))
		 :dimensions (mapcar #'- (unit-dimensions arg)))
      (/ arg)))

(defun divide-unit (unit1 &rest units)
  "数値に対しての/と同様"
  (if (null units)
      (invert-unit unit1)
      (apply #'multiply-unit unit1 (mapcar #'invert-unit units))))

(defun add-unit (&rest args)
  "+と同じ。次元が合わないとエラーを吐く"
  (if (remove-if #'unit-p args)
      (apply #'+ args)
      (if (null (cdr args))
	  (car args)
	  (if (equal (unit-dimensions (car args)) (unit-dimensions (cadr args)))
	      (let
		  ((next (apply #'add-unit (cdr args))))
		(make-unit :value 
			   (+ (unit-value (car args)) (unit-value next))
			   :dimensions
			   (unit-dimensions (car args))))
	      (error "dimension does't match")))))

(defun subtract-unit (&rest args)
  "-と同じ。次元が合わないとエラーを吐く"
  (if (numberp (car args))
      (apply #'- args)
      (if (null (cdr args))
	  (multiply-unit -1 (car args))
	  (if (equal (unit-dimensions (car args)) (unit-dimensions (cadr args)))
	      (let
		  ((rest (apply #'add-unit (cdr args))))
		(make-unit :value
			   (- (unit-value (car args)) (unit-value rest))
			   :dimensions
			   (unit-dimensions (car args))))
	      (error "dimension doesn't match")))))

(defun expt-unit (base power)
  "baseのpower乗。powerは数値でないといけない"
  (if (unit-p base)
      (make-unit :value 
		 (expt (unit-value base) power)
		 :dimensions 
		 (mapcar #'(lambda (dim) (* dim power)) (unit-dimensions base)))
      (expt base power)))

(defun sqrt-unit (u)
  "平方根"
  (expt-unit u 1/2))

(defun isqrt-unit (u)
  "平方根の整数部分。引数は整数(単位を持っていても良い)"
  (if (unit-p u)
      (make-unit
       :value (isqrt (unit-value u))
       :dimensions (mapcar #'(lambda (dim) (* dim 1/2)) (unit-dimensions u))
      (isqrt u))))

(defvar *unit-calculation-function-alist*
  `((+ . add-unit)
    (- . subtract-unit)
    (* . multiply-unit)
    (/ . divide-unit)
    (|expt| . expt-unit)
    (|pow| . expt-unit)
    (|sqrt| . sqrt-unit)
    (isqrt . isqrt-unit)
    (|isqrt| . isqrt-unit)
    (expt . expt-unit) 	;下の定数定義部分で使おうとするとこんなことになる
    (pow . expt-unit)
    (sqrt . sqrt-unit)
  ))

(defmacro defun-unit-with-1-arg (func args)	
  "数値を1引数かつ数値だけに影響し次元を変えない関数を定義し*unit-calculation-function-alist*に登録する"
  (let ((func_ (symb (string-downcase (symbol-name func)))))
  `(progn
     (push (cons ',func_ ',(symb (symbol-name func) '-unit)) ;どうせ小文字
	   *unit-calculation-function-alist*)
     (defun ,(symb (symbol-name func) '-unit) ,args
       (if (unit-p ,@args)
	 (make-unit
	  :value (,func_ (unit-value ,@args))
	  :dimensions (unit-dimensions ,@args))
	 (,func_ ,@args))))))

(defmacro defun-unit-with-1-arg-nodim (func args)	
  "1引数かつ数値だけに影響し無次元量を返す関数を定義し*unit-calculation-function-alist*に登録する"
  (let ((func_ (symb (string-downcase (symbol-name func)))))
    `(progn
       (push (cons ',func_ ',(symb (symbol-name func) '-unit)) ;どうせ小文字
	     *unit-calculation-function-alist*)
       (defun ,(symb (symbol-name func) '-unit) ,args
	 (,func_ (if (unit-p ,@args) (unit-value ,@args) args))))))

(defun-unit-with-1-arg abs (x))
(defun-unit-with-1-arg floor (x))
(defun-unit-with-1-arg ceil (x))
(defun-unit-with-1-arg ceiling (x))
(defun-unit-with-1-arg trunc (x))
(defun-unit-with-1-arg truncate (x))
(defun-unit-with-1-arg round (x))
(defun-unit-with-1-arg conjugate (x))
(defun-unit-with-1-arg ffloor (x))
(defun-unit-with-1-arg fceil (x))
(defun-unit-with-1-arg fceiling (x))
(defun-unit-with-1-arg ftrunc (x))
(defun-unit-with-1-arg ftruncate (x))
(defun-unit-with-1-arg fround (x))
(defun-unit-with-1-arg float (x))
(defun-unit-with-1-arg int (x))
(defun-unit-with-1-arg ration (x))
(defun-unit-with-1-arg complex (x))
(defun-unit-with-1-arg int (x))
(defun-unit-with-1-arg realpart (complex))
(defun-unit-with-1-arg real (complex))
(defun-unit-with-1-arg imagpart (complex))
(defun-unit-with-1-arg imag (complex))

(defun-unit-with-1-arg-nodim phase (complex))
(defun-unit-with-1-arg-nodim signum (x))

(defun-unit-with-1-arg numerator (ratio))	;次元は分子側に押し付けてみた
(defun-unit-with-1-arg num (ratio))
(defun-unit-with-1-arg-nodim denominator (ratio))
(defun-unit-with-1-arg-nodim den (ratio))

(defun si-unit-to-unit (s)
  "未束縛のシンボルが 'si接頭辞+定義済みの単位' の形の名前なら元の単位での値に直して返す"
  (if (not (or (numberp s) (unit-p s)))
      (if (> (length (symbol-name s)) 1)
	  (let
	      ((prefix (intern (subseq (symbol-name s) 0 1) :calc))
	       (unit (intern (subseq (symbol-name s) 1) :calc)))
	    (if (and (not (boundp s))
		     (boundp unit)
		     (assoc prefix *si-prefix-alist*))
		(multiply-unit
		 (cdr (assoc prefix *si-prefix-alist*))
		 (eval unit))
		s))
	  s)
      s))

(defun print-unit (arg &optional (bits) (radix 10) (unsigned) (ratio '|ratio|))
  "単位構造体or数値を与えるとストリームに印字"
  (flet ((prnt (x)
	   (let ((result x))
	     (when (rationalp result)
	       (cond ((eq ratio '|int|)
		      (setf result (truncate result)))
		     ((eq ratio '|float|)
		      (setf result (float result 1d0)))))
	     (when (and (integerp result) bits)
	       (setq result (logand result (- (ash 1 bits) 1)))
	       (when (and (null unsigned)
			  (logbitp (- bits 1) result))
		 (setq result (- result (ash 1 bits)))))
	     (if (integerp result)
		 (cond ((eql radix 2)
			(format t "0b~b" result))
		       ((eql radix 8)
			(format t "0~o" result))
		       ((eql radix 16)
			(format t "0x~x" result))
		       (t
			(format t "~d" result)))
		 (format t "~a" result)))))
    (if (unit-p arg)
	(progn
	  (prnt (unit-value arg))
	  (dolist (u *fundamental-units*)
	    (let
		((deg (nth (position u *fundamental-units*) (unit-dimensions arg))))
	      (unless (zerop deg)
	      (if (= 1 deg)
		  (format t " ~A" (symbol-name u))
		  (format t " ~A^~D" (symbol-name u) deg))))))
	(prnt arg))))

(defun with-unit-calculation% (tree)
  "数値に対して書くようなS式( (* 3 K) とか)を与えると単位を持った量を扱うよう変換する"
  (if tree
    (if (listp tree)
	(cons
	 (with-unit-calculation% (car tree))
	 (with-unit-calculation% (cdr tree)))
	(if
	 (assoc tree *unit-calculation-function-alist*)
	 (cdr (assoc tree *unit-calculation-function-alist*))
	 (if (symbolp tree)
	     (si-unit-to-unit tree)
	     tree)))))

(defun zero-of-unit (u)
  (if (numberp u)
      0
      (make-unit :value 0 :dimensions (unit-dimensions u))))

(defmacro with-unit-calculation (sexp)
  `(eval (with-unit-calculation% ',sexp)))



(defmacro define-unit (name definition)
  "単位を定義"
  `(defvar ,name (with-unit-calculation ,definition)))

(defmacro define-const (name value &optional (unit 1))
  "定数を定義。"
  `(defvar ,name (with-unit-calculation (* ,value ,unit))))


;;;;

;; mathematical constants
(define-const |pi| pi)
(define-const |Pi| pi)
(define-const |E| (exp 1.0d0))

; combined units and constants

;; aux.
;(define-unit |mol| 1)
(define-unit |sec| 1)
(define-unit |meter| 1)
(define-unit |rad| 1)
(define-unit |deg| (/ pi 180))

;; combined units
(define-unit |N| (* |kg| (/ |m| |s| |s|)))
(define-unit |J| (* N |m|))
(define-unit |A| (/ C |s|))
(define-unit |V| (/ J C))
(define-unit |H| (* V (/ A) |s|))
(define-unit |F| (/ C V))
(define-unit |Ohm| (/ V A))
(define-unit |W| (/ J |s|))
(define-unit |Weber| (* N |m| (/ A)))
(define-unit |Wb| |Weber|)
(define-unit |Tesla| (/ |Weber| |m| |m|))
(define-unit |Pa| (/ N |m| |m|))

;; scientific constants in MKSA
(define-const |e|    1.60217646d-19 C)
(define-const |c|    299792458 (/ |m| |s|))
(define-const |h|    6.626068d-34 (* J |s|))
(define-const |hbar| (/ |h| (* 2 pi)))
(define-const |G|    6.67300d-11 (* (expt |m| 3) (expt |kg| -1) (expt |s| -2)))
(define-const |g_|    9.80665 (/ |m| |s| |s|))
(define-const |kB|   1.3806503d-23 (/ J K))
(define-const |NA|   6.0221415d23 (/ |mol|))

(define-const |e0|   8.85418782d-12 (/ F |m|))
(define-const |u0|   (* 4 pi 1d-7) (/ H |m|))

(define-const |m_e|  9.1093826d-31 |kg|)
(define-const |m_p|  1.672621637d-27 |kg|)
(define-const |m_n|  1.674927211d-27 |kg|)

(define-unit |u| (* 1.66053886d-27 |kg|)) ;atomic mass unit

;; misc.
(define-unit |eV| (* |e| V))
(define-unit |erg| (* 1d-7 J))

(define-unit |g| (* 1/1000 |kg|))
(define-unit |gram| (* 1/1000 |kg|))
(define-unit |ton|  (* 1000 |kg|))

(define-unit |Hz| (/ |s|))

(define-unit |Bq| (/ |s|))

;; (define-unit |km| (* 1000 |m|)) ;このへんがなくても動く
;; (define-unit |cm| (* 1/100 |m|))
;; (define-unit |mm| (* 1/1000 |m|))
;; (define-unit |nm| (* 1/1000000 |m|))
(define-unit |angstrom| (* (expt 10 -10) |m|))

(define-unit |litre| (* 1/1000 (expt |m| 3)))
(define-unit |l| |litre|)
(define-unit |cc| (expt |cm| 3))

;(defun yard-pound () 
  (define-unit |inch| (* 0.0254 |m|))
  (define-unit |feet| (* 12 |inch|))
  ;;(define-unit |ft| |feet|)
  (define-unit |yard| (* 36 |inch|))
  ;;(define-unit |yd| |yard|)
  (define-unit |mile| (* 1760 |yard|))
  ;;(define-unit |ml| |mile|)
  (define-unit |league| (* 3 |mile|))
;  'yard-pound
;  )

(define-unit |barn| (* 100 |fm| |fm|))

(define-unit |Gauss| (* 1/10000 |Tesla|))

(define-unit |atm| (* 101325 |Pa|))
(define-unit |torr| (* 101325/760 |Pa|))
(define-unit |mmHg| (* 133.3224 |Pa|))

(define-unit |ly| (* 9.46073d15 |m|))
(define-unit |AU| (* 1.49598d11 |m|))
(define-unit |pc| (* 3.08568d16 |m|))

(define-unit |minute| (* 60 |s|))
(define-unit |hour| (* 60 |minute|))
(define-unit |day| (* 60 |hour|))
(define-unit |year| (* 365.25 |day|))


(define-const |Celcius| 273.15 K)

;; plack units
(define-unit |t_pl| (sqrt (/ (* |hbar| G) (expt |c| 5))))
(define-unit |l_pl| (* |c| |t_pl|))
(define-unit |m_pl| (sqrt (/ (* |hbar| |c|) G)))
;(define-unit |q_pl| (* ))
;(define-unit |T_pl| (* ))

*1:まあ元もcommon lisp(準拠)なので移植というべきかだいぶ微妙ではある

*2:kg,m,s,C,K,molだけは大元にしているので一番上で定義してある