乗算器のない8bit CPUで高速にべき乗剰余演算するコード

警告
現在、WZetaは8bit CPUから16bit CPUへバージョンアップしています。非常に良く似た命令セットですが互換性はありません。この べき乗剰余演算のWZetaのアセンブラコードは参考程度となっています。

目的

軽量で非常に良くできた8bit CPUのオープンソースWZetaを普及させるためにRSA暗号や楕円暗号などの公開鍵暗号で使われるべき乗剰余演算のアセンブラコードをCC0ライセンス(パブリックドメインと同じ)で公開したことのお知らせ。

警告
CC0であるのは、べき乗剰余演算のWZetaのアセンブラコードのみ

WZetaとは

WZetaはトランジスタ数当たりの性能、命令セット内パリティ、ハードマクロ命令など8bit CPUにして新技術を多数持ち、常識を破る変則的な命令セットですが、高い効率であることが確認されつつあるCPU。

1024bitのべき乗剰余演算するコード(事前定数有版)

###################################################################
### WZeta powmod 1024bit 2022/05/27 by Naoki Hirayama 
### 
### These codes are licensed under CC0.(Public Domain)
### http://creativecommons.org/publicdomain/zero/1.0/deed.ja
### 
### Montgomery multiplication method with radix 2,
### which is famous for being used in ICF3(1999).
###
### This code is for Open source 8bit CPU WZeta.
### https://wzeta.idletime.tokyo/
###
### download WZeta simulator
### https://subnote.icf3.net/wz660/download/
### 
### How to do simulation.
### > wzsim (this file)
### 
###################################################################
### ---------------------------------------------------------------
###  Subroutine powmod1024  y = g^x mod p
###  INOUT g,r(=y) 128,128 !0x100
###  IN    x,p     128,128 !0x100
###  r = R^2 mod p
### ---------------------------------------------------------------
	MEMORY 4
	MODEL TINY
###     MEMORY 8
### 	MODEL HALF
	
	.PROG
	LD A,&powmod_gr.H
	ST %_b_powmod_gr,A
	LD A,&powmod_px.H
	ST %_b_powmod_px,A

	LD A, ^_b_powmod1024.H
	LD B, ^_b_powmod1024.L
	BAL A:B
	
	LD A,32
	ST %_b_powmod_i,A
	LD A, &ANS.H
	LD B,0x80
        LD C,0x00
copy_ans:			# ans <- powmod_y
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	DEC %_b_powmod_i
	JRZ0 ^copy_ans
	
        $PRINT &ANS 128
	NOP
        $PRINT &powmod_exp 128
        EXIT	

### ---------------------------------------------------------------
###  [%h:%l] =[%h:%l] * g R-1mod p
###  IN %h:%l : pointer
###  CONST p,g
### ---------------------------------------------------------------
_b_powmod1024_mm:
	ST %_b_powmod_retA,A
	ST %_b_powmod_retB,B

	LD A,&_b_powmod_a.H	# a = 0
	LD B,&_b_powmod_a.L
	LD C,128		# repeat 129
	LOOPZERO

	ZERO %_b_powmod_m
_b_powmod1024_mm_loop_m:		# for m=0 to 127
	LD A,8
	ST %_b_powmod_n,A

        LD A, %_b_powmod_m
	ADD A, %_b_powmod_l
	LD B,A
        LD A, %_b_powmod_h
	LD A,[A:B]
	ST %_b_powmod_x,A
_b_powmod1024_mm_loop_n:		# for n=0 to 7
### loop body start --------------------------
	LD A, %_b_powmod_x
	LD B, 0
	AND A,[&powmod_g:B]
	XOR A,[&_b_powmod_a:B]
	ST %_b_powmod_u,A
	
	SHRC %_b_powmod_x
	JRC0 ^_b_powmod1024_mm2	# yi = 0 skip a += g

	LD B,0			# a += g
	CLC
	LD C,15
_b_powmod1024_mm1:
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	ADDC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm1
	INCC [&_b_powmod_a:B]
	
_b_powmod1024_mm2:
	SHRC %_b_powmod_u
	JRC0 ^_b_powmod1024_mm4	# u=0 skip +p

	LD B,0			# a += p
	CLC
	LD C,15
_b_powmod1024_mm3:
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	ADDC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm3
	INCC [&_b_powmod_a:B]

_b_powmod1024_mm4: 		# a >>= 1
	LD A, &_b_powmod_a.H
	LD B,128
	LD C,128
	CLC
	LOOPSHRC
### loop body end   --------------------------
	DEC %_b_powmod_n
 	JRZ0 ^_b_powmod1024_mm_loop_n

	INCX %_b_powmod_m	# A=m++
	XOR A,0x7F
 	JRZ0 ^_b_powmod1024_mm_loop_m

### Final Substruct
	INC %_b_powmod_n	# n=1
_b_powmod1024_mm5:	
	LD A,32			# y=a(before sub)
	ST %_b_powmod_m,A	# 32 = 128 / 4
	LD A,%_b_powmod_h
	LD C,%_b_powmod_l
	LD B,0
_b_powmod1024_mm6:
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	DEC %_b_powmod_m
	JRZ0 ^_b_powmod1024_mm6

	LD C,%_b_powmod_n
	LOOPINC ^_b_powmod1024_mm8 # if C!=0 jump

_b_powmod1024_mm7:		# return
	LD A,%_b_powmod_retA
	LD B,%_b_powmod_retB
	JMP A:B

_b_powmod1024_mm8:
	LD B,0			# a -= p
	LD A,0
	SUB A,C			# CF=1
	LD C,31
_b_powmod1024_mm9:
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	SUBC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm9
	DECC [&_b_powmod_a:B]
	JRC0 ^_b_powmod1024_mm7
	ZERO %_b_powmod_n
	JR ^_b_powmod1024_mm5
### ---------------------------------------------------------------

### ---------------------------------------------------------------
###  Subroutine powmod1024
###        INOUT g,r(=y) 128,128 !0x100
###  const IN    x,p     128,128 !0x100
###  work _b_powmod_a    129     !0x100
###  r : R2modP
### ---------------------------------------------------------------
_b_powmod1024:
	ST %_b_powmod_RETA,A
	ST %_b_powmod_RETB,B
	
### CALL powmod_1024mm
	LD A, %_b_powmod_gr	# y = y(r) * g R-1mod p
	ST %_b_powmod_h , A
	LD A, 0x80
	ST %_b_powmod_l , A
	BR ^_b_powmod1024_mm

### COPY y -> g
	LD A, %_b_powmod_gr
        LD C,0x7F
	LD B,0xFF
copy_y_g:	
	STAC [A:B]
	LOOPDEC ^copy_y_g 		# if c--==0 jump -1
	
# powmod_y = 1 start	
	LD A, %_b_powmod_gr
	LD B,0x80		# y
	LD C, 1
	!ST [A:B],C		# y[0] = 1
	LD C,126		# repeat 127
	LOOPZERO
	
	ZERO %_b_powmod_i			# for i=0 to 127
powmod_loop_i:
	LD A,8
	ST %_b_powmod_j,A			# for j=0 to 7
	
	LD A, &powmod_x.L
	ADD A, %_b_powmod_i
	LD B,A
	LD A,[&powmod_x:B]
	ST %_b_powmod_y,A
powmod_loop_j:
### 
### loop main
### 
	SHRC %_b_powmod_y
	JRC0 ^powmod1024_lbl3	# skip yi=0
	
### CALL powmod_1024mm (y*g)
	LD A, &powmod_y.H
	ST %_b_powmod_h , A
	LD A, &powmod_y.L
	ST %_b_powmod_l , A
	LD B, ^_b_powmod1024_mm.L
	BAL ^_b_powmod1024_mm:B

powmod1024_lbl3:
### CALL powmod_1024mm (g*g)
	LD A, &powmod_g.H
	ST %_b_powmod_h , A
	LD A, &powmod_g.L
	ST %_b_powmod_l , A
	LD B, ^_b_powmod1024_mm.L
	BAL ^_b_powmod1024_mm:B

	DEC %_b_powmod_j
 	JRZ0 ^powmod_loop_j
	
	INCX %_b_powmod_i		  	# A=[i++]
	XOR A,0x7F
 	JRZ0 ^powmod_loop_i

	LD A,%_b_powmod_RETA
	LD B,%_b_powmod_RETB
	JMP A:B
	
	.DATA
	REG 16 _reg02
	ALIAS _b_powmod_gr _reg02 0
	ALIAS _b_powmod_px _reg02 1
	ALIAS _b_powmod_x  _reg02 2
	ALIAS _b_powmod_y  _reg02 3
	ALIAS _b_powmod_i  _reg02 4
	ALIAS _b_powmod_j  _reg02 5
	ALIAS _b_powmod_m  _reg02 6
	ALIAS _b_powmod_n  _reg02 7
	ALIAS _b_powmod_h  _reg02 8
	ALIAS _b_powmod_l  _reg02 9
	ALIAS _b_powmod_u  _reg02 10
	ALIAS _b_powmod_retA  _reg02 11 # child sub
	ALIAS _b_powmod_retB  _reg02 12 # child sub
	ALIAS _b_powmod_RETA  _reg02 13
	ALIAS _b_powmod_RETB  _reg02 14

	MEM 129 _system_work !0x100
	ALIAS _b_powmod_a _system_work 0
	
	MEM 256 powmod_gr !0x100 \
	    0x3A 0x43 0xCE 0xEB 0xEC 0x52 0x88 0x94 0x4B 0x42 0x09 0xF6 0x91 0x70 0x7B 0x7F
	    0x25 0x1F 0xB8 0x44 0x32 0x0F 0x05 0x7F 0xF9 0xDD 0xA5 0x94 0x58 0x32 0x70 0x1E
	    0xCD 0x91 0x16 0xAC 0x21 0x16 0x90 0x79 0xCC 0xC2 0x66 0xC1 0x88 0x85 0x06 0x3B
	    0xF5 0x66 0xF2 0x31 0x89 0x74 0x2C 0x87 0x1D 0xBF 0x93 0x85 0x80 0x40 0x20 0xD3
	    0x18 0x47 0x16 0xE1 0xC5 0xB7 0xEB 0x70 0xF6 0xF0 0x2C 0x34 0x1A 0xBA 0xD5 0x9C
	    0x8C 0xA4 0xEE 0xDD 0x3C 0xD8 0xD9 0x7B 0xC4 0xA6 0x15 0xB3 0xAC 0x26 0x3D 0x45
	    0x21 0x7B 0x67 0x43 0xA7 0x63 0x16 0x7B 0x83 0x3C 0x1B 0xAE 0x6A 0x8F 0x9E 0x49
	    0xBE 0x99 0x01 0x92 0x1A 0xDE 0x6A 0xC5 0x7B 0x48 0xBA 0xF9 0x0A 0x29 0x82 0x46
	    0x3E 0x9C 0xD5 0x81 0xDB 0xC7 0x70 0xFD 0xEB 0x6C 0x52 0x37 0x07 0xF8 0x61 0x3F
	    0x79 0x31 0x57 0x32 0xC4 0x44 0x51 0x54 0x10 0x14 0x98 0x82 0x6B 0x0D 0xE5 0x50
	    0xA8 0xC0 0xC2 0x11 0x48 0xF1 0x38 0x4B 0x86 0xAC 0xF0 0xA5 0x34 0x37 0x52 0x67
	    0x89 0x88 0x93 0x4A 0x08 0x1E 0x92 0x17 0xE5 0x15 0xE5 0x14 0xB9 0x0C 0xCE 0x89
	    0xC7 0x0F 0x92 0xC5 0x1C 0xA3 0xF6 0x77 0xB8 0x88 0x50 0xF0 0x57 0xC8 0x62 0xED
	    0x34 0x89 0xE8 0x33 0xD2 0x22 0x6D 0x1B 0x45 0xC7 0xEB 0x84 0x32 0x13 0x18 0x84
	    0x84 0x69 0x24 0xD5 0x87 0xC0 0x5D 0x82 0xE8 0x55 0x2A 0x1A 0x77 0xEC 0x04 0xA5
	    0x21 0x86 0x16 0x04 0x25 0x9D 0x6D 0x66 0x11 0x9A 0x0E 0xBD 0xA9 0xB1 0xD1 0x4A
	ALIAS powmod_g powmod_gr 0x00
	ALIAS powmod_y powmod_gr 0x80 # = R2modP

	MEM 256 powmod_px !0x100 \
	    0xA1 0x63 0xD2 0x90 0x0B 0xD0 0x3A 0xAC 0xB4 0x15 0x75 0xEA 0x5E 0x05 0x57 0x86
	    0x55 0x12 0x36 0xE2 0xF8 0x07 0xCD 0x3E 0xE5 0x65 0x6E 0x70 0xE0 0x57 0xAE 0x5E
	    0x2D 0x21 0x53 0xE1 0xF0 0xC4 0xE5 0x96 0xB9 0x83 0xB3 0x5C 0x37 0xE2 0x3C 0x54
	    0xBF 0xF7 0x0F 0x4A 0xD8 0x19 0x2C 0xCA 0x3C 0xCA 0xD9 0x28 0xC3 0x76 0x75 0x5B
	    0x13 0x93 0xA3 0xC4 0x98 0x4E 0x6B 0x99 0x33 0x79 0xDC 0xCE 0x9C 0xB0 0xE8 0xBA
	    0x71 0xCD 0x27 0x5E 0xFA 0x8D 0xB0 0x7C 0x08 0xCF 0xF6 0x24 0x1B 0xD4 0xDA 0xCC
	    0x09 0x27 0x72 0x9C 0x5F 0x8D 0x08 0x6D 0x95 0x33 0x28 0xEE 0xAF 0x5D 0xED 0xB2
	    0x6E 0xB0 0xC3 0xE4 0xA3 0x70 0x4E 0x6B 0xF2 0xF9 0x6E 0x9C 0x4F 0xAB 0xFB 0xD8
	    0xC0 0xFA 0x69 0x4E 0x7C 0xBF 0x6D 0x94 0x18 0xE7 0x5C 0xB7 0xE8 0xE2 0xF4 0xCB
	    0x4B 0x75 0x3B 0x93 0x78 0x55 0xFC 0x0A 0x0D 0x00 0xEE 0xFE 0x16 0x0A 0x1D 0xAF
	    0x28 0xD0 0x9B 0xF2 0x09 0x4B 0x31 0xDE 0xF5 0x54 0x0E 0xCD 0xFE 0x4E 0x98 0x73
	    0x31 0x07 0x7E 0x90 0x1E 0xC6 0x98 0xD0 0xD8 0x1F 0xE8 0x0B 0x34 0xB7 0xFF 0x73
	    0x4A 0x86 0x82 0x48 0x20 0x52 0x9C 0x57 0xD4 0x38 0x5F 0x8F 0x43 0x33 0xB6 0x91
	    0xC1 0x14 0x5D 0xA6 0xAF 0x7B 0x13 0xA8 0xFB 0x89 0xAC 0xC8 0xFA 0x40 0x1F 0x21
	    0x14 0x7C 0xE7 0x64 0xCE 0xFF 0x48 0x49 0x7F 0xFD 0x9B 0x1A 0xA5 0x97 0x42 0x4A
	    0x94 0xB9 0xD6 0x70 0x83 0x38 0xAE 0xB3 0x78 0xBF 0xA6 0x48 0x24 0x4D 0x81 0x8C
	ALIAS powmod_p powmod_px 0x00
	ALIAS powmod_x powmod_px 0x80
	
	MEM 128 ANS &0xF00
	MEM 128 powmod_exp \
	    0xB4 0xBB 0x5B 0xB6 0x58 0x7A 0x9D 0x7A 0xEB 0x1C 0xE9 0x79 0xBF 0x94 0xF6 0x60
	    0x57 0x37 0x87 0x60 0x20 0x85 0xB3 0xB6 0xDD 0x9B 0xF4 0xA5 0xD8 0xFA 0x4C 0xC6
	    0xAD 0x84 0x39 0xCC 0xEC 0xD3 0xDD 0x23 0x86 0x4A 0x2C 0x50 0x93 0xEC 0xD6 0xAD
	    0x0D 0xD6 0xC4 0xBC 0xA1 0xD7 0x36 0x4F 0x2A 0x80 0x3F 0xC3 0x0C 0x0F 0xCD 0x4F
	    0x6C 0x6D 0xF1 0x5F 0x3D 0x58 0xAF 0xB5 0x1A 0x5B 0x85 0xF4 0xB8 0x77 0x11 0xD0
	    0x48 0x94 0xD7 0xC4 0xFF 0xFD 0x21 0xC1 0x41 0xA4 0x34 0xBB 0xC3 0x01 0x3A 0x36
	    0x4A 0x30 0xC3 0x3E 0x91 0xB2 0x86 0x2E 0x82 0xB0 0x87 0x63 0x81 0xCF 0x9E 0x33
	    0xC0 0x22 0x73 0x3B 0x05 0x95 0x8D 0xEA 0x80 0xA1 0xEE 0x33 0x47 0xD7 0xF2 0x84

解説(事前定数有版)

y = g ^ x mod p ( r: R^2 mod p )

rはpの値によって決まるモンゴメリ乗算の事前定数。モンゴメリ乗算の基数2はICF3(1999年)でも使われていました。1024bitべき乗剰余演算のサブルーチンの入力パラメータは256バイトのブロック2個。powmod_grとpowmod_pxです。grは前半128バイトがg、後半128バイトがr。　pxは前半128バイトがp、後半128バイトがx。サブルーチンが終了するとrの場所に演算結果yが入っています。

1024bitのべき乗剰余演算するコード(事前定数無版)

###################################################################
### WZeta powmod 1024bit 2022/06/10 by Naoki Hirayama 
### 
### These codes are licensed under CC0.(Public Domain)
### http://creativecommons.org/publicdomain/zero/1.0/deed.ja
### 
### Montgomery multiplication method with radix 2,
### which is famous for being used in ICF3(1999).
###
### This code is for Open source 8bit CPU WZeta.
### https://wzeta.idletime.tokyo/
###
### download WZeta simulator
### https://subnote.icf3.net/wz660/download/
### 
### How to do simulation.
### > wzsim (this file)
### 
###################################################################
### ---------------------------------------------------------------
###  Subroutine powmod1024  y = g^x mod p
###  INOUT g,y     128,128 !0x100
###  IN    x,p     128,128 !0x100
### ---------------------------------------------------------------
	MEMORY 4
	MODEL TINY
###     MEMORY 8
### 	MODEL HALF
	
	.PROG
	LD A,&powmod_gy.H
	ST %_b_powmod_gy,A
	LD A,&powmod_px.H
	ST %_b_powmod_px,A

	LD A, ^_b_powmod1024.H
	LD B, ^_b_powmod1024.L
	BAL A:B
	
	LD A,32
	ST %_b_powmod_i,A
	LD A, &ANS.H
	LD B,0x80
        LD C,0x00
copy_ans:			# ans <- powmod_y
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	STACP [&powmod_y:B]
	DEC %_b_powmod_i
	JRZ0 ^copy_ans
	
        $PRINT &ANS 128
	NOP
        $PRINT &powmod_exp 128
        EXIT	

### ---------------------------------------------------------------
###  [%h:%l] =[%h:%l] * g R-1mod p
###  IN %h:%l : pointer
###  CONST p,g
### ---------------------------------------------------------------
_b_powmod1024_mm:
	ST %_b_powmod_retA,A
	ST %_b_powmod_retB,B

	LD A,&_b_powmod_a.H	# a = 0
	LD B,&_b_powmod_a.L
	LD C,128		# repeat 129
	LOOPZERO

	ZERO %_b_powmod_m
_b_powmod1024_mm_loop_m:		# for m=0 to 127
	LD A,8
	ST %_b_powmod_n,A

        LD A, %_b_powmod_m
	ADD A, %_b_powmod_l
	LD B,A
        LD A, %_b_powmod_h
	LD A,[A:B]
	ST %_b_powmod_x,A
_b_powmod1024_mm_loop_n:		# for n=0 to 7
### loop body start --------------------------
	LD A, %_b_powmod_x
	LD B, 0
	AND A,[&powmod_g:B]
	XOR A,[&_b_powmod_a:B]
	ST %_b_powmod_u,A
	
	SHRC %_b_powmod_x
	JRC0 ^_b_powmod1024_mm2	# yi = 0 skip a += g

	LD B,0			# a += g
	CLC
	LD C,15
_b_powmod1024_mm1:
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_g:B]
	ADDC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm1
	INCC [&_b_powmod_a:B]
	
_b_powmod1024_mm2:
	SHRC %_b_powmod_u
	JRC0 ^_b_powmod1024_mm4	# u=0 skip +p

	LD B,0			# a += p
	CLC
	LD C,15
_b_powmod1024_mm3:
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	!ADDC [&_b_powmod_a:B],A
	LD  A,[&powmod_p:B]
	ADDC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm3
	INCC [&_b_powmod_a:B]

_b_powmod1024_mm4: 		# a >>= 1
	LD A, &_b_powmod_a.H
	LD B,128
	LD C,128
	CLC
	LOOPSHRC
### loop body end   --------------------------
	DEC %_b_powmod_n
 	JRZ0 ^_b_powmod1024_mm_loop_n

	INCX %_b_powmod_m	# A=m++
	XOR A,0x7F
 	JRZ0 ^_b_powmod1024_mm_loop_m

### Final Substruct
	INC %_b_powmod_n	# n=1
_b_powmod1024_mm5:	
	LD A,32			# y=a(before sub)
	ST %_b_powmod_m,A	# 32 = 128 / 4
	LD A,%_b_powmod_h
	LD C,%_b_powmod_l
	LD B,0
_b_powmod1024_mm6:
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	STACP [&_b_powmod_a:B]
	DEC %_b_powmod_m
	JRZ0 ^_b_powmod1024_mm6

	LD C,%_b_powmod_n
	LOOPINC ^_b_powmod1024_mm8 # if C!=0 jump

_b_powmod1024_mm7:		# return
	LD A,%_b_powmod_retA
	LD B,%_b_powmod_retB
	JMP A:B

_b_powmod1024_mm8:
	LD B,0			# a -= p
	LD A,0
	SUB A,C			# CF=1
	LD C,31
_b_powmod1024_mm9:
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	!SUBC [&_b_powmod_a:B],A
	LD A,[&powmod_p:B]
	SUBC [&_b_powmod_a:B],A
	LOOPINC ^_b_powmod1024_mm9
	DECC [&_b_powmod_a:B]
	JRC0 ^_b_powmod1024_mm7
	ZERO %_b_powmod_n
	JR ^_b_powmod1024_mm5
### ---------------------------------------------------------------

### ---------------------------------------------------------------
###  Subroutine powmod1024
###        INOUT g,y     128,128 !0x100
###  const IN    x,p     128,128 !0x100
###  work _b_powmod_a    129     !0x100
### ---------------------------------------------------------------
_b_powmod1024:
	ST %_b_powmod_RETA,A
	ST %_b_powmod_RETB,B
	
### g = gR mod p
	LD A,4
	ST %_b_powmod_i,A	# for i=4 to 1
	ZERO %_b_powmod_j 	# for j=0 to 255
powmod_loop_gi:
	### always %_b_powmod_j = 0
powmod_loop_gj:

	LD A,&powmod_g.H	# g<<=1
	LD B,0
	LD C,0x7F
	CLC
	LOOPSHLC
	GETFLAG [0]
	JRC1 ^_b_powmod1024_lbl1 # skip copy

	LD A,&_b_powmod_a.H	 # copy a <- g
	LD B,0
	LD C,0x7F
_b_powmod1024_loop1:
	STAB [&powmod_g:B]
	LOOPINC ^_b_powmod1024_loop1

_b_powmod1024_lbl1:
	LD B,0
	LD A,[&powmod_p:B]
	SUB [&powmod_g:B],A
	!LD C,126
_b_powmod1024_loop2:
	LD A,[&powmod_p:B]
	SUBC [&powmod_g:B],A
	LOOPINC ^_b_powmod1024_loop2
	JRC1 ^_b_powmod1024_lbl2
	SHL [0]
	JRC1 ^_b_powmod1024_lbl2
	
	LD A,&powmod_g.H	 # copy g <- a
	LD B,0
	LD C,0x7F
_b_powmod1024_loop3:
	STAB [&_b_powmod_a:B]
	LOOPINC ^_b_powmod1024_loop3

_b_powmod1024_lbl2:
	DEC %_b_powmod_j
 	JRZ0 ^powmod_loop_gj
	DEC %_b_powmod_i
 	JRZ0 ^powmod_loop_gi
	
# powmod_y = 1 start	
	LD A, %_b_powmod_gy
	LD B,0x80		# y
	LD C, 1
	!ST [A:B],C		# y[0] = 1
	LD C,126		# repeat 127
	LOOPZERO
	
	ZERO %_b_powmod_i			# for i=0 to 127
powmod_loop_i:
	LD A,8
	ST %_b_powmod_j,A			# for j=0 to 7
	
	LD A, &powmod_x.L
	ADD A, %_b_powmod_i
	LD B,A
	LD A,[&powmod_x:B]
	ST %_b_powmod_y,A
powmod_loop_j:
### 
### loop main
### 
	SHRC %_b_powmod_y
	JRC0 ^powmod1024_lbl3	# skip yi=0
	
### CALL powmod_1024mm (y*g)
	LD A, &powmod_y.H
	ST %_b_powmod_h , A
	LD A, &powmod_y.L
	ST %_b_powmod_l , A
	LD B, ^_b_powmod1024_mm.L
	BAL ^_b_powmod1024_mm:B

powmod1024_lbl3:
### CALL powmod_1024mm (g*g)
	LD A, &powmod_g.H
	ST %_b_powmod_h , A
	LD A, &powmod_g.L
	ST %_b_powmod_l , A
	LD B, ^_b_powmod1024_mm.L
	BAL ^_b_powmod1024_mm:B

	DEC %_b_powmod_j
 	JRZ0 ^powmod_loop_j
	
	INCX %_b_powmod_i		  	# A=[i++]
	XOR A,0x7F
 	JRZ0 ^powmod_loop_i

	LD A,%_b_powmod_RETA
	LD B,%_b_powmod_RETB
	JMP A:B
	
	.DATA
	REG 16 _reg02
	ALIAS _b_powmod_gy _reg02 0
	ALIAS _b_powmod_px _reg02 1
	ALIAS _b_powmod_x  _reg02 2
	ALIAS _b_powmod_y  _reg02 3
	ALIAS _b_powmod_i  _reg02 4
	ALIAS _b_powmod_j  _reg02 5
	ALIAS _b_powmod_m  _reg02 6
	ALIAS _b_powmod_n  _reg02 7
	ALIAS _b_powmod_h  _reg02 8
	ALIAS _b_powmod_l  _reg02 9
	ALIAS _b_powmod_u  _reg02 10
	ALIAS _b_powmod_retA  _reg02 11 # child sub
	ALIAS _b_powmod_retB  _reg02 12 # child sub
	ALIAS _b_powmod_RETA  _reg02 13
	ALIAS _b_powmod_RETB  _reg02 14

	MEM 129 _system_work !0x100
	ALIAS _b_powmod_a _system_work 0
	
	MEM 256 powmod_gy &0xE00 \
	    0x3A 0x43 0xCE 0xEB 0xEC 0x52 0x88 0x94 0x4B 0x42 0x09 0xF6 0x91 0x70 0x7B 0x7F
	    0x25 0x1F 0xB8 0x44 0x32 0x0F 0x05 0x7F 0xF9 0xDD 0xA5 0x94 0x58 0x32 0x70 0x1E
	    0xCD 0x91 0x16 0xAC 0x21 0x16 0x90 0x79 0xCC 0xC2 0x66 0xC1 0x88 0x85 0x06 0x3B
	    0xF5 0x66 0xF2 0x31 0x89 0x74 0x2C 0x87 0x1D 0xBF 0x93 0x85 0x80 0x40 0x20 0xD3
	    0x18 0x47 0x16 0xE1 0xC5 0xB7 0xEB 0x70 0xF6 0xF0 0x2C 0x34 0x1A 0xBA 0xD5 0x9C
	    0x8C 0xA4 0xEE 0xDD 0x3C 0xD8 0xD9 0x7B 0xC4 0xA6 0x15 0xB3 0xAC 0x26 0x3D 0x45
	    0x21 0x7B 0x67 0x43 0xA7 0x63 0x16 0x7B 0x83 0x3C 0x1B 0xAE 0x6A 0x8F 0x9E 0x49
	    0xBE 0x99 0x01 0x92 0x1A 0xDE 0x6A 0xC5 0x7B 0x48 0xBA 0xF9 0x0A 0x29 0x82 0x46
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	    0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00 0x00
	ALIAS powmod_g powmod_gy 0x00
	ALIAS powmod_y powmod_gy 0x80

	MEM 256 powmod_px &0xD00 \
	    0xA1 0x63 0xD2 0x90 0x0B 0xD0 0x3A 0xAC 0xB4 0x15 0x75 0xEA 0x5E 0x05 0x57 0x86
	    0x55 0x12 0x36 0xE2 0xF8 0x07 0xCD 0x3E 0xE5 0x65 0x6E 0x70 0xE0 0x57 0xAE 0x5E
	    0x2D 0x21 0x53 0xE1 0xF0 0xC4 0xE5 0x96 0xB9 0x83 0xB3 0x5C 0x37 0xE2 0x3C 0x54
	    0xBF 0xF7 0x0F 0x4A 0xD8 0x19 0x2C 0xCA 0x3C 0xCA 0xD9 0x28 0xC3 0x76 0x75 0x5B
	    0x13 0x93 0xA3 0xC4 0x98 0x4E 0x6B 0x99 0x33 0x79 0xDC 0xCE 0x9C 0xB0 0xE8 0xBA
	    0x71 0xCD 0x27 0x5E 0xFA 0x8D 0xB0 0x7C 0x08 0xCF 0xF6 0x24 0x1B 0xD4 0xDA 0xCC
	    0x09 0x27 0x72 0x9C 0x5F 0x8D 0x08 0x6D 0x95 0x33 0x28 0xEE 0xAF 0x5D 0xED 0xB2
	    0x6E 0xB0 0xC3 0xE4 0xA3 0x70 0x4E 0x6B 0xF2 0xF9 0x6E 0x9C 0x4F 0xAB 0xFB 0xD8
	    0xC0 0xFA 0x69 0x4E 0x7C 0xBF 0x6D 0x94 0x18 0xE7 0x5C 0xB7 0xE8 0xE2 0xF4 0xCB
	    0x4B 0x75 0x3B 0x93 0x78 0x55 0xFC 0x0A 0x0D 0x00 0xEE 0xFE 0x16 0x0A 0x1D 0xAF
	    0x28 0xD0 0x9B 0xF2 0x09 0x4B 0x31 0xDE 0xF5 0x54 0x0E 0xCD 0xFE 0x4E 0x98 0x73
	    0x31 0x07 0x7E 0x90 0x1E 0xC6 0x98 0xD0 0xD8 0x1F 0xE8 0x0B 0x34 0xB7 0xFF 0x73
	    0x4A 0x86 0x82 0x48 0x20 0x52 0x9C 0x57 0xD4 0x38 0x5F 0x8F 0x43 0x33 0xB6 0x91
	    0xC1 0x14 0x5D 0xA6 0xAF 0x7B 0x13 0xA8 0xFB 0x89 0xAC 0xC8 0xFA 0x40 0x1F 0x21
	    0x14 0x7C 0xE7 0x64 0xCE 0xFF 0x48 0x49 0x7F 0xFD 0x9B 0x1A 0xA5 0x97 0x42 0x4A
	    0x94 0xB9 0xD6 0x70 0x83 0x38 0xAE 0xB3 0x78 0xBF 0xA6 0x48 0x24 0x4D 0x81 0x8C
	ALIAS powmod_p powmod_px 0x00
	ALIAS powmod_x powmod_px 0x80
	
	MEM 128 ANS &0xF00
	MEM 128 powmod_exp \
	    0xB4 0xBB 0x5B 0xB6 0x58 0x7A 0x9D 0x7A 0xEB 0x1C 0xE9 0x79 0xBF 0x94 0xF6 0x60
	    0x57 0x37 0x87 0x60 0x20 0x85 0xB3 0xB6 0xDD 0x9B 0xF4 0xA5 0xD8 0xFA 0x4C 0xC6
	    0xAD 0x84 0x39 0xCC 0xEC 0xD3 0xDD 0x23 0x86 0x4A 0x2C 0x50 0x93 0xEC 0xD6 0xAD
	    0x0D 0xD6 0xC4 0xBC 0xA1 0xD7 0x36 0x4F 0x2A 0x80 0x3F 0xC3 0x0C 0x0F 0xCD 0x4F
	    0x6C 0x6D 0xF1 0x5F 0x3D 0x58 0xAF 0xB5 0x1A 0x5B 0x85 0xF4 0xB8 0x77 0x11 0xD0
	    0x48 0x94 0xD7 0xC4 0xFF 0xFD 0x21 0xC1 0x41 0xA4 0x34 0xBB 0xC3 0x01 0x3A 0x36
	    0x4A 0x30 0xC3 0x3E 0x91 0xB2 0x86 0x2E 0x82 0xB0 0x87 0x63 0x81 0xCF 0x9E 0x33
	    0xC0 0x22 0x73 0x3B 0x05 0x95 0x8D 0xEA 0x80 0xA1 0xEE 0x33 0x47 0xD7 0xF2 0x84

解説(事前定数無版)

y = g ^ x mod p

モンゴメリ乗算の基数2はICF3(1999年)でも使われていました。1024bitべき乗剰余演算のサブルーチンの入力パラメータは256バイトのブロック2個。powmod_gyとpowmod_pxです。gyは前半128バイトがg、後半128バイトは演算結果。　pxは前半128バイトがp(奇数)、後半128バイトがx。

実行環境

WZetaシミュレータで実際に演算させてみることが可能です。

> wzsim (アセンブラのファイル)

警告
シミュレータにはC言語実装のシミュレータとverilogのバイナリがあります。
verilogでは1024bitのべき乗剰余演算に１週間以上かかることもあるので、ご注意ください。

実行結果

実行後、演算が終了すると演算結果に続いて期待値が表示されます。期待値はコードの最後にあるpowmod_expの値が、そのまま表示されます。

性能について

乗算命令が無いため8bitの加算器1個で演算しています。絶対時間では非常に遅いため軽量な楕円暗号か、数時間の演算時間が許容されるIoTシステムなどでの利用に向いています。BIOSを通してアプリを作成すればアプリを改変することなく超高速な暗号プロセッサSnakeCubeに置き換えられることを考えているため、性能が必要になったところで、暗号プロセッサSnakeCubeに置き換えることが可能になる予想です。シミュレータの実行終了後、命令数が表示されます。WZetaのSDogコアでは1命令4サイクルなので命令数を4倍して周波数を計算すると1024bitのべき乗剰余演算１回の時間が求まります。値に依存しますが、上記、サンプルコードでは666530315命令。CPUの周波数に8MHzを想定した場合、1サイクル125[ns]なので
666530315 × 4 × 125 ÷ 10^9 ≒ 333秒 ( 512bitのべき乗剰余演算だと約42秒 )

脆弱性対策

今回のコードは脆弱性対策の無い単純に最も高速なプログラムになっています。また最上位ビットを高速モードとして使っています。サイドチャネル攻撃を対策をするには値に依存しないようにダミーの計算時間を入れます。WZetaの命令セットはオペコードに依存せずメモリアクセスをするという低消費電力でない特性がありますが、これが電力解析などのサイドチャネル攻撃への耐性を高めることにもなっています。トランジスタ数が少ないのでオペコードに依存しない動作で0、1の切り替えが多くなったとしても、全体としての消費電力は少ないという予想です。