ia64.S

		add		r23=r23,r22	};;
{ .mfb;	getf.sig	r26=f94			}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p0=r23,r22
		add		r23=r23,carry1	};;
{ .mfb;	getf.sig	r27=f85			}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p8=r23,carry1};;
{ .mii;	getf.sig	r28=f76
	add		r25=r25,r24
	mov		carry1=0		}
{ .mii;		st8		[r32]=r23,16
	(p7)	add		carry2=1,carry3
	(p8)	add		carry2=0,carry3	};;

{ .mfb;	nop.m	0x0				}
{ .mii;	getf.sig	r29=f67
	cmp.ltu		p6,p0=r25,r24
	add		r26=r26,r25		};;
{ .mfb;	getf.sig	r30=f58			}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,r25
	add		r27=r27,r26		};;
{ .mfb;		getf.sig	r16=f113	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r27,r26
	add		r28=r28,r27		};;
{ .mfb;		getf.sig	r17=f104	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r28,r27
	add		r29=r29,r28		};;
{ .mfb;		getf.sig	r18=f95		}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r29,r28
	add		r30=r30,r29		};;
{ .mii;		getf.sig	r19=f86
		add		r17=r17,r16
		mov		carry3=0	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r30,r29
	add		r30=r30,carry2		};;
{ .mii;		getf.sig	r20=f77
		cmp.ltu		p7,p0=r17,r16
		add		r18=r18,r17	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r30,carry2	};;
{ .mfb;		getf.sig	r21=f68		}
{ .mii;	st8		[r33]=r30,16
(p6)	add		carry1=1,carry1		};;

{ .mfb;	getf.sig	r24=f114		}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p0=r18,r17
		add		r19=r19,r18	};;
{ .mfb;	getf.sig	r25=f105		}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p0=r19,r18
		add		r20=r20,r19	};;
{ .mfb;	getf.sig	r26=f96			}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p0=r20,r19
		add		r21=r21,r20	};;
{ .mfb;	getf.sig	r27=f87			}
{ .mii;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p0=r21,r20
		add		r21=r21,carry1	};;
{ .mib;	getf.sig	r28=f78			
	add		r25=r25,r24		}
{ .mib;	(p7)	add		carry3=1,carry3
		cmp.ltu		p7,p8=r21,carry1};;
{ .mii;		st8		[r32]=r21,16
	(p7)	add		carry2=1,carry3
	(p8)	add		carry2=0,carry3	}

{ .mii;	mov		carry1=0
	cmp.ltu		p6,p0=r25,r24
	add		r26=r26,r25		};;
{ .mfb;		getf.sig	r16=f115	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,r25
	add		r27=r27,r26		};;
{ .mfb;		getf.sig	r17=f106	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r27,r26
	add		r28=r28,r27		};;
{ .mfb;		getf.sig	r18=f97		}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r28,r27
	add		r28=r28,carry2		};;
{ .mib;		getf.sig	r19=f88
		add		r17=r17,r16	}
{ .mib;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r28,carry2	};;
{ .mii;	st8		[r33]=r28,16
(p6)	add		carry1=1,carry1		}

{ .mii;		mov		carry2=0
		cmp.ltu		p7,p0=r17,r16
		add		r18=r18,r17	};;
{ .mfb;	getf.sig	r24=f116		}
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r18,r17
		add		r19=r19,r18	};;
{ .mfb;	getf.sig	r25=f107		}
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r19,r18
		add		r19=r19,carry1	};;
{ .mfb;	getf.sig	r26=f98			}
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r19,carry1};;
{ .mii;		st8		[r32]=r19,16
	(p7)	add		carry2=1,carry2	}

{ .mfb;	add		r25=r25,r24		};;

{ .mfb;		getf.sig	r16=f117	}
{ .mii;	mov		carry1=0
	cmp.ltu		p6,p0=r25,r24
	add		r26=r26,r25		};;
{ .mfb;		getf.sig	r17=f108	}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,r25
	add		r26=r26,carry2		};;
{ .mfb;	nop.m	0x0				}
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,carry2	};;
{ .mii;	st8		[r33]=r26,16
(p6)	add		carry1=1,carry1		}

{ .mfb;		add		r17=r17,r16	};;
{ .mfb;	getf.sig	r24=f118		}
{ .mii;		mov		carry2=0
		cmp.ltu		p7,p0=r17,r16
		add		r17=r17,carry1	};;
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r17,carry1};;
{ .mii;		st8		[r32]=r17
	(p7)	add		carry2=1,carry2	};;
{ .mfb;	add		r24=r24,carry2		};;
{ .mib;	st8		[r33]=r24		}

{ .mib;	rum		1<<5		// clear um.mfh
	br.ret.sptk.many	b0	};;
.endp	bn_mul_comba8#
#undef	carry3
#undef	carry2
#undef	carry1
#endif

#if 1
// It's possible to make it faster (see comment to bn_sqr_comba8), but
// I reckon it doesn't worth the effort. Basically because the routine
// (actually both of them) practically never called... So I just play
// same trick as with bn_sqr_comba8.
//
// void bn_sqr_comba4(BN_ULONG *r, BN_ULONG *a)
//
.global	bn_sqr_comba4#
.proc	bn_sqr_comba4#
.align	64
bn_sqr_comba4:
	.prologue
	.save	ar.pfs,r2
#if defined(_HPUX_SOURCE) && !defined(_LP64)
{ .mii;	alloc   r2=ar.pfs,2,1,0,0
	addp4	r32=0,r32
	addp4	r33=0,r33		};;
{ .mii;
#else
{ .mii;	alloc	r2=ar.pfs,2,1,0,0
#endif
	mov	r34=r33
	add	r14=8,r33		};;
	.body
{ .mii;	add	r17=8,r34
	add	r15=16,r33
	add	r18=16,r34		}
{ .mfb;	add	r16=24,r33
	br	.L_cheat_entry_point4	};;
.endp	bn_sqr_comba4#
#endif

#if 1
// Runs in ~115 cycles and ~4.5 times faster than C. Well, whatever...
//
// void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
//
#define	carry1	r14
#define	carry2	r15
.global	bn_mul_comba4#
.proc	bn_mul_comba4#
.align	64
bn_mul_comba4:
	.prologue
	.save	ar.pfs,r2
#if defined(_HPUX_SOURCE) && !defined(_LP64)
{ .mii;	alloc   r2=ar.pfs,3,0,0,0
	addp4	r33=0,r33
	addp4	r34=0,r34		};;
{ .mii;	addp4	r32=0,r32
#else
{ .mii;	alloc	r2=ar.pfs,3,0,0,0
#endif
	add	r14=8,r33
	add	r17=8,r34		}
	.body
{ .mii;	add	r15=16,r33
	add	r18=16,r34
	add	r16=24,r33		};;
.L_cheat_entry_point4:
{ .mmi;	add	r19=24,r34

	ldf8	f32=[r33]		}

{ .mmi;	ldf8	f120=[r34]
	ldf8	f121=[r17]		};;
{ .mmi;	ldf8	f122=[r18]
	ldf8	f123=[r19]		}

{ .mmi;	ldf8	f33=[r14]
	ldf8	f34=[r15]		}
{ .mfi;	ldf8	f35=[r16]

		xma.hu	f41=f32,f120,f0		}
{ .mfi;		xma.lu	f40=f32,f120,f0		};;
{ .mfi;		xma.hu	f51=f32,f121,f0		}
{ .mfi;		xma.lu	f50=f32,f121,f0		};;
{ .mfi;		xma.hu	f61=f32,f122,f0		}
{ .mfi;		xma.lu	f60=f32,f122,f0		};;
{ .mfi;		xma.hu	f71=f32,f123,f0		}
{ .mfi;		xma.lu	f70=f32,f123,f0		};;//
// Major stall takes place here, and 3 more places below. Result from
// first xma is not available for another 3 ticks.
{ .mfi;	getf.sig	r16=f40
		xma.hu	f42=f33,f120,f41
	add		r33=8,r32		}
{ .mfi;		xma.lu	f41=f33,f120,f41	};;
{ .mfi;	getf.sig	r24=f50
		xma.hu	f52=f33,f121,f51	}
{ .mfi;		xma.lu	f51=f33,f121,f51	};;
{ .mfi;	st8		[r32]=r16,16
		xma.hu	f62=f33,f122,f61	}
{ .mfi;		xma.lu	f61=f33,f122,f61	};;
{ .mfi;		xma.hu	f72=f33,f123,f71	}
{ .mfi;		xma.lu	f71=f33,f123,f71	};;//
//-------------------------------------------------//
{ .mfi;	getf.sig	r25=f41
		xma.hu	f43=f34,f120,f42	}
{ .mfi;		xma.lu	f42=f34,f120,f42	};;
{ .mfi;	getf.sig	r16=f60
		xma.hu	f53=f34,f121,f52	}
{ .mfi;		xma.lu	f52=f34,f121,f52	};;
{ .mfi;	getf.sig	r17=f51
		xma.hu	f63=f34,f122,f62
	add		r25=r25,r24		}
{ .mfi;	mov		carry1=0
		xma.lu	f62=f34,f122,f62	};;
{ .mfi;	st8		[r33]=r25,16
		xma.hu	f73=f34,f123,f72
	cmp.ltu		p6,p0=r25,r24		}
{ .mfi;		xma.lu	f72=f34,f123,f72	};;//
//-------------------------------------------------//
{ .mfi;	getf.sig	r18=f42
		xma.hu	f44=f35,f120,f43
(p6)	add		carry1=1,carry1		}
{ .mfi;	add		r17=r17,r16
		xma.lu	f43=f35,f120,f43
	mov		carry2=0		};;
{ .mfi;	getf.sig	r24=f70
		xma.hu	f54=f35,f121,f53
	cmp.ltu		p7,p0=r17,r16		}
{ .mfi;		xma.lu	f53=f35,f121,f53	};;
{ .mfi;	getf.sig	r25=f61
		xma.hu	f64=f35,f122,f63
	add		r18=r18,r17		}
{ .mfi;		xma.lu	f63=f35,f122,f63
(p7)	add		carry2=1,carry2		};;
{ .mfi;	getf.sig	r26=f52
		xma.hu	f74=f35,f123,f73
	cmp.ltu		p7,p0=r18,r17		}
{ .mfi;		xma.lu	f73=f35,f123,f73
	add		r18=r18,carry1		};;
//-------------------------------------------------//
{ .mii;	st8		[r32]=r18,16
(p7)	add		carry2=1,carry2
	cmp.ltu		p7,p0=r18,carry1	};;

{ .mfi;	getf.sig	r27=f43	// last major stall
(p7)	add		carry2=1,carry2		};;
{ .mii;		getf.sig	r16=f71
	add		r25=r25,r24
	mov		carry1=0		};;
{ .mii;		getf.sig	r17=f62	
	cmp.ltu		p6,p0=r25,r24
	add		r26=r26,r25		};;
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,r25
	add		r27=r27,r26		};;
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r27,r26
	add		r27=r27,carry2		};;
{ .mii;		getf.sig	r18=f53
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r27,carry2	};;
{ .mfi;	st8		[r33]=r27,16
(p6)	add		carry1=1,carry1		}

{ .mii;		getf.sig	r19=f44
		add		r17=r17,r16
		mov		carry2=0	};;
{ .mii;	getf.sig	r24=f72
		cmp.ltu		p7,p0=r17,r16
		add		r18=r18,r17	};;
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r18,r17
		add		r19=r19,r18	};;
{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r19,r18
		add		r19=r19,carry1	};;
{ .mii;	getf.sig	r25=f63
	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r19,carry1};;
{ .mii;		st8		[r32]=r19,16
	(p7)	add		carry2=1,carry2	}

{ .mii;	getf.sig	r26=f54
	add		r25=r25,r24
	mov		carry1=0		};;
{ .mii;		getf.sig	r16=f73
	cmp.ltu		p6,p0=r25,r24
	add		r26=r26,r25		};;
{ .mii;
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,r25
	add		r26=r26,carry2		};;
{ .mii;		getf.sig	r17=f64
(p6)	add		carry1=1,carry1
	cmp.ltu		p6,p0=r26,carry2	};;
{ .mii;	st8		[r33]=r26,16
(p6)	add		carry1=1,carry1		}

{ .mii;	getf.sig	r24=f74
		add		r17=r17,r16	
		mov		carry2=0	};;
{ .mii;		cmp.ltu		p7,p0=r17,r16
		add		r17=r17,carry1	};;

{ .mii;	(p7)	add		carry2=1,carry2
		cmp.ltu		p7,p0=r17,carry1};;
{ .mii;		st8		[r32]=r17,16
	(p7)	add		carry2=1,carry2	};;

{ .mii;	add		r24=r24,carry2		};;
{ .mii;	st8		[r33]=r24		}

{ .mib;	rum		1<<5		// clear um.mfh
	br.ret.sptk.many	b0	};;
.endp	bn_mul_comba4#
#undef	carry2
#undef	carry1
#endif

#if 1
//
// BN_ULONG bn_div_words(BN_ULONG h, BN_ULONG l, BN_ULONG d)
//
// In the nutshell it's a port of my MIPS III/IV implementation.
//
#define	AT	r14
#define	H	r16
#define	HH	r20
#define	L	r17
#define	D	r18
#define	DH	r22
#define	I	r21

#if 0
// Some preprocessors (most notably HP-UX) appear to be allergic to
// macros enclosed to parenthesis [as these three were].
#define	cont	p16
#define	break	p0	// p20
#define	equ	p24
#else
cont=p16
break=p0
equ=p24
#endif

.global	abort#
.global	bn_div_words#
.proc	bn_div_words#
.align	64
bn_div_words:
	.prologue
	.save	ar.pfs,r2
{ .mii;	alloc		r2=ar.pfs,3,5,0,8
	.save	b0,r3
	mov		r3=b0
	.save	pr,r10
	mov		r10=pr		};;
{ .mmb;	cmp.eq		p6,p0=r34,r0
	mov		r8=-1
(p6)	br.ret.spnt.many	b0	};;

	.body
{ .mii;	mov		H=r32		// save h
	mov		ar.ec=0		// don't rotate at exit
	mov		pr.rot=0	}
{ .mii;	mov		L=r33		// save l
	mov		r36=r0		};;

.L_divw_shift:	// -vv- note signed comparison
{ .mfi;	(p0)	cmp.lt		p16,p0=r0,r34	// d
	(p0)	shladd		r33=r34,1,r0	}
{ .mfb;	(p0)	add		r35=1,r36
	(p0)	nop.f		0x0
(p16)	br.wtop.dpnt		.L_divw_shift	};;

{ .mii;	mov		D=r34
	shr.u		DH=r34,32
	sub		r35=64,r36		};;
{ .mii;	setf.sig	f7=DH
	shr.u		AT=H,r35
	mov		I=r36			};;
{ .mib;	cmp.ne		p6,p0=r0,AT
	shl		H=H,r36
(p6)	br.call.spnt.clr	b0=abort	};;	// overflow, die...

{ .mfi;	fcvt.xuf.s1	f7=f7
	shr.u		AT=L,r35		};;
{ .mii;	shl		L=L,r36
	or		H=H,AT			};;

{ .mii;	nop.m		0x0
	cmp.leu		p6,p0=D,H;;
(p6)	sub		H=H,D			}

{ .mlx;	setf.sig	f14=D
	movl		AT=0xffffffff		};;
///////////////////////////////////////////////////////////
{ .mii;	setf.sig	f6=H
	shr.u		HH=H,32;;
	cmp.eq		p6,p7=HH,DH		};;
{ .mfb;
(p6)	setf.sig	f8=AT
(p7)	fcvt.xuf.s1	f6=f6
(p7)	br.call.sptk	b6=.L_udiv64_32_b6	};;

{ .mfi;	getf.sig	r33=f8				// q
	xmpy.lu		f9=f8,f14		}
{ .mfi;	xmpy.hu		f10=f8,f14
	shrp		H=H,L,32		};;

{ .mmi;	getf.sig	r35=f9				// tl
	getf.sig	r31=f10			};;	// th

.L_divw_1st_iter:
{ .mii;	(p0)	add		r32=-1,r33
	(p0)	cmp.eq		equ,cont=HH,r31		};;
{ .mii;	(p0)	cmp.ltu		p8,p0=r35,D
	(p0)	sub		r34=r35,D
	(equ)	cmp.leu		break,cont=r35,H	};;
{ .mib;	(cont)	cmp.leu		cont,break=HH,r31
	(p8)	add		r31=-1,r31
(cont)	br.wtop.spnt		.L_divw_1st_iter	};;
///////////////////////////////////////////////////////////
{ .mii;	sub		H=H,r35
	shl		r8=r33,32
	shl		L=L,32			};;
///////////////////////////////////////////////////////////
{ .mii;	setf.sig	f6=H
	shr.u		HH=H,32;;
	cmp.eq		p6,p7=HH,DH		};;
{ .mfb;
(p6)	setf.sig	f8=AT
(p7)	fcvt.xuf.s1	f6=f6
(p7)	br.call.sptk	b6=.L_udiv64_32_b6	};;

{ .mfi;	getf.sig	r33=f8				// q
	xmpy.lu		f9=f8,f14		}
{ .mfi;	xmpy.hu		f10=f8,f14
	shrp		H=H,L,32		};;

{ .mmi;	getf.sig	r35=f9				// tl
	getf.sig	r31=f10			};;	// th

.L_divw_2nd_iter:
{ .mii;	(p0)	add		r32=-1,r33
	(p0)	cmp.eq		equ,cont=HH,r31		};;
{ .mii;	(p0)	cmp.ltu		p8,p0=r35,D
	(p0)	sub		r34=r35,D
	(equ)	cmp.leu		break,cont=r35,H	};;
{ .mib;	(cont)	cmp.leu		cont,break=HH,r31
	(p8)	add		r31=-1,r31
(cont)	br.wtop.spnt		.L_divw_2nd_iter	};;
///////////////////////////////////////////////////////////
{ .mii;	sub	H=H,r35
	or	r8=r8,r33
	mov	ar.pfs=r2		};;
{ .mii;	shr.u	r9=H,I			// remainder if anybody wants it
	mov	pr=r10,0x1ffff		}
{ .mfb;	br.ret.sptk.many	b0	};;

// Unsigned 64 by 32 (well, by 64 for the moment) bit integer division
// procedure.
//
// inputs:	f6 = (double)a, f7 = (double)b
// output:	f8 = (int)(a/b)
// clobbered:	f8,f9,f10,f11,pred
pred=p15
// One can argue that this snippet is copyrighted to Intel
// Corporation, as it's essentially identical to one of those
// found in "Divide, Square Root and Remainder" section at
// http://www.intel.com/software/products/opensource/libraries/num.htm.
// Yes, I admit that the referred code was used as template,
// but after I realized that there hardly is any other instruction
// sequence which would perform this operation. I mean I figure that
// any independent attempt to implement high-performance division
// will result in code virtually identical to the Intel code. It
// should be noted though that below division kernel is 1 cycle
// faster than Intel one (note commented splits:-), not to mention
// original prologue (rather lack of one) and epilogue.
.align	32
.skip	16
.L_udiv64_32_b6:
	frcpa.s1	f8,pred=f6,f7;;		// [0]  y0 = 1 / b

(pred)	fnma.s1		f9=f7,f8,f1		// [5]  e0 = 1 - b * y0
(pred)	fmpy.s1		f10=f6,f8;;		// [5]  q0 = a * y0
(pred)	fmpy.s1		f11=f9,f9		// [10] e1 = e0 * e0
(pred)	fma.s1		f10=f9,f10,f10;;	// [10] q1 = q0 + e0 * q0
(pred)	fma.s1		f8=f9,f8,f8	//;;	// [15] y1 = y0 + e0 * y0
(pred)	fma.s1		f9=f11,f10,f10;;	// [15] q2 = q1 + e1 * q1
(pred)	fma.s1		f8=f11,f8,f8	//;;	// [20] y2 = y1 + e1 * y1
(pred)	fnma.s1		f10=f7,f9,f6;;		// [20] r2 = a - b * q2
(pred)	fma.s1		f8=f10,f8,f9;;		// [25] q3 = q2 + r2 * y2

	fcvt.fxu.trunc.s1	f8=f8		// [30] q = trunc(q3)
	br.ret.sptk.many	b6;;
.endp	bn_div_words#
#endif