// file kernel/n/alpha/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 n > 0
   #
   # sortie :
   # b <- (a - b) mod (BASE^n - 1), non normalis

#ifdef assembly_sn_ssub
#undef L
#define L(x) .Lsn_ssub_##x

        .align 5
        .globl sn_ssub
        .ent   sn_ssub
sn_ssub:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):

	bis     $18, $18, $5	# r5 <- &b
	
	# si la <= n, b <- a + BASE^n - 1 - b
	cmpule  $17, $19, $0
	beq     $0,  L(big)
	addq    $31, 1,   $0	# r0 <- 1 (retenue)
	beq     $17, 1f		# si la = 0, rien  ajouter
	subq    $31, $17, $2	# r2 <- -la
	and     $2,  31,  $3	# r3 <- (-la) mod 32
	bic     $2,  31,  $2	# r2 <- -32*ceil(la/32)
	sll     $3,  3,   $3	# r3 <- 8*((-la) mod 32)
	lda     $27, sn_subloop
	subq    $16, $3,  $16	# cadre les pointeurs
	subq    $18, $3,  $18
	bis     $18, $18, $20
	s4addq  $3,  $27, $27
	jsr     $27, ($27)	# b[0..la-1] <- a - b[0..la-1] - 1
1:
	subq    $19, $17, $2	# r2 <- n - la
2:
	ldq     $1,  0($18)	# b[la..n-1] <- - b[la..n-1] - ret
	addq    $1,  $0,  $1
	cmpult  $1,  $0,  $0
	subq    $31, $1,  $1
	stq     $1,  0($18)
	cmpult  $31, $1,  $1
	addq    $1,  $0,  $0
	lda     $2,  -1($2)
	lda     $18, 8($18)
	bne     $2, 2b

	# recycle la retenue
	xor     $0,  1,   $0
	br      $31, L(cycle)

	# cas la > n : cumule par blocs de n chiffres
	.align 5
L(big):
	subq    $31, $19, $2
	and     $2,  31,  $4
	bic     $2,  31,  $2	# r2 <- -32*ceil(n/32)
	sll     $4,  3,   $4	# r4 <- 8*((-n) mod 32
	lda     $27, sn_subloop
	lda     $28, sn_addloop
	subq    $16, $4,  $16	# cadre les pointeurs
	subq    $18, $4,  $18
	bis     $18, $18, $20
	s4addq  $4,  $27, $27	# r27 <- adresse de saut pour subloop
	s4addq  $4,  $28, $28	# r28 <- adresse de saut pour addloop
	bis     $2,  $2,  $6	# r6 <- -32*ceil(n/32)
	subq    $17, $19, $17	# la -= 2n
	subq    $17, $19, $17
	addq    $31, 1,   $0	# r0 <- 1 (retenue)
	jsr     $27, ($27)	# b <- a[0..n-1] - b - 1
	xor     $0,  1,   $0	# complmente la retenue
	blt     $17, 2f
1:
	subq    $17, $19, $17	# la -= n
	subq    $16, $4,  $16	# recadre les pointeurs
	subq    $5,  $4,  $18
	subq    $5,  $4,  $20
	bis     $6,  $6,  $2	# r2 <- -32*ceil(n/32)
	jsr     $27, ($28)	# b += a[i..i+n-1] + ret
	bge     $17, 1b
2:

	# ajoute le dernier bloc
	addq    $17, $19, $2	# r2 <- -la
	subq    $31, $2,  $2
	beq     $2, L(cycle)
	and     $2,  31,  $4
	sll     $4,  3,   $4	# r4 <- 8*((-la) mod 32
	bic     $2,  31,  $2	# r2 <- -32*ceil(la/32)
	lda     $28, sn_addloop
	s4addq  $4,  $28, $28	# r28 <- adresse de saut pour addloop
	subq    $16, $4,  $16	# cadre les pointeurs
	subq    $5,  $4,  $18
	subq    $5,  $4,  $20
	jsr     $27, ($28)	# b += a[i..la-1] + ret
	subq    $31, $17, $2	# r2 <- n-la
	bsr     $27, sn_incloop # propage la retenue

	# recycle la dernire retenue
L(cycle):
	bis     $5,  $5,  $20	# r20 <- &b
	bis     $19, $19, $2	# r2  <- n
	bsr     $27, sn_incloop # b += ret
	beq     $0,  1f
	stq     $0,  0($5)	# si la retenue ressort, b <- 1
1:
	ret     $31, ($26),1

        .end   sn_ssub
#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

#define L(x) .Lsn_sred_##x

        .align 5
        .globl sn_sred
        .ent   sn_sred
sn_sred:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):

	# si la <= n, copie a dans b et complte par des zros
	cmpule  $17, $19, $0
	beq     $0,  L(big)
	bis     $18, $18, $20	# r20 <- &b
	beq     $17, 1f		# si la = 0, rien  copier
	subq    $31, $17, $2	# r2 <- -la
	subq    $19, $17, $19	# n <- n - la
	and     $2,  31,  $3	# r3 <- (-la) mod 32
	bic     $2,  31,  $2	# r2 <- -32*ceil(la/32)
	sll     $3,  3,   $3	# r3 <- 8*((-la) mod 32)
	lda     $27, sn_cpuploop
	subq    $16, $3,  $16	# cadre les pointeurs
	subq    $20, $3,  $20
	addq    $27, $3,  $27
	jsr     $27, ($27)	# copie a
	beq     $19, 2f
1:
	lda     $19, -1($19)
	stq     $31, 0($20)	# complte par des zros
	lda     $20, 8($20)
	bne     $19, 1b
2:
	ret     $31, ($26),1

	# cas la > n : cumule par blocs de n chiffres
	.align 5
L(big):
	subq    $31, $19, $2
	and     $2,  31,  $4
	bic     $2,  31,  $2	# r2 <- -32*ceil(n/32)
	sll     $4,  3,   $4	# r4 <- 8*((-n) mod 32
	lda     $27, sn_cpuploop
	lda     $28, sn_addloop
	subq    $16, $4,  $16	# cadre les pointeurs
	subq    $18, $4,  $20
	addq    $27, $4,  $27	# r27 <- adresse de saut pour cpuploop
	s4addq  $4,  $28, $28	# r28 <- adresse de saut pour addloop
	bis     $18, $18, $5	# r5 <- &b
	bis     $2,  $2,  $6	# r6 <- -32*ceil(n/32)
	bis     $31, $31, $0	# r0 <- 0 (retenue  recycler)
	subq    $17, $19, $17	# la -= 2n
	subq    $17, $19, $17
	jsr     $27, ($27)	# b <- a[0..n-1]
	blt     $17, 2f
1:
	subq    $17, $19, $17	# la -= n
	subq    $16, $4,  $16	# recadre les pointeurs
	subq    $5,  $4,  $18
	subq    $5,  $4,  $20
	bis     $6,  $6,  $2	# r2 <- -32*ceil(n/32)
	jsr     $27, ($28)	# b += a[i..i+n-1] + ret
	bge     $17, 1b
2:

	# ajoute le dernier bloc
	addq    $17, $19, $2	# r2 <- -la
	subq    $31, $2,  $2
	beq     $2, L(cycle)
	and     $2,  31,  $4
	sll     $4,  3,   $4	# r4 <- 8*((-la) mod 32
	bic     $2,  31,  $2	# r2 <- -32*ceil(la/32)
	lda     $28, sn_addloop
	s4addq  $4,  $28, $28	# r28 <- adresse de saut pour addloop
	subq    $16, $4,  $16	# cadre les pointeurs
	subq    $5,  $4,  $18
	subq    $5,  $4,  $20
	jsr     $27, ($28)	# b += a[i..la-1] + ret
	subq    $31, $17, $2	# r2 <- n-la
	bsr     $27, sn_incloop # propage la retenue

	# recycle la dernire retenue
L(cycle):
	bis     $5,  $5,  $20	# r20 <- &b
	bis     $19, $19, $2	# r2  <- n
	bsr     $27, sn_incloop # b += ret
	beq     $0,  1f
	stq     $0,  0($5)	# si la retenue ressort, b <- 1
1:
	ret     $31, ($26),1

        .end   sn_sred
#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

#define L(x) .Lsn_split_even_##x

        .align 5
        .globl sn_split_even
        .ent   sn_split_even
sn_split_even:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):

	subq   $31, $19, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2	# r2 <- -32*ceil(p/32)
	sll    $3,  3,   $3	# r3 <- 8*((-p) mod 32)
	subq   $16, $3,  $16	# cadre les pointeurs
	subq   $17, $3,  $20
	subq   $18, $3,  $21
	s8addq $19, $16, $18
	lda    $27, sn_addsubloop
	s8subq $3,  $3,  $3	# r3 <- 56*((-p) mod 32)
	addq   $27, $3,  $27	# r27 <- adresse de saut dans addsubloop
	bis    $31, $31, $0	# init retenues
	bis    $31, $31, $1
	jsr    $27, ($27)	# b <- a0+a1, c <- a0-a1
	bis    $1,  $1,  $3	# r3 <- retenue sortant de c
	sll    $19, 3,   $19	# r19 <- 8*p
	subq   $20, $19, $20	# r20 <- &b
	subq   $31, 1,   $2	# recycle la retenue sortant de b
	bsr    $27, sn_incloop	# (ne peut pas dborder)
	subq   $21, $19, $20	# r20 <- &c
	srl    $19, 3,   $2	# r2  <- p
	bis    $3,  $3,  $0	# r0 <- retenue sortant de c
	bsr    $27, sn_incloop	# recycle la retenue sortant de c
	stq    $0,  0($21)	# c[p] <- retenue sortante
	ret    $31, ($26),1	

        .end   sn_split_even
#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, normalis
   #   c  <- a mod BASE^(p+1/2) + 1, normalis

#define L(x) .Lsn_split_odd_##x

        .align 5
        .globl sn_split_odd
        .ent   sn_split_odd
sn_split_odd:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):

	s8addq $19, $16, $20	# r20 <- &a0[p]
	bis    $19, $19, $2	# r2 <- p
	bis    $31, $31, $0	# init retenues
	bis    $31, $31, $1
	ldq    $5,  0($20)
	bis    $5,  $5,  $8	# r8 <- high(a0):low(a1)
	bis    $17, $17, $21	# r21 <- &b
	bis    $18, $18, $22	# r22 <- &c

	.align 5
1:
	lda    $20, 8($20)
	srl    $5,  32, $3	# r3 <- low(a1[i])
	lda    $2,  -1($2)
	ldq    $4,  0($16)	# r4 <- a0[i]
	ldq    $5,  0($20)	# r5 <- high(a1[i]):low(a1[i+1])
	sll    $5,  32, $6
	addq   $6,  $3, $3	# r3 <- a1[i]
	addq   $4,  $0, $0	# r0 <- a0[i] + ret(a0+a1)
	addq   $3,  $1, $1	# r1 <- a1[i] + ret(a0-a1)
	cmpult $0,  $4, $6
	cmpult $1,  $3, $7
	addq   $0,  $3, $0	# r0 <- a0[i] + a1[i] + ret(a0+a1)
	subq   $4,  $1, $1	# r1 <- a0[i] - a1[i] - ret(a0-a1)
	stq    $0,  0($17)	# sauve b[i]
	stq    $1,  0($18)	# sauve c[i]
	cmpult $0,  $3, $0
	cmpult $4,  $1, $1
	addq   $0,  $6, $0	# r0 <- ret(a0+a1)
	addq   $1,  $7, $1	# r1 <- ret(a0-a1)
	lda    $16, 8($16)	# avance les pointeurs
	lda    $17, 8($17)
	lda    $18, 8($18)
	bne    $2,  1b

	# dernier demi-chiffre
	sll    $8,  32,  $8	# r8 <- high(a0)
	srl    $8,  32,  $8
	srl    $5,  32,  $5	# r5 <- high(a1)
	addq   $8,  $0,  $0	# r0 <- high(a0) + ret(a0+a1)
	addq   $5,  $1,  $1	# r1 <- high(a1) + ret(a0-a1)
	addq   $0,  $5,  $3	# r3 <- high(a0) + high(a1) + ret(a0+a1)
	subq   $8,  $1,  $4	# r4 <- high(a0) - high(a1) - ret(a0-a1)
	sll    $3,  32,  $0	# r0 <- r3 mod 2^32
	srl    $0,  32,  $0
	sll    $4,  32,  $1	# r1 <- r4 mod 2^32
	srl    $1,  32,  $1
	stq    $0,  0($17)	# sauve b[p]
	stq    $1,  0($18)	# sauve c[p]

	# recycle la retenue sortant de b
	srl    $3,  32,  $0	# r0 <- retenue
1:
	bis    $21, $21, $20	# r20 <- &b
	bis    $19, $19, $2	# r2  <- p
	bsr    $27, sn_incloop
	bne    $0,  1b
	
	# recycle la retenue sortant de c
	srl    $4,  63,  $0	# r0 <- retenue
	bis    $22, $22, $20	# r20 <- &c
	subq   $31, 1,   $2
	bis    $26, $26, $27	# retour terminal
	br     $31, sn_incloop

        .end   sn_split_odd
#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

        .align 5
#ifdef debug_smul
        .globl sn_smul_buggy
        .ent   sn_smul_buggy
sn_smul_buggy:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
#else
        .globl sn_smul
        .ent   sn_smul
sn_smul:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):
#endif

	#define _d_   48($30)
	#define _b_    0($30)
	#define _lb_   8($30)
	#define _c_   16($30)
	#define _n_   24($30)
	#define _ra_  32($30)

	# rserve 2n+6. chiffres dans la pile
	s4addq $21, 0,   $0
	s4addq $0,  48,  $0
	subq   $30, $0,  $30
	stq    $18, _b_
	stq    $19, _lb_
	stq    $20, _c_
	stq    $21, _n_
	stq    $26, _ra_
	
	# d0 <- a mod (BASE^n - 1)
	lda    $18, _d_
	bis    $21, $21, $19
	bsr    $26, .Lsn_sred_nogp
	
	# d1 <- b mod (BASE^n - 1)
	ldq    $16, _b_
	ldq    $17, _lb_
	lda    $18, _d_
	ldq    $19, _n_
	s8addq $19, $18, $18
	bsr    $26, .Lsn_sred_nogp

	# c <- a*b mod (BASE^n - 1)
	lda    $16, _d_
	ldq    $17, _c_
	ldq    $18, _n_
	bsr    $27, L(recurse)
	ldq    $26, _ra_	# rcupre l adresse de retour
	lda    $16, _d_		# nettoie la pile
	ldq    $17, _n_
	s8addq $17, $16, $16
	s8addq $17, $16, $30
	ret    $31, ($26),1

	#undef _d_
	#undef _b_
	#undef _lb_
	#undef _c_
	#undef _n_
	#undef _ra_
	
   # entre rcursive
   # r16 = &a (b = a+n)
   # r17 = &c
   # r18 = n
   # r27 = ra
L(recurse):
	
	# aiguillage selon la longueur et la parit
	blbs   $18, 1f
	cmpule $18, smul_lim_even, $0
	beq    $0,  L(big_even)
	br     $31, L(small)
	.align 5
1:
	cmpule $18, smul_lim_odd, $0
	beq    $0,  L(big_odd)
	
	# petite multiplication => Toom puis rduction
L(small):
	#define _d_   32($30)
	#define _a_    0($30)
	#define _c_    8($30)
	#define _n_   16($30)
	#define _ra_  24($30)

	# rserve 2n+4. chiffres dans la pile
	s4addq $18, 0,   $0
	s4addq $0,  32,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	stq    $18, _n_
	stq    $27, _ra_

	# d1:d2 <- a*b mod (BASE^n - 1)
	bis    $18, $18, $17
	bis    $18, $18, $19
	s8addq $18, $16, $18
	lda    $20, _d_
	bsr    $26, .Lsn_toommul_nogp

	# point d entre pour ssqr
L(small_aux):
	
	# c <- d1 + d2
	lda    $16, _d_
	ldq    $0,  _n_
	ldq    $20, _c_
	subq   $31, $0,  $2
	and    $2,  31,  $3
	bic    $2,  31,  $2	# r2 <- 32*ceil(n/32)
	sll    $3,  3,   $3	# r3 <- 8*((-n) mod 32)
	lda    $27, sn_addloop	# cadre les pointeurs
	subq   $16, $3,  $16
	subq   $20, $3,  $20
	s4addq $3,  $27, $27
	s8addq $0,  $16, $18
	bis    $31, $31, $0	# init retenue
	jsr    $27, ($27)	# c <- d1 + d2
	ldq    $27, _ra_	# rcupre l adresse de retour
	ldq    $20, _c_
	bis    $18, $18, $30	# nettoie la pile
	subq   $31, 1,   $2
	br     $31, sn_incloop	# recycle la retenue

	#undef _d_
	#undef _a_
	#undef _c_
	#undef _n_
	#undef _ra_

	# cas n grand pair
	.align 5
L(big_even):

	#define _d_   32($30)
	#define _a_    0($30)
	#define _r_    0($30)
	#define _c_    8($30)
	#define _p_   16($30)
	#define _ra_  24($30)

	# rserve 3p+5. chiffres dans la pile
	s4subq $18, 0,   $0
	s4addq $0,  48,  $0
	bic    $0,  15,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	srl    $18, 1,   $19	# r19 <- p
	stq    $19, _p_
	stq    $27, _ra_

	# dcompose a et b modulo BASE^p +/- 1
	bis    $17, $17, $18	# r18 <- &c
	lda    $17, _d_
	bsr    $26, .Lsn_split_even_nogp
	ldq    $16, _a_
	lda    $17, _d_
	ldq    $19, _p_
	s8addq $19, $16, $16	# r16 <- &b
	s8addq $19, $16, $16
	s8addq $19, $17, $17	# r17 <- &d1
	s8addq $19, $17, $18	# r18 <- &d2
	bsr    $26, .Lsn_split_even_nogp

	# c[0..p] <- a*b mod BASE^p + 1
	ldq    $16, _c_
	lda    $17, _d_
	ldq    $18, _p_
	s8addq $18, $17, $17
	s8addq $18, $17, $17	# r17 <- &d2
	bsr    $26, .Lsn_mmul_nogp

	# c[p..2p-1] <- a*b mod BASE^p - 1
	lda    $16, _d_
	ldq    $17, _c_
	ldq    $18, _p_
	s8addq $18, $17, $17	# r17 <- &c1
	ldq    $0,  0($17)	# sauve c0[p]
	stq    $0,  _r_
	bsr    $27, L(recurse)

	# point d entre pour ssqr
L(even_aux):
	
	# force c0+c1 pair en ajoutant BASE^p+1  c0:r si ncessaire
	# la version assembleur de sn_mmul produit un chiffre de
	# poids fort valant 0 ou 1, donc il ne peut pas y avoir
	# dbordement sur r.
	ldq    $17, _c_		# r17 <- &c0
	ldq    $20, _p_
	s8addq $20,  $17, $16	# r16 <- &c1
	ldq    $0, 0($16)
	ldq    $1, 0($17)
	addq   $0, $1, $0
	blbc   $0, 3f
	ldq    $2, _r_
	addq   $1, 1,  $1	# ajoute 1  c0 et  r
	addq   $2, 1,  $2
	stq    $1, 0($17)	# sauve c0[0]
	cmpult $1, 1,  $0
	beq    $0, 2f		# propage la retenue
	lda    $18, 8($17)
	lda    $19, -1($20)
1:
	ldq    $0,  0($18)
	addq   $0,  1,   $0
	stq    $0,  0($18)
	bne    $0,  2f
	lda    $18, 8($18)
	lda    $19, -1($19)
	bne    $19, 1b
	addq   $2,  1,  $2
2:
	stq    $2, _r_
3:
	# c0 <- (c0+c1)/2, c1 <- (c1-c0)/2
	bis    $17, $17, $18	# r18 <- &c0
	bis    $16, $16, $19	# r19 <- &c1
	bsr    $26, .Lsn_half_addsub_nogp

	# retenues  propager : r = u + 2*v
	# r*BASE^p/2
	# -r*BASE^(2p)/2
	# -BASE^2p si c1[p-1] >= BASE/2 (ie. c1 < c0)

	# ajoute r*B^p/2
	ldq    $3,  _r_
	srl    $3,  1,   $4	# r4 <- v
	sll    $3,  63,  $3	# r3 <- u*BASE/2
	ldq    $5,  _p_
	ldq    $16, _c_
	s8addq $5,  $16, $20	# r20 <- &c1
	ldq    $0,  -8($20)	# ajoute u  c0[p-1] = c1[-1]
	addq   $0,  $3,  $0
	stq    $0,  -8($20)
	cmpult $0,  $3,  $0
	addq   $0,  $4,  $0	# r0 <- v + ret
	bis    $5,  $5,  $2
	bsr    $27, sn_incloop	# propage sur c1

	# corrige le dernier chiffre de c1
	s8addq $5,  $16, $20
	s8addq $5,  $20, $20
	ldq    $0,  -8($20)
	srl    $0,  63,  $1
	addq   $1,  $4,  $4	# r4 <- v + signe(c1[p-1])
	subq   $0,  $3,  $3
	stq    $3,  -8($20)	# c1[p-1] -= u
	cmpult $0,  $3,  $0
	addq   $0,  $4,  $0	# r0 <- retenue ngative sortant de c1

	# recycle la retenue sortant de c1
1:
	bis    $16, $16, $20	# r20 <- &c
	addq   $5,  $5,  $2	# r2 <- 2p
	bsr    $27, sn_decloop	# c -= ret
	bne    $0,  1b

	# nettoie la pile
	ldq    $27, _ra_
	s4subq $5,  0,   $0
	s8addq $0,  48,  $0
	bic    $0,  15,  $0
	addq   $30, $0,  $30
	ret    $31, ($27),1
	
	#undef _d_
	#undef _a_
	#undef _r_
	#undef _c_
	#undef _p_
	#undef _ra_

	# cas n grand impair
	.align 5
L(big_odd):

	#define _d_   32($30)
	#define _a_    0($30)
	#define _c_    8($30)
	#define _p_   16($30)
	#define _ra_  24($30)

	# rserve 6p+10. chiffres dans la pile
	s4subq $18, $18, $0
	s8addq $0,  56,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	srl    $18, 1,   $19	# r19 <- p
	stq    $19, _p_
	stq    $27, _ra_

	# dcompose a et b modulo BASE^(p+1/2) +/- 1
	lda    $17, _d_
	s8addq $18, $16, $16	# r16 <- &b
	s8addq $18, $17, $18	# r18 <- &d2
	lda    $18, 8($18)
	bsr    $26, .Lsn_split_odd_nogp
	ldq    $16, _a_
	lda    $17, _d_
	ldq    $19, _p_
	addq   $19, 1,   $0	# r0  <- p+1
	s8addq $0,  $17, $17	# r17 <- &d1
	s8addq $0,  $17, $18
	s8addq $0,  $18, $18	# r18 <- &d3
	bsr    $26, .Lsn_split_odd_nogp

	# d4:d5 <- a*b mod BASE^(p+1/2) + 1
	lda    $16, _d_
	ldq    $17, _p_
	addq   $17, 1,   $17	# r17 <- p+1
	s8addq $17, $16, $16	# r16 <- &d2
	s8addq $17, $16, $16
	s8addq $17, $16, $18	# r18 <- &d3
	bis    $17, $17, $19	# r19 <- p+1
	s8addq $17, $18, $20	# r20 <- &d4
	bsr    $26, .Lsn_toommul_nogp

	# d2:d3 <- a*b mod BASE^(p+1/2) - 1
	lda    $16, _d_
	ldq    $17, _p_
	addq   $17, 1,   $17	# r17 <- p+1
	s8addq $17, $16, $18	# r18 <- &d1
	bis    $17, $17, $19	# r19 <- p+1
	s8addq $17, $18, $20	# r20 <- &d2
	bsr    $26, .Lsn_toommul_nogp
	
	# point d entre pour ssqr
L(odd_aux):
	
	# c <- d2:d3 + d4:d5 mod (BASE^n - 1)
	# d <- d2:d3 - d4:d5 mod (BASE^n - 1)
	ldq    $17, _p_
	addq   $17, 1,   $17
	sll    $17, 1,   $17	# r17 <- 2*p+2
	lda    $21, _d_
	s4subq $17, $17, $19
	s8addq $19, $21, $16	# r16 <- &d6
	ldq    $0,  -8($16)	# r0 <- d5[2p+1] (retenue pour addsub)
	bis    $0,  $0,  $1	# r1 <- d2[2p+1]
	ldq    $20, _c_
	addq   $31, 1,   $2
	subq   $2,  $17, $2	# r2 <- -(2p+1)
	and    $2,  31,  $3
	bic    $2,  31,  $2	# r2 <- 32*ceil((2p+1)/32)
	sll    $3,  3,   $3	# r3 <- 8*(-(2p+1) mod 32)
	lda    $27, sn_addsubloop
	subq   $20, $3,  $20	# cadre les pointeurs
	subq   $21, $3,  $21
	s8subq $3,  $3,  $3
	addq   $27, $3,  $27
	s8addq $17, $21, $16
	s8addq $17, $16, $18
	jsr    $27, ($27)	# effectue l addition/soustraction
	bis    $1,  $1,  $3	# r3 <- retenue sortant de d
1:
	ldq    $20, _c_		# recycle la retenue sortant de c
	subq   $17, 1,   $2
	bsr    $27, sn_incloop
	bne    $0,  1b
	bis    $3,  $3,  $0
1:
	lda    $20, _d_		# recycle la retenue sortant de d
	subq   $17, 1,   $2
	bsr    $27, sn_decloop
	bne    $0,  1b
	
	# c <- (c + d*BASE^(p+1/2))/2 mod BASE^(2p+1)-1
	ldq    $16, _c_
	bis    $16, $16, $21	# r21 <- &c
	lda    $18, _d_
	s4addq $17, $18, $19	# r19 <- &d[p+1]
	subq   $17, 4,   $2	# r2  <- 2p-2

	# chiffre 0
	ldq    $1,  -8($19)
	srl    $1,  32,  $1	# r1 <- low(d1[0])
	ldq    $4,  0($16)	# r4 <- c[0]
	ldq    $5,  0($19)
	sll    $5,  32,  $6
	addq   $6,  $1,  $1	# r1 <- d1[0]
	addq   $4,  $1,  $4	# r4 <- (c+d1)[i]
	cmpult $4,  $1,  $0	# r0 <- retenue
	and    $4,  1,   $8	# r8 <- (c+d1)[0] mod 2
	srl    $4,  1,   $3	# r3 <- ((c+d1)/2)[0] mod BASE/2

	# chiffres 1..p-1
	.align 5
1:
	lda    $19, 8($19)	# d1++
	lda    $2,  -2($2)
	srl    $5,  32,  $1	# r1 <- low(d1[i])
	ldq    $4,  8($16)	# r4 <- c[i]
	ldq    $5,  0($19)
	sll    $5,  32,  $6
	addq   $6,  $1,  $1	# r1 <- d1[i]
	addq   $4,  $0,  $4
	cmpult $4,  $0,  $0
	addq   $4,  $1,  $4	# r4 <- (c+d1)[i]
	cmpult $4,  $1,  $1
	addq   $1,  $0,  $0	# r0 <- retenue
	sll    $4,  63,  $6
	addq   $6,  $3,  $3	# r3 <- ((c+d1)/2)[i-1]
	stq    $3,  0($16)	# sauve c[i-1]
	lda    $16, 8($16)	# c++
	srl    $4,  1,   $3	# r3 <- ((c+d1)/2)[i] mod BASE/2
	bne    $2,  1b
	
	# chiffres p..2p
	bis    $17, $17, $2	# r2  <- 2p+2
	.align 5
1:
	lda    $2,  -2($2)
	srl    $5,  32,  $1	# r1 <- low(d0[i])
	ldq    $4,  8($16)	# r4 <- c[i]
	ldq    $5,  0($18)
	sll    $5,  32,  $6
	addq   $6,  $1,  $1	# r1 <- d0[i]
	addq   $4,  $0,  $4
	cmpult $4,  $0,  $0
	addq   $4,  $1,  $4	# r4 <- (c+d0)[i]
	cmpult $4,  $1,  $1
	addq   $1,  $0,  $0	# r0 <- retenue
	sll    $4,  63,  $6
	addq   $6,  $3,  $3	# r3 <- ((c+d0)/2)[i-1]
	stq    $3,  0($16)	# sauve c[i-1]
	lda    $18, 8($18)	# d0++
	lda    $16, 8($16)	# c++
	srl    $4,  1,   $3	# r3 <- ((c+d0)/2)[i] mod BASE/2
	bne    $2,  1b

	# dernier chiffre
	addq   $0,  $8,  $0	# incorpore le bit des units
	sll    $0,  63,  $1
	addq   $1,  $3,  $3
	stq    $3,  0($16)	# sauve c[2p]
	srl    $0,  1,   $0	# r0 <- retenue
	subq   $17, 2,   $17	# r17 <- 2p
1:
	bis    $21, $21, $20	# r20 <- &c
	srl    $17, 1,   $2	# r2  <- p
	bsr    $27, sn_incloop	# recycle la retenue
	bne    $0,  1b

	# nettoie la pile
	ldq    $27, _ra_
	ldq    $0,  _p_
	s4subq $0,  $0,  $0
	sll    $0,  1,   $0
	addq   $0,  10,  $0
	s8addq $0,  $30, $30
	ret    $31, ($27),1
	
	#undef _d_
	#undef _a_
	#undef _c_
	#undef _p_
	#undef _ra_

#ifdef debug_smul
        .end sn_smul_buggy
#else
        .end sn_smul
#endif
#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

        .align 5
#ifdef debug_smul
        .globl sn_ssqr_buggy
        .ent   sn_ssqr_buggy
sn_ssqr_buggy:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
#else
        .globl sn_ssqr
        .ent   sn_ssqr
sn_ssqr:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):
#endif

	#define _d_   24($30)
	#define _c_    0($30)
	#define _n_    8($30)
	#define _ra_  16($30)

	# rserve n+3. chiffres dans la pile
	s8addq $19, 32,  $0
	bic    $0,  15,  $0
	subq   $30, $0,  $30
	stq    $18, _c_
	stq    $19, _n_
	stq    $26, _ra_
	
	# d <- a mod (BASE^n - 1)
	lda    $18, _d_
	bsr    $26, .Lsn_sred_nogp
	
	# c <- a^2 mod (BASE^n - 1)
	lda    $16, _d_
	ldq    $17, _c_
	ldq    $18, _n_
	bsr    $27, L(recurse)
	ldq    $26, _ra_	# rcupre l adresse de retour
	ldq    $0,  _n_		# nettoie la pile
	s8addq $0,  32,  $0
	bic    $0,  15,  $0
	addq   $30, $0,  $30
	ret    $31, ($26),1

	#undef _d_
	#undef _c_
	#undef _n_
	#undef _ra_
	
   # entre rcursive
   # r16 = &a
   # r17 = &c
   # r18 = n
   # r27 = ra
L(recurse):
	
	# aiguillage selon la longueur et la parit
	blbs   $18, 1f
	cmpule $18, ssqr_lim_even, $0
	beq    $0,  L(big_even)
	br     $31, L(small)
	.align 5
1:
	cmpule $18, ssqr_lim_odd, $0
	beq    $0,  L(big_odd)
	
	# petit carr => Toom puis rduction
L(small):
	#define _d_   32($30)
	#define _a_    0($30)
	#define _c_    8($30)
	#define _n_   16($30)
	#define _ra_  24($30)

	# rserve 2n+4. chiffres dans la pile
	s4addq $18, 0,   $0
	s4addq $0,  32,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	stq    $18, _n_
	stq    $27, _ra_

	# d1:d2 <- a^2 mod (BASE^n - 1)
	bis    $18, $18, $17
	lda    $18, _d_
	bsr    .Lsn_toomsqr_nogp
	br     $31, .Lsn_smul_small_aux # continue auvec smul

	#undef _d_
	#undef _a_
	#undef _c_
	#undef _n_
	#undef _ra_

	# cas n grand pair
	.align 5
L(big_even):

	#define _d_   32($30)
	#define _a_    0($30)
	#define _r_    0($30)
	#define _c_    8($30)
	#define _p_   16($30)
	#define _ra_  24($30)

	# rserve 3p+5. chiffres dans la pile
	s4subq $18, 0,   $0
	s4addq $0,  48,  $0
	bic    $0,  15,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	srl    $18, 1,   $19	# r19 <- p
	stq    $19, _p_
	stq    $27, _ra_

	# dcompose a modulo BASE^p +/- 1
	bis    $17, $17, $18	# r18 <- &c
	lda    $17, _d_
	bsr    $26, .Lsn_split_even_nogp

	# c[0..p] <- a^2 mod BASE^p + 1
	ldq    $16, _c_
	ldq    $17, _p_
	bsr    $26, .Lsn_msqr_nogp

	# c[p..2p-1] <- a^2 mod BASE^p - 1
	lda    $16, _d_
	ldq    $17, _c_
	ldq    $18, _p_
	s8addq $18, $17, $17	# r17 <- &c1
	ldq    $0,  0($17)	# sauve c0[p]
	stq    $0,  _r_
	bsr    $27, L(recurse)
	br     $31, .Lsn_smul_even_aux # continue avec smul
	
	#undef _d_
	#undef _a_
	#undef _r_
	#undef _c_
	#undef _p_
	#undef _ra_

	# cas n grand impair
	.align 5
L(big_odd):

	#define _d_   32($30)
	#define _a_    0($30)
	#define _c_    8($30)
	#define _p_   16($30)
	#define _ra_  24($30)

	# rserve 6p+10. chiffres dans la pile
	s4subq $18, $18, $0
	s8addq $0,  56,  $0
	subq   $30, $0,  $30
	stq    $16, _a_
	stq    $17, _c_
	srl    $18, 1,   $19	# r19 <- p
	stq    $19, _p_
	stq    $27, _ra_

	# dcompose a modulo BASE^(p+1/2) +/- 1
	lda    $17, _d_
	s8addq $18, $17, $18	# r18 <- &d2
	lda    $18, 8($18)
	bsr    $26, .Lsn_split_odd_nogp

	# d4:d5 <- a^2 mod BASE^(p+1/2) + 1
	lda    $16, _d_
	ldq    $17, _p_
	addq   $17, 1,   $17	# r17 <- p+1
	s8addq $17, $16, $16	# r16 <- &d2
	s8addq $17, $16, $16
	s8addq $17, $16, $18
	s8addq $17, $18, $18	# r18 <- &d4
	bsr    $26, .Lsn_toomsqr_nogp

	# d2:d3 <- a^2 mod BASE^(p+1/2) - 1
	lda    $16, _d_
	ldq    $17, _p_
	addq   $17, 1,   $17	# r17 <- p+1
	s8addq $17, $16, $18
	s8addq $17, $18, $18	# r18 <- &d2
	bsr    $26, .Lsn_toomsqr_nogp
	br     $31, .Lsn_smul_odd_aux # continue avec smul
	
	#undef _d_
	#undef _a_
	#undef _c_
	#undef _p_
	#undef _ra_

#ifdef debug_smul
        .end sn_ssqr_buggy
#else
        .end sn_sqr
#endif
#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
#define L(x) .Lsn_sjoin3_##x

        .align 5
#ifdef debug_sjoin
        .globl sn_sjoin3_buggy
        .ent   sn_sjoin3_buggy
sn_sjoin3_buggy:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
#else
        .globl sn_sjoin3
        .ent   sn_sjoin3
sn_sjoin3:
        .frame $30,0,$26,0
        .prologue 1
        ldgp   $gp,  0($27)
L(nogp):
#endif

	#define _a_ $21
	#define _b_ $22
	#define _c_ $23
	#define _k_ $24
	#define _q_ $25

	# sauvegarde les paramtres
	bis    $16, $16, _a_
	bis    $18, $18, _k_
	mulq   $17, $18, _q_
	sll    _q_, 1,   _q_
	s8addq _q_, _a_, _b_
	s8addq _k_, _b_, _b_
	s8addq _k_, _b_, _b_
	s8addq _q_, _b_, _c_
	s8addq _k_, _c_, _c_
	
	# normalise a
	addq   _q_, _k_, $2
	addq   $2,  _k_, $2
1:
	ldq    $1,  0($16)
	addq   $1,  1,   $1
	lda    $2,  -1($2)
	lda    $16, 8($16)
	bne    $1,  3f
	bne    $2,  1b
	addq   _q_, _k_, $2
	addq   $2,  _k_, $2
2:
	lda    $16, -8($16)
	lda    $2,  -1($2)
	stq    $31, 0($16)
	bne    $2,  2b
3:

	# b <- (a - b) mod (BASE^p - 1), normalis vers le haut
	# soustrait les p premiers chiffres
	subq   $31, _q_, $2
	subq   $2,  _k_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_subloop
	subq   _a_, $3,  $16
	subq   _b_, $3,  $18
	subq   _b_, $3,  $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)
	bis    $0,  $0,  $4	# r4 <- retenue ngative

	# ajoute les k derniers chiffres de a
	subq   $31, _k_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_addloop
	subq   $16, $3,  $16
	subq   _b_, $3,  $18
	subq   _b_, $3,  $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)

	# propage la retenue positive
	bis    _q_, _q_, $2
1:
	ldq    $1,  0($20)
	addq   $1,  $0,  $1
	cmpult $1,  $0,  $0
	stq    $1,  0($20)
	beq    $0,  2f
	lda    $2,  -1($2)
	lda    $20, 8($20)
	bne    $2,  1b
2:
	# propage la retenue ngative
	subq   $4,  $0,  $0
	beq    $0,  3f
1:
	bis    _b_, _b_, $20
	addq   _q_, _k_, $2
2:
	ldq    $1,  0($20)
	subq   $1,  1,   $0
	stq    $0,  0($20)
	lda    $2, -1($2)
	lda    $20, 8($20)
	bne    $1,  3f
	bne    $2,  2b
	br     $31, 1b
	.align 5

	# normalise b vers le haut
3:
	bis    _b_, _b_, $20
	addq   _q_, _k_, $2
1:
	ldq    $0,  0($20)
	bne    $0,  3f
	lda    $2,  -1($2)
	lda    $20, 8($20)
	bne    $2,  1b
	addq   _q_, _k_, $2
	subq   $31, 1,   $0
2:
	lda    $20, -8($20)
	lda    $2,  -1($2)
	stq    $0,  0($20)
	bne    $2,  2b
3:

	# c <- (c - a) + (BASE^k + 1)*b - (BASE^(2*k) - 1) mod (BASE^q - 1)
	# chiffres 0..k-1
	bis    _c_, _c_, $16	# r16 <- &c
	bis    _a_, _a_, $17	# r17 <- &a
	s8addq _q_, _a_, $18	# r18 <- &a[q]
	bis    _b_, _b_, $19	# r19 <- &b
	s8addq _q_, _b_, $20	# r20 <- &b[q]
	sll    _k_, 3,   $0
	subq   $20, $0,  $28	# r21 <- &b[q-k]
	bis    _k_, _k_, $2
	addq   $31, 1,   $0	# r0 <- 1 (retenue)
1:
	ldq    $1,  0($16)	# c[i]
	ldq    $3,  0($17)	# a[i]
	ldq    $4,  0($18)	# a[i+q]
	ldq    $5,  0($19)	# b[i]
	ldq    $6,  0($20)	# b[i+q]
	ldq    $7,  0($28)	# b[i+q-k]
	addq   $1,  $0,  $1	# c[i] + ret
	cmpult $1,  $0,  $8
	sra    $0,  63,  $0
	addq   $8,  $0,  $0
	cmpult $1,  $3,  $8	# c[i] - a[i] + ret
	subq   $1,  $3,  $1
	subq   $0,  $8,  $0
	cmpult $1,  $4,  $8	# c[i] - a[i] - a[i+q] + ret
	subq   $1,  $4,  $1
	subq   $0,  $8,  $0
	addq   $1,  $5,  $1	# c[i] - a[i] - a[i+q] + b[i] + ret
	cmpult $1,  $5,  $8
	addq   $0,  $8,  $0
	addq   $1,  $6,  $1	# c[i] - a[i] - a[i+q] + b[i] + b[i+q] + ret
	cmpult $1,  $6,  $8
	addq   $0,  $8,  $0
	addq   $1,  $7,  $1	# c[i] - a[i] - a[i+q] + b[i] + b[i+q] + b[i+q-k] + ret
	cmpult $1,  $7,  $8
	addq   $0,  $8,  $0
	stq    $1,  0($16)	# sauve c[i]
	lda    $2,  -1($2)	# avance les pointeurs
	lda    $16, 8($16)
	lda    $17, 8($17)
	lda    $18, 8($18)
	lda    $19, 8($19)
	lda    $20, 8($20)
	lda    $28, 8($28)
	bne    $2, 1b

	# chiffres k..2k-1
	bis    _b_, _b_, $20
	bis    _k_, _k_, $2
1:
	ldq    $1,  0($16)	# c[i]
	ldq    $3,  0($17)	# a[i]
	ldq    $4,  0($18)	# a[i+q]
	ldq    $5,  0($19)	# b[i]
	ldq    $6,  0($20)	# b[i-k]
	ldq    $7,  0($28)	# b[i+q-k]
	addq   $1,  $0,  $1	# c[i] + ret
	cmpult $1,  $0,  $8
	sra    $0,  63,  $0
	addq   $8,  $0,  $0
	cmpult $1,  $3,  $8	# c[i] - a[i] + ret
	subq   $1,  $3,  $1
	subq   $0,  $8,  $0
	cmpult $1,  $4,  $8	# c[i] - a[i] - a[i+q] + ret
	subq   $1,  $4,  $1
	subq   $0,  $8,  $0
	addq   $1,  $5,  $1	# c[i] - a[i] - a[i+q] + b[i] + ret
	cmpult $1,  $5,  $8
	addq   $0,  $8,  $0
	addq   $1,  $6,  $1	# c[i] - a[i] - a[i+q] + b[i] + b[i-k] + ret
	cmpult $1,  $6,  $8
	addq   $0,  $8,  $0
	addq   $1,  $7,  $1	# c[i] - a[i] - a[i+q] + b[i] + b[i-k] + b[i+q-k] + ret
	cmpult $1,  $7,  $8
	addq   $0,  $8,  $0
	stq    $1,  0($16)	# sauve c[i]
	lda    $2,  -1($2)	# avance les pointeurs
	lda    $16, 8($16)
	lda    $17, 8($17)
	lda    $18, 8($18)
	lda    $19, 8($19)
	lda    $20, 8($20)
	lda    $28, 8($28)
	bne    $2, 1b

	# chiffres 2k..q-1
	subq   _q_, _k_, $2
	subq   $2,  _k_, $2
	subq   $0,  1,   $0
1:
	ldq    $1,  0($16)	# c[i]
	ldq    $3,  0($17)	# a[i]
	ldq    $5,  0($19)	# b[i]
	ldq    $6,  0($20)	# b[i-k]
	addq   $1,  $0,  $1	# c[i] + ret
	cmpult $1,  $0,  $8
	sra    $0,  63,  $0
	addq   $8,  $0,  $0
	cmpult $1,  $3,  $8	# c[i] - a[i] + ret
	subq   $1,  $3,  $1
	subq   $0,  $8,  $0
	addq   $1,  $5,  $1	# c[i] - a[i] + b[i] + ret
	cmpult $1,  $5,  $8
	addq   $0,  $8,  $0
	addq   $1,  $6,  $1	# c[i] - a[i] + b[i] + b[i-k] + ret
	cmpult $1,  $6,  $8
	addq   $0,  $8,  $0
	stq    $1,  0($16)	# sauve c[i]
	lda    $2,  -1($2)	# avance les pointeurs
	lda    $16, 8($16)
	lda    $17, 8($17)
	lda    $19, 8($19)
	lda    $20, 8($20)
	bne    $2, 1b

	# recycle la retenue
1:
	bis    _c_, _c_, $16
	bis    _q_, _q_, $2
2:
	ldq    $1,  0($16)
	addq   $1,  $0,  $1
	stq    $1,  0($16)
	cmpult $1,  $0,  $1
	sra    $0,  63,  $0
	addq   $1,  $0,  $0
	beq    $0,  3f
	lda    $2,  -1($2)
	lda    $16, 8($16)
	bne    $2,  2b
	br     $31, 1b
	.align 5
3:

	# c = 0 mod (BASE^q - 1) ?
	ldq    $0,  0(_c_)
	addq   $0,  1,   $1
	cmovne $1,  $0,  $1
	bne    $1,  L(c_non_nul)
	lda    $16, 8(_c_)
	subq   _q_, 1,   $2
1:
	ldq    $1,  0($16)
	xor    $1,  $0,  $1
	bne    $1,  L(c_non_nul)
	lda    $2,  -1($2)
	lda    $16, 8($16)
	bne    $2,  1b

	# si c = 0, alors c <- BASE^q - 1
	bne    $0,  2f
	bis    _q_, _q_, $2
	subq   $31, 1,   $0
1:
	lda    $2,  -1($2)
	lda    $16, -8($16)
	stq    $0,  0($16)
	bne    $2,  1b
2:
	# c:b += 1
	bis    _b_, _b_, $20
	addq   _q_, _k_, $2
	addq   _q_, $2,  $2
	br     $31, L(inc_b)
	.align 5

	# sinon, c:b <- b + (BASE^p - 1)*(c/(1-BASE^(2k)) - 1) + BASE^q
L(c_non_nul):
	# c <- c/(1-BASE^(2k)) - 1
	addq   _k_, _k_, $2
	subq   $2,  _q_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_addloop
	subq   _c_, $3,  $16
	s8addq _k_, $16, $18
	s8addq _k_, $18, $18
	bis    $18, $18, $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)
	bis    _c_, _c_, $20
	addq   $31, 1,   $0
	bis    _q_, _q_, $2
	bsr    $27, sn_decloop

	# c:b <- c:b - c
	subq   $31, _q_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_subloop
	subq   _b_, $3,  $16
	subq   _c_, $3,  $18
	subq   _b_, $3,  $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)
	addq   _q_, _k_, $2
	bsr    $27, sn_decloop

	# c:b <- c:b + BASE^q
	s8addq _q_, _b_, $20
	addq   _k_, _q_, $2
L(inc_b):
	addq   $31, 1,   $0
	bsr    $27, sn_incloop
	
	# c:b:a <- a + (BASE^n - 1)*(c:b)/(1-BASE^k)
	subq   $31, _q_, $2
	subq   $2,  _q_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_addloop
	subq   _b_, $3,  $16
	s8addq _k_, $16  $18
	s8addq _k_, $16  $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)
	subq   $31, _q_, $2
	subq   $2,  _q_, $2
	subq   $2,  _k_, $2
	and    $2,  31,  $3
	bic    $2,  31,  $2
	sll    $3,  3,   $3
	lda    $27, sn_subloop
	subq   _a_, $3,  $16
	subq   _b_, $3,  $18
	subq   _a_, $3,  $20
	s4addq $3,  $27, $27
	bis    $31, $31, $0
	jsr    $27, ($27)
	addq   _q_, _k_, $2
	addq   $2,  _k_, $2
	bsr    $27, sn_decloop
	
	# termin
	ret   $31, ($26),1
	
	#undef _a_
	#undef _b_
	#undef _c_
	#undef _k_
	#undef _q_	

#ifdef debug_sjoin
        .end sn_sjoin3_buggy
#else
        .end sn_sqr
#endif
#undef L
#endif /* assembly_sn_sjoin3 */