Added math.scm to the repository, which contains bignum functions.
[picobit/chj.git] / math.scm
blob04636b388548230ae7b03d644ede1b50746e93de
1 ;------------------------------------------------------------------------------
3 (define-type num
4   lo16
5   shr16
8 (define min-fix -2)
9 (define max-fix  2)
11 (define fixnums (make-vector (+ (- max-fix min-fix) 1)))
13 (define fixnum
14   (lambda (n)
15     (vector-ref fixnums (- n min-fix))))
17 (let loop ((n min-fix))
18   (if (<= n max-fix)
19       (begin
20         (vector-set! fixnums (- n min-fix) (make-num (modulo n #x10000) #f))
21         (loop (+ n 1)))))
23 (let loop ((n min-fix))
24   (if (<= n max-fix)
25       (begin
26         (num-shr16-set!
27          (fixnum n)
28          (if (< n 0) (fixnum -1) (fixnum 0)))
29         (loop (+ n 1)))))
31 (define minus1 (fixnum -1))
32 (define zero   (fixnum 0))
33 (define one    (fixnum 1))
35 (define make-num
36   (let ((mn make-num))
37     (lambda (lo16 shr16)
38       (if (or (< lo16 0) (> lo16 #xffff))
39           (error "lo16 out of range"))
40       (mn lo16 shr16))))
42 (define check-norm
43   (lambda (name result)
44     (cond ((and (eq? (num-shr16 result) zero)
45                 (<= (num-lo16 result) max-fix))
46            (if (not (eq? result (fixnum (num-lo16 result))))
47                (pp (list name "nonneg number is not normalized" (dec result) (num-lo16 result)))))
48           ((and (eq? (num-shr16 result) minus1)
49                 (>= (num-lo16 result) (+ #x10000 min-fix)))
50            (if (not (eq? result (fixnum (- (num-lo16 result) #x10000))))
51                (pp (list name "neg number is not normalized" (dec result) (- (num-lo16 result) #x10000))))))))
53 (define check-normalized
54   (lambda (name f)
55     (lambda args
56       (let ((result (apply f args)))
57         (check-norm name result)
58         result))))
60 (define invert
61   (lambda (x)
62     (if (eq? x zero) minus1 zero)))
64 (define sign?
65   (lambda (x)
66     (check-norm 'sign?.x x)
67     (or (eq? x zero) (eq? x minus1))))
69 (define neg?
70   (lambda (x)
71     (check-norm 'neg?.x x)
72     (cond ((eq? x zero)   #f)
73           ((eq? x minus1) #t)
74           (else           (neg? (num-shr16 x))))))
76 (define lt?
77   (lambda (x y)
78     (check-norm 'lt?.x x)
79     (check-norm 'lt?.y y)
80     (< (cmp x y) 0)))
82 (define le?
83   (lambda (x y)
84     (check-norm 'le?.x x)
85     (check-norm 'le?.y y)
86     (<= (cmp x y) 0)))
88 (define cmp
89   (lambda (x y)
90     (check-norm 'cmp.x x)
91     (check-norm 'cmp.y y)
92     (let loop ((x x) (y y) (c 0))
93       (cond ((sign? x)
94              (cond ((eq? x y) c)
95                    ((neg? y)  1)
96                    (else      -1)))
97             ((sign? y)
98              (cond ;((eq? x y) c)
99                    ((neg? x)  -1)
100                    (else      1)))
101             (else
102              (loop (num-shr16 x)
103                    (num-shr16 y)
104                    (let ((x-lo16 (num-lo16 x))
105                          (y-lo16 (num-lo16 y)))
106                      (cond ((= x-lo16 y-lo16) c)
107                            ((< x-lo16 y-lo16) -1)
108                            (else              1)))))))))
110 (define int-len-nonneg
111   (lambda (x)
112     (check-norm 'int-len-nonneg.x x)
113     (let loop1 ((x x) (n 0))
114       (let ((shr16 (num-shr16 x)))
115         (if (eq? shr16 zero)
116             (let loop2 ((lo16 (num-lo16 x)) (n n))
117               (if (= lo16 0)
118                   n
119                   (loop2 (arithmetic-shift lo16 -1)
120                          (+ n 1))))
121             (loop1 shr16 (+ 16 n)))))))
123 (define shift-left
124   (check-normalized 'shift-left
125    (lambda (x n)
126      (check-norm 'shift-left.x x)
127      (if (< n 0) (error "xxxxxxxxxxx"))
128      (if (eq? x zero)
129          zero
130          (let loop1 ((x x) (n n))
131            (if (= (bitwise-and n 15) 0)
132                (let loop2 ((x x) (n n))
133                  (if (= n 0)
134                      x
135                      (loop2 (make-num #x0000 x)
136                             (- n 16))))
137                (loop1 (shl x)
138                       (- n 1))))))))
140 (define shl
141   (check-normalized 'shl
142    (lambda (x)
143      (check-norm 'shl.x x)
144      (let loop ((x x) (minus-c zero) (r #f))
145        (cond ((eq? x minus-c)
146               (norm r x))
147              (else
148               (let ((z (+ (* 2 (num-lo16 x)) (if (eq? minus-c zero) 0 1))))
149                 (loop (num-shr16 x)
150                       (if (>= z #x10000) minus1 zero)
151                       (make-num (modulo z #x10000) r)))))))))
153 (define shr
154   (check-normalized 'shr
155    (lambda (x)
156      (check-norm 'shr.x x)
157      (let loop ((x x) (r #f))
158        (cond ((eq? x zero)
159               (norm r x))
160              ((eq? x minus1)
161               (norm r x))
162              (else
163               (let ((z (+ (arithmetic-shift (num-lo16 x) -1)
164                           (if (even? (num-lo16 (num-shr16 x))) #x0000 #x8000))))
165                 (loop (num-shr16 x)
166                       (make-num z r)))))))))
168 (define add
169   (check-normalized 'add
170    (lambda (x y)
171      (check-norm 'add.x x)
172      (check-norm 'add.y y)
173      (let loop ((x x) (y y) (minus-c zero) (r #f))
174        (cond ((eq? x minus-c)
175               (norm r y))
176              ((eq? y minus-c)
177               (norm r x))
178              (else
179               (let ((z (+ (num-lo16 x) (num-lo16 y) (if (eq? minus-c zero) 0 1))))
180                 (loop (num-shr16 x)
181                       (num-shr16 y)
182                       (if (>= z #x10000) minus1 zero)
183                       (make-num (modulo z #x10000) r)))))))))
185 (define sub
186   (check-normalized 'sub
187    (lambda (x y)
188      (check-norm 'sub.x x)
189      (check-norm 'sub.y y)
190      (let loop ((x x) (y y) (minus-c minus1) (r #f))
191        (cond ((and (eq? x minus-c) (sign? y))
192               (norm r (invert y)))
193              ((eq? y (invert minus-c))
194               (norm r x))
195              (else
196               (let ((z (+ (num-lo16 x) (- #xffff (num-lo16 y)) (if (eq? minus-c zero) 0 1))))
197                 (loop (num-shr16 x)
198                       (num-shr16 y)
199                       (if (>= z #x10000) minus1 zero)
200                       (make-num (modulo z #x10000) r)))))))))
202 (define scale
203   (check-normalized 'scale
204    (lambda (n16 x)
205      (check-norm 'scale.x x)
206      (cond ((or (= n16 0)
207                 (eq? x zero))
208             zero)
209            ((eq? n16 1)
210             x)
211            (else
212             (let loop ((x x) (carry 0) (r #f))
213               (cond ((eq? x zero)
214                      (norm r (build-nonneg carry)))
215                     ((eq? x minus1)
216                      (norm r (build-neg (+ (modulo (- n16) #x10000) carry))))
217                     (else
218                      (let ((m (+ (* n16 (num-lo16 x)) carry)))
219                        (loop (num-shr16 x)
220                              (arithmetic-shift m -16)
221                              (make-num (modulo m #x10000) r)))))))))))
223 (define neg
224   (check-normalized 'neg
225    (lambda (x)
226      (check-norm 'neg.x x)
227      (sub zero x))))
229 (define build-nonneg
230   (check-normalized 'build-nonneg
231    (lambda (lo16)
232      (if (<= lo16 max-fix) (fixnum lo16) (make-num lo16 zero)))))
234 (define build-neg
235   (check-normalized 'build-neg
236    (lambda (lo16)
237      (if (>= lo16 (+ #x10000 min-fix)) (fixnum (- lo16 #x10000)) (make-num lo16 minus1)))))
239 (define mul
240   (check-normalized 'mul
241    (lambda (x y)
242      (check-norm 'mul.x x)
243      (check-norm 'mul.y y)
244      (if (neg? x)
245          (neg (mulnonneg (neg x) y))
246          (mulnonneg x y)))))
248 (define mulnonneg
249   (lambda (x y) ; x is nonnegative
250     (check-norm 'mulnonneg.x x)
251     (check-norm 'mulnonneg.y y)
252     (let loop ((x x) (s zero) (r #f))
253       (if (eq? x zero)
254           (norm r s)
255           (let ((a (add s (scale (num-lo16 x) y))))
256             (loop (num-shr16 x)
257                   (num-shr16 a)
258                   (make-num (num-lo16 a) r)))))))
260 (define rem
261   (check-normalized 'rem
262    (lambda (x y)
263      (check-norm 'rem.x x)
264      (check-norm 'rem.y y)
265      (if (neg? x)
266          (if (neg? y)
267              (neg (cdr (divnonneg (neg x) (neg y))))
268              (neg (cdr (divnonneg (neg x) y))))
269          (if (neg? y)
270              (cdr (divnonneg x (neg y)))
271              (cdr (divnonneg x y)))))))
273 (define div
274   (check-normalized 'div
275    (lambda (x y)
276      (check-norm 'div.x x)
277      (check-norm 'div.y y)
278      (if (neg? x)
279          (if (neg? y)
280              (car (divnonneg (neg x) (neg y)))
281              (neg (car (divnonneg (neg x) y))))
282          (if (neg? y)
283              (neg (car (divnonneg x (neg y))))
284              (car (divnonneg x y)))))))
286 (define divnonneg2
287   (lambda (x y)
288     (let ((x (decode x)) (y (decode y)))
289     (let ((lx (integer-length x))
290           (ly (integer-length y)))
291       (let ((b (- lx ly)))
292         (let loop ((i b) (a x) (b (arithmetic-shift y b)) (r 0))
293 ;          (pp (list i a b r))
294           (if (< i 0)
295               (cons (encode r) (encode a))
296               (begin
297                 (if (<= b a)
298                     (loop (- i 1)
299                           (- a b)
300                           (arithmetic-shift b -1)
301                           (+ 1 (* 2 r)))
302                     (loop (- i 1)
303                           a
304                           (arithmetic-shift b -1)
305                           (* 2 r)))))))))))
307 (define divnonneg
308   (lambda (x y)
309     (check-norm 'divnonneg.x x)
310     (check-norm 'divnonneg.y y)
311     (let ((lx (int-len-nonneg x))
312           (ly (int-len-nonneg y)))
313       (if (not (and (= lx (integer-length (decode x)))
314                     (= ly (integer-length (decode y)))))
315           (error "!!!!"))
316       (let ((b (- lx ly)))
317         (if (< b 0)
318             (cons zero x)
319             (let loop ((i b) (a x) (b (shift-left y b)) (r zero))
320 ;              (pp (list i (decode a) (decode b) (decode r)))
321               (if (< i 0)
322                   (cons r a)
323                   (if (le? b a)
324                       (loop (- i 1)
325                             (sub a b)
326                             (shr b)
327                             (add one (shl r)))
328                       (loop (- i 1)
329                             a
330                             (shr b)
331                             (shl r))))))))))
333 (define norm2
334   (check-normalized 'norm2
335    (lambda (shr16 r)
336      (if r
337          (let ((lo16 (num-lo16 r))
338                (newr (num-shr16 r)))
339            (norm newr
340                  (if (cond ((eq? shr16 zero)   (= lo16 #x0000))
341                            ((eq? shr16 minus1) (= lo16 #xffff))
342                            (else               #f))
343                      shr16
344                      (begin
345                        (num-shr16-set! r shr16)
346                        r))))
347          shr16))))
349 (define norm
350   (check-normalized 'norm
351    (lambda (r shr16)
352      (check-norm 'norm.shr16 shr16)
353      (if r
354          (let ((lo16 (num-lo16 r))
355                (newr (num-shr16 r)))
356            (cond ((eq? shr16 zero)
357                   (cond ((<= lo16 max-fix)
358                          (norm newr (fixnum lo16)))
359                         (else
360                          (num-shr16-set! r shr16)
361                          (norm newr r))))
362                  ((eq? shr16 minus1)
363                   (cond ((>= lo16 (+ #x10000 min-fix))
364                          (norm newr (fixnum (- lo16 #x10000))))
365                         (else
366                          (num-shr16-set! r shr16)
367                          (norm newr r))))
368                  (else
369                   (num-shr16-set! r shr16)
370                   (norm newr r))))
371          shr16))))
373 (define #%list->string
374   (lambda (lst)
375     (list->string (map integer->char (map decode lst)))))
377 (define #%cons cons)
379 (define #%< lt?)
380 (define #%+ add)
381 (define #%neg neg)
382 (define #%quotient div)
383 (define #%remainder rem)
385 (define #%number->string
386   (lambda (n)
387     (#%list->string
388      (if (#%< n zero)
389          (#%cons (make-num 45 zero) (#%number->string-aux (#%neg n) '()))
390          (#%number->string-aux n '())))))
392 (define #%number->string-aux
393   (lambda (n lst)
394     (let ((rest (#%cons (#%+ (make-num 48 zero) (#%remainder n (make-num 10 zero))) lst)))
395       (if (#%< n (make-num 10 zero))
396           rest
397           (#%number->string-aux (#%quotient n (make-num 10 zero)) rest)))))
399 ;------------------------------------------------------------------------------
401 (define fn
402  (lambda (i)
403    (list (- (expt 2 (+ i -1)) 3)
404          (- (expt 2 (+ i -1)) 2)
405          (- (expt 2 (+ i -1)) 1)
406          (expt 2 (+ i -1))
407          (+ (expt 2 (+ i -1)) 1)
408          (+ (expt 2 (+ i -1)) 2)
409          (+ (expt 2 (+ i -1)) 3)
410          
411          (- (expt 2 (+ i 0)) 3)
412          (- (expt 2 (+ i 0)) 2)
413          (- (expt 2 (+ i 0)) 1)
414          (expt 2 (+ i 0))
415          (+ (expt 2 (+ i 0)) 1)
416          (+ (expt 2 (+ i 0)) 2)
417          (+ (expt 2 (+ i 0)) 3)
418          
419          (- (expt 2 (+ i 1)) 3)
420          (- (expt 2 (+ i 1)) 2)
421          (- (expt 2 (+ i 1)) 1)
422          (expt 2 (+ i 1))
423          (+ (expt 2 (+ i 1)) 1)
424          (+ (expt 2 (+ i 1)) 2)
425          (+ (expt 2 (+ i 1)) 3)
426          )))
428 (define nums1
429   (append (fn (* 16 1))
430           (fn (* 16 2))
431           (fn (* 16 3))))
433 (define nums
434   (append (map - (reverse nums1))
435           (list -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10)
436           nums1))
438 (define dec
439   (lambda (n)
440     (cond ((eq? n zero) '(0))
441           ((eq? n minus1) '(-1))
442           (else (cons (num-lo16 n) (dec (num-shr16 n)))))))
444 (define decode
445   (lambda (n)
446     (cond ((eq? n zero) 0)
447           ((eq? n minus1) -1)
448           (else (+ (* #x10000 (decode (num-shr16 n))) (num-lo16 n))))))
450 (define encode
451   (lambda (n)
452     (cond ((and (>= n min-fix)
453                 (<= n max-fix))
454            (fixnum n))
455           (else
456            (make-num (modulo n #x10000)
457                      (encode (arithmetic-shift n -16)))))))
459 ;------------------------------------------------------------------------------
461 (pp 'cmp)
462 (for-each
463  (lambda (x)
464    (for-each
465     (lambda (y)
466       (let ((a (encode x))
467             (b (encode y))
468             (c (cond ((< x y) -1) ((= x y) 0) (else 1))))
469         (if (not (equal? c (cmp a b)))
470             (error "incorrect" x y (dec a) (dec b) c (cmp a b)))))
471     nums))
472  nums)
474 (pp 'int-len-nonneg)
475 (for-each
476  (lambda (x)
477    (if (>= x 0)
478        (let ((a (encode x))
479              (c (integer-length x)))
480          (if (not (equal? c (int-len-nonneg a)))
481              (error "incorrect" x (dec a) c (int-len-nonneg a))))))
482  nums)
484 (pp 'shl)
485 (for-each
486  (lambda (x)
487    (let ((a (encode x))
488          (c (encode (arithmetic-shift x 1))))
489      (if (not (equal? (dec c) (dec (shl a))))
490          (error "incorrect" x (dec a) (dec c) (dec (shl a))))))
491  nums)
493 (pp 'shr)
494 (for-each
495  (lambda (x)
496    (let ((a (encode x))
497          (c (encode (arithmetic-shift x -1))))
498      (if (not (equal? (dec c) (dec (shr a))))
499          (error "incorrect" x (dec a) (dec c) (dec (shr a))))))
500  nums)
502 (pp 'add)
503 (for-each
504  (lambda (x)
505    (for-each
506     (lambda (y)
507       (let ((a (encode x))
508             (b (encode y))
509             (c (encode (+ x y))))
510         (if (not (equal? (dec c) (dec (add a b))))
511             (error "incorrect" x y (dec a) (dec b) (dec c) (dec (add a b))))))
512     nums))
513  nums)
515 (pp 'sub)
516 (for-each
517  (lambda (x)
518    (for-each
519     (lambda (y)
520       (let ((a (encode x))
521             (b (encode y))
522             (c (encode (- x y))))
523         (if (not (equal? (dec c) (dec (sub a b))))
524             (error "incorrect" x y (dec a) (dec b) (dec c) (dec (sub a b))))))
525     nums))
526  nums)
528 (pp 'scale)
529 (for-each
530  (lambda (n16)
531    (for-each
532     (lambda (x)
533       (let ((a (encode x))
534             (c (encode (* n16 x))))
535         (if (not (equal? (dec c) (dec (scale n16 a))))
536             (error "incorrect" x (dec a) (dec c) (dec (scale n16 a))))))
537     nums))
538  '(0 1 2 3 4 5 6 7 8 9 10 32767 32768 65535))
540 (pp 'mul)
541 (for-each
542  (lambda (x)
543    (for-each
544     (lambda (y)
545       (let ((a (encode x))
546             (b (encode y))
547             (c (encode (* x y))))
548         (if (not (equal? (dec c) (dec (mul a b))))
549             (error "incorrect" x y (dec a) (dec b) (dec c) (dec (mul a b))))))
550     nums))
551  nums)
553 (pp 'div)
554 (for-each
555  (lambda (x)
556    (for-each
557     (lambda (y)
558       (if (not (= y 0))
559           (let ((a (encode x))
560                 (b (encode y))
561                 (c (encode (quotient x y))))
562             (if (not (equal? (dec c) (dec (div a b))))
563                 (error "incorrect" x y (dec a) (dec b) (dec c) (dec (div a b)))))))
564     nums))
565  nums)
567 (pp 'rem)
568 (for-each
569  (lambda (x)
570    (for-each
571     (lambda (y)
572       (if (not (= y 0))
573           (let ((a (encode x))
574                 (b (encode y))
575                 (c (encode (remainder x y))))
576             (if (not (equal? (dec c) (dec (rem a b))))
577                 (error "incorrect" x y (dec a) (dec b) (dec c) (dec (rem a b)))))))
578     nums))
579  nums)
581 (pp 'fact)
582 (define fact
583   (lambda (n)
584     (if (eq? n zero)
585         one
586         (mul n (fact (sub n one))))))
588 (define n (time (fact (encode 50))))
589 (pp (time (#%number->string n)))