Ademの関係（Schemeについて教えて下さい）

Ademの関係について、
If a＜pb then

P^a・P^b ＝ Σ[t=0～(a/p)]｛(-1)^(a+t)・((p-1)(b-t)-1 a-pt)・P^(a+b-t)・P^t｝


If a＜pb then

P^a・β・P^b ＝ Σ[t=0～(a/p)]｛(-1)^(a+t)・((p-1)(b-t) a-pt)・β・P^(a+b-t)・P^t｝ + Σ[t=0～((a-1)/p)]｛(-1)^(a+t-1)・((p-1)(b-t)-1 a-pt-1)・P^(a+b-t)・β・P^t｝



とした時、この規則を適用して
1.有限回の書き換えでadmissible(既約形)な単次式の和となる
2.出てくる答えはただ一つ
となるような計算をするプログラムをSchemeでかきたいのですが、分かる方がいましたらヒントをいただけないでしょうか?
p＝2の場合のsampleは手元にあるのですが、現在pが奇素数の場合を考えています。
どのように変更&追加すればよいかお願いします。





【追伸】
p(奇素数)を法としたSteenrod代数Apは、Zpを係数とし、記号β,p^0,p^1,p^2,… で生成される。
degree(β)＝1
degree(p^k)＝2k(p-1)

Apの単項式 : β^ε0・p^s1・β^ε1・…・p^sk・β^sk
ただし ε0,ε1,…,εk＝0又は1
β^2＝0

s1>=ε1+p・s2
s2>=ε2+p・s3
・
・
なる式はadmissibleである。


以上、分かっていることを全て書きました。少しの情報でも(曖昧でも)かまわないのでお願い致します。

No.6 ベストアンサー 回答者： hpsk 回答日時： 2005/01/31 00:56

No.5で書いた私の勘違いの修正と、

2個以上のβの削除に対応させました。
No.3の関数から、下のものだけ差し替えてください。

できれば、正しい変形例などがあれば確認がしやすくて嬉しいのですが。。。

> 余分な補足になってしまいましたらご了承下さい。
いえ、その部分がまさに私が確認したかった点です。

;;
(define (adem2-aux a b p)
(let ((na (car a)) (nb (car b))
(ea (cdr a)) (eb (cdr b)))
(letrec
((iter
(lambda (stk i cond-end getkeisuu beta1 beta2)
(if (cond-end i) ;; cond-end
(reverse stk)
;; keisuu
(let ((keisuu (mod (getkeisuu i) p)))
(if (zero? keisuu)
(iter stk (+ i 1) cond-end getkeisuu beta1 beta2)
;; even-arg
(let ((beta1-cons
(if (zero? beta1)
'()
(list (cons 0 1)) )))
(if (zero? i)
(iter
(cons `(,keisuu
,@beta1-cons
(,(+ na nb) . ,(+ beta2 eb)))
stk)
(+ i 1) cond-end getkeisuu beta1 beta2)
(iter
(cons `(,keisuu
,@beta1-cons
(,(- (+ na nb) i) . ,ea)
(,i . ,eb))
stk)
(+ i 1) cond-end getkeisuu beta1 beta2))))))))
(cond-end-a (lambda (i) (> (* i p) na)))
(cond-end-b (lambda (i) (> (* i p) (1- na))))
(binom-arg1-a (lambda (i) (1- (* (1- p) (- nb i))) ))
(binom-arg1-b (lambda (i) (* (1- p) (- nb i)) ))
(binom-arg2-a (lambda (i) (- na (* p i)) ))
(binom-arg2-b (lambda (i) (1- (- na (* p i))) ))
(even-arg-a (lambda (i) (+ na i) ))
(even-arg-b (lambda (i) (1- (+ na i)) )))
(if (>= na (* p nb))
(list (list a b))
(if (= 1 ea)
(append
(iter '() 0 cond-end-a
(lambda (i)
(* (expt-1 (even-arg-a i))
(binom (binom-arg1-b i) (binom-arg2-a i))))
1 0)
(iter '() 0 cond-end-b
(lambda (i)
(* (expt-1 (even-arg-b i))
(binom (binom-arg1-a i) (binom-arg2-b i))))
0 1))
(iter '() 0 cond-end-a
(lambda (i)
(* (expt-1 (even-arg-a i))
(binom (binom-arg1-a i) (binom-arg2-a i))))
0 0) )))))

;; P^^a P^^b==> ? (modulo p)
;; P^ P^8 == (cons 10 1) (cons 8 0)
(define (adem2 a b p)
(cond
((zero? (car b))
(list (cons 1 (list (cons (car a) (+ (cdr a) (cdr b)))))) )
((zero? (car a))
(list (cons 1 (list a b))) )
((>= (car a) (* p (car b)))
(list (cons 1 (list a b))) )
(else
(adem2-aux a b p))))

(define (sq-adem cal-stack p)
(display cal-stack)(newline)

(cond ((null? cal-stack) ; no calc status remained
(full-pop-heap-mod-p p))
((null? (cadar cal-stack )) ; no <rhead> remained
(let ((k (caar cal-stack))
(tl (caddar cal-stack)))
(if (and (not (= 0 k))
(not (null? tl)))
(push-heap (cons k tl))))
(sq-adem (cdr cal-stack) p))
(else
(let ((k (caar cal-stack))
(rhd (cadar cal-stack))
(tl (caddar cal-stack)))
(cond
((and (= 0 (caar rhd))
(not (null? (cdr rhd))))
(sq-adem (cons (list k
(cons
(cons (car (cadr rhd))
(+ (cdr (car rhd))
(cdr (cadr rhd))))
(cddr rhd))
tl)
(cdr cal-stack))
p))
((or (null? tl)
(= 0 (caar rhd))
(>= (caar rhd) (* p (caar tl))))
(let* ((car-rhd (cons (caar rhd)
(if (= 1 (cdar rhd)) 1 0)))
(next-tl (if (and (= 0 (car car-rhd))
(= 0 (cdr car-rhd)))
tl
(cons car-rhd tl))))
(sq-adem (cons (list k
(cdr rhd)
next-tl)
(cdr cal-stack))
p) ))
(else
(let ((res (adem2 (car rhd) (car tl) p)))
(if (null? res)
(sq-adem (cdr cal-stack) p)
(sq-adem
(append
(map
(lambda (x) ;; x := (<k> . <cons-list>)
(list (mod (* (car x) k) p)
(pop-push (cdr x) (cdr rhd))
(cdr tl)))
res)
(cdr cal-stack))
p) ))))))))

;; convert Ap into reversed cons-list
(define (ap2cons-list ap)
(letrec
((iter (lambda (rem acc)
(cond
((null? rem) acc)
((eq? 'b (car rem))
(iter (cdr rem) (cons (cons 0 1) acc)))
(else
(iter (cdr rem) (cons (cons (car rem) 0) acc)) )))))
(iter ap '() ) ))

この回答への補足

回答ありがとうございます。

この問題の正しい変換例としましては、
mod 3 の場合、
P^2・P^1 ==> 0
P^5・P^3 ==> 2・P^8
P^3・β・P^1 ==> β・P^3・P^1
P^3・β・P^2・β ==> β・P^4・P^1・β + 2・β・P^5・β

mod 5 の場合、
P^2・P^3 ==> 0
P^5・P^3・P^5 ==> P^12・P^1 + 4・P^13
P^6・β・P^2 ==> β・P^7・P^1 + P^7・β・P^1 + 3・β・P^8 + 4・P^8・β

などがあります。

補足日時：2005/02/03 13:28

この回答へのお礼

大変ありがとうございます。
おかげ様でようやくAdemの関係式に従って式変形を行うプログラムができあがりました! (^^)

このプログラムでも本当に十分なのですが、計算量が多くなると値が返ってこなくなります。
例えば、mod 3で (P^27)^11 などとしたときには、パソコンの性能にもよると思いますが、難しいようです。

前に示したp=2の場合のサンプルプログラムの一部では、２項係数を計算するところでベクトルを用いて計算結果を表に蓄えておき、同じ計算を行うときはそこから取り出すという形で手間を省いて、効率を上げているように思います。

あと、Ademの関係で式変形を行う部分でも、同様に表として変換した結果を記憶しておき、同じ変形の計算を二度行わないという方法もあるように思われます。

計算効率をあげるにはどのようにすればよろしいでしょうか？

お礼日時：2005/02/08 10:26

No.7 回答者： hpsk 回答日時： 2005/02/11 11:02

変換例をいただいておきながら、私自身まだ確認できていないのですが。

うまく動きましたか！

> 計算効率をあげるにはどのようにすればよろしいでしょうか？
一般に、
(adem2 m n p)
の計算は、一般に他のadem2の計算の部分計算としても何度も使われることが多いので、
一度計算したものは覚えておきましょうというのが基本的なアイデアです。（p=2のときの元のプログラムもそうでした）

ちょっとプログラムを書いている時間がないので、とりあえず参考になるURLを示しておきます。

参考URL：http://www.geocities.co.jp/SiliconValley-PaloAlt …

この回答へのお礼

計算結果をキャッシュする関数、参考にさせていただきました！
p=2のときのプログラムを参考に、一度計算したものを覚えておく関数を加えると、若干速度が上がったような気がします。

本当にありがとうございました。

お礼日時：2005/02/23 15:30

No.5 回答者： hpsk 回答日時： 2005/01/20 10:31

> (1) β・P^(a+b)

> (2) P^(a+b)・β
> の二通り形があるように思われ
すみません。完全に勘違いしていまして、(2)の場合しか考えていませんでした。
というか、そもそもP^0って消滅すると思い込んでいたのですが、質問文の定義を見直すとそういうわけではないんですね。
そうなるとadem2-aux中の(zero? i)の条件分岐は必要ないです。

> βは2回以上連続して現れると全て0になることがわかりました
変換途中にこれが現れる場合を考慮にいれると、それなりに修正する必要がありそうです。
考え直してみます。

> na,nbはPの指数、ea,ebはβの指数になるように作られていると考えてよろしいでしょうか?
> あと、ここでの(adem2-aux a b p)や(adem2 a b p)のa,bというのも(8 ．0)や(9．1)のようなドットペアを想定しているのでしょうか?
この２点についてはその通りです。

この回答への補足

P^(na)・β^(ea)・P^(nb)・β^(eb)
というPから始まるペアの形を常に想定しているのですね。
ありがとうございます。

＞そもそもP^0って消滅すると思い込んでいたのですが
の部分について補足なのですが、
P^0=1で考えていただいて大丈夫です。
もし私がこの言葉の意味をとらえ違えて、余分な補足になってしまいましたらご了承下さい。

補足日時：2005/01/20 11:25

No.4 回答者： hpsk 回答日時： 2005/01/19 01:07

#3のお礼のところのご質問に関してですが、

私が示したプログラムでは、ユーザはsq-ademではなく、ademを利用することを想定しています。
ademの中で ap2cons-list を呼び出し、
(10 9 b 8) -> ((8 . 0) (9 . 1) (10 . 0))
のように表現方法の変換＋並び順を逆順にしてから、これを<rhead>として、
sq-ademを呼び出しています。
つまり、<rhead>と<tail>はこのようなドットペアのリストです。<k>は整数（係数）です。
（人間は左の「整数またはbのリスト」の表現が見やすいですが、内部的には右の「(sk,εk)の組のリスト」のほうが扱いやすかったからです）

sq-adem内部の計算についても補足したほうがよろしいでしょうか？
とりあえず、
(trace sq-adem)
(trace adem2)
としてから
(adem (list '(1 (10 9 b 8))) 3)
などを実行してみると様子がわかるかと思います。

この回答へのお礼

ご回答ありがとうございます。
おがげ様で実行してみることができました! Pの指数とβの指数をひとつの組みとしてみているのですね。

最初にβ・β＝0とお伝えしましたが、
β・β・β＝0、β・β・β・β＝0、 …
のようにβは2回以上連続して現れると全て0になることがわかりました。つまりそういった場合その単項式自体が消えてしまうみたいなのです。
具体的には,
(adem (list '(1 (b b))) 3) ;;(=β・β)
などと実行した場合に、
()
のような空リストが返ってくるようにしたいです。いまさら申し訳ございません。。これでもし変更部分があるようでしたら教えていただきたいと思いまして・・・

この部分と、adem2-aux,adem2の式変形を行っている部分に関してなのですが、
iが0のとき、
(iter (cons (cons keisuu (list (cons (+ na nb) (mod2+ ea eb)))) stk) (+ i 1) cond-end getkeisuu)
とありますが、iが0の場合、係数部分を除くと
(1) β・P^(a+b)
(2) P^(a+b)・β
の二通り形があるように思われ、ちょっと理解できませんでした。βの位置が異なれば式変形も変わってくるのでは?と感じたもので。
na,nbはPの指数、ea,ebはβの指数になるように作られていると考えてよろしいでしょうか?
あと、ここでの(adem2-aux a b p)や(adem2 a b p)のa,bというのも(8 ．0)や(9．1)のようなドットペアを想定しているのでしょうか?

以上、お礼を申し上げなければならないにもかかわらず合わせて質問を繰り返し恐縮ですがお願い致します。

お礼日時：2005/01/19 19:18

No.3 回答者： hpsk 回答日時： 2005/01/13 07:41

遅くなりました。

とりあえず、残りの部分も変更し、答らしきものが出てくるようにしました。
だた・・・変換順序によって解が一意にならないので、どこか間違っているような気がします。

示していただいたコードのうち、呼び出されない関数がたくさんあるので、下に示したもので全部です。

(define **orig-heapsize** 5000)

(define (fact n)
(letrec ((fact-t (lambda (nn acc)
(if (= nn 0) acc (fact-t (- nn 1) (* nn acc))))))
(fact-t n 1)))

(define (1- n) (- n 1))
(define (1+ n) (+ n 1))

;; binomial coefficient
(define (binom n r)
(if (or (< n 0) (< r 0) (< n r))
0
(/ (fact n) (fact r) (fact (- n r)))))

(define (mod m p)
(let ((rem (remainder m p)))
(if (>= rem 0) rem (+ p rem))))

(define (mod2+ m n)
(mod (+ m n) 2))

(define (expt-1 n)
(if (even? n) 1 -1))

;;
(define (adem2-aux a b p)
(let ((na (car a)) (nb (car b))
(ea (cdr a)) (eb (cdr b)))
(letrec
((iter
(lambda (stk i cond-end getkeisuu)
(if (cond-end i) ;; cond-end
(reverse stk)
;; keisuu
(let ((keisuu (mod (getkeisuu i) p)))
(if (zero? keisuu)
(iter stk (+ i 1) cond-end getkeisuu)
;; even-arg
(if (zero? i)
(iter
(cons (cons keisuu
(list (cons (+ na nb) (mod2+ ea eb))))
stk)
(+ i 1) cond-end getkeisuu)
(iter
(cons (cons keisuu
(list (cons (- (+ na nb) i) ea)
(cons i eb)))
stk)
(+ i 1) cond-end getkeisuu)))))))
(cond-end-a (lambda (i) (> (* i p) na)))
(cond-end-b (lambda (i) (> (* i p) (1- na))))
(binom-arg1-a (lambda (i) (1- (* (1- p) (- nb i))) ))
(binom-arg1-b (lambda (i) (* (1- p) (- nb i)) ))
(binom-arg2-a (lambda (i) (- na (* p i)) ))
(binom-arg2-b (lambda (i) (1- (- na (* p i))) ))
(even-arg-a (lambda (i) (+ na i) ))
(even-arg-b (lambda (i) (1- (+ na i)) )))
(if (>= na (* p nb))
(list (list a b))
(if (zero? ea)
(iter '() 0 cond-end-a
(lambda (i)
(* (expt-1 (even-arg-a i))
(binom (binom-arg1-a i) (binom-arg2-a i)))))
(iter '() 0 cond-end-a
(lambda (i)
(+ (* (expt-1 (even-arg-a i))
(binom (binom-arg1-b i) (binom-arg2-a i)))
(if (cond-end-b i)
0
(* (expt-1 (even-arg-b i))
(binom (binom-arg1-a i) (binom-arg2-b i)))))) ))))))

;; P^a β^a P^b β^b==> ? (modulo p)
;; P^10 β P^8 == (cons 10 1) (cons 8 0)
(define (adem2 a b p)
(cond
((>= (car a) (* p (car b)))
(list (cons 1 (list a b))) )
((zero? (car a)) (list (cons 1 (list a b))) )
((zero? (car b))
(list (list (cons (car a) (mod2+ (cdr a) (cdr b))))))
(else
(adem2-aux a b p))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; (append (reverse a) b)
(define (pop-push a b)
(if (null? a)
b
(pop-push (cdr a) (cons (car a) b))))

;; global variables for heap sort
(define *orig-heapsize* **orig-heapsize**)
(define *heapsize* *orig-heapsize*)
(define *heap* (make-vector (+ *heapsize* 1)))
(define *cnt* 0)
(define *heavier?* <)

;; Heap sort
(define (upward-check! i)
(let ((up (quotient i 2))
(curr (vector-ref *heap* i)))
(if (> i 1)
(if (*heavier?*
(vector-ref *heap* up)
curr)
(begin
(vector-set! *heap*
i
(vector-ref *heap* up))
(vector-set! *heap* up curr)
(upward-check! up))))))

(define (downward-check! i)
(let ((left (* i 2))
(right (+ (* i 2) 1))
(curr (vector-ref *heap* i)))
(cond ((<= right *cnt*)
(if (*heavier?* (vector-ref *heap* right)
(vector-ref *heap* left))
(if (*heavier?* curr
(vector-ref *heap* left))
(begin
(vector-set! *heap*
i
(vector-ref *heap* left))
(vector-set! *heap* left curr)
(downward-check! left)))
(if (*heavier?* curr
(vector-ref *heap* right))
(begin
(vector-set! *heap*
i
(vector-ref *heap* right))
(vector-set! *heap* right curr)
(downward-check! right)))))
((= *cnt* left)
(if (*heavier?* curr (vector-ref *heap* left))
(begin
(vector-set! *heap*
i
(vector-ref *heap* left))
(vector-set! *heap* left curr)))))))

(define (push-heap obj)
(if (<= *heapsize* *cnt*)
(begin
(let* ((new-heapsize (* *heapsize* 2))
(new-heap (make-vector (+ new-heapsize 1))))
(do ((i 0 (+ i 1)))
((> i *heapsize*))
(vector-set! new-heap
i
(vector-ref *heap* i)))
(set! *heapsize* new-heapsize)
(set! *heap* new-heap))))
(set! *cnt* (1+ *cnt*))
(vector-set! *heap* *cnt* obj)
(upward-check! *cnt*))

(define (pop-heap)
(if (> *cnt* 0)
(let ((lightest
(vector-ref *heap* 1)))
(vector-set! *heap*
1
(vector-ref *heap* *cnt*))
(set! *cnt* (- *cnt* 1))
(downward-check! 1)
lightest)
'*empty-heap*))

(define (full-pop-heap-mod-p p)
(define (iter res)
(cond ((zero? *cnt*) res)
((null? res)
(iter (cons (pop-heap) res)))
(else (let ((top (pop-heap)))
(if (equal? (cdr top) (cdar res))
(let ((k-sum (mod (+ (car top) (caar res)) p)))
(if (= 0 k-sum)
(iter (cdr res))
(iter (cons (cons k-sum (cdr top))
res))))
(iter (cons top res)))))))
(iter '()))

;; Adem relations in Steenrod algebra : Sq operations

;; Sq monomial = list of integers e.g. (8 3 1) for
;; Sq^8 Sq^3 Sq^1
;; Sq expressions = list of Sq monomials e.g. ((12)(11 1)(8 3 1)) for
;; Sq^12 + Sq^11 Sq^1 + Sq^8 Sq^3 Sq^1

;;;
;; (define (sqm>? a b)
;; (define (same-length-iter aa bb)
;; (cond ((null? aa) #f)
;; ((< (car aa) (car bb)) #t)
;; ((> (car aa) (car bb)) #f)
;; (else (same-length-iter (cdr aa) (cdr bb)))))
;; (let ((la (length a)) (lb (length b)))
;; (cond ((> la lb) #t)
;; ((< la lb) #f)
;; (else (same-length-iter a b)))))

(define (apm>? a b)
(letrec ((same-length-iter (lambda (aa bb)
(if (null? aa)
#f
(let ((aan (+ (* 2 (caar aa)) (cdar aa)))
(bbn (+ (* 2 (caar bb)) (cdar bb))))
(cond
((< aan bbn) #t)
((> aan bbn) #f)
(else (same-length-iter (cdr aa) (cdr bb)))))))))
(let ((at (cdr a)) (bt (cdr b)))
(let ((la (length at)) (lb (length bt)))
(cond ((> la lb) #t)
((< la lb) #f)
(else (same-length-iter at bt)))))))

;; cal-stack is a list of calculation status:
;; a calculation status is
;; (<k> <rhead> <tail>) where <tail> is an admissible
;; sq monomial and <rhead> is reversed head part.
;; <k> is coefficient.
;; acc is a binary tree of sqm's.
;;
(set! *heavier?* apm>?)
(set! *cnt* 0)

(define (sq-adem cal-stack p)
(cond ((null? cal-stack) ; no calc status remained
(full-pop-heap-mod-p p))
((null? (cadar cal-stack )) ; no <rhead> remained
(let ((k (caar cal-stack)))
(if (not (= 0 k))
(push-heap (cons k (caddar cal-stack)))))
(sq-adem (cdr cal-stack) p))
(else
(let ((k (caar cal-stack))
(rhd (cadar cal-stack))
(tl (caddar cal-stack)))
(if (or (null? tl)
(= 0 (caar rhd))
(>= (caar rhd) (* p (caar tl))))
(sq-adem (cons (list k
(cdr rhd)
(cons (car rhd) tl))
(cdr cal-stack))
p)
(let ((res (adem2 (car rhd) (car tl) p)))
(if (null? res)
(sq-adem (cdr cal-stack) p)
(sq-adem
(append
(map
(lambda (x) ;; x := (<k> . <cons-list>)
(list (mod (* (car x) k) p)
(pop-push (cdr x) (cdr rhd))
(cdr tl)))
res)
(cdr cal-stack))
p) )
))))))

;; convert Ap into reversed cons-list
(define (ap2cons-list ap)
(letrec
((iter (lambda (rem acc)
(cond
((null? rem) acc)
((eq? 'b (car rem))
(iter (cons 0 rem) acc))
((or (null? (cdr rem))
(not (eq? 'b (cadr rem))))
(iter (cdr rem)
(cons (cons (car rem) 0) acc)))
(else
(iter (cddr rem)
(cons (cons (car rem) 1) acc)))))))
(iter (if (eq? 'b (car ap)) (cons 0 ap) ap)
'() )))

;; convert cons-list into Ap
(define (cons-list2ap cons-list)
(apply append
(map (lambda (x)
(if (= (car x) 0)
(if (= 0 (cdr x))
'()
'(b))
(if (= (cdr x) 0)
(list (car x))
(list (car x) 'b))))
cons-list)))

;; aps : list of Ap
;; Ap : (k <integer-or-beta list>)
;; * 3 P^5 β P^4 + β P^10 P^9 => ( (3 (5 b 4)) (1 (b 10 9)) )
;; p : prime number ( 0 < k < p )
;; ex. (adem '((1 (9 b 8))) 3) ==> ((1 (14 b 3)) (2 (17 b)))
(define (adem aps p)
(let ((res
(sq-adem
(map (lambda (ap)
(list (car ap) (ap2cons-list (cadr ap)) '()))
aps)
p)))
(map (lambda (res1)
(list (car res1) (cons-list2ap (cdr res1))))
res)))

この回答へのお礼

本当にありがとうございます! (^ ^)
プログラムを見せていただいて感心するばかりなのですが、(sq-adem)のところで<cal-stack>は具体的にどのようなデータ型を想定しているのでしょうか?
この部分で簡略化の操作を最終的に既約形になるまで繰り返し適用していると思われるため、重要な部分だと感じているのですが、どうもうまく実行することができません。
それと<rhd>や<tl>が(9 b 8)のような係数以外の部分のリストだと思ったので、それのcaarをとって条件分岐しているところが理解できませんでした。。
これだけのものを作っていただいて初歩的な質問を繰り返し心苦しいばかりなのですがお願い致します。

お礼日時：2005/01/18 18:24

No.2 回答者： hpsk 回答日時： 2004/12/07 19:26

#1です。

すみません。時間がとれなくて、まだ適当に作ったものしかできてません。
ですが、一応忘れてないですよということを表明するために（笑）現時点のものを送っておきます。
（変換例がないので、これで正しいのかどうかもわかりません・・・）
とりあえず、変換後、同類項をまとめないといけないのですが、そのへんが今のところ未実装です。
元のソースにあった大域変数は、実行結果をキャッシュして効率を上げるためのもののようですので、バッサリ切り捨てました（残そうと思えば残せますが、プログラムが長くなるので）。

(define (fact n)
(letrec ((fact-t (lambda (nn acc)
(if (= nn 0) acc (fact-t (- nn 1) (* nn acc))))))
(fact-t n 1)))

;; binomial coefficient
(define (binom n r)
(/ (fact n) (fact r) (fact (- n r))))

(define (mod m p)
(let ((rem (remainder m p)))
(if (>= rem 0) rem (+ p rem))))

(define (mod2+ m n)
(mod (+ m n) 2))

;;
(define (adem2-aux a b p)
(let ((na (car a)) (nb (car b))
(ea (cdr a)) (eb (cdr b)))
(letrec ((iter
(lambda (stk i cond-end binom-arg1 binom-arg2 even-arg)
(if (cond-end i) ;; cond-end
stk
;; binom-arg1, binom-arg2
(let ((bn (remainder (binom (binom-arg1 i) (binom-arg2 i))
p)))
(if (zero? bn)
(iter stk (+ i 1)
cond-end binom-arg1 binom-arg2 even-arg)
;; even-arg
(let ((k (if (even? (even-arg i)) bn (- p bn))))
(if (zero? i)
(iter
(cons (cons k
(list (cons (+ na nb) (mod2+ ea eb))))
stk)
(+ i 1)
cond-end binom-arg1 binom-arg2 even-arg)
(iter
(cons (cons k (list
(cons (- (+ na nb) i) ea)
(cons i eb)))
stk)
(+ i 1)
cond-end binom-arg1 binom-arg2 even-arg))))))))
(cond-end-a (lambda (i) (> (* i p) na)))
(cond-end-b (lambda (i) (> (* i p) (1- na))))
(binom-arg1-a (lambda (i) (1- (* (1- p) (- nb i))) ))
(binom-arg1-b (lambda (i) (* (1- p) (- nb i)) ))
(binom-arg2-a (lambda (i) (- na (* p i)) ))
(binom-arg2-b (lambda (i) (1- (- na (* p i))) ))
(even-arg-a (lambda (i) (+ na i) ))
(even-arg-b (lambda (i) (1- (+ na i)) )))
(if (>= na (* p nb))
(list (list a b))
(if (zero? ea)
(reverse (iter '() 0 cond-end-a
binom-arg1-a binom-arg2-a even-arg-a))
(append-steenrod
(reverse (iter '() 0 cond-end-a
binom-arg1-b binom-arg2-a even-arg-a))
(reverse (iter '() 0 cond-end-b
binom-arg1-a binom-arg2-b even-arg-b))
) )))))

;; P^a β^a P^b β^b==> ? (modulo p)
;; P^10 β P^8 == (cons 10 1) (cons 8 0)
(define (adem2 a b p)
(cond
((>= (car a) (* p (car b)))
(list (list a b)))
((zero? (car a)) (list (list a b)))
((zero? (car b))
(list (list (cons (car a) (mod2+ (cdr a) (cdr b))))))
(else
(adem2-aux a b p))))

この回答への補足

;; Heap sort
(define (upward-check! i)
(let ((up (quotient i 2))
(curr (vector-ref *heap* i)))
(if (> i 1)
(if (*heavier?*
(vector-ref *heap* up)
curr)
(begin
(vector-set! *heap*
i
(vector-ref *heap* up))
(vector-set! *heap* up curr)
(upward-check! up))))))

(define (downward-check! i)
(let ((left (* i 2))
(right (+ (* i 2) 1))
(curr (vector-ref *heap* i)))
(cond ((<= right *cnt*)
(if (*heavier?* (vector-ref *heap* right)
(vector-ref *heap* left))
(if (*heavier?* curr
(vector-ref *heap* left))
(begin
(vector-set! *heap*
i
(vector-ref *heap* left))
(vector-set! *heap* left curr)
(downward-check! left)))
(if (*heavier?* curr
(vector-ref *heap* right))
(begin
(vector-set! *heap*
i
(vector-ref *heap* right))
(vector-set! *heap* right curr)
(downward-check! right)))))
((= *cnt* left)
(if (*heavier?* curr (vector-ref *heap* left))
(begin
(vector-set! *heap*
i
(vector-ref *heap* left))
(vector-set! *heap* left curr)))))))

(define (push-heap obj)
(if (<= *heapsize* *cnt*)
(begin
(let* ((new-heapsize (* *heapsize* 2))
(new-heap (make-vector (+ new-heapsize 1))))
(do ((i 0 (+ i 1)))
((> i *heapsize*))
(vector-set! new-heap
i
(vector-ref *heap* i)))
(set! *heapsize* new-heapsize)
(set! *heap* new-heap))))
(set! *cnt* (+ *cnt* 1))
(vector-set! *heap* *cnt* obj)
(upward-check! *cnt*))

(define (pop-heap)
(if (> *cnt* 0)
(let ((lightest
(vector-ref *heap* 1)))
(vector-set! *heap*
1
(vector-ref *heap* *cnt*))
(set! *cnt* (- *cnt* 1))
(downward-check! 1)
lightest)
'*empty-heap*))

(define (push-list-to-heap lst)
(do ((rem lst (cdr rem)))
((null? rem))
(push-heap (car rem))))

(define (full-pop-heap)
(do ((res '() (cons (pop-heap) res)))
((zero? *cnt*)
(begin (set! *cnt* 0)
res))))

(define (full-pop-heap-mod2)
(define (iter res)
(cond ((zero? *cnt*) res)
((null? res)
(iter (cons (pop-heap) res)))
(else (let ((top (pop-heap)))
(if (equal? top (car res))
(iter (cdr res))
(iter (cons top res)))))))
(iter '()))

;; Adem relations in Steenrod algebra : Sq operations

;; Sq monomial = list of integers e.g. (8 3 1) for
;; Sq^8 Sq^3 Sq^1
;; Sq expressions = list of Sq monomials e.g. ((12)(11 1)(8 3 1)) for
;; Sq^12 + Sq^11 Sq^1 + Sq^8 Sq^3 Sq^1

(define *debug* #t)

;; If the same two items appear in the list
;; they are removed. If the same item appears
;; even times, they are all removed, whereas
;; if it appears odd times, only one item remains.
;; e.g (remove-same-pairs '(1 1 3 2 1 3 2 2 1 3)) => (2 3)

(define (remove-same-pairs lis)
(define (iter res rem)
(if (or (null? rem) (null? (cdr rem)))
(pop-push res rem)
(let ((i (find-pos (car rem) (cdr rem))))
(if (null? i)
(iter (cons (car rem) res)
(cdr rem))
(iter res
(delete-nth i (cdr rem)))))))
(if (null? lis)
lis
(iter '() lis)))

(define (find-pos item lst)
(define (iter i inp)
(if (null? inp)
'()
(if (equal? item (car inp))
i
(iter (+ i 1) (cdr inp)))))
(iter 0 lst))

(define (sqm>? a b)
(define (same-length-iter aa bb)
(cond ((null? aa) #f)
((< (car aa) (car bb)) #t)
((> (car aa) (car bb)) #f)
(else (same-length-iter (cdr aa) (cdr bb)))))
(let ((la (length a)) (lb (length b)))
(cond ((> la lb) #t)
((< la lb) #f)
(else (same-length-iter a b)))))

(define (add-asqe asqe1 asqe2)
(define (iter rem1 rem2 stk)
(cond ((null? rem1)(pop-push stk rem2))
((null? rem2)(pop-push stk rem1))
((sqm>? (car rem1) (car rem2))
(iter (cdr rem1) rem2 (cons (car rem1) stk)))
((sqm>? (car rem2) (car rem1))
(iter rem1 (cdr rem2) (cons (car rem2) stk)))
(else
(iter (cdr rem1) (cdr rem2) stk))))
(cond ((null? asqe1) asqe2)
((null? asqe2) asqe1)
(else (iter asqe1 asqe2 '()))))

(define (add-asqe-list asqe-list)
(define (iter rem res)
(if (null? rem)
res
(iter (cdr rem)
(add-asqe (car rem) res))))
(iter asqe-list '()))

補足日時：2004/12/10 11:37

この回答へのお礼

ありがとうございます。前に記載した重要と思われる部分以外にもいろいろ変更しなければならないところがあるみたいなのですが、ちょっと分からななかったので記載した部分以外のプログラムもかいておきたいと思います。本当に役立たずですみません・・　先の回答は現時点のものということですが、すごいですね！これだけでも本当に感謝しております！(^^)
長すぎてここにはかききれないので回答No2のお礼、補足、回答No1のお礼に分けて記載させていただきます。分かりにくくてすみません。。。

(define (pop-push a b)
(if (null? a)
b
(pop-push (cdr a) (cons (car a) b))))

(define (delete-nth n lst)
(define (iter stk k rem)
(if (null? rem)
lst
(if (zero? k)
(pop-push stk (cdr rem))
(iter (cons (car rem) stk)
(- k 1)
(cdr rem)))))
(iter '() n lst))

(define (replace-nth n lst new-obj)
(define (iter stk k rem)
(if (zero? k)
(pop-push stk (cons new-obj (cdr rem)))
(iter (cons (car rem) stk) (- k 1) (cdr rem)))) (if (>= n (length lst))
lst
(iter '() n lst)))

(define (remove-if pred? lst)
(define (iter rem res)
(if (null? rem)
(reverse res)
(if (pred? (car rem))
(iter (cdr rem) res)
(iter (cdr rem) (cons (car rem) res)))))
(iter lst '()))

(define (preserve-if pred? lst)
(define (iter rem res)
(if (null? rem)
(reverse res)
(if (pred? (car rem))
(iter (cdr rem) (cons (car rem) res))
(iter (cdr rem) res))))
(iter lst '()))

(define (max-list-length lstlst)
(define (iter rem m)
(if (null? rem)
m
(let ((m1 (length (car rem))))
(iter (cdr rem)
(if (> m1 m)
m1
m)))))
(iter lstlst 0))

(define (sorted-member obj lst leq?)
(cond ((null? lst) #f)
((equal? obj (car lst)) lst)
((leq? (car lst) obj)
(sorted-member obj (cdr lst) leq?))
(else #f)))

;; global variables for heap sort
(define *orig-heapsize* **orig-heapsize**)
(define *heapsize* *orig-heapsize*)
(define *heap* (make-vector (+ *heapsize* 1)))
(define *cnt* 0)
(define *heavier?* <)

お礼日時：2004/12/10 12:13

No.1 回答者： hpsk 回答日時： 2004/12/01 22:47

なんというか、、、専門的過ぎてとても回答がつきそうにない御質問ですね。

。。

Schemeはわかるのですが、Adem Relationsというものを全く知らないので、この一週間くらい時間のあるときにネットでSteenrod代数とか少しかじってみました。
まだあまりよくわかってないのですが、要するに、書き換え規則に従って、例えば
P^8 β P^9 P^10
のような項をadmissibleな項の和に書き換えられれば良いわけですよね。（admissibleの定義は理解したつもりです）

なんとかなりそうな気がしたりしなかったり(^^;)するのですが、以下の４点を確認させていただけないでしょうか？

（１）
degreeの意味がまだわかっていませんが、今の問題を解くためだけなら必要ないでしょうか？

（２）
ネットではp=2のとき以外の説明はあまり見つけられないので確認しておきたいのですが、各項の係数である
(-1)^(a+t)・((p-1)(b-t)-1 a-pt)
などは、mod pで計算すればいいのですね？

（３）
> ただし ε0,ε1,…,εk＝0又は1
なのに
> β^2＝0
というのがよくわからないのですが、、、

（４）
Adem Relations（Steenrod代数）を理解している人にとっては、単に書き換え規則を適用するプログラムを書けばいいだけのようにも見えますが、障害となっている部分はどこでしょうか？

この回答への補足

回答ありがとうございます！
プログラミング自体が苦手で、質問者でありながらうまくお答えできるかわからず申し訳ないのですが・・・

（１）おそらく私も直接関わりはないのではないかなと思います。はっきりお答えすることができずすみません・・

（２）はい！mod pで計算してください。(^-^)
ちなみに、(a b) mod pの求め方は
a=[a(r) a(r-1) …　a(1) a(0)]=Σ[i=0～r]a(i)p^i
b=[b(r) b(r-1) …　b(1) b(0)]=Σ[i=0～r]b(i)p^i
とすると、
(a b)=(a(0) b(0))・(a(1) b(1))・…・(a(r) b(r))
となるそうです。

（３）β^2=0であるから、ε0,ε1,…,εkが0又は1となるため、 P^a・P^b の場合と、 P^a・β・P^b の場合の２通りの書き換え規則があるのだと思います。

（４）
(define **max-dim** 220)
(define **adem-max-dim** 110)
(define **orig-heapsize** 5000)

(define *derived-zero-functions* (make-vector (+ **max-dim** 1) #f))
(define *zero-functions* (make-vector (+ **max-dim** 1) #f))
(define *generating-zero-functions* (make-vector (+ **max-dim** 1) #f))
(define *sq-table* (make-vector (+ **max-dim** 1) #t))

(define (init max-dim)
;; (require 'sort)
(set! **max-dim** max-dim)
(set! **adem-max-dim**
(if (even? **max-dim**)
(/ **max-dim** 2)
(/ (+ **max-dim** 1) 2)))
(set! *derived-zero-functions* (make-vector (+ **max-dim** 1) #f))
(set! *zero-functions* (make-vector (+ **max-dim** 1) #f))
(set! *generating-zero-functions* (make-vector (+ **max-dim** 1) #f))
(set! *sq-table* (make-vector (+ **max-dim** 1) #t))
(do ((i 0 (+ i 1)))
((> i 8))
(vector-set! *zero-functions* i '())
(vector-set! *derived-zero-functions* i '())
(vector-set! *generating-zero-functions* i '()))
(set! *adem-max-dim* **adem-max-dim**)
(set! *adem-table*
(make-vector (+ *adem-max-dim* 1)))
(do ((i 1 (+ i 1)))
((> i *adem-max-dim*))
(vector-set! *adem-table* i (make-vector (* 2 i) #f))))

;;(vector-set! *zero-functions* 9 '(((6 (2 1)) (2 (4 2 1)))))

;; Global table for adem2
(define *adem-max-dim* **adem-max-dim**)
(define *adem-table*
(make-vector (+ *adem-max-dim* 1)))
(do ((i 1 (+ i 1)))
((> i *adem-max-dim*))
(vector-set! *adem-table* i (make-vector (* 2 i) #f)))

;; Binomial coefficient reduced modulo 2
(define (binom2 n m)
(cond ((or (< n 0) (< m 0) (< n m)) 0)
((or (zero? m) (= n m)) 1)
((and (even? n) (odd? m)) 0)
(else (binom2 (quotient n 2) (quotient m 2)))))

(define (adem2-aux a b)
(define (iter stk i)
(if (> (* i 2) a)
(reverse stk)
(if (zero? (binom2 (- b 1 i) (- a (* i 2))))
(iter stk (+ i 1))
(if (zero? i)
(iter
(cons (list (+ a b)) stk)
(+ i 1))
(iter
(cons (list (- (+ a b) i) i) stk)
(+ i 1))))))
(if (>= a (+ b b))
(list (list a b))
(reverse (iter '() 0))))

(define (adem2 a b)
(cond ((>= a (+ b b)) (list (list a b)))
((zero? a) (list (list b)))
((zero? b) (list (list a)))
((and (<= b *adem-max-dim*)
(<= a (+ b b)))
(let ((val (vector-ref
(vector-ref *adem-table* b)
a)))
(if val
val
(let ((res (adem2-aux a b)))
(vector-set!
(vector-ref *adem-table* b)
a
res)
res))))
(else (adem2-aux a b))))

;; Ex. (adem2 4 4) => ((6 2) (7 1))
;; Sq^4 Sq^4 = sq^6 Sq&2 + Sq^7 Sq^1

以上、p=2の場合の与えられたサンプルプログラムの重要だと思われる部分を記述したのですが、最終的にpが奇素数の場合に上の例のような簡略化を実行するようなプログラムに書き換えたいです。
Schemeどころかプログラミング知識も乏しい(*_*)ために、どのように変更すればよいかさっぱりわからず困っています・・・

不完全な補足かと思いますが、他にも不明な点や知りたいことがありましたらわかっている範囲でできる限りお答えいたしますので、よろしくお願いします！(*^_^*)

補足日時：2004/12/03 05:53

この回答へのお礼

;; cal-stack is a list of calculation status:
;; a calculation status is
;; (<rhead> <tail>) where <tail> is an admissible
;; sq monomial and <rhead> is reversed head part.
;; acc is a binary tree of sqm's.
;;
(set! *heavier?* sqm>?)
(set! *cnt* 0)

(define (sq-adem cal-stack)
(cond ((null? cal-stack)
(full-pop-heap-mod2))
((null? (caar cal-stack))
(push-heap (cadar cal-stack))
(sq-adem (cdr cal-stack)))
(else
(let ((rhd (caar cal-stack))
(tl (cadar cal-stack)))
(if (or (null? tl)
(>= (car rhd) (* 2 (car tl))))
(sq-adem (cons (list (cdr rhd)
(cons (car rhd) tl))
(cdr cal-stack)))
(let ((res (adem2 (car rhd) (car tl))))
(if (null? res)
(sq-adem (cdr cal-stack))
(sq-adem
(append
(map
(lambda (x)
(list (pop-push x (cdr rhd))
(cdr tl)))
res)
(cdr cal-stack))))))))))

;; compatibility with previous version
;;
(define (adem-monom monom)
(sq-adem (list (list (reverse monom) '()))))

(define (sqleq a b)
(let ((la (length a)) (lb (length b)))
(cond ((null? a) #t)
((null? b) #f)
((> la lb) #t)
((< la lb) #f)
((< (car a) (car b)) #t)
((> (car a) (car b)) #f)
(else (sqleq (cdr a) (cdr b))))))

(define sq
(lambda x (sq-adem (list (list (reverse x) '())))))

(define sqadd add-asqe)

(define (sqmulmonos sqm1 sqm2)
(sq-adem (list (list (reverse (append sqm1 sqm2)) '()))))

(define (sqmulone sqm sqe)
(let ((rh (reverse sqm)))
(sq-adem
(map (lambda (x) (list rh x)) sqe))))

(define (sqmul sqe1 sqe2)
(define (iter rem res)
(if (null? rem)
res
(let ((rhd (reverse (car rem))))
(iter
(cdr rem)
(add-asqe
(sq-adem
(map (lambda (tl)
(list rhd tl))
sqe2))
res)))))
(iter sqe1 '()))

(define (sqbra sqe1 sqe2)
(sqadd (sqmul sqe1 sqe2) (sqmul sqe2 sqe1)))

お礼日時：2004/12/10 12:13