// file kernel/n/ppc32/smod.S: operations on residues modulo BASE^n - 1
/*-----------------------------------------------------------------------+
 |  Copyright 2005-2006, Michel Quercia (michel.quercia@prepas.org)      |
 |                                                                       |
 |  This file is part of Numerix. Numerix is free software; you can      |
 |  redistribute it and/or modify it under the terms of the GNU Lesser   |
 |  General Public License as published by the Free Software Foundation; |
 |  either version 2.1 of the License, or (at your option) any later     |
 |  version.                                                             |
 |                                                                       |
 |  The Numerix Library is distributed in the hope that it will be       |
 |  useful, but WITHOUT ANY WARRANTY; without even the implied warranty  |
 |  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU  |
 |  Lesser General Public License for more details.                      |
 |                                                                       |
 |  You should have received a copy of the GNU Lesser General Public     |
 |  License along with the GNU MP Library; see the file COPYING. If not, |
 |  write to the Free Software Foundation, Inc., 59 Temple Place -       |
 |  Suite 330, Boston, MA 02111-1307, USA.                               |
 +-----------------------------------------------------------------------+
 |                                                                       |
 |                     Arithmtique modulo BASE^n - 1                    |
 |                                                                       |
 +-----------------------------------------------------------------------*/

                    ; +----------------------------------+
                    ; |  Soustraction modulo BASE^n - 1  |
                    ; +----------------------------------+

; void xn(ssub)(chiffre *a, long la, chiffre *b, long lb)
;
; entre :
; a = naturel de longueur la
; b = naturel de longueur lb > 0
;
; sortie :
; b <- (a - b) mod (BASE^lb - 1), non normalis

#ifdef assembly_sn_ssub
.globl _sn_ssub
_sn_ssub:

#define L(x) Lsn_ssub_##x

        # variables locales
	#define _d_    r2
	#define _x_    r3
	#define _y_    r4
        #define _b_    r5
        #define _lb_   r6
        #define _add_  r7
        #define _ctr_  r8
        #define _la_   r9

	subf.  _la_, _lb_, r4	; la <- la - lb
	bge    L(big_a)
	
        ; cas a petit (ne pas optimiser, ce cas est inutilis)
	subi   r10,  r3,   4	; r10  <- &a[-1]
	subi   _b_,  _b_,  4	; b    <- &b[-1]
	mr     r11,  _b_	; r11  <- &b[-1]
	subfc  _x_,  _x_,  _x_	; CA <- 1
	addc.  _la_, _lb_, _la_	; restaure la, CA <- 1
	beq    2f
	mtctr  _la_
1:
	lwzu   _x_,  4(r10)	; soustrait les chiffres communs
	lwzu   _y_,  4(r11)
	subfe  _x_,  _y_,  _x_
	stw    _x_,  0(r11)
	bdnz   1b
2:
	subf   _x_,  _la_, _lb_	; ctr <- lb - la
	mtctr  _x_
1:
	lwzu   _x_,  4(r11)	; soustrait le reste de b de 0
	subfze _x_,  _x_
	stw    _x_,  0(r11)
	bdnz   1b

	subfe. _x_,  _x_,  _x_	; recycle la retenue ngative
	beq    2f
	mtctr  _lb_
1:
	lwzu   _x_,  4(_b_)
	and.   _x_,  _x_,  _x_
	subi   _x_,  _x_,  1
	stw    _x_,  0(_b_)
	beq    1b
2:
	blr

	; cas la >= lb
L(big_a):

	; prpare le droulement de l'addition et de la soustraction
        mflr   r0		; r0  <- adresse de retour
	mr     r10,  r3		; r10 <- a
	neg    _d_,  _lb_
        clrlslwi _d_,_d_, 27,2  ; d <- 4*((-lb) % 32)
        slwi   _x_,  _d_,  2	; add  <- point d entre pour l addition
	bcl    20,31, L(here)	; lr   <- point d entre pour la soustraction 
L(here):
	mflr   _add_
	add    _add_, _x_, _add_
	addi   _x_,  _add_, lo16(Lsn_subloop - L(here))
/*	addis  _x_,  _x_,   ha16(Lsn_subloop - L(here)) */
	mtlr   _x_
	addi   _add_, _add_, lo16(Lsn_addloop - L(here))
/*	addis  _add_, _add_, ha16(Lsn_addloop - L(here)) */
	
	addi   _ctr_, _lb_,   31
	srawi  _ctr_, _ctr_,  5	; ctr <- ceil(lb/32), CA <- 0
	
        ; additionne les lb chiffres suivants et recycle la retenue
L(loop):
	subf   r10,  _d_,  r10	; recadre les pointeurs
	subf   r11,  _d_,  _b_
	mtctr  _ctr_
	mr     r12,  r11
	blrl			; b += a[i..i+lb-1] + ret
	mtlr   _add_
	subf.  _la_, _lb_, _la_
	bge    L(loop)

	; dernire tranche, incomplte
	add.   _la_, _lb_, _la_
	beq    L(cycle)
	neg    _x_,  _la_
        clrlslwi _x_,_x_,  27,0 ; x <- (-la) % 32
	add    _y_,  _x_,  _la_	; y <- 32*ceil(la/32)
	slwi   _x_,  _x_,  2	; x <- 4*((-la) % 32)
	subf   r10,  _x_,  r10  ; cadre a,b sur le multiple de 32 prcdent
	subf   r11,  _x_,  _b_
	mr     r12,  r11
	subf   _x_,  _d_,  _x_	; add  <- point d entre pour l addition
        slwi   _x_,  _x_,  2
	add    _add_, _x_, _add_
	mtlr   _add_
	srwi   _y_,  _y_, 5	; ctr <- ceil(la/32)
	mtctr  _y_
	blrl			; b += reste(a) + ret
	subfe. _x_,  _x_, _x_	; propage la retenue
	bne    L(done)
	subf   _x_,  _la_, _lb_
	mtctr  _x_
	subi   r11,  r11,  4
	
	; recycle la retenue
1:
	lwzu   _x_,  4(r11)
	addze. _x_,  _x_
	stw    _x_,  0(r11)
	bdnzt  eq,   1b
L(cycle):
	subi   r11,  _b_,  4
	mtctr  _lb_
	subfe. _x_,  _x_, _x_
	beq    1b

L(done):
	mtlr   r0		; termin
	blrl
	
	#undef _d_
	#undef _x_
	#undef _y_
        #undef _b_
        #undef _lb_
        #undef _add_
        #undef _ctr_
        #undef _la_
	
#undef L
#endif /* assembly_sn_ssub */

#if defined(assembly_sn_smul) || defined(assembly_sn_ssqr)

                      ; +-------------------------------+
                      ; |  Rduction modulo BASE^n - 1  |
                      ; +-------------------------------+

; void xn(sred)(chiffre *a, long la, chiffre *b, long n)
;
; entre :
;   a = naturel de longueur la
;   b = naturel de longueur n
;
; contraintes : n > 0, la >= 0
;
; sortie :
;   b <- a mod BASE^n - 1
;
; remarque :
;   fonction non implmente en C
;
; registres modifis : r0, r2-r12 <- ind.

#define L(x) .Lsn_sred_##x
Lsn_sred:
	
        ; variables locales
	#define _d_    r2
	#define _x_    r3
	#define _y_    r4
        #define _b_    r5
        #define _n_    r6
        #define _add_  r7
        #define _ctr_  r8
        #define _la_   r9

	subf.  _la_, _n_, r4	; la <- la - n
	bgt    L(medium_a)

	; si la <= n, copie a dans b et complte par des zros
	subi   r10,  r3,   4	; recule les pointeurs
	subi   _b_,  _b_,  4
	add.   _la_, _n_,  _la_	; restaure la
	beq    2f
	mtctr  _la_
1:
	lwzu   _x_,  4(r10)	; b <- a
	stwu   _x_,  4(_b_)
	bdnz   1b
2:
	subf.  _n_,  _la_,  _n_	; complte par des zros
	beq    2f
	mtctr  _n_
	li     _x_,  0
1:	
	stwu   _x_,  4(_b_)
	bdnz   1b
2:
	blr

	; cas la > n
L(medium_a):
        mflr   r0		; r0  <- adresse de retour
	bcl    20,31, L(here)
L(here):
	mflr   _add_		; add <- adresse dbut de la boucle d addition
	mr     r11,  r3		; r11 <- &a[0]
	slwi   _x_,  _n_,  2	; r10 <- &a[n]
	add    r10,  _x_,  r11
	subf.  _la_, _n_, _la_
	bge    L(big_a)
	
	; si n < la < 2n, b <- a[0..n-1] + a[n..la-1]
	add    _la_, _n_,  _la_
	neg    _d_,  _la_
        clrlslwi _d_,_d_, 27,2  ; d <- 4*((n-la) % 32)
        slwi   _x_,  _d_,  2	; add  <- point d entre pour l addition
	add    _add_, _x_, _add_
	addi   _add_, _add_, lo16(Lsn_addloop - L(here))
/*	addis  _add_, _add_, ha16(Lsn_addloop - L(here)) */
	addi   _ctr_, _la_,  31
	srawi  _ctr_, _ctr_,  5	; ctr <- ceil((la-n)/32), CA <- 0
	
	mtlr   _add_
	subf   r10,  _d_,  r10	; cadre les pointeurs
	subf   r11,  _d_,  r11
	subf   r12,  _d_,  _b_
	mtctr  _ctr_
	blrl			; b <- a[0..n-1] + a[n..la-1]

	subi   r11,  r11,  4	; propage la retenue
	subi   r12,  r12,  4
	subf   _x_,  _la_, _n_
	mtctr  _x_
1:
	lwzu   _x_,  4(r11)
	addze  _x_,  _x_
	stwu   _x_,  4(r12)
	bdnz   1b
	b      L(cycle)

	; cas la >= 2n
L(big_a):
	neg    _d_,  _n_
        clrlslwi _d_,_d_, 27,2  ; d <- 4*((-n) % 32)
        slwi   _x_,  _d_,  2	; add  <- point d entre pour l addition
	add    _add_, _x_, _add_
	addi   _add_, _add_, lo16(Lsn_addloop - L(here))
/*	addis  _add_, _add_, ha16(Lsn_addloop - L(here)) */
	addi   _ctr_, _n_,   31
	srawi  _ctr_, _ctr_,  5	; ctr <- ceil(n/32), CA <- 0
	
        ; additionne les n chiffres suivants et recycle la retenue
L(loop):
	mtlr   _add_
	subf   r10,  _d_,  r10	; recadre les pointeurs
	subf   r11,  _d_,  r11
	subf   r12,  _d_,  _b_
	mtctr  _ctr_
	blrl			; b += a[i..i+n-1] + ret
	subf.  _la_, _n_, _la_
	mr     r11,  _b_
	bge    L(loop)

	; dernire tranche, incomplte
	add.   _la_, _n_, _la_
	beq    L(cycle)
	neg    _x_,  _la_
        clrlslwi _x_,_x_,  27,0 ; x <- (-la) % 32
	add    _y_,  _x_,  _la_	; y <- 32*ceil(la/32)
	slwi   _x_,  _x_,  2	; x <- 4*((-la) % 32)
	subf   r10,  _x_,  r10  ; cadre a,b sur le multiple de 32 prcdent
	subf   r11,  _x_,  _b_
	mr     r12,  r11
	subf   _x_,  _d_,  _x_	; add  <- point d entre pour l addition
        slwi   _x_,  _x_,  2
	add    _add_, _x_, _add_
	mtlr   _add_
	srwi   _y_,  _y_, 5	; ctr <- ceil(la/32)
	mtctr  _y_
	blrl			; b += reste(a) + ret
	subfe. _x_,  _x_, _x_	; propage la retenue
	bne    L(done)
	subf   _x_,  _la_, _n_
	mtctr  _x_
	subi   r11,  r11,  4
	
	; recycle la retenue
1:
	lwzu   _x_,  4(r11)
	addze. _x_,  _x_
	stw    _x_,  0(r11)
	bdnzt  eq,   1b
L(cycle):
	subi   r11,  _b_,  4
	mtctr  _n_
	subfe. _x_,  _x_, _x_
	beq    1b

L(done):
	mtlr   r0		; termin
	blrl
	
	#undef _d_
	#undef _x_
	#undef _y_
        #undef _b_
        #undef _n_
        #undef _add_
        #undef _ctr_
        #undef _la_
	
#undef L
             ; +-------------------------------------------------+
             ; |  Dcomposition modulo BASE^p - 1 et BASE^p + 1  |
             ; +-------------------------------------------------+

; void xn(split_even)(chiffre *a, chiffre *b, chiffre *c, long p)
; entre :
;  a = naturel de longueur 2p
;  b = naturel de longueur p
;  c = naturel de longueur p+1
;
; contrainte :  p > 0, a,b,c non confondus
;
; sortie :
;   b  <- a mod BASE^p - 1
;   c  <- a mod BASE^p + 1
;
; remarque :
;   fonction non implmente en C
;
; registres modifis : r0, r2-r7, r10-r12 <- ind.

#define L(x) .Lsn_split_even_##x
Lsn_split_even:

        ; variables locales
	#define _d_   r2
	#define _x_   r3
	#define _y_   r4
	#define _c_   r5
	#define _p_   r6
	#define _sub_ r7

	; prpare le droulement de l'addition et de la soustraction
        mflr   r0		; r0  <- adresse de retour
	neg    _d_,  _p_
        clrlslwi _d_,_d_, 27,2  ; d <- 4*((-p) % 32)
	subf   r10,  _d_,  r3	; cadre les pointeurs
	subf   r12,  _d_,  r4
        slwi   _x_,  _d_,  2	; lr  <- point d entre pour l addition
	bcl    20,31, L(here)	; sub <- point d entre pour la soustraction 
L(here):
	mflr   _sub_
	add    _sub_, _x_, _sub_
	addi   _x_,  _sub_, lo16(Lsn_addloop - L(here))
/*	addis  _x_,  _x_,   ha16(Lsn_addloop - L(here)) */
	mtlr   _x_
	addi   _sub_, _sub_, lo16(Lsn_subloop - L(here))
/*	addis  _sub_, _sub_, ha16(Lsn_subloop - L(here)) */

	; b <- a0 + a1 mod BASE^p - 1
	slwi   _x_, _p_,  2	; cadre le pointeur sur a1
	add    r11, _x_,  r10
	addi   _x_, _p_,  31	; ctr <- ceil(p/32), CA <- 0
	srawi  _x_, _x_,  5
	mtctr  _x_
	blrl
	
	subf   r11, _d_,  r10	; recadre les pointeurs sur a0 et a1
	slwi   _y_, _p_,  2
	subf   r10, _y_,  r11
	subfe. _x_,  _x_, _x_	; recycle la retenue
	bne    L(sub)
	subi   r12, r12,  4
	subf   r12, _y_, r12
1:
	lwzu   _x_,  4(r12)
	addze. _x_,  _x_
	stw    _x_,  0(r12)
	bdnzt  eq,   1b

	; c <- a0 - a1 mod BASE^p + 1
L(sub):
	mtlr   _sub_
	subf   r12,  _d_,  _c_
	addi   _x_, _p_,  31	; ctr <- ceil(p/32)
	srwi   _x_, _x_,  5
	mtctr  _x_
	subfc  _x_,  _x_,  _x_	; CA <- 1
	blrl

	li     _x_, 0		; recycle la retenue
	stw    _x_, 0(r12)
	subfe. _x_,  _x_, _x_
	beq    L(done)
	subi   r12, _c_, 4
	subfc  _x_, _x_, _x_
1:
	lwzu   _x_,  4(r12)
	addze. _x_,  _x_
	stw    _x_,  0(r12)
	bdnzt  eq,   1b
	
L(done):
	mtlr   r0		; termin
	blrl
	
	#undef _d_
	#undef _x_
	#undef _y_
	#undef _c_
	#undef _p_
	#undef _sub_
	
#undef L

       ; +-------------------------------------------------------------+
       ; |  Dcomposition modulo BASE^(p+1/2) - 1 et BASE^(p+1/2) + 1  |
       ; +-------------------------------------------------------------+

; void xn(split_odd)(chiffre *a, chiffre *b, chiffre *c, long p)
; entre :
;  a = naturel de longueur 2p+1
;  b = naturel de longueur p+1
;  c = naturel de longueur p+1
;
; contrainte :  p > 0, a,b,c non confondus
;
; sortie :
;   b  <- a mod BASE^(p+1/2) - 1
;   c  <- a mod BASE^(p+1/2) + 1
;
; remarque :
;   fonction non implmente en C
;
; registres modifis : r0, r2-r12 <- ind.

#define L(x) .Lsn_split_odd_##x
Lsn_split_odd:

        ; variables locales
	#define _d_   r2
	#define _x_   r3
	#define _y_   r4
	#define _c_   r5
	#define _p_   r6
	#define _add_ r7
	#define _a_   r8
	#define _b_   r9

	; b <- a1 << 16
	slwi   r7,  _p_, 2	; r10 <- &a[p]
	add    r10, r7, r3
	subi   r11, r4,  4	; r11 <- &b[-1]
	mtctr  _p_
	li     r8,  0		; init retenue
1:
	lwzu   r7,  4(r10)
	insrwi r8,  r7, 16, 0
	stwu   r8,  4(r11)
	srwi   r8,  r7, 16
	bdnz   1b
	stw    r8,  4(r11)

	; prpare le droulement de l'addition et de la soustraction
        mflr   r0		; r0  <- adresse de retour
	addi   _p_,  _p_,  1	; p <- p+1
	neg    _d_,  _p_
        clrlslwi _d_,_d_, 27,2  ; d <- 4*((-p-1) % 32)
	subf   _a_,  _d_,  r3	; cadre les pointeurs
	subf   _b_,  _d_,  r4
	subf   r12,  _d_,  _c_
        slwi   _x_,  _d_,  2	; lr  <- point d entre pour la soustraction 
	bcl    20,31, L(here)	; add <- point d entre pour l addition
L(here):
	mflr   _add_
	add    _add_, _x_, _add_
	addi   _x_,  _add_, lo16(Lsn_subloop - L(here))
/*	addis  _x_,  _x_,   ha16(Lsn_subloop - L(here)) */
	mtlr   _x_
	addi   _add_, _add_, lo16(Lsn_addloop - L(here))
/*	addis  _add_, _add_, ha16(Lsn_addloop - L(here)) */

	; c <- a0 - (a1 << 16)
	mr     r10, _a_
	mr     r11, _b_
	addi   _x_, _p_,  31	; ctr <- ceil((p+1)/32)
	srwi   _x_, _x_,  5
	mtctr  _x_
	subfc  _x_,  _x_,  _x_	; CA <- 1
	blrl
	subfe. _y_,  _y_,  _y_	; recycle la retenue
	beq    L(add)
	andi.  _x_,  _x_,  0xffff
	stw    _x_,  -4(r12)
	subi   r12,  _c_,  4
	subfc  _x_,  _x_,  _x_
1:
	lwzu   _x_,  4(r12)
	addze. _x_,  _x_
	stw    _x_,  0(r12)
	beq    1b

	; b <- a0 + (a1 << 16)
L(add):
	mtlr   _add_
	mr     r10, _a_
	mr     r11, _b_
	mr     r12, _b_
	addi   _x_, _p_,  31	; ctr <- ceil((p+1)/32), CA <- 0
	srawi   _x_, _x_,  5
	mtctr  _x_
	blrl
	subfe. _y_,  _y_,  _y_	; recycle la retenue
	bne    L(done)
	lwzux  _x_,  _b_,  _d_
	lis    _y_,  1
	addc   _x_,  _y_,  _x_
	stw    _x_,  0(_b_)
	subfe. _x_,  _x_,  _x_
	bne    L(done)
1:
	lwzu   _x_,  4(_b_)
	addze. _x_,  _x_
	stw    _x_,  0(_b_)
	beq    1b
	
L(done):
	mtlr   r0		; termin
	blrl
	
	#undef _d_
	#undef _x_
	#undef _y_
	#undef _c_
	#undef _p_
	#undef _add_
	#undef _a_
	#undef _b_
	
#undef L

#endif /* defined(assembly_sn_smul) || defined(assembly_sn_ssqr) */

                   ; +------------------------------------+
                   ; |  Multiplication modulo BASE^n - 1  |
                   ; +------------------------------------+

; void xn(smul) (chiffre *a, long la, chiffre *b, long lb, chiffre *c, long n)
;
; entre :
; a = naturel de longueur la
; b = naturel de longueur lb
; c = naturel de longueur n
;
; contraintes : n > 0, 0 <= lb <= la
;
; sortie :
; c <- a*b mod (BASE^n - 1)

#ifdef assembly_sn_smul
#define L(x) Lsn_smul_##x

#ifdef debug_smul
.globl _sn_smul_buggy
_sn_smul_buggy:
#else
.globl _sn_smul
_sn_smul:
#endif

        ; variables locales
	#define _b_  r31
	#define _lb_ r30
	#define _c_  r29
	#define _n_  r28
	#define _a_  32(r1)

	; rserve 2n chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r28,  8(r1)	; sauvegarde r28-r31
	slwi   r2,   r8,   3
	addi   r2,   r2,  44
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2

	; sauve les paramtres
	mr     _b_,  r5
	mr     _lb_, r6
	mr     _c_,  r7
	mr     _n_,  r8

	; rduit a et b modulo BASE^n - 1
	la     r5,   _a_
	mr     r6,   _n_
	bl     Lsn_sred
	mr     r3,  _b_
	mr     r4,  _lb_
	slwi   r6,  _n_,  2
	la     r5,  _a_
	add    r5,  r6,  r5
	mr     r6,  _n_
	bl     Lsn_sred

	; effectue la multiplication
	la     r3,  _a_
	slwi   r4,  _n_,  2
	add    r4,  r4,  r3
	mr     r5,  _c_
	mr     r6,  _n_
	bl     Lsn_smul

	; termin
	lwz   r1,    0(r1)	; nettoie la pile
	lwz   r0,    4(r1)	; rcupre l adresse de retour
	mtlr  r0
	lmw   r28,   8(r1)	; rcupre r28-r31
	blr

	#undef _b_
	#undef _lb_
	#undef _c_
	#undef _n_
	#undef _a_

	; point d entre rcursif, paramtres = a,b,c,n
	; aiguillage selon la longueur et la parit
Lsn_smul:
	andi.  r2,  r6,  1
	beq    1f
	cmpwi  cr0, r6, smul_lim_odd
	bgt    L(big_odd)
	b      L(small)
1:
	cmpwi  cr0, r6, smul_lim_even
	bgt    L(big_even)

	; petite multiplication => Karatsuba
L(small):

	#define _c_ r31
	#define _n_ r30
	#define _x_ 32(r1)

	; rserve 2n chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r30,  8(r1)	; sauvegarde r30-r31
	slwi   r2,   r6,   3
	addi   r2,   r2,  44
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _c_,  r5		; sauve les paramtres
	mr     _n_,  r6

	; x <- a*b
	mr     r5,   r4
	mr     r4,   r6
	la     r7,   _x_
	bl     Lsn_karamul

	; point d entre pour ssqr
Lsn_smul_aux_small:

	; c <- x mod BASE^n - 1
	la     r3,   _x_
	slwi   r4,   _n_,  1
	mr     r5,   _c_
	mr     r6,   _n_
	bl     Lsn_sred

	; termin
	lwz   r1,    0(r1)	; nettoie la pile
	lwz   r0,    4(r1)	; rcupre l adresse de retour
	mtlr  r0
	lmw   r30,   8(r1)	; rcupre r30-r31
	blr

	#undef _c_
	#undef _n_
	#undef _x_

        ; cas n grand pair : dcompose en deux produits modulaires
L(big_even):

	#define _r_ r31
	#define _c_ r30
	#define _p_ r29
	#define _x_ 32(r1)

	; rserve 3n/2+1 chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r29,  8(r1)	; sauvegarde r29-r31
	mulli  r2,   r6,   6
	addi   r2,   r2,  48
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _r_,  r4		; sauve les paramtres
	mr     _c_,  r5
	srwi   _p_,  r6,  1

	; dcompose a et b
	la     r4,   _x_
	mr     r6,   _p_
	bl     Lsn_split_even	; x0 <- a mod BASE^p-1, c0 <- a mod BASE^p+1
	mr     r3,   _r_
	la     r4,   _x_
	slwi   r5,   _p_,  2
	add    r4,   r5,   r4	; r4 <- &x1
	add    r5,   r5,   r4	; r5 <- &x2
	mr     r6,   _p_
	bl     Lsn_split_even	; x1 <- b mod BASE^p-1, x2 <- a mod BASE^p+1

	; x2 <- a*b mod BASE^p + 1
	la     r3,   _x_
	slwi   r5,   _p_,  3
	add    r3,   r5,   r3	; r3 <- &x2
	mr     r4,   _c_
	mr     r5,   _p_
	bl     Lsn_mmul

	; c1 <- a*b mod BASE^p - 1
	la     r3,   _x_
	slwi   r4,   _p_,  2
	add    r5,   _c_,  r4	; r5 <- &c1
	add    r4,   r3,   r4	; r4 <- &x1
	mr     r6,   _p_
	bl     Lsn_smul
	
	; point d entre pour ssqr
Lsn_smul_aux_big_even:
	
	; c0 <- floor((x2+c1)/2), c1 <- floor((c1-x2)/2) mod BASE^p
	la     r4,   _x_	; r4 <- &x2
	slwi   r3,   _p_,  2
	add    r4,   r4,   r3
	add    r4,   r4,   r3
	lwzx   _r_,  r3,   r4	; r  <- x2[p]
	add    r3,   _c_,  r3	; r3 <- &c1
	mr     r5,   _c_	; r5 <- &c0
	mr     r6,   r3		; r6 <- &c1
	mr     r7,   _p_
	bl     Lsn_half_add_sub

	; retenue = (CA - r - 1) + (r6 + r)*BASE^p/2 + (r6 + r)*BASE^(2p)/2
	add    r6,   _r_,  r6	; r6 <- retenue sur BASE^p/2 et BASE^(2p)/2
	slwi   r5,   r6,   31	; r5 <- retenue sur BASE^(p-1) et BASE^(2p-1)
	srwi   r6,   r6,    1	; r6 <- retenue sur BASE^p
	subfe  _r_,  _r_,  r6	; r  <- retenue  recycler
	slwi   r3,   _p_,   2
	add    r3,   _c_,  r3	; r3 <- &c1
	lwz    r4,   -4(r3)	; ajoute la retenue sur BASE^(p-1)
	addc   r4,   r5,   r4
	stw    r4,   -4(r3)
	lwzu   r4,   0(r3)	; ajoute la retenue sur BASE^p
	adde   r4,   r6,   r4
	stw    r4,   0(r3)
	subfe. r4,   r4,   r4
	bne    2f
	subi   r4,   _p_,   2	; propage jusqu a BASE^(2p-1)
	mtctr  r4
1:
	lwzu   r4,   4(r3)
	addze. r4,   r4
	stw    r4,   0(r3)
	bdnzt  eq,   1b
2:
	slwi   r3,   _p_,   3
	add    r3,   _c_,  r3	; r3 <- &c1[p]
	lwz    r4,   -4(r3)	; ajoute la retenue sur BASE^(2p-1)
	adde   r4,   r5,   r4
	stw    r4,   -4(r3)

	addze  _r_,  _r_	; r <- retenue  recycler (entre -2 et 2)
	srawi  r5,   _r_,  2	; r5 <- signe(r)
1:
	slwi   r4,   _p_,  1
	mtctr  r4
	subi   r3,   _c_,  4
2:
	lwzu   r4,    4(r3)
	add    r4,   _r_,  r4
	stw    r4,   0(r3)
	addze. _r_,  r5
	bdnzf  eq,   2b
	bne    1b

	; termin
	lwz   r1,    0(r1)	; nettoie la pile
	lwz   r0,    4(r1)	; rcupre l adresse de retour
	mtlr  r0
	lmw   r29,   8(r1)	; rcupre r29-r31
	blr

	#undef _r_
	#undef _c_
	#undef _p_
	#undef _x_

        ; cas n grand impair : dcompose en deux produits
L(big_odd):

	#define _d_ r31
	#define _c_ r30
	#define _p_ r29
	#define _x_ 32(r1)

	; rserve 6p+6 chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r29,  8(r1)	; sauvegarde r29-r31
	mulli  r2,   r6,  12
	addi   r2,   r2,  56
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _d_,  r4		; sauve les paramtres
	mr     _c_,  r5
	srwi   _p_,  r6,  1

	; dcompose a et b
	la     r4,   _x_	; r4 <- &x5
	mulli  r6,   _p_,  20
	add    r4,   r6,   r4
	addi   r4,   r4,   20
	mr     r6,   _p_
	bl     Lsn_split_odd	; x5 <- a mod BASE^(p+1/2)-1, c0 <- a mod BASE^(p+1/2)+1

	mr     r3,   _d_
	la     r4,   _x_
	slwi   r5,   _p_,  4
	addi   r5,   r5,   16
	add    r4,   r5,   r4	; r4 <- &x4
	srwi   r5,   r5,   2
	subf   r5,   r5,   r4	; r5 <- &x3
	mr     r6,   _p_
	bl     Lsn_split_odd	; x4 <- b mod BASE^(p+1/2)-1, x3 <- b mod BASE^(p+1/2)+1

        ; x1:x0 <- a*b mod BASE^(p+1/2) + 1
	mr     r3,   _c_
	addi   r4,   _p_,  1	; r4 <- p+1
	la     r7,   _x_
	mr     r6,   r4		; r6 <- p+1
	mulli  r5,   r4,   12
	add    r5,   r5,   r7	; r5 <- &x3
	bl     Lsn_toommul
	
	; rinjecte le chiffre de rang 2p+1
	slwi   r3,   _p_,  3
	add    r3,   r1,   r3
	lwz    r3,   36(r3)	; r3 <- x0[2p+1]
1:
	la     r5,   28(r1)	; r5 <- &x0[-1]
	addi   r4,   _p_,  1
	add    r4,   _p_,  r4	; r4 <- 2p+1
	mtctr  r4
	li     r6,   0
2:
	lwzu   r4,   4(r5)
	addc   r4,   r3,   r4
	stw    r4,   0(r5)
	addze. r3,   r6
	bdnzf  eq,   2b
	bne    1b

        ; x3:x2 <- a*b mod BASE^(p+1/2) - 1
	addi   r4,   _p_,  1	; r4 <- p+1
	la     r7,   _x_
	mr     r6,   r4		; r6 <- p+1
	slwi   r5,   r4,   3
	add    r7,   r5,   r7	; r7 <- &x2
	add    r3,   r5,   r7	; r3 <- &x4
	srwi   r5,   r5,   1
	add    r5,   r5,   r3	; r5 <- &x5
	bl     Lsn_toommul

	; point d entre pour ssqr
Lsn_smul_aux_big_odd:
	
	; rinjecte le chiffre de rang 2p+1
	slwi   r3,   _p_,  4
	add    r3,   r1,   r3
	lwz    r3,   44(r3)	; r3 <- x2[2p+1]
1:
	addi   r4,   _p_,  1
	add    r4,   _p_,  r4	; r4 <- 2p+1
	slwi   r5,   r4,   2
	add    r5,   r5,   r1
	addi   r5,   r5,   32	; r5 <- &x2[-1]
	mtctr  r4
	li     r6,   0
2:
	lwzu   r4,   4(r5)
	addc   r4,   r3,   r4
	stw    r4,   0(r5)
	addze. r3,   r6
	bdnzf  eq,   2b
	bne    1b

	; c <- floor((x3:x2 + x1:x0)/2)
	; x <- floor((x3:x2 - x1:x0)/2) mod BASE^(2p+1) - 1
	la     r4,   _x_
	slwi   r3,   _p_,   3
	add    r3,   r4,   r3
	addi   r3,   r3,    8	; r3 <- &x2
	mr     r5,   _c_
	mr     r6,   r4		; r6 <- &x0
	addi   r7,   _p_,   1
	add    r7,   _p_,  r7	; r7 <- 2p+1
	bl     Lsn_half_add_sub

	subfe. r3,   r3,   r3
	beq    2f
	la     r3,   28(r1)	; reycle la retenue
1:
	lwzu   r4,   4(r3)
	and.   r4,   r4,   r4
	subi   r4,   r4,   1
	stw    r4,   0(r3)
	beq    1b
2:

	; x <- x << 16 + retenue de la division par 2
	la     r3,   28(r1)
	addi   r4,   _p_,   1	; r4 <- 2p+1
	add    r4,   _p_,  r4
	mtctr  r4
	slwi   r5,   r6,   15	; r5 <- retenue
1:
	lwzu   r4,   4(r3)
	insrwi r5,   r4, 16, 0
	stw    r5,   0(r3)
	srwi   r5,   r4,   16
	bdnz   1b
	stw    r5,   4(r3)

	; c <- c + r6*BASE^(2p+1)/2
	slwi.  r6,   r6,  31
	beq    2f
	slwi   r3,   _p_,  3
	lwzx   r4,   _c_,  r3	; r4 <- c[2p]
	addc.  r4,   r6,   r4
	stwx   r4,   _c_,  r3	; c[2p] <- c[2p] + BASE/2
	blt    2f
	subi   r3,   _c_,  4	; recycle la retenue
1:
	lwzu   r4,   4(r3)
	addze. r4,   r4
	stw    r4,   0(r3)
	beq    1b
2:

	; prpare le droulement des additions
	la     r10,  _x_
	slwi   r11,  _p_,  2
	add    r11,  _c_,  r11	; r11 <- &c[p]
	addi   _d_,  _p_,  1
	neg    _d_,  _d_
        clrlslwi _d_, _d_, 27,2 ; d <- 4*((-p-1) % 32)
	subf   r10,  _d_,  r10	; cadre les pointeurs
	subf   r11,  _d_,  r11
        slwi   r3,   _d_,  2	; r2  <- point d entre pour la soustraction 
	bcl    20,31, L(here)
L(here):
	mflr   r2
	add    r2,   r3,   r2
	addi   r2,   r2, lo16(Lsn_addloop - L(here))
/*	addis  r2,   r2, ha16(Lsn_addloop - L(here)) */
	addic  r5,  _p_,  32	; r5 <- ceil((p+1)/32), CA <- 0
	srwi   r5,  r5,   5
	
	mtlr   r2		; c1 <- c1 + x0
	mtctr  r5
	mr     r12,  r11
	blrl

	mtlr   r2		; c0 <- c0 + x1 + ret
	subf   r10,  _d_,  r10
	subf   r11,  _d_,  _c_
	mr     r12,  r11
	mtctr  r5
	blrl
	
	; propage et recycle la retenue
	subfe. r3,   r3,   r3
	bne    L(done)
	subi   r12,  r12,  4
	mtctr  _p_
1:
	lwzu   r3,   4(r12)
	addze. r3,   r3
	stw    r3,   0(r12)
	bdnzt  eq,   1b
	bne    L(done)
	subi   r12,  _c_,  4
1:
	lwzu   r3,   4(r12)
	addze. r3,   r3
	stw    r3,   0(r12)
	beq    1b
	
	; termin
L(done):
	lwz   r1,    0(r1)	; nettoie la pile
	lwz   r0,    4(r1)	; rcupre l adresse de retour
	mtlr  r0
	lmw   r29,   8(r1)	; rcupre r29-r31
	blr

	#undef _d_
	#undef _c_
	#undef _p_
	#undef _x_

#undef L
#endif /* assembly_sn_smul */

                        ; +---------------------------+
                        ; |  Carr modulo BASE^n - 1  |
                        ; +---------------------------+

; void xn(ssqr) (chiffre *a, long la, chiffre *c, long n)
;
; entre :
; a = naturel de longueur la
; c = naturel de longueur n
;
; contraintes : n > 0, 0 <= la
;
; sortie :
; c <- a^2 mod (BASE^n - 1)

#ifdef assembly_sn_ssqr
#define L(x) Lsn_ssqr_##x

#ifdef debug_smul
.globl _sn_ssqr_buggy
_sn_ssqr_buggy:
#else
.globl _sn_ssqr
_sn_ssqr:
#endif

        ; variables locales
	#define _c_  r31
	#define _n_  r30
	#define _a_  32(r1)

	; rserve n chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r30,  8(r1)	; sauvegarde r30-r31
	slwi   r2,   r6,   2
	addi   r2,   r2,  44
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2

	; sauve les paramtres
	mr     _c_,  r5
	mr     _n_,  r6

	; rduit a modulo BASE^n - 1
	la     r5,   _a_
	bl     Lsn_sred

	; effectue l lvation au carr
	la     r3,  _a_
	mr     r5,  _c_
	mr     r6,  _n_
	bl     Lsn_ssqr

	; termin
	lwz   r1,    0(r1)	; nettoie la pile
	lwz   r0,    4(r1)	; rcupre l adresse de retour
	mtlr  r0
	lmw   r30,   8(r1)	; rcupre r30-r31
	blr

	#undef _c_
	#undef _n_
	#undef _a_

	; point d entre rcursif, paramtres = a,.,c,n
	; aiguillage selon la longueur et la parit
Lsn_ssqr:
	andi.  r2,  r6,  1
	beq    1f
	cmpwi  cr0, r6, ssqr_lim_odd
	bgt    L(big_odd)
	b      L(small)
1:
	cmpwi  cr0, r6, ssqr_lim_even
	bgt    L(big_even)

	; petit carr => Karatsuba
L(small):

	#define _c_ r31
	#define _n_ r30
	#define _x_ 32(r1)

	; rserve 2n chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r30,  8(r1)	; sauvegarde r30-r31
	slwi   r2,   r6,   3
	addi   r2,   r2,  44
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _c_,  r5		; sauve les paramtres
	mr     _n_,  r6

	; x <- a^2
	mr     r4,   r6
	la     r5,   _x_
	bl     Lsn_karasqr

	; continue avec smul
	b      Lsn_smul_aux_small

	#undef _c_
	#undef _n_
	#undef _x_

        ; cas n grand pair : dcompose en deux carrs modulaires
L(big_even):

	#define _r_ r31
	#define _c_ r30
	#define _p_ r29
	#define _x_ 32(r1)

	; rserve 3n/2+1 chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r29,  8(r1)	; sauvegarde r29-r31
	mulli  r2,   r6,   6
	addi   r2,   r2,  48
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _c_,  r5		; sauve les paramtres
	srwi   _p_,  r6,  1

	; dcompose a
	la     r4,   _x_
	slwi   r5,   _p_,  3
	add    r5,   r5,   r4	; r5 <- &x2
	mr     r6,   _p_
	bl     Lsn_split_even	; x0 <- a mod BASE^p-1, x2 <- a mod BASE^p+1

	; x2 <- a^2 mod BASE^p + 1
	la     r3,   _x_
	slwi   r5,   _p_,  3
	add    r3,   r5,   r3	; r3 <- &x2
	mr     r4,   _p_
	bl     Lsn_msqr

	; c1 <- a^2 mod BASE^p - 1
	la     r3,   _x_
	slwi   r4,   _p_,  2
	add    r5,   _c_,  r4	; r5 <- &c1
	mr     r6,   _p_
	bl     Lsn_ssqr

	; continue avec smul
	b      Lsn_smul_aux_big_even

	#undef _r_
	#undef _c_
	#undef _p_
	#undef _x_

        ; cas n grand impair : dcompose en deux carrs
L(big_odd):

	#define _d_ r31
	#define _c_ r30
	#define _p_ r29
	#define _x_ 32(r1)

	; rserve 6p+6 chiffres + cadre de pile
	mflr   r0
	stw    r0,   4(r1)	; sauvegarde l adresse de retour
	stmw   r29,  8(r1)	; sauvegarde r29-r31
	mulli  r2,   r6,  12
	addi   r2,   r2,  56
	clrrwi r2,   r2,   4	; arrondi  un multiple de 16 octets
	neg    r2,   r2
	stwux  r1,   r1,  r2
	mr     _c_,  r5		; sauve les paramtres
	srwi   _p_,  r6,  1

	; dcompose a
	la     r4,   _x_	; r4 <- &x4
	slwi   r6,   _p_,  4
	add    r4,   r6,   r4
	addi   r4,   r4,   16
	mr     r6,   _p_
	bl     Lsn_split_odd	; x4 <- a mod BASE^(p+1/2)-1, c0 <- a mod BASE^(p+1/2)+1

        ; x1:x0 <- a^2 mod BASE^(p+1/2) + 1
	mr     r3,   _c_
	addi   r4,   _p_,  1	; r4 <- p+1
	la     r5,   _x_
	bl     Lsn_toomsqr
	
	; rinjecte le chiffre de rang 2p+1
	slwi   r3,   _p_,  3
	add    r3,   r1,   r3
	lwz    r3,   36(r3)	; r3 <- x0[2p+1]
1:
	la     r5,   28(r1)	; r5 <- &x0[-1]
	addi   r4,   _p_,  1
	add    r4,   _p_,  r4	; r4 <- 2p+1
	mtctr  r4
	li     r6,   0
2:
	lwzu   r4,   4(r5)
	addc   r4,   r3,   r4
	stw    r4,   0(r5)
	addze. r3,   r6
	bdnzf  eq,   2b
	bne    1b

        ; x3:x2 <- a^2 mod BASE^(p+1/2) - 1
	addi   r4,   _p_,  1	; r4 <- p+1
	la     r5,   _x_
	slwi   r3,   r4,   3
	add    r5,   r3,   r5	; r5 <- &x2
	add    r3,   r5,   r3	; r3 <- &x4
	bl     Lsn_toomsqr

	; continue avec smul
	b      Lsn_smul_aux_big_odd

	#undef _d_
	#undef _c_
	#undef _p_
	#undef _x_

#undef L
#endif /* assembly_sn_ssqr */

                       ; +----------------------------+
                       ; |  Combinaison de 3 rsidus  |
                       ; +----------------------------+

; void xn(sjoin3)(chiffre *a, long h, long k)
;
;  entre :
;  a = naturel de longueur n+p+q
;  n = (2h+2)k, p = (2h+1)k, q = (2h)k
;
;  contraintes : h >= 2, k >= 2
;
;  sortie :
;  a <- x mod ppcm(BASE^n - 1, BASE^p - 1, BASE^q - 1) normalis
;  avec
;    a[0..n-1]       = x mod (BASE^n - 1),
;    a[n..n+p-1]     = x mod (BASE^p - 1),
;    a[n+p..n+p+q-1] = x mod (BASE^q - 1)
;
;  remarque : ppcm = produit/(BASE^k - 1)/(BASE^(2k) - 1)

#ifdef assembly_sn_sjoin3
#undef L
#define L(x) .Lsn_sjoin3_##x
#ifdef debug_sjoin
.globl _sn_sjoin3_buggy
_sn_sjoin3_buggy:
#else
.globl _sn_sjoin3
_sn_sjoin3:
#endif

	; variables locales
	#define _a_   r2
	#define _p_   r6
	#define _k_   r5
	#define _add_ r7
	#define _sub_ r8
	#define _r_   r9

	mflr   r0
	mr     _a_,  r3
	add    _p_,  r4,  r4
	addi   _p_,  _p_, 1
	mullw  _p_,  _p_, _k_
	bcl    20,31, L(here)
L(here):
	mflr   _add_
	addi   _sub_, _add_, lo16(Lsn_subloop - L(here))
/*	addis  _sub_, _sub_, ha16(Lsn_subloop - L(here)) */
	addi   _add_, _add_, lo16(Lsn_addloop - L(here))
/*	addis  _add_, _add_, ha16(Lsn_addloop - L(here)) */

	; normalise a vers le bas
	subi   r10,  _a_,   4
	add    r4,   _p_,   _k_	; r4 <- n
	mtctr  r4		; ctr <- n
	subfc  r3,   r3,   r3	; CA <- 1
1:
	lwzu   r3,   4(r10)
	addze. r3,   r3
	bdnzt  eq,   1b
	bne    2f		; a = BASE^n - 1 ?
	subi   r10,  _a_,   4	; si oui, a <- 0
	mtctr  r4		; ctr <- n
1:
	stwu   r3,   4(r10)
	bdnz   1b
2:

	; b <- (a-b) mod BASE^p - 1
	neg    r3,   _p_
        clrlslwi r3, r3, 27,2
	subf   r10,  r3,   _a_	; cadre r10 sur a
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _sub_, r3
	mtlr   r3
	add    r3,   _p_,  _k_	; r3 <- 4n
	slwi   r3,   r3,   2
	add    r11,  r10,  r3	; cadre r11, r12 sur b
	mr     r12,  r11
	addi   r3,   _p_,  31	; ctr <- ceil(p/32), CA <- 0
	srawi  r3,   r3,   5
	mtctr  r3
	blrl			; b <- a[0..p-1] - b - 1

	neg    r3,   _k_
        clrlslwi r3, r3, 27,2
	subf   r10,  r3,   r10	; cadre r10 sur a[p]
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _k_,   2
	add    r11,  r10,  r3	; cadre r11, r12 sur b
	mr     r12,  r11
	addi   r3,   _k_,  31	; ctr <- ceil(k/32)
	srwi   r3,   r3,   5
	mtctr  r3
	blrl			; b += a[p..n-1] + retenue

	subfe. r3,   r3,   r3	; propage la retenue
	bne    2f
	subf   r3,   _k_,  _p_
	mtctr  r3
	subi   r12,  r12,  4
1:
	lwzu   r3,   4(r12)
	addze. r3,   r3
	stw    r3,   0(r12)
	bdnzt  eq,   1b
	bne    2f
	mtctr  _p_		; si elle ressort, le rinjecte
	subi   r12,  r10,  4
	b      1b
2:

	; normalise b vers le haut
	subi   r10,  r10,  4
	mtctr  _p_
1:
	lwzu   r3,   4(r10)	; b = 0 ?
	and.   r3,  r3,  r3
	bdnzt  eq,  1b
	bne    2f
	li     r3,  -1		;  alors b <- BASE^p - 1
	mtctr  _p_
	addi   r10,  r10,  4
1:
	stwu   r3,  -4(r10)
	bdnz   1b
2:

        ; c <- (c - a) + (BASE^k + 1)*b - (BASE^(2*k) - 1) mod (BASE^q - 1)
	subf   r3,  _p_,  _k_
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   _a_	; cadre r11 sur a
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _sub_, r3
	mtlr   r3
	add    r3,   _p_,  _k_	; r3 <- 4(n+p)
	add    r3,   r3,   _p_
	slwi   r3,   r3,   2
	add    r10,  r11,  r3	; cadre r10, r12 sur c
	mr     r12,  r10
	subf   r3,   _k_,  _p_	; ctr <- ceil(q/32), CA <- 0
	addi   r3,   r3,   31
	srawi  r3,   r3,   5
	mtctr  r3
	blrl			; c <- c - a[0..q-1] - 1

	subf   r3,  _p_,   _k_
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur b
	slwi   r4,  _k_,   3
	add    r11,  r4,   r11
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _p_,  2	; cadre r10 et r12 sur c
	add    r10,  r3,   r11
	mr     r12,  r10
	subf   r3,   _k_,  _p_	; ctr <- ceil(q/32)
	addi   r3,   r3,   31
	srwi   r3,   r3,   5
	mtctr  r3
	blrl			; c <- c + b[0..q-1] + retenue

	neg    r3,   _k_
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur b[q]
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _k_,  2	; cadre r10 et r12 sur c
	add    r10,  r3,   r11
	mr     r12,  r10
	addi   r3,   _k_,  31	; ctr <- ceil(k/32)
	srwi   r3,   r3,   5
	mtctr  r3
	blrl			; c <- c + b[q..p-1] + retenue

	subf   r4,   _k_,  _p_	; r4 <- q-k
	subf   r4,   _k_,  r4
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r10,  r3,   r10	; cadre r10 et r12 sur c
	mr     r12,  r10
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	add    r3,   _p_,  _k_	; cadre r11 sur b
	slwi   r3,   r3,   2
	subf   r11,  r3,   r10
	addi   r4,   r4,   31	; ctr <- ceil((q-k)/32)
	srwi   r4,   r4,   5
	mtctr  r4
	blrl			; c <- c + BASE^k*b[0..q-k-1] + retenue

	add    r4,   _k_,  _k_
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur b[q-k]
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _k_,  3	; cadre r10 et r12 sur c
	add    r10,  r3,   r11
	mr     r12,  r10
	addi   r4,   r4,  31	; ctr <- ceil((2k)/32)
	srwi   r4,   r4,   5
	mtctr  r4
	blrl			; c <- c + b[q-k..p-1] + retenue
	subfe  _r_,  _r_,  _r_	; r <- retenue - 1

	add    r4,   _k_,  _k_
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r10,  r3,   r11	; cadre r10, r12 sur c
	mr     r12,  r10
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _sub_, r3
	mtlr   r3
	add    r3,   r4,   _p_	; cadre r11 sur a[q]
	slwi   r3,   r3,   2
	subf   r11,  r3,   r10
	addi   r4,   r4,   31	; ctr <- ceil((2k)/32)
	srwi   r4,   r4,   5
	mtctr  r4
	subfc  r3,   r3,   r3	; CA <- 1
	blrl			; c <- c - a[q..n-1]
	addme. _r_,  _r_	; r -= retenue
	slwi   r3,   _p_,  2	; r10 <- &c
	add    r10,  r3,  r11
	beq    2f

	subi   r12,  r12,  4	; propage la retenue
	subf   r3,   _k_,  _p_	; ctr <- q-2k
	subf   r3,   _k_,  r3
	subf   r3,   _k_,  r3
	mtctr  r3
1:
	lwzu   r3,   4(r12)
	addc   r3,   _r_,  r3
	subfe. _r_,  _r_,  _r_
	stw    r3,   0(r12)
	bdnzf  eq,   1b
	bne    4f		; si la retenue sort, c est fini
2:
	subi   r12,  r10,  4	; sinon, ajoute 1
	subfc  r3,   _k_, _p_	; ctr <- q, CA <- 1
	mtctr  r3
3:
	lwzu   r3,   4(r12)
	addze. r3,   r3
	stw    r3,   0(r12)
	bdnzt  eq,   3b
	beq    2b
4:

	; c = 0 mod BASE^q - 1 ?
	mr     r11,  r10	; r11 <- &c
	subf   r3,   _k_,  _p_
	subi   r3,   r3,   1	; ctr <- q-1
	mtctr  r3
	lwz    r3,   0(r10)	; teste si c[0] = 0 ou BASE-1
	addi   r4,   r3,   1
	srwi.  r4,   r4,   1
	bne    L(c_non_nul)
1:
	lwzu   r4,   4(r10)	; si oui, compare les autres chiffres  c[0]
	cmpw   cr0,  r3,   r4
	bdnzt  eq,  1b
	bne    L(c_non_nul)
	and.   r3,   r3,   r3
	bne    2f
	li     r3,   -1		;  si c = 0, c <- BASE^q - 1
	subf   r4,   _k_,  _p_
	mtctr  r4
	addi   r10,  r10,  4
1:
	stwu   r3,  -4(r10)
	bdnz   1b
2:
	add    r3,   _p_, _k_	; r12 <- &b
	slwi   r4,   r3,   2
	add    r12,  r4,  _a_
	subfc  r3,  _k_,  _p_	; ctr <- p+q, CA <- 1
	add    r3,  r3,   _p_
	mtctr  r3
	b      L(inc_b)		; c:b += 1

        ; si c <> 0 mod (BASE^q - 1),
        ; alors c:b <- b + (BASE^p - 1)*(c/(1-BASE^(2*k)) - 1) + BASE^q
L(c_non_nul):
	subf   r4,   _k_,  _p_	; r4 <- q-2k
	subf   r4,   _k_,  r4
	subf   r4,   _k_,  r4
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur c
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _k_,  3	; cadre r10 et r12 sur c[2k]
	add    r10,  r3,   r11
	mr     r12,  r10
	addi   r4,   r4,  31	; ctr <- ceil((q-2k)/32), CA <- 0
	srawi  r4,   r4,   5
	mtctr  r4
	blrl			; c <- c/(1 - BASE^(2k))

	subf   r3,   _k_,  _p_	; r11 <- &c
	slwi   r3,   r3,   2
	subf   r11,  r3,   r10
	subi   r10,  r11,  4	; c <- c - 1
1:
	lwzu   r3,   4(r10)
	and.   r3,   r3,   r3
	subi   r3,   r3,   1
	stw    r3,   0(r10)
	beq    1b

	subf   r3,  _p_,  _k_
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur c
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _sub_, r3
	mtlr   r3
	slwi   r3,   _p_,  2
	subf   r10,  r3,   r11	; cadre r10, r12 sur b
	mr     r12,  r10
	subfc  r3,   _k_,  _p_	; ctr <- ceil(q/32), CA <- 1
	addi   r3,   r3,   31
	srwi   r3,   r3,   5
	mtctr  r3
	blrl			; b <- b - c

	subfe. r3,   r3,   r3	; s il y a retenue, n ajoute pas BASE^q
	bne    2f
# passage non test
#ifdef debug_sjoin
	TRACE("sjoin3, checkpoint 1\n")
#endif
	mtctr  _p_
L(inc_b):
	subi   r12,  r12,  4
1:
	lwzu   r3,   4(r12)
	addze. r3,   r3
	stw    r3,   0(r12)
	bdnzt  eq,   1b
2:
        ; c:b:a <- a + (BASE^n - 1)*(c:b)/(1-BASE^k)
	add    r3,   _p_,  _k_	; r11 <- &b
	slwi   r3,   r3,   2
	add    r11,  r3,   _a_
	subf   r4,   _k_,  _p_	; r4 <- 2q
	add    r4,   r4,   r4
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r11,  r3,   r11	; cadre r11 sur b
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _add_, r3
	mtlr   r3
	slwi   r3,   _k_,  2	; cadre r10, r12 sur b[k]
	add    r10,  r3,   r11
	mr     r12,  r10
	addi   r4,   r4,   31	; ctr <- ceil((2q)/32), CA <- 0
	srawi  r4,   r4,   5
	mtctr  r4
	blrl			; c:b <- c:b/(1 - BASE^k)

	subfc  r4,   _k_,  _p_	; r4 <- p+q, CA <- 1
	add    r4,   r4,   _p_
	neg    r3,   r4
        clrlslwi r3, r3, 27,2
	subf   r10,  r3,   _a_	; cadre r10, r12 sur a
	mr     r12,  r10
	slwi   r3,   r3,   2	; lr <- point d entre dans la boucle
	add    r3,  _sub_, r3
	mtlr   r3
	add    r3,   _p_,  _k_	; cadre r11 sur b
	slwi   r3,   r3,   2
	add    r11,  r3,   r10
	addi   r4,   r4,   31	; ctr <- ceil((2q)/32)
	srwi   r4,   r4,    5
	mtctr  r4
	blrl			; c:b:a <- c:b:a - c:b

	subfe. r3,   r3,   r3	;  propage la retenue
	beq    2f
# passage non test
#ifdef debug_sjoin
	TRACE("sjoin3, checkpoint 2\n")
#endif
	subi   r12,  r12,  4
1:
	lwzu   r3,   4(r12)
	and.   r3,   r3,   r3
	subi   r3,   r3,   1
	stw    r3,   0(r12)
	beq    1b
2:

	; termin
	mtlr   r0
	blr

	#undef _a_
	#undef _p_
	#undef _k_
	#undef _add_
	#undef _sub_
	#undef _r_
	
#undef L
#endif /* assembly_sn_sjoin3 */