share/contrib/lll.lisp

   1 ; file lll.lisp
   2 ; Maxima Interface
   3
   4 (defun maxima::$latticereduce (L) ; L is maxima list of maxima lists
   5   (if (or
   6        (not (listp L))
   7        (< (length L) 2)
   8        (not (equal (apply 'max (mapcar 'length (cdr L)))
   9                    (apply 'min (mapcar 'length (cdr L))) ))  )
  10       (maxima::merror "Latticereduce needs a list lists as input"))
  11   (let* ((vs (mapcar 'cdr (cdr L)))
  12          (n (length vs))
  13          (v (make-array (list n))))
  14     (loop for i from 0 to (- n 1) do
  15           (setf (aref v i) (make-array (list n) :initial-contents (nth i vs))))
  16     (setf v (lll v))
  17     (cons '(maxima::mlist)
  18           (mapcar #'(lambda (u) (cons '(maxima::mlist) (coerce u 'list)))
  19                   (coerce v 'list)))
  20 ))
  21
  22
  23 (defun maxima::$integerrelations (L)
  24   (let ((res (integerrelations (cdr L))))
  25     (cons '(maxima::mlist) (car res))
  26     ))
  27
  28 (defun maxima::$floatrelations (L)
  29   (let ((res (floatrelations (cdr L))))
  30     (cons '(maxima::mlist) (car res))
  31     ))
  32
  33 (defun maxima::$recognize (x)
  34   (convert2AlgNum x))
  35
  36 ; linear algebra helpers
  37
  38 (defun norm2 (a)
  39   (declare (type (vector number *) a))
  40   (loop for i from 0 to (- (length a) 1)  sum
  41         (* (aref a i) (aref a i) )))
  42
  43
  44 (defun vadd (a b)
  45   (declare (type (vector number *) a))
  46   (declare (type (vector number *) b))
  47   (if (not (= (length a) (length b)))
  48       (error "vadd: length must coincide")
  49     (let ((res (make-array (list (length a)) ))
  50         (n (length a)))
  51       (loop for i from 0 to (- n 1) do
  52             (setf (aref res i) (+ (aref a i) (aref b i))))
  53       res)))
  54
  55 (defun vsub (a b)
  56   (declare (type (vector number *) a))
  57   (declare (type (vector number *) b))
  58   (if (not (= (length a) (length b)))
  59       (error "vadd: length must coincide")
  60     (let ((res (make-array (list (length a)) ))
  61         (n (length a)))
  62       (loop for i from 0 to (- n 1) do
  63             (setf (aref res i) (- (aref a i) (aref b i))))
  64       res)))
  65
  66 (defun smul (s a)
  67   (declare (type (vector number *) a))
  68   (declare (type (vector number *) b))
  69   (let ((res (make-array (list (length a)) ))
  70         (n (length a)))
  71       (loop for i from 0 to (- n 1) do
  72             (setf (aref res i) (* s (aref a i))))
  73       res))
  74
  75
  76 (defun dot (a b)
  77   (declare (type (vector number *) a))
  78   (declare (type (vector number *) b))
  79   (if (not (= (length a) (length b)))
  80       (error "vadd: length must coincide")
  81     (loop for i from 0 to (- (length a) 1) sum
  82           (* (aref a i) (aref b i)))))
  83
  84 (defvar *gsbasis* nil)
  85
  86
  87 (defun gso (f mu) ; returns Gram-Schmidt OGS, updates mu
  88   (declare (type vector a))
  89   (declare (type (array number (* *) mu)))
  90   (let* ((n (length f))
  91          (g (make-array (list n) )))
  92     (loop for j from 0 to (- n 1) do (setf (aref g j) (aref f j)))
  93     (loop for j from 0 to (- n 1) do
  94           (loop for k from 0 to (- j 1) do
  95                 (if (> (norm2 (aref g k)) 0)
  96                     (progn
  97                       (setf (aref mu j k)
  98                             (/ (dot (aref f j) (aref g k))
  99                                (dot (aref g k) (aref g k))))
 100                       (setf (aref g j) (vsub (aref g j) (smul (aref mu j k) (aref g k))))
 101                       ))
 102                 ))
 103     g
 104     ))
 105
 106 (defun lll (f) ; f is vetor of n vectors
 107   (declare (type (vector vector *) f))
 108   (let* ((n (length f))
 109         (g (make-array (list n) ))
 110         (gg nil)
 111         (i 1)
 112         (temp nil)
 113         (mu (make-array (list n n) :initial-element 0)))
 114     (loop for j from 0 to (- n 1) do
 115           (setf (aref g j) (aref f j))
 116           (setf (aref mu j j) 1))
 117     (setf gg (gso g mu))
 118     (do () ((>= i n) )
 119       (loop for j from (- i 1) downto 0 do
 120             (setf (aref g i)
 121                   (vsub (aref g i) (smul (nint (aref mu i j)) (aref g j))))
 122             (setf gg (gso g mu))
 123             )
 124       (if (and (> i 0) (> (norm2 (aref gg (- i 1))) (* 2 (norm2 (aref gg i)))))
 125           (progn
 126             (setf temp (aref g (- i 1)))
 127             (setf (aref g (- i 1)) (aref g i))
 128             (setf (aref g i) temp)
 129             (setf gg (gso g mu))
 130             (setf i (- i 1))
 131             )
 132         (setf i (+ i 1)))
 133       )
 134     g
 135     )
 136   )
 137
 138 (defun nint (q)
 139        (car (list (floor (+ q (/ 1 2))))))
 140
 141 (defun test1 ()
 142   (lll #( #(12 2) #(13 4))))
 143
 144 (defun test2 ()
 145   (lll #( #(12 2) #(1 2))))
 146
 147 (defun test3 ()
 148   (lll #( #(1 2) #(12 2) )))
 149
 150 (defun test4 ()
 151   (lll #( #(1 2) #(9 -4))))
 152
 153 (defun integerRelations (L) ; L ist list of integers
 154 ; returns (rel def) such that dotproduct def of rel and L is small
 155   (let* ((n (length L))
 156         (z nil) (g nil) (res nil) (optdef 0) (optres nil) (def 0)
 157         (f (make-array (list (+ n 1)))))
 158     (setf (aref f 0) (make-array (list (+ n 1)) :initial-element 0))
 159     (loop for i from 1 to n do
 160           (setf z  (make-array (list (+ n 1)) :initial-element 0))
 161           (setf (aref z 0) (nth (- i 1) L))
 162           (setf (aref z i) 1)
 163           (setf (aref f i) z))
 164     (setf g (lll f))
 165     (setf optdef most-POSITIVE-fixnum)
 166     (loop for k from 1 to n do
 167           (setf def 0)
 168           (setf res nil)
 169           (loop for i from n downto 1 do
 170           (setf res (cons (aref (aref g k) i) res))
 171           (setf def (+ def (* (nth (- i 1) L) (aref (aref g k) i))))
 172           )
 173           (if (< (abs def) (abs optdef))
 174               (progn
 175                 (setf optdef def)
 176                 (setf optres res)
 177                 ))
 178           )
 179     (list optres optdef)
 180     ))
 181
 182 (defun test5 ()
 183   (print "Should give ((3 -7 4) 0)")
 184   (integerRelations '(5707963267 4142135623 2967764890)))
 185
 186
 187 (defun nkomma (a)
 188   (- a (car (list (floor a)))))
 189
 190
 191 (defun FloatRelations (x &optional (digits 10)) ; x is a list of floats or doubles
 192   (let ((n (length x)) (digs nil) (i 0) (nkomma 0)
 193         (nkommaStellen0) (xx 0) (expo nil) (rel nil) (def 0))
 194     (setf nkommaStellen
 195           (lambda (x) ; x ist float
 196             (let ((s 0) (i 0))
 197               (do () ((< (abs (nkomma (* x (expt 10d0 s)))) (expt 1d-1 digits)))
 198                 (setf s (+ s 1)))
 199               s)))
 200   (setf digs (mapcar nkommaStellen x))
 201   (setf expo (apply 'max digs))
 202   (setf xx (mapcar #'(lambda (u) (car (list (floor (* u (expt 10d0 expo))))))  x))
 203   (setf rel (IntegerRelations xx))
 204   (list (car rel) (/ (cadr rel) (expt 10.0d0 expo)))
 205   ))
 206
 207
 208 (defun test6 ()
 209   (FloatRelations (list .5707963267d0 .4142135623d0 .2967764890d0)))
 210
 211
 212 (defun mkdif (a b) (list '(MAXIMA::MPLUS) a (list '(MAXIMA::MMINUS) b)) )
 213 (defun mksum (a b) (list '(MAXIMA::MPLUS) a b) )
 214 (defun mkprod (a b) (list '(MAXIMA::MTIMES) a b) )
 215 (defun mkexpt (a b) (list '(MAXIMA::MEXPT) a b))
 216 (defun mkeq (a b) (list '(MAXIMA::MEQUAL) a b))
 217
 218 (defun max2str (expr)
 219   (maxima::mfuncall 'maxima::$string expr)
 220   )
 221
 222 (defun convert2AlgNum (x &optional (digits 10)) ; x is float or double
 223   (let ((n 0) (nn 6) ;; nn: maximal search degree
 224         (xs nil) (res nil) (def 0) (sol nil)
 225         (opt nil) (mini 0)
 226         (var ($gensym "v")))
 227     (loop for n from 2 to nn do
 228           (if (not (null opt)) (return opt))
 229           (setf xs (loop for i from 0 to n collect (expt x i)))
 230           (setf res  (FloatRelations xs digits))
 231           (setf def (cadr res))
 232           (if (< (abs def) (expt 10d0 (- digits)))
 233               (progn
 234                 (setf pol 0)
 235                 (loop for i from 0 to n do
 236                       (setf pol (mksum pol (mkprod (nth i (car res)) (mkexpt var i)))))
 237                 (setf sol (maxima::meval (maxima::$solve (mkeq pol 0) var)))
 238                 (if (not (null (cdr sol))) ; a solution has been found
 239                     (progn (setq *sol* sol)
 240                       (setf sol (mapcar 'third (cdr sol)))
 241                       (setf sol
 242                             (mapcar #'(lambda (u)
 243                                         (let ((r (maxima::$realpart u))
 244                                               (i (maxima::$float (maxima::$imagpart u))))
 245                                           (if (or (equal i 0) (equal i 0.0) (equal i 0.0d0))
 246                                               (list u (abs (- x (maxima::$float u))))
 247                                               (list u most-POSITIVE-DOUBLE-FLOAT)
 248                                               )))
 249                             sol))
 250                       (setf mini most-POSITIVE-DOUBLE-FLOAT)
 251                       (loop for p in sol do (if (< (second p) mini) (progn (setf opt (car p)) (setq mini (second p)))))
 252                 ))))
 253           ); loop
 254     opt
 255     ) ; let
 256   ) ;defun
 257
 258 (defun test7 () ; fails
 259        (convert2AlgNum 1.414213562373d0))
 260
 261
 262 (defun test8 () ; works
 263        (convert2AlgNum (+ 2 (* 3 1.414213562373d0))))
 264
 265
 266 (defun test9 () :fails
 267        (convert2AlgNum (+ 1d0 (sqrt 2d0) (sqrt 3d0)) ))