xref: /openssl/crypto/bn/asm/ia64.S (revision df443918)
1.explicit
2.text
3.ident	"ia64.S, Version 2.1"
4.ident	"IA-64 ISA artwork by Andy Polyakov <appro@openssl.org>"
5
6// Copyright 2001-2018 The OpenSSL Project Authors. All Rights Reserved.
7//
8// Licensed under the Apache License 2.0 (the "License").  You may not use
9// this file except in compliance with the License.  You can obtain a copy
10// in the file LICENSE in the source distribution or at
11// https://www.openssl.org/source/license.html
12
13//
14// ====================================================================
15// Written by Andy Polyakov <appro@openssl.org> for the OpenSSL
16// project.
17//
18// Rights for redistribution and usage in source and binary forms are
19// granted according to the License. Warranty of any kind is disclaimed.
20// ====================================================================
21//
22// Version 2.x is Itanium2 re-tune. Few words about how Itanium2 is
23// different from Itanium to this module viewpoint. Most notably, is it
24// "wider" than Itanium? Can you experience loop scalability as
25// discussed in commentary sections? Not really:-( Itanium2 has 6
26// integer ALU ports, i.e. it's 2 ports wider, but it's not enough to
27// spin twice as fast, as I need 8 IALU ports. Amount of floating point
28// ports is the same, i.e. 2, while I need 4. In other words, to this
29// module Itanium2 remains effectively as "wide" as Itanium. Yet it's
30// essentially different in respect to this module, and a re-tune was
31// required. Well, because some instruction latencies has changed. Most
32// noticeably those intensively used:
33//
34//			Itanium	Itanium2
35//	ldf8		9	6		L2 hit
36//	ld8		2	1		L1 hit
37//	getf		2	5
38//	xma[->getf]	7[+1]	4[+0]
39//	add[->st8]	1[+1]	1[+0]
40//
41// What does it mean? You might ratiocinate that the original code
42// should run just faster... Because sum of latencies is smaller...
43// Wrong! Note that getf latency increased. This means that if a loop is
44// scheduled for lower latency (as they were), then it will suffer from
45// stall condition and the code will therefore turn anti-scalable, e.g.
46// original bn_mul_words spun at 5*n or 2.5 times slower than expected
47// on Itanium2! What to do? Reschedule loops for Itanium2? But then
48// Itanium would exhibit anti-scalability. So I've chosen to reschedule
49// for worst latency for every instruction aiming for best *all-round*
50// performance.
51
52// Q.	How much faster does it get?
53// A.	Here is the output from 'openssl speed rsa dsa' for vanilla
54//	0.9.6a compiled with gcc version 2.96 20000731 (Red Hat
55//	Linux 7.1 2.96-81):
56//
57//	                  sign    verify    sign/s verify/s
58//	rsa  512 bits   0.0036s   0.0003s    275.3   2999.2
59//	rsa 1024 bits   0.0203s   0.0011s     49.3    894.1
60//	rsa 2048 bits   0.1331s   0.0040s      7.5    250.9
61//	rsa 4096 bits   0.9270s   0.0147s      1.1     68.1
62//	                  sign    verify    sign/s verify/s
63//	dsa  512 bits   0.0035s   0.0043s    288.3    234.8
64//	dsa 1024 bits   0.0111s   0.0135s     90.0     74.2
65//
66//	And here is similar output but for this assembler
67//	implementation:-)
68//
69//	                  sign    verify    sign/s verify/s
70//	rsa  512 bits   0.0021s   0.0001s    549.4   9638.5
71//	rsa 1024 bits   0.0055s   0.0002s    183.8   4481.1
72//	rsa 2048 bits   0.0244s   0.0006s     41.4   1726.3
73//	rsa 4096 bits   0.1295s   0.0018s      7.7    561.5
74//	                  sign    verify    sign/s verify/s
75//	dsa  512 bits   0.0012s   0.0013s    891.9    756.6
76//	dsa 1024 bits   0.0023s   0.0028s    440.4    376.2
77//
78//	Yes, you may argue that it's not fair comparison as it's
79//	possible to craft the C implementation with BN_UMULT_HIGH
80//	inline assembler macro. But of course! Here is the output
81//	with the macro:
82//
83//	                  sign    verify    sign/s verify/s
84//	rsa  512 bits   0.0020s   0.0002s    495.0   6561.0
85//	rsa 1024 bits   0.0086s   0.0004s    116.2   2235.7
86//	rsa 2048 bits   0.0519s   0.0015s     19.3    667.3
87//	rsa 4096 bits   0.3464s   0.0053s      2.9    187.7
88//	                  sign    verify    sign/s verify/s
89//	dsa  512 bits   0.0016s   0.0020s    613.1    510.5
90//	dsa 1024 bits   0.0045s   0.0054s    221.0    183.9
91//
92//	My code is still way faster, huh:-) And I believe that even
93//	higher performance can be achieved. Note that as keys get
94//	longer, performance gain is larger. Why? According to the
95//	profiler there is another player in the field, namely
96//	BN_from_montgomery consuming larger and larger portion of CPU
97//	time as keysize decreases. I therefore consider putting effort
98//	to assembler implementation of the following routine:
99//
100//	void bn_mul_add_mont (BN_ULONG *rp,BN_ULONG *np,int nl,BN_ULONG n0)
101//	{
102//	int      i,j;
103//	BN_ULONG v;
104//
105//	for (i=0; i<nl; i++)
106//		{
107//		v=bn_mul_add_words(rp,np,nl,(rp[0]*n0)&BN_MASK2);
108//		nrp++;
109//		rp++;
110//		if (((nrp[-1]+=v)&BN_MASK2) < v)
111//			for (j=0; ((++nrp[j])&BN_MASK2) == 0; j++) ;
112//		}
113//	}
114//
115//	It might as well be beneficial to implement even combaX
116//	variants, as it appears as it can literally unleash the
117//	performance (see comment section to bn_mul_comba8 below).
118//
119//	And finally for your reference the output for 0.9.6a compiled
120//	with SGIcc version 0.01.0-12 (keep in mind that for the moment
121//	of this writing it's not possible to convince SGIcc to use
122//	BN_UMULT_HIGH inline assembler macro, yet the code is fast,
123//	i.e. for a compiler generated one:-):
124//
125//	                  sign    verify    sign/s verify/s
126//	rsa  512 bits   0.0022s   0.0002s    452.7   5894.3
127//	rsa 1024 bits   0.0097s   0.0005s    102.7   2002.9
128//	rsa 2048 bits   0.0578s   0.0017s     17.3    600.2
129//	rsa 4096 bits   0.3838s   0.0061s      2.6    164.5
130//	                  sign    verify    sign/s verify/s
131//	dsa  512 bits   0.0018s   0.0022s    547.3    459.6
132//	dsa 1024 bits   0.0051s   0.0062s    196.6    161.3
133//
134//	Oh! Benchmarks were performed on 733MHz Lion-class Itanium
135//	system running Redhat Linux 7.1 (very special thanks to Ray
136//	McCaffity of Williams Communications for providing an account).
137//
138// Q.	What's the heck with 'rum 1<<5' at the end of every function?
139// A.	Well, by clearing the "upper FP registers written" bit of the
140//	User Mask I want to excuse the kernel from preserving upper
141//	(f32-f128) FP register bank over process context switch, thus
142//	minimizing bus bandwidth consumption during the switch (i.e.
143//	after PKI operation completes and the program is off doing
144//	something else like bulk symmetric encryption). Having said
145//	this, I also want to point out that it might be good idea
146//	to compile the whole toolkit (as well as majority of the
147//	programs for that matter) with -mfixed-range=f32-f127 command
148//	line option. No, it doesn't prevent the compiler from writing
149//	to upper bank, but at least discourages to do so. If you don't
150//	like the idea you have the option to compile the module with
151//	-Drum=nop.m in command line.
152//
153
154#if defined(_HPUX_SOURCE) && !defined(_LP64)
155#define	ADDP	addp4
156#else
157#define	ADDP	add
158#endif
159#ifdef __VMS
160.alias abort, "decc$abort"
161#endif
162
163#if 1
164//
165// bn_[add|sub]_words routines.
166//
167// Loops are spinning in 2*(n+5) ticks on Itanium (provided that the
168// data reside in L1 cache, i.e. 2 ticks away). It's possible to
169// compress the epilogue and get down to 2*n+6, but at the cost of
170// scalability (the neat feature of this implementation is that it
171// shall automagically spin in n+5 on "wider" IA-64 implementations:-)
172// I consider that the epilogue is short enough as it is to trade tiny
173// performance loss on Itanium for scalability.
174//
175// BN_ULONG bn_add_words(BN_ULONG *rp, BN_ULONG *ap, BN_ULONG *bp,int num)
176//
177.global	bn_add_words#
178.proc	bn_add_words#
179.align	64
180.skip	32	// makes the loop body aligned at 64-byte boundary
181bn_add_words:
182	.prologue
183	.save	ar.pfs,r2
184{ .mii;	alloc		r2=ar.pfs,4,12,0,16
185	cmp4.le		p6,p0=r35,r0	};;
186{ .mfb;	mov		r8=r0			// return value
187(p6)	br.ret.spnt.many	b0	};;
188
189{ .mib;	sub		r10=r35,r0,1
190	.save	ar.lc,r3
191	mov		r3=ar.lc
192	brp.loop.imp	.L_bn_add_words_ctop,.L_bn_add_words_cend-16
193					}
194{ .mib;	ADDP		r14=0,r32		// rp
195	.save	pr,r9
196	mov		r9=pr		};;
197	.body
198{ .mii;	ADDP		r15=0,r33		// ap
199	mov		ar.lc=r10
200	mov		ar.ec=6		}
201{ .mib;	ADDP		r16=0,r34		// bp
202	mov		pr.rot=1<<16	};;
203
204.L_bn_add_words_ctop:
205{ .mii;	(p16)	ld8		r32=[r16],8	  // b=*(bp++)
206	(p18)	add		r39=r37,r34
207	(p19)	cmp.ltu.unc	p56,p0=r40,r38	}
208{ .mfb;	(p0)	nop.m		0x0
209	(p0)	nop.f		0x0
210	(p0)	nop.b		0x0		}
211{ .mii;	(p16)	ld8		r35=[r15],8	  // a=*(ap++)
212	(p58)	cmp.eq.or	p57,p0=-1,r41	  // (p20)
213	(p58)	add		r41=1,r41	} // (p20)
214{ .mfb;	(p21)	st8		[r14]=r42,8	  // *(rp++)=r
215	(p0)	nop.f		0x0
216	br.ctop.sptk	.L_bn_add_words_ctop	};;
217.L_bn_add_words_cend:
218
219{ .mii;
220(p59)	add		r8=1,r8		// return value
221	mov		pr=r9,0x1ffff
222	mov		ar.lc=r3	}
223{ .mbb;	nop.b		0x0
224	br.ret.sptk.many	b0	};;
225.endp	bn_add_words#
226
227//
228// BN_ULONG bn_sub_words(BN_ULONG *rp, BN_ULONG *ap, BN_ULONG *bp,int num)
229//
230.global	bn_sub_words#
231.proc	bn_sub_words#
232.align	64
233.skip	32	// makes the loop body aligned at 64-byte boundary
234bn_sub_words:
235	.prologue
236	.save	ar.pfs,r2
237{ .mii;	alloc		r2=ar.pfs,4,12,0,16
238	cmp4.le		p6,p0=r35,r0	};;
239{ .mfb;	mov		r8=r0			// return value
240(p6)	br.ret.spnt.many	b0	};;
241
242{ .mib;	sub		r10=r35,r0,1
243	.save	ar.lc,r3
244	mov		r3=ar.lc
245	brp.loop.imp	.L_bn_sub_words_ctop,.L_bn_sub_words_cend-16
246					}
247{ .mib;	ADDP		r14=0,r32		// rp
248	.save	pr,r9
249	mov		r9=pr		};;
250	.body
251{ .mii;	ADDP		r15=0,r33		// ap
252	mov		ar.lc=r10
253	mov		ar.ec=6		}
254{ .mib;	ADDP		r16=0,r34		// bp
255	mov		pr.rot=1<<16	};;
256
257.L_bn_sub_words_ctop:
258{ .mii;	(p16)	ld8		r32=[r16],8	  // b=*(bp++)
259	(p18)	sub		r39=r37,r34
260	(p19)	cmp.gtu.unc	p56,p0=r40,r38	}
261{ .mfb;	(p0)	nop.m		0x0
262	(p0)	nop.f		0x0
263	(p0)	nop.b		0x0		}
264{ .mii;	(p16)	ld8		r35=[r15],8	  // a=*(ap++)
265	(p58)	cmp.eq.or	p57,p0=0,r41	  // (p20)
266	(p58)	add		r41=-1,r41	} // (p20)
267{ .mbb;	(p21)	st8		[r14]=r42,8	  // *(rp++)=r
268	(p0)	nop.b		0x0
269	br.ctop.sptk	.L_bn_sub_words_ctop	};;
270.L_bn_sub_words_cend:
271
272{ .mii;
273(p59)	add		r8=1,r8		// return value
274	mov		pr=r9,0x1ffff
275	mov		ar.lc=r3	}
276{ .mbb;	nop.b		0x0
277	br.ret.sptk.many	b0	};;
278.endp	bn_sub_words#
279#endif
280
281#if 0
282#define XMA_TEMPTATION
283#endif
284
285#if 1
286//
287// BN_ULONG bn_mul_words(BN_ULONG *rp, BN_ULONG *ap, int num, BN_ULONG w)
288//
289.global	bn_mul_words#
290.proc	bn_mul_words#
291.align	64
292.skip	32	// makes the loop body aligned at 64-byte boundary
293bn_mul_words:
294	.prologue
295	.save	ar.pfs,r2
296#ifdef XMA_TEMPTATION
297{ .mfi;	alloc		r2=ar.pfs,4,0,0,0	};;
298#else
299{ .mfi;	alloc		r2=ar.pfs,4,12,0,16	};;
300#endif
301{ .mib;	mov		r8=r0			// return value
302	cmp4.le		p6,p0=r34,r0
303(p6)	br.ret.spnt.many	b0		};;
304
305{ .mii;	sub	r10=r34,r0,1
306	.save	ar.lc,r3
307	mov	r3=ar.lc
308	.save	pr,r9
309	mov	r9=pr			};;
310
311	.body
312{ .mib;	setf.sig	f8=r35	// w
313	mov		pr.rot=0x800001<<16
314			// ------^----- serves as (p50) at first (p27)
315	brp.loop.imp	.L_bn_mul_words_ctop,.L_bn_mul_words_cend-16
316					}
317
318#ifndef XMA_TEMPTATION
319
320{ .mmi;	ADDP		r14=0,r32	// rp
321	ADDP		r15=0,r33	// ap
322	mov		ar.lc=r10	}
323{ .mmi;	mov		r40=0		// serves as r35 at first (p27)
324	mov		ar.ec=13	};;
325
326// This loop spins in 2*(n+12) ticks. It's scheduled for data in Itanium
327// L2 cache (i.e. 9 ticks away) as floating point load/store instructions
328// bypass L1 cache and L2 latency is actually best-case scenario for
329// ldf8. The loop is not scalable and shall run in 2*(n+12) even on
330// "wider" IA-64 implementations. It's a trade-off here. n+24 loop
331// would give us ~5% in *overall* performance improvement on "wider"
332// IA-64, but would hurt Itanium for about same because of longer
333// epilogue. As it's a matter of few percents in either case I've
334// chosen to trade the scalability for development time (you can see
335// this very instruction sequence in bn_mul_add_words loop which in
336// turn is scalable).
337.L_bn_mul_words_ctop:
338{ .mfi;	(p25)	getf.sig	r36=f52			// low
339	(p21)	xmpy.lu		f48=f37,f8
340	(p28)	cmp.ltu		p54,p50=r41,r39	}
341{ .mfi;	(p16)	ldf8		f32=[r15],8
342	(p21)	xmpy.hu		f40=f37,f8
343	(p0)	nop.i		0x0		};;
344{ .mii;	(p25)	getf.sig	r32=f44			// high
345	.pred.rel	"mutex",p50,p54
346	(p50)	add		r40=r38,r35		// (p27)
347	(p54)	add		r40=r38,r35,1	}	// (p27)
348{ .mfb;	(p28)	st8		[r14]=r41,8
349	(p0)	nop.f		0x0
350	br.ctop.sptk	.L_bn_mul_words_ctop	};;
351.L_bn_mul_words_cend:
352
353{ .mii;	nop.m		0x0
354.pred.rel	"mutex",p51,p55
355(p51)	add		r8=r36,r0
356(p55)	add		r8=r36,r0,1	}
357{ .mfb;	nop.m	0x0
358	nop.f	0x0
359	nop.b	0x0			}
360
361#else	// XMA_TEMPTATION
362
363	setf.sig	f37=r0	// serves as carry at (p18) tick
364	mov		ar.lc=r10
365	mov		ar.ec=5;;
366
367// Most of you examining this code very likely wonder why in the name
368// of Intel the following loop is commented out? Indeed, it looks so
369// neat that you find it hard to believe that it's something wrong
370// with it, right? The catch is that every iteration depends on the
371// result from previous one and the latter isn't available instantly.
372// The loop therefore spins at the latency of xma minus 1, or in other
373// words at 6*(n+4) ticks:-( Compare to the "production" loop above
374// that runs in 2*(n+11) where the low latency problem is worked around
375// by moving the dependency to one-tick latent integer ALU. Note that
376// "distance" between ldf8 and xma is not latency of ldf8, but the
377// *difference* between xma and ldf8 latencies.
378.L_bn_mul_words_ctop:
379{ .mfi;	(p16)	ldf8		f32=[r33],8
380	(p18)	xma.hu		f38=f34,f8,f39	}
381{ .mfb;	(p20)	stf8		[r32]=f37,8
382	(p18)	xma.lu		f35=f34,f8,f39
383	br.ctop.sptk	.L_bn_mul_words_ctop	};;
384.L_bn_mul_words_cend:
385
386	getf.sig	r8=f41		// the return value
387
388#endif	// XMA_TEMPTATION
389
390{ .mii;	nop.m		0x0
391	mov		pr=r9,0x1ffff
392	mov		ar.lc=r3	}
393{ .mfb;	rum		1<<5		// clear um.mfh
394	nop.f		0x0
395	br.ret.sptk.many	b0	};;
396.endp	bn_mul_words#
397#endif
398
399#if 1
400//
401// BN_ULONG bn_mul_add_words(BN_ULONG *rp, BN_ULONG *ap, int num, BN_ULONG w)
402//
403.global	bn_mul_add_words#
404.proc	bn_mul_add_words#
405.align	64
406.skip	48	// makes the loop body aligned at 64-byte boundary
407bn_mul_add_words:
408	.prologue
409	.save	ar.pfs,r2
410{ .mmi;	alloc		r2=ar.pfs,4,4,0,8
411	cmp4.le		p6,p0=r34,r0
412	.save	ar.lc,r3
413	mov		r3=ar.lc	};;
414{ .mib;	mov		r8=r0		// return value
415	sub		r10=r34,r0,1
416(p6)	br.ret.spnt.many	b0	};;
417
418{ .mib;	setf.sig	f8=r35		// w
419	.save	pr,r9
420	mov		r9=pr
421	brp.loop.imp	.L_bn_mul_add_words_ctop,.L_bn_mul_add_words_cend-16
422					}
423	.body
424{ .mmi;	ADDP		r14=0,r32	// rp
425	ADDP		r15=0,r33	// ap
426	mov		ar.lc=r10	}
427{ .mii;	ADDP		r16=0,r32	// rp copy
428	mov		pr.rot=0x2001<<16
429			// ------^----- serves as (p40) at first (p27)
430	mov		ar.ec=11	};;
431
432// This loop spins in 3*(n+10) ticks on Itanium and in 2*(n+10) on
433// Itanium 2. Yes, unlike previous versions it scales:-) Previous
434// version was performing *all* additions in IALU and was starving
435// for those even on Itanium 2. In this version one addition is
436// moved to FPU and is folded with multiplication. This is at cost
437// of propagating the result from previous call to this subroutine
438// to L2 cache... In other words negligible even for shorter keys.
439// *Overall* performance improvement [over previous version] varies
440// from 11 to 22 percent depending on key length.
441.L_bn_mul_add_words_ctop:
442.pred.rel	"mutex",p40,p42
443{ .mfi;	(p23)	getf.sig	r36=f45			// low
444	(p20)	xma.lu		f42=f36,f8,f50		// low
445	(p40)	add		r39=r39,r35	}	// (p27)
446{ .mfi;	(p16)	ldf8		f32=[r15],8		// *(ap++)
447	(p20)	xma.hu		f36=f36,f8,f50		// high
448	(p42)	add		r39=r39,r35,1	};;	// (p27)
449{ .mmi;	(p24)	getf.sig	r32=f40			// high
450	(p16)	ldf8		f46=[r16],8		// *(rp1++)
451	(p40)	cmp.ltu		p41,p39=r39,r35	}	// (p27)
452{ .mib;	(p26)	st8		[r14]=r39,8		// *(rp2++)
453	(p42)	cmp.leu		p41,p39=r39,r35		// (p27)
454	br.ctop.sptk	.L_bn_mul_add_words_ctop};;
455.L_bn_mul_add_words_cend:
456
457{ .mmi;	.pred.rel	"mutex",p40,p42
458(p40)	add		r8=r35,r0
459(p42)	add		r8=r35,r0,1
460	mov		pr=r9,0x1ffff	}
461{ .mib;	rum		1<<5		// clear um.mfh
462	mov		ar.lc=r3
463	br.ret.sptk.many	b0	};;
464.endp	bn_mul_add_words#
465#endif
466
467#if 1
468//
469// void bn_sqr_words(BN_ULONG *rp, BN_ULONG *ap, int num)
470//
471.global	bn_sqr_words#
472.proc	bn_sqr_words#
473.align	64
474.skip	32	// makes the loop body aligned at 64-byte boundary
475bn_sqr_words:
476	.prologue
477	.save	ar.pfs,r2
478{ .mii;	alloc		r2=ar.pfs,3,0,0,0
479	sxt4		r34=r34		};;
480{ .mii;	cmp.le		p6,p0=r34,r0
481	mov		r8=r0		}	// return value
482{ .mfb;	ADDP		r32=0,r32
483	nop.f		0x0
484(p6)	br.ret.spnt.many	b0	};;
485
486{ .mii;	sub	r10=r34,r0,1
487	.save	ar.lc,r3
488	mov	r3=ar.lc
489	.save	pr,r9
490	mov	r9=pr			};;
491
492	.body
493{ .mib;	ADDP		r33=0,r33
494	mov		pr.rot=1<<16
495	brp.loop.imp	.L_bn_sqr_words_ctop,.L_bn_sqr_words_cend-16
496					}
497{ .mii;	add		r34=8,r32
498	mov		ar.lc=r10
499	mov		ar.ec=18	};;
500
501// 2*(n+17) on Itanium, (n+17) on "wider" IA-64 implementations. It's
502// possible to compress the epilogue (I'm getting tired to write this
503// comment over and over) and get down to 2*n+16 at the cost of
504// scalability. The decision will very likely be reconsidered after the
505// benchmark program is profiled. I.e. if performance gain on Itanium
506// will appear larger than loss on "wider" IA-64, then the loop should
507// be explicitly split and the epilogue compressed.
508.L_bn_sqr_words_ctop:
509{ .mfi;	(p16)	ldf8		f32=[r33],8
510	(p25)	xmpy.lu		f42=f41,f41
511	(p0)	nop.i		0x0		}
512{ .mib;	(p33)	stf8		[r32]=f50,16
513	(p0)	nop.i		0x0
514	(p0)	nop.b		0x0		}
515{ .mfi;	(p0)	nop.m		0x0
516	(p25)	xmpy.hu		f52=f41,f41
517	(p0)	nop.i		0x0		}
518{ .mib;	(p33)	stf8		[r34]=f60,16
519	(p0)	nop.i		0x0
520	br.ctop.sptk	.L_bn_sqr_words_ctop	};;
521.L_bn_sqr_words_cend:
522
523{ .mii;	nop.m		0x0
524	mov		pr=r9,0x1ffff
525	mov		ar.lc=r3	}
526{ .mfb;	rum		1<<5		// clear um.mfh
527	nop.f		0x0
528	br.ret.sptk.many	b0	};;
529.endp	bn_sqr_words#
530#endif
531
532#if 1
533// Apparently we win nothing by implementing special bn_sqr_comba8.
534// Yes, it is possible to reduce the number of multiplications by
535// almost factor of two, but then the amount of additions would
536// increase by factor of two (as we would have to perform those
537// otherwise performed by xma ourselves). Normally we would trade
538// anyway as multiplications are way more expensive, but not this
539// time... Multiplication kernel is fully pipelined and as we drain
540// one 128-bit multiplication result per clock cycle multiplications
541// are effectively as inexpensive as additions. Special implementation
542// might become of interest for "wider" IA-64 implementation as you'll
543// be able to get through the multiplication phase faster (there won't
544// be any stall issues as discussed in the commentary section below and
545// you therefore will be able to employ all 4 FP units)... But these
546// Itanium days it's simply too hard to justify the effort so I just
547// drop down to bn_mul_comba8 code:-)
548//
549// void bn_sqr_comba8(BN_ULONG *r, BN_ULONG *a)
550//
551.global	bn_sqr_comba8#
552.proc	bn_sqr_comba8#
553.align	64
554bn_sqr_comba8:
555	.prologue
556	.save	ar.pfs,r2
557#if defined(_HPUX_SOURCE) && !defined(_LP64)
558{ .mii;	alloc	r2=ar.pfs,2,1,0,0
559	addp4	r33=0,r33
560	addp4	r32=0,r32		};;
561{ .mii;
562#else
563{ .mii;	alloc	r2=ar.pfs,2,1,0,0
564#endif
565	mov	r34=r33
566	add	r14=8,r33		};;
567	.body
568{ .mii;	add	r17=8,r34
569	add	r15=16,r33
570	add	r18=16,r34		}
571{ .mfb;	add	r16=24,r33
572	br	.L_cheat_entry_point8	};;
573.endp	bn_sqr_comba8#
574#endif
575
576#if 1
577// I've estimated this routine to run in ~120 ticks, but in reality
578// (i.e. according to ar.itc) it takes ~160 ticks. Are those extra
579// cycles consumed for instructions fetch? Or did I misinterpret some
580// clause in Itanium µ-architecture manual? Comments are welcomed and
581// highly appreciated.
582//
583// On Itanium 2 it takes ~190 ticks. This is because of stalls on
584// result from getf.sig. I do nothing about it at this point for
585// reasons depicted below.
586//
587// However! It should be noted that even 160 ticks is darn good result
588// as it's over 10 (yes, ten, spelled as t-e-n) times faster than the
589// C version (compiled with gcc with inline assembler). I really
590// kicked compiler's butt here, didn't I? Yeah! This brings us to the
591// following statement. It's damn shame that this routine isn't called
592// very often nowadays! According to the profiler most CPU time is
593// consumed by bn_mul_add_words called from BN_from_montgomery. In
594// order to estimate what we're missing, I've compared the performance
595// of this routine against "traditional" implementation, i.e. against
596// following routine:
597//
598// void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
599// {	r[ 8]=bn_mul_words(    &(r[0]),a,8,b[0]);
600//	r[ 9]=bn_mul_add_words(&(r[1]),a,8,b[1]);
601//	r[10]=bn_mul_add_words(&(r[2]),a,8,b[2]);
602//	r[11]=bn_mul_add_words(&(r[3]),a,8,b[3]);
603//	r[12]=bn_mul_add_words(&(r[4]),a,8,b[4]);
604//	r[13]=bn_mul_add_words(&(r[5]),a,8,b[5]);
605//	r[14]=bn_mul_add_words(&(r[6]),a,8,b[6]);
606//	r[15]=bn_mul_add_words(&(r[7]),a,8,b[7]);
607// }
608//
609// The one below is over 8 times faster than the one above:-( Even
610// more reasons to "combafy" bn_mul_add_mont...
611//
612// And yes, this routine really made me wish there were an optimizing
613// assembler! It also feels like it deserves a dedication.
614//
615//	To my wife for being there and to my kids...
616//
617// void bn_mul_comba8(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
618//
619#define	carry1	r14
620#define	carry2	r15
621#define	carry3	r34
622.global	bn_mul_comba8#
623.proc	bn_mul_comba8#
624.align	64
625bn_mul_comba8:
626	.prologue
627	.save	ar.pfs,r2
628#if defined(_HPUX_SOURCE) && !defined(_LP64)
629{ .mii;	alloc	r2=ar.pfs,3,0,0,0
630	addp4	r33=0,r33
631	addp4	r34=0,r34		};;
632{ .mii;	addp4	r32=0,r32
633#else
634{ .mii;	alloc   r2=ar.pfs,3,0,0,0
635#endif
636	add	r14=8,r33
637	add	r17=8,r34		}
638	.body
639{ .mii;	add	r15=16,r33
640	add	r18=16,r34
641	add	r16=24,r33		}
642.L_cheat_entry_point8:
643{ .mmi;	add	r19=24,r34
644
645	ldf8	f32=[r33],32		};;
646
647{ .mmi;	ldf8	f120=[r34],32
648	ldf8	f121=[r17],32		}
649{ .mmi;	ldf8	f122=[r18],32
650	ldf8	f123=[r19],32		};;
651{ .mmi;	ldf8	f124=[r34]
652	ldf8	f125=[r17]		}
653{ .mmi;	ldf8	f126=[r18]
654	ldf8	f127=[r19]		}
655
656{ .mmi;	ldf8	f33=[r14],32
657	ldf8	f34=[r15],32		}
658{ .mmi;	ldf8	f35=[r16],32;;
659	ldf8	f36=[r33]		}
660{ .mmi;	ldf8	f37=[r14]
661	ldf8	f38=[r15]		}
662{ .mfi;	ldf8	f39=[r16]
663// -------\ Entering multiplier's heaven /-------
664// ------------\                    /------------
665// -----------------\          /-----------------
666// ----------------------\/----------------------
667		xma.hu	f41=f32,f120,f0		}
668{ .mfi;		xma.lu	f40=f32,f120,f0		};; // (*)
669{ .mfi;		xma.hu	f51=f32,f121,f0		}
670{ .mfi;		xma.lu	f50=f32,f121,f0		};;
671{ .mfi;		xma.hu	f61=f32,f122,f0		}
672{ .mfi;		xma.lu	f60=f32,f122,f0		};;
673{ .mfi;		xma.hu	f71=f32,f123,f0		}
674{ .mfi;		xma.lu	f70=f32,f123,f0		};;
675{ .mfi;		xma.hu	f81=f32,f124,f0		}
676{ .mfi;		xma.lu	f80=f32,f124,f0		};;
677{ .mfi;		xma.hu	f91=f32,f125,f0		}
678{ .mfi;		xma.lu	f90=f32,f125,f0		};;
679{ .mfi;		xma.hu	f101=f32,f126,f0	}
680{ .mfi;		xma.lu	f100=f32,f126,f0	};;
681{ .mfi;		xma.hu	f111=f32,f127,f0	}
682{ .mfi;		xma.lu	f110=f32,f127,f0	};;//
683// (*)	You can argue that splitting at every second bundle would
684//	prevent "wider" IA-64 implementations from achieving the peak
685//	performance. Well, not really... The catch is that if you
686//	intend to keep 4 FP units busy by splitting at every fourth
687//	bundle and thus perform these 16 multiplications in 4 ticks,
688//	the first bundle *below* would stall because the result from
689//	the first xma bundle *above* won't be available for another 3
690//	ticks (if not more, being an optimist, I assume that "wider"
691//	implementation will have same latency:-). This stall will hold
692//	you back and the performance would be as if every second bundle
693//	were split *anyway*...
694{ .mfi;	getf.sig	r16=f40
695		xma.hu	f42=f33,f120,f41
696	add		r33=8,r32		}
697{ .mfi;		xma.lu	f41=f33,f120,f41	};;
698{ .mfi;	getf.sig	r24=f50
699		xma.hu	f52=f33,f121,f51	}
700{ .mfi;		xma.lu	f51=f33,f121,f51	};;
701{ .mfi;	st8		[r32]=r16,16
702		xma.hu	f62=f33,f122,f61	}
703{ .mfi;		xma.lu	f61=f33,f122,f61	};;
704{ .mfi;		xma.hu	f72=f33,f123,f71	}
705{ .mfi;		xma.lu	f71=f33,f123,f71	};;
706{ .mfi;		xma.hu	f82=f33,f124,f81	}
707{ .mfi;		xma.lu	f81=f33,f124,f81	};;
708{ .mfi;		xma.hu	f92=f33,f125,f91	}
709{ .mfi;		xma.lu	f91=f33,f125,f91	};;
710{ .mfi;		xma.hu	f102=f33,f126,f101	}
711{ .mfi;		xma.lu	f101=f33,f126,f101	};;
712{ .mfi;		xma.hu	f112=f33,f127,f111	}
713{ .mfi;		xma.lu	f111=f33,f127,f111	};;//
714//-------------------------------------------------//
715{ .mfi;	getf.sig	r25=f41
716		xma.hu	f43=f34,f120,f42	}
717{ .mfi;		xma.lu	f42=f34,f120,f42	};;
718{ .mfi;	getf.sig	r16=f60
719		xma.hu	f53=f34,f121,f52	}
720{ .mfi;		xma.lu	f52=f34,f121,f52	};;
721{ .mfi;	getf.sig	r17=f51
722		xma.hu	f63=f34,f122,f62
723	add		r25=r25,r24		}
724{ .mfi;		xma.lu	f62=f34,f122,f62
725	mov		carry1=0		};;
726{ .mfi;	cmp.ltu		p6,p0=r25,r24
727		xma.hu	f73=f34,f123,f72	}
728{ .mfi;		xma.lu	f72=f34,f123,f72	};;
729{ .mfi;	st8		[r33]=r25,16
730		xma.hu	f83=f34,f124,f82
731(p6)	add		carry1=1,carry1		}
732{ .mfi;		xma.lu	f82=f34,f124,f82	};;
733{ .mfi;		xma.hu	f93=f34,f125,f92	}
734{ .mfi;		xma.lu	f92=f34,f125,f92	};;
735{ .mfi;		xma.hu	f103=f34,f126,f102	}
736{ .mfi;		xma.lu	f102=f34,f126,f102	};;
737{ .mfi;		xma.hu	f113=f34,f127,f112	}
738{ .mfi;		xma.lu	f112=f34,f127,f112	};;//
739//-------------------------------------------------//
740{ .mfi;	getf.sig	r18=f42
741		xma.hu	f44=f35,f120,f43
742	add		r17=r17,r16		}
743{ .mfi;		xma.lu	f43=f35,f120,f43	};;
744{ .mfi;	getf.sig	r24=f70
745		xma.hu	f54=f35,f121,f53	}
746{ .mfi;	mov		carry2=0
747		xma.lu	f53=f35,f121,f53	};;
748{ .mfi;	getf.sig	r25=f61
749		xma.hu	f64=f35,f122,f63
750	cmp.ltu		p7,p0=r17,r16		}
751{ .mfi;	add		r18=r18,r17
752		xma.lu	f63=f35,f122,f63	};;
753{ .mfi;	getf.sig	r26=f52
754		xma.hu	f74=f35,f123,f73
755(p7)	add		carry2=1,carry2		}
756{ .mfi;	cmp.ltu		p7,p0=r18,r17
757		xma.lu	f73=f35,f123,f73
758	add		r18=r18,carry1		};;
759{ .mfi;
760		xma.hu	f84=f35,f124,f83
761(p7)	add		carry2=1,carry2		}
762{ .mfi;	cmp.ltu		p7,p0=r18,carry1
763		xma.lu	f83=f35,f124,f83	};;
764{ .mfi;	st8		[r32]=r18,16
765		xma.hu	f94=f35,f125,f93
766(p7)	add		carry2=1,carry2		}
767{ .mfi;		xma.lu	f93=f35,f125,f93	};;
768{ .mfi;		xma.hu	f104=f35,f126,f103	}
769{ .mfi;		xma.lu	f103=f35,f126,f103	};;
770{ .mfi;		xma.hu	f114=f35,f127,f113	}
771{ .mfi;	mov		carry1=0
772		xma.lu	f113=f35,f127,f113
773	add		r25=r25,r24		};;//
774//-------------------------------------------------//
775{ .mfi;	getf.sig	r27=f43
776		xma.hu	f45=f36,f120,f44
777	cmp.ltu		p6,p0=r25,r24		}
778{ .mfi;		xma.lu	f44=f36,f120,f44
779	add		r26=r26,r25		};;
780{ .mfi;	getf.sig	r16=f80
781		xma.hu	f55=f36,f121,f54
782(p6)	add		carry1=1,carry1		}
783{ .mfi;		xma.lu	f54=f36,f121,f54	};;
784{ .mfi;	getf.sig	r17=f71
785		xma.hu	f65=f36,f122,f64
786	cmp.ltu		p6,p0=r26,r25		}
787{ .mfi;		xma.lu	f64=f36,f122,f64
788	add		r27=r27,r26		};;
789{ .mfi;	getf.sig	r18=f62
790		xma.hu	f75=f36,f123,f74
791(p6)	add		carry1=1,carry1		}
792{ .mfi;	cmp.ltu		p6,p0=r27,r26
793		xma.lu	f74=f36,f123,f74
794	add		r27=r27,carry2		};;
795{ .mfi;	getf.sig	r19=f53
796		xma.hu	f85=f36,f124,f84
797(p6)	add		carry1=1,carry1		}
798{ .mfi;		xma.lu	f84=f36,f124,f84
799	cmp.ltu		p6,p0=r27,carry2	};;
800{ .mfi;	st8		[r33]=r27,16
801		xma.hu	f95=f36,f125,f94
802(p6)	add		carry1=1,carry1		}
803{ .mfi;		xma.lu	f94=f36,f125,f94	};;
804{ .mfi;		xma.hu	f105=f36,f126,f104	}
805{ .mfi;	mov		carry2=0
806		xma.lu	f104=f36,f126,f104
807	add		r17=r17,r16		};;
808{ .mfi;		xma.hu	f115=f36,f127,f114
809	cmp.ltu		p7,p0=r17,r16		}
810{ .mfi;		xma.lu	f114=f36,f127,f114
811	add		r18=r18,r17		};;//
812//-------------------------------------------------//
813{ .mfi;	getf.sig	r20=f44
814		xma.hu	f46=f37,f120,f45
815(p7)	add		carry2=1,carry2		}
816{ .mfi;	cmp.ltu		p7,p0=r18,r17
817		xma.lu	f45=f37,f120,f45
818	add		r19=r19,r18		};;
819{ .mfi;	getf.sig	r24=f90
820		xma.hu	f56=f37,f121,f55	}
821{ .mfi;		xma.lu	f55=f37,f121,f55	};;
822{ .mfi;	getf.sig	r25=f81
823		xma.hu	f66=f37,f122,f65
824(p7)	add		carry2=1,carry2		}
825{ .mfi;	cmp.ltu		p7,p0=r19,r18
826		xma.lu	f65=f37,f122,f65
827	add		r20=r20,r19		};;
828{ .mfi;	getf.sig	r26=f72
829		xma.hu	f76=f37,f123,f75
830(p7)	add		carry2=1,carry2		}
831{ .mfi;	cmp.ltu		p7,p0=r20,r19
832		xma.lu	f75=f37,f123,f75
833	add		r20=r20,carry1		};;
834{ .mfi;	getf.sig	r27=f63
835		xma.hu	f86=f37,f124,f85
836(p7)	add		carry2=1,carry2		}
837{ .mfi;		xma.lu	f85=f37,f124,f85
838	cmp.ltu		p7,p0=r20,carry1	};;
839{ .mfi;	getf.sig	r28=f54
840		xma.hu	f96=f37,f125,f95
841(p7)	add		carry2=1,carry2		}
842{ .mfi;	st8		[r32]=r20,16
843		xma.lu	f95=f37,f125,f95	};;
844{ .mfi;		xma.hu	f106=f37,f126,f105	}
845{ .mfi;	mov		carry1=0
846		xma.lu	f105=f37,f126,f105
847	add		r25=r25,r24		};;
848{ .mfi;		xma.hu	f116=f37,f127,f115
849	cmp.ltu		p6,p0=r25,r24		}
850{ .mfi;		xma.lu	f115=f37,f127,f115
851	add		r26=r26,r25		};;//
852//-------------------------------------------------//
853{ .mfi;	getf.sig	r29=f45
854		xma.hu	f47=f38,f120,f46
855(p6)	add		carry1=1,carry1		}
856{ .mfi;	cmp.ltu		p6,p0=r26,r25
857		xma.lu	f46=f38,f120,f46
858	add		r27=r27,r26		};;
859{ .mfi;	getf.sig	r16=f100
860		xma.hu	f57=f38,f121,f56
861(p6)	add		carry1=1,carry1		}
862{ .mfi;	cmp.ltu		p6,p0=r27,r26
863		xma.lu	f56=f38,f121,f56
864	add		r28=r28,r27		};;
865{ .mfi;	getf.sig	r17=f91
866		xma.hu	f67=f38,f122,f66
867(p6)	add		carry1=1,carry1		}
868{ .mfi;	cmp.ltu		p6,p0=r28,r27
869		xma.lu	f66=f38,f122,f66
870	add		r29=r29,r28		};;
871{ .mfi;	getf.sig	r18=f82
872		xma.hu	f77=f38,f123,f76
873(p6)	add		carry1=1,carry1		}
874{ .mfi;	cmp.ltu		p6,p0=r29,r28
875		xma.lu	f76=f38,f123,f76
876	add		r29=r29,carry2		};;
877{ .mfi;	getf.sig	r19=f73
878		xma.hu	f87=f38,f124,f86
879(p6)	add		carry1=1,carry1		}
880{ .mfi;		xma.lu	f86=f38,f124,f86
881	cmp.ltu		p6,p0=r29,carry2	};;
882{ .mfi;	getf.sig	r20=f64
883		xma.hu	f97=f38,f125,f96
884(p6)	add		carry1=1,carry1		}
885{ .mfi;	st8		[r33]=r29,16
886		xma.lu	f96=f38,f125,f96	};;
887{ .mfi;	getf.sig	r21=f55
888		xma.hu	f107=f38,f126,f106	}
889{ .mfi;	mov		carry2=0
890		xma.lu	f106=f38,f126,f106
891	add		r17=r17,r16		};;
892{ .mfi;		xma.hu	f117=f38,f127,f116
893	cmp.ltu		p7,p0=r17,r16		}
894{ .mfi;		xma.lu	f116=f38,f127,f116
895	add		r18=r18,r17		};;//
896//-------------------------------------------------//
897{ .mfi;	getf.sig	r22=f46
898		xma.hu	f48=f39,f120,f47
899(p7)	add		carry2=1,carry2		}
900{ .mfi;	cmp.ltu		p7,p0=r18,r17
901		xma.lu	f47=f39,f120,f47
902	add		r19=r19,r18		};;
903{ .mfi;	getf.sig	r24=f110
904		xma.hu	f58=f39,f121,f57
905(p7)	add		carry2=1,carry2		}
906{ .mfi;	cmp.ltu		p7,p0=r19,r18
907		xma.lu	f57=f39,f121,f57
908	add		r20=r20,r19		};;
909{ .mfi;	getf.sig	r25=f101
910		xma.hu	f68=f39,f122,f67
911(p7)	add		carry2=1,carry2		}
912{ .mfi;	cmp.ltu		p7,p0=r20,r19
913		xma.lu	f67=f39,f122,f67
914	add		r21=r21,r20		};;
915{ .mfi;	getf.sig	r26=f92
916		xma.hu	f78=f39,f123,f77
917(p7)	add		carry2=1,carry2		}
918{ .mfi;	cmp.ltu		p7,p0=r21,r20
919		xma.lu	f77=f39,f123,f77
920	add		r22=r22,r21		};;
921{ .mfi;	getf.sig	r27=f83
922		xma.hu	f88=f39,f124,f87
923(p7)	add		carry2=1,carry2		}
924{ .mfi;	cmp.ltu		p7,p0=r22,r21
925		xma.lu	f87=f39,f124,f87
926	add		r22=r22,carry1		};;
927{ .mfi;	getf.sig	r28=f74
928		xma.hu	f98=f39,f125,f97
929(p7)	add		carry2=1,carry2		}
930{ .mfi;		xma.lu	f97=f39,f125,f97
931	cmp.ltu		p7,p0=r22,carry1	};;
932{ .mfi;	getf.sig	r29=f65
933		xma.hu	f108=f39,f126,f107
934(p7)	add		carry2=1,carry2		}
935{ .mfi;	st8		[r32]=r22,16
936		xma.lu	f107=f39,f126,f107	};;
937{ .mfi;	getf.sig	r30=f56
938		xma.hu	f118=f39,f127,f117	}
939{ .mfi;		xma.lu	f117=f39,f127,f117	};;//
940//-------------------------------------------------//
941// Leaving multiplier's heaven... Quite a ride, huh?
942
943{ .mii;	getf.sig	r31=f47
944	add		r25=r25,r24
945	mov		carry1=0		};;
946{ .mii;		getf.sig	r16=f111
947	cmp.ltu		p6,p0=r25,r24
948	add		r26=r26,r25		};;
949{ .mfb;		getf.sig	r17=f102	}
950{ .mii;
951(p6)	add		carry1=1,carry1
952	cmp.ltu		p6,p0=r26,r25
953	add		r27=r27,r26		};;
954{ .mfb;	nop.m	0x0				}
955{ .mii;
956(p6)	add		carry1=1,carry1
957	cmp.ltu		p6,p0=r27,r26
958	add		r28=r28,r27		};;
959{ .mii;		getf.sig	r18=f93
960		add		r17=r17,r16
961		mov		carry3=0	}
962{ .mii;
963(p6)	add		carry1=1,carry1
964	cmp.ltu		p6,p0=r28,r27
965	add		r29=r29,r28		};;
966{ .mii;		getf.sig	r19=f84
967		cmp.ltu		p7,p0=r17,r16	}
968{ .mii;
969(p6)	add		carry1=1,carry1
970	cmp.ltu		p6,p0=r29,r28
971	add		r30=r30,r29		};;
972{ .mii;		getf.sig	r20=f75
973		add		r18=r18,r17	}
974{ .mii;
975(p6)	add		carry1=1,carry1
976	cmp.ltu		p6,p0=r30,r29
977	add		r31=r31,r30		};;
978{ .mfb;		getf.sig	r21=f66		}
979{ .mii;	(p7)	add		carry3=1,carry3
980		cmp.ltu		p7,p0=r18,r17
981		add		r19=r19,r18	}
982{ .mfb;	nop.m	0x0				}
983{ .mii;
984(p6)	add		carry1=1,carry1
985	cmp.ltu		p6,p0=r31,r30
986	add		r31=r31,carry2		};;
987{ .mfb;		getf.sig	r22=f57		}
988{ .mii;	(p7)	add		carry3=1,carry3
989		cmp.ltu		p7,p0=r19,r18
990		add		r20=r20,r19	}
991{ .mfb;	nop.m	0x0				}
992{ .mii;
993(p6)	add		carry1=1,carry1
994	cmp.ltu		p6,p0=r31,carry2	};;
995{ .mfb;		getf.sig	r23=f48		}
996{ .mii;	(p7)	add		carry3=1,carry3
997		cmp.ltu		p7,p0=r20,r19
998		add		r21=r21,r20	}
999{ .mii;
1000(p6)	add		carry1=1,carry1		}
1001{ .mfb;	st8		[r33]=r31,16		};;
1002
1003{ .mfb;	getf.sig	r24=f112		}
1004{ .mii;	(p7)	add		carry3=1,carry3
1005		cmp.ltu		p7,p0=r21,r20
1006		add		r22=r22,r21	};;
1007{ .mfb;	getf.sig	r25=f103		}
1008{ .mii;	(p7)	add		carry3=1,carry3
1009		cmp.ltu		p7,p0=r22,r21
1010		add		r23=r23,r22	};;
1011{ .mfb;	getf.sig	r26=f94			}
1012{ .mii;	(p7)	add		carry3=1,carry3
1013		cmp.ltu		p7,p0=r23,r22
1014		add		r23=r23,carry1	};;
1015{ .mfb;	getf.sig	r27=f85			}
1016{ .mii;	(p7)	add		carry3=1,carry3
1017		cmp.ltu		p7,p8=r23,carry1};;
1018{ .mii;	getf.sig	r28=f76
1019	add		r25=r25,r24
1020	mov		carry1=0		}
1021{ .mii;		st8		[r32]=r23,16
1022	(p7)	add		carry2=1,carry3
1023	(p8)	add		carry2=0,carry3	};;
1024
1025{ .mfb;	nop.m	0x0				}
1026{ .mii;	getf.sig	r29=f67
1027	cmp.ltu		p6,p0=r25,r24
1028	add		r26=r26,r25		};;
1029{ .mfb;	getf.sig	r30=f58			}
1030{ .mii;
1031(p6)	add		carry1=1,carry1
1032	cmp.ltu		p6,p0=r26,r25
1033	add		r27=r27,r26		};;
1034{ .mfb;		getf.sig	r16=f113	}
1035{ .mii;
1036(p6)	add		carry1=1,carry1
1037	cmp.ltu		p6,p0=r27,r26
1038	add		r28=r28,r27		};;
1039{ .mfb;		getf.sig	r17=f104	}
1040{ .mii;
1041(p6)	add		carry1=1,carry1
1042	cmp.ltu		p6,p0=r28,r27
1043	add		r29=r29,r28		};;
1044{ .mfb;		getf.sig	r18=f95		}
1045{ .mii;
1046(p6)	add		carry1=1,carry1
1047	cmp.ltu		p6,p0=r29,r28
1048	add		r30=r30,r29		};;
1049{ .mii;		getf.sig	r19=f86
1050		add		r17=r17,r16
1051		mov		carry3=0	}
1052{ .mii;
1053(p6)	add		carry1=1,carry1
1054	cmp.ltu		p6,p0=r30,r29
1055	add		r30=r30,carry2		};;
1056{ .mii;		getf.sig	r20=f77
1057		cmp.ltu		p7,p0=r17,r16
1058		add		r18=r18,r17	}
1059{ .mii;
1060(p6)	add		carry1=1,carry1
1061	cmp.ltu		p6,p0=r30,carry2	};;
1062{ .mfb;		getf.sig	r21=f68		}
1063{ .mii;	st8		[r33]=r30,16
1064(p6)	add		carry1=1,carry1		};;
1065
1066{ .mfb;	getf.sig	r24=f114		}
1067{ .mii;	(p7)	add		carry3=1,carry3
1068		cmp.ltu		p7,p0=r18,r17
1069		add		r19=r19,r18	};;
1070{ .mfb;	getf.sig	r25=f105		}
1071{ .mii;	(p7)	add		carry3=1,carry3
1072		cmp.ltu		p7,p0=r19,r18
1073		add		r20=r20,r19	};;
1074{ .mfb;	getf.sig	r26=f96			}
1075{ .mii;	(p7)	add		carry3=1,carry3
1076		cmp.ltu		p7,p0=r20,r19
1077		add		r21=r21,r20	};;
1078{ .mfb;	getf.sig	r27=f87			}
1079{ .mii;	(p7)	add		carry3=1,carry3
1080		cmp.ltu		p7,p0=r21,r20
1081		add		r21=r21,carry1	};;
1082{ .mib;	getf.sig	r28=f78
1083	add		r25=r25,r24		}
1084{ .mib;	(p7)	add		carry3=1,carry3
1085		cmp.ltu		p7,p8=r21,carry1};;
1086{ .mii;		st8		[r32]=r21,16
1087	(p7)	add		carry2=1,carry3
1088	(p8)	add		carry2=0,carry3	}
1089
1090{ .mii;	mov		carry1=0
1091	cmp.ltu		p6,p0=r25,r24
1092	add		r26=r26,r25		};;
1093{ .mfb;		getf.sig	r16=f115	}
1094{ .mii;
1095(p6)	add		carry1=1,carry1
1096	cmp.ltu		p6,p0=r26,r25
1097	add		r27=r27,r26		};;
1098{ .mfb;		getf.sig	r17=f106	}
1099{ .mii;
1100(p6)	add		carry1=1,carry1
1101	cmp.ltu		p6,p0=r27,r26
1102	add		r28=r28,r27		};;
1103{ .mfb;		getf.sig	r18=f97		}
1104{ .mii;
1105(p6)	add		carry1=1,carry1
1106	cmp.ltu		p6,p0=r28,r27
1107	add		r28=r28,carry2		};;
1108{ .mib;		getf.sig	r19=f88
1109		add		r17=r17,r16	}
1110{ .mib;
1111(p6)	add		carry1=1,carry1
1112	cmp.ltu		p6,p0=r28,carry2	};;
1113{ .mii;	st8		[r33]=r28,16
1114(p6)	add		carry1=1,carry1		}
1115
1116{ .mii;		mov		carry2=0
1117		cmp.ltu		p7,p0=r17,r16
1118		add		r18=r18,r17	};;
1119{ .mfb;	getf.sig	r24=f116		}
1120{ .mii;	(p7)	add		carry2=1,carry2
1121		cmp.ltu		p7,p0=r18,r17
1122		add		r19=r19,r18	};;
1123{ .mfb;	getf.sig	r25=f107		}
1124{ .mii;	(p7)	add		carry2=1,carry2
1125		cmp.ltu		p7,p0=r19,r18
1126		add		r19=r19,carry1	};;
1127{ .mfb;	getf.sig	r26=f98			}
1128{ .mii;	(p7)	add		carry2=1,carry2
1129		cmp.ltu		p7,p0=r19,carry1};;
1130{ .mii;		st8		[r32]=r19,16
1131	(p7)	add		carry2=1,carry2	}
1132
1133{ .mfb;	add		r25=r25,r24		};;
1134
1135{ .mfb;		getf.sig	r16=f117	}
1136{ .mii;	mov		carry1=0
1137	cmp.ltu		p6,p0=r25,r24
1138	add		r26=r26,r25		};;
1139{ .mfb;		getf.sig	r17=f108	}
1140{ .mii;
1141(p6)	add		carry1=1,carry1
1142	cmp.ltu		p6,p0=r26,r25
1143	add		r26=r26,carry2		};;
1144{ .mfb;	nop.m	0x0				}
1145{ .mii;
1146(p6)	add		carry1=1,carry1
1147	cmp.ltu		p6,p0=r26,carry2	};;
1148{ .mii;	st8		[r33]=r26,16
1149(p6)	add		carry1=1,carry1		}
1150
1151{ .mfb;		add		r17=r17,r16	};;
1152{ .mfb;	getf.sig	r24=f118		}
1153{ .mii;		mov		carry2=0
1154		cmp.ltu		p7,p0=r17,r16
1155		add		r17=r17,carry1	};;
1156{ .mii;	(p7)	add		carry2=1,carry2
1157		cmp.ltu		p7,p0=r17,carry1};;
1158{ .mii;		st8		[r32]=r17
1159	(p7)	add		carry2=1,carry2	};;
1160{ .mfb;	add		r24=r24,carry2		};;
1161{ .mib;	st8		[r33]=r24		}
1162
1163{ .mib;	rum		1<<5		// clear um.mfh
1164	br.ret.sptk.many	b0	};;
1165.endp	bn_mul_comba8#
1166#undef	carry3
1167#undef	carry2
1168#undef	carry1
1169#endif
1170
1171#if 1
1172// It's possible to make it faster (see comment to bn_sqr_comba8), but
1173// I reckon it doesn't worth the effort. Basically because the routine
1174// (actually both of them) practically never called... So I just play
1175// same trick as with bn_sqr_comba8.
1176//
1177// void bn_sqr_comba4(BN_ULONG *r, BN_ULONG *a)
1178//
1179.global	bn_sqr_comba4#
1180.proc	bn_sqr_comba4#
1181.align	64
1182bn_sqr_comba4:
1183	.prologue
1184	.save	ar.pfs,r2
1185#if defined(_HPUX_SOURCE) && !defined(_LP64)
1186{ .mii;	alloc   r2=ar.pfs,2,1,0,0
1187	addp4	r32=0,r32
1188	addp4	r33=0,r33		};;
1189{ .mii;
1190#else
1191{ .mii;	alloc	r2=ar.pfs,2,1,0,0
1192#endif
1193	mov	r34=r33
1194	add	r14=8,r33		};;
1195	.body
1196{ .mii;	add	r17=8,r34
1197	add	r15=16,r33
1198	add	r18=16,r34		}
1199{ .mfb;	add	r16=24,r33
1200	br	.L_cheat_entry_point4	};;
1201.endp	bn_sqr_comba4#
1202#endif
1203
1204#if 1
1205// Runs in ~115 cycles and ~4.5 times faster than C. Well, whatever...
1206//
1207// void bn_mul_comba4(BN_ULONG *r, BN_ULONG *a, BN_ULONG *b)
1208//
1209#define	carry1	r14
1210#define	carry2	r15
1211.global	bn_mul_comba4#
1212.proc	bn_mul_comba4#
1213.align	64
1214bn_mul_comba4:
1215	.prologue
1216	.save	ar.pfs,r2
1217#if defined(_HPUX_SOURCE) && !defined(_LP64)
1218{ .mii;	alloc   r2=ar.pfs,3,0,0,0
1219	addp4	r33=0,r33
1220	addp4	r34=0,r34		};;
1221{ .mii;	addp4	r32=0,r32
1222#else
1223{ .mii;	alloc	r2=ar.pfs,3,0,0,0
1224#endif
1225	add	r14=8,r33
1226	add	r17=8,r34		}
1227	.body
1228{ .mii;	add	r15=16,r33
1229	add	r18=16,r34
1230	add	r16=24,r33		};;
1231.L_cheat_entry_point4:
1232{ .mmi;	add	r19=24,r34
1233
1234	ldf8	f32=[r33]		}
1235
1236{ .mmi;	ldf8	f120=[r34]
1237	ldf8	f121=[r17]		};;
1238{ .mmi;	ldf8	f122=[r18]
1239	ldf8	f123=[r19]		}
1240
1241{ .mmi;	ldf8	f33=[r14]
1242	ldf8	f34=[r15]		}
1243{ .mfi;	ldf8	f35=[r16]
1244
1245		xma.hu	f41=f32,f120,f0		}
1246{ .mfi;		xma.lu	f40=f32,f120,f0		};;
1247{ .mfi;		xma.hu	f51=f32,f121,f0		}
1248{ .mfi;		xma.lu	f50=f32,f121,f0		};;
1249{ .mfi;		xma.hu	f61=f32,f122,f0		}
1250{ .mfi;		xma.lu	f60=f32,f122,f0		};;
1251{ .mfi;		xma.hu	f71=f32,f123,f0		}
1252{ .mfi;		xma.lu	f70=f32,f123,f0		};;//
1253// Major stall takes place here, and 3 more places below. Result from
1254// first xma is not available for another 3 ticks.
1255{ .mfi;	getf.sig	r16=f40
1256		xma.hu	f42=f33,f120,f41
1257	add		r33=8,r32		}
1258{ .mfi;		xma.lu	f41=f33,f120,f41	};;
1259{ .mfi;	getf.sig	r24=f50
1260		xma.hu	f52=f33,f121,f51	}
1261{ .mfi;		xma.lu	f51=f33,f121,f51	};;
1262{ .mfi;	st8		[r32]=r16,16
1263		xma.hu	f62=f33,f122,f61	}
1264{ .mfi;		xma.lu	f61=f33,f122,f61	};;
1265{ .mfi;		xma.hu	f72=f33,f123,f71	}
1266{ .mfi;		xma.lu	f71=f33,f123,f71	};;//
1267//-------------------------------------------------//
1268{ .mfi;	getf.sig	r25=f41
1269		xma.hu	f43=f34,f120,f42	}
1270{ .mfi;		xma.lu	f42=f34,f120,f42	};;
1271{ .mfi;	getf.sig	r16=f60
1272		xma.hu	f53=f34,f121,f52	}
1273{ .mfi;		xma.lu	f52=f34,f121,f52	};;
1274{ .mfi;	getf.sig	r17=f51
1275		xma.hu	f63=f34,f122,f62
1276	add		r25=r25,r24		}
1277{ .mfi;	mov		carry1=0
1278		xma.lu	f62=f34,f122,f62	};;
1279{ .mfi;	st8		[r33]=r25,16
1280		xma.hu	f73=f34,f123,f72
1281	cmp.ltu		p6,p0=r25,r24		}
1282{ .mfi;		xma.lu	f72=f34,f123,f72	};;//
1283//-------------------------------------------------//
1284{ .mfi;	getf.sig	r18=f42
1285		xma.hu	f44=f35,f120,f43
1286(p6)	add		carry1=1,carry1		}
1287{ .mfi;	add		r17=r17,r16
1288		xma.lu	f43=f35,f120,f43
1289	mov		carry2=0		};;
1290{ .mfi;	getf.sig	r24=f70
1291		xma.hu	f54=f35,f121,f53
1292	cmp.ltu		p7,p0=r17,r16		}
1293{ .mfi;		xma.lu	f53=f35,f121,f53	};;
1294{ .mfi;	getf.sig	r25=f61
1295		xma.hu	f64=f35,f122,f63
1296	add		r18=r18,r17		}
1297{ .mfi;		xma.lu	f63=f35,f122,f63
1298(p7)	add		carry2=1,carry2		};;
1299{ .mfi;	getf.sig	r26=f52
1300		xma.hu	f74=f35,f123,f73
1301	cmp.ltu		p7,p0=r18,r17		}
1302{ .mfi;		xma.lu	f73=f35,f123,f73
1303	add		r18=r18,carry1		};;
1304//-------------------------------------------------//
1305{ .mii;	st8		[r32]=r18,16
1306(p7)	add		carry2=1,carry2
1307	cmp.ltu		p7,p0=r18,carry1	};;
1308
1309{ .mfi;	getf.sig	r27=f43	// last major stall
1310(p7)	add		carry2=1,carry2		};;
1311{ .mii;		getf.sig	r16=f71
1312	add		r25=r25,r24
1313	mov		carry1=0		};;
1314{ .mii;		getf.sig	r17=f62
1315	cmp.ltu		p6,p0=r25,r24
1316	add		r26=r26,r25		};;
1317{ .mii;
1318(p6)	add		carry1=1,carry1
1319	cmp.ltu		p6,p0=r26,r25
1320	add		r27=r27,r26		};;
1321{ .mii;
1322(p6)	add		carry1=1,carry1
1323	cmp.ltu		p6,p0=r27,r26
1324	add		r27=r27,carry2		};;
1325{ .mii;		getf.sig	r18=f53
1326(p6)	add		carry1=1,carry1
1327	cmp.ltu		p6,p0=r27,carry2	};;
1328{ .mfi;	st8		[r33]=r27,16
1329(p6)	add		carry1=1,carry1		}
1330
1331{ .mii;		getf.sig	r19=f44
1332		add		r17=r17,r16
1333		mov		carry2=0	};;
1334{ .mii;	getf.sig	r24=f72
1335		cmp.ltu		p7,p0=r17,r16
1336		add		r18=r18,r17	};;
1337{ .mii;	(p7)	add		carry2=1,carry2
1338		cmp.ltu		p7,p0=r18,r17
1339		add		r19=r19,r18	};;
1340{ .mii;	(p7)	add		carry2=1,carry2
1341		cmp.ltu		p7,p0=r19,r18
1342		add		r19=r19,carry1	};;
1343{ .mii;	getf.sig	r25=f63
1344	(p7)	add		carry2=1,carry2
1345		cmp.ltu		p7,p0=r19,carry1};;
1346{ .mii;		st8		[r32]=r19,16
1347	(p7)	add		carry2=1,carry2	}
1348
1349{ .mii;	getf.sig	r26=f54
1350	add		r25=r25,r24
1351	mov		carry1=0		};;
1352{ .mii;		getf.sig	r16=f73
1353	cmp.ltu		p6,p0=r25,r24
1354	add		r26=r26,r25		};;
1355{ .mii;
1356(p6)	add		carry1=1,carry1
1357	cmp.ltu		p6,p0=r26,r25
1358	add		r26=r26,carry2		};;
1359{ .mii;		getf.sig	r17=f64
1360(p6)	add		carry1=1,carry1
1361	cmp.ltu		p6,p0=r26,carry2	};;
1362{ .mii;	st8		[r33]=r26,16
1363(p6)	add		carry1=1,carry1		}
1364
1365{ .mii;	getf.sig	r24=f74
1366		add		r17=r17,r16
1367		mov		carry2=0	};;
1368{ .mii;		cmp.ltu		p7,p0=r17,r16
1369		add		r17=r17,carry1	};;
1370
1371{ .mii;	(p7)	add		carry2=1,carry2
1372		cmp.ltu		p7,p0=r17,carry1};;
1373{ .mii;		st8		[r32]=r17,16
1374	(p7)	add		carry2=1,carry2	};;
1375
1376{ .mii;	add		r24=r24,carry2		};;
1377{ .mii;	st8		[r33]=r24		}
1378
1379{ .mib;	rum		1<<5		// clear um.mfh
1380	br.ret.sptk.many	b0	};;
1381.endp	bn_mul_comba4#
1382#undef	carry2
1383#undef	carry1
1384#endif
1385
1386#if 1
1387//
1388// BN_ULONG bn_div_words(BN_ULONG h, BN_ULONG l, BN_ULONG d)
1389//
1390// In the nutshell it's a port of my MIPS III/IV implementation.
1391//
1392#define	AT	r14
1393#define	H	r16
1394#define	HH	r20
1395#define	L	r17
1396#define	D	r18
1397#define	DH	r22
1398#define	I	r21
1399
1400#if 0
1401// Some preprocessors (most notably HP-UX) appear to be allergic to
1402// macros enclosed to parenthesis [as these three were].
1403#define	cont	p16
1404#define	break	p0	// p20
1405#define	equ	p24
1406#else
1407cont=p16
1408break=p0
1409equ=p24
1410#endif
1411
1412.global	abort#
1413.global	bn_div_words#
1414.proc	bn_div_words#
1415.align	64
1416bn_div_words:
1417	.prologue
1418	.save	ar.pfs,r2
1419{ .mii;	alloc		r2=ar.pfs,3,5,0,8
1420	.save	b0,r3
1421	mov		r3=b0
1422	.save	pr,r10
1423	mov		r10=pr		};;
1424{ .mmb;	cmp.eq		p6,p0=r34,r0
1425	mov		r8=-1
1426(p6)	br.ret.spnt.many	b0	};;
1427
1428	.body
1429{ .mii;	mov		H=r32		// save h
1430	mov		ar.ec=0		// don't rotate at exit
1431	mov		pr.rot=0	}
1432{ .mii;	mov		L=r33		// save l
1433	mov		r25=r0		// needed if abort is called on VMS
1434	mov		r36=r0		};;
1435
1436.L_divw_shift:	// -vv- note signed comparison
1437{ .mfi;	(p0)	cmp.lt		p16,p0=r0,r34	// d
1438	(p0)	shladd		r33=r34,1,r0	}
1439{ .mfb;	(p0)	add		r35=1,r36
1440	(p0)	nop.f		0x0
1441(p16)	br.wtop.dpnt		.L_divw_shift	};;
1442
1443{ .mii;	mov		D=r34
1444	shr.u		DH=r34,32
1445	sub		r35=64,r36		};;
1446{ .mii;	setf.sig	f7=DH
1447	shr.u		AT=H,r35
1448	mov		I=r36			};;
1449{ .mib;	cmp.ne		p6,p0=r0,AT
1450	shl		H=H,r36
1451(p6)	br.call.spnt.clr	b0=abort	};;	// overflow, die...
1452
1453{ .mfi;	fcvt.xuf.s1	f7=f7
1454	shr.u		AT=L,r35		};;
1455{ .mii;	shl		L=L,r36
1456	or		H=H,AT			};;
1457
1458{ .mii;	nop.m		0x0
1459	cmp.leu		p6,p0=D,H;;
1460(p6)	sub		H=H,D			}
1461
1462{ .mlx;	setf.sig	f14=D
1463	movl		AT=0xffffffff		};;
1464///////////////////////////////////////////////////////////
1465{ .mii;	setf.sig	f6=H
1466	shr.u		HH=H,32;;
1467	cmp.eq		p6,p7=HH,DH		};;
1468{ .mfb;
1469(p6)	setf.sig	f8=AT
1470(p7)	fcvt.xuf.s1	f6=f6
1471(p7)	br.call.sptk	b6=.L_udiv64_32_b6	};;
1472
1473{ .mfi;	getf.sig	r33=f8				// q
1474	xmpy.lu		f9=f8,f14		}
1475{ .mfi;	xmpy.hu		f10=f8,f14
1476	shrp		H=H,L,32		};;
1477
1478{ .mmi;	getf.sig	r35=f9				// tl
1479	getf.sig	r31=f10			};;	// th
1480
1481.L_divw_1st_iter:
1482{ .mii;	(p0)	add		r32=-1,r33
1483	(p0)	cmp.eq		equ,cont=HH,r31		};;
1484{ .mii;	(p0)	cmp.ltu		p8,p0=r35,D
1485	(p0)	sub		r34=r35,D
1486	(equ)	cmp.leu		break,cont=r35,H	};;
1487{ .mib;	(cont)	cmp.leu		cont,break=HH,r31
1488	(p8)	add		r31=-1,r31
1489(cont)	br.wtop.spnt		.L_divw_1st_iter	};;
1490///////////////////////////////////////////////////////////
1491{ .mii;	sub		H=H,r35
1492	shl		r8=r33,32
1493	shl		L=L,32			};;
1494///////////////////////////////////////////////////////////
1495{ .mii;	setf.sig	f6=H
1496	shr.u		HH=H,32;;
1497	cmp.eq		p6,p7=HH,DH		};;
1498{ .mfb;
1499(p6)	setf.sig	f8=AT
1500(p7)	fcvt.xuf.s1	f6=f6
1501(p7)	br.call.sptk	b6=.L_udiv64_32_b6	};;
1502
1503{ .mfi;	getf.sig	r33=f8				// q
1504	xmpy.lu		f9=f8,f14		}
1505{ .mfi;	xmpy.hu		f10=f8,f14
1506	shrp		H=H,L,32		};;
1507
1508{ .mmi;	getf.sig	r35=f9				// tl
1509	getf.sig	r31=f10			};;	// th
1510
1511.L_divw_2nd_iter:
1512{ .mii;	(p0)	add		r32=-1,r33
1513	(p0)	cmp.eq		equ,cont=HH,r31		};;
1514{ .mii;	(p0)	cmp.ltu		p8,p0=r35,D
1515	(p0)	sub		r34=r35,D
1516	(equ)	cmp.leu		break,cont=r35,H	};;
1517{ .mib;	(cont)	cmp.leu		cont,break=HH,r31
1518	(p8)	add		r31=-1,r31
1519(cont)	br.wtop.spnt		.L_divw_2nd_iter	};;
1520///////////////////////////////////////////////////////////
1521{ .mii;	sub	H=H,r35
1522	or	r8=r8,r33
1523	mov	ar.pfs=r2		};;
1524{ .mii;	shr.u	r9=H,I			// remainder if anybody wants it
1525	mov	pr=r10,0x1ffff		}
1526{ .mfb;	br.ret.sptk.many	b0	};;
1527
1528// Unsigned 64 by 32 (well, by 64 for the moment) bit integer division
1529// procedure.
1530//
1531// inputs:	f6 = (double)a, f7 = (double)b
1532// output:	f8 = (int)(a/b)
1533// clobbered:	f8,f9,f10,f11,pred
1534pred=p15
1535// This snippet is based on text found in the "Divide, Square
1536// Root and Remainder" section at
1537// http://www.intel.com/software/products/opensource/libraries/num.htm.
1538// Yes, I admit that the referred code was used as template,
1539// but after I realized that there hardly is any other instruction
1540// sequence which would perform this operation. I mean I figure that
1541// any independent attempt to implement high-performance division
1542// will result in code virtually identical to the Intel code. It
1543// should be noted though that below division kernel is 1 cycle
1544// faster than Intel one (note commented splits:-), not to mention
1545// original prologue (rather lack of one) and epilogue.
1546.align	32
1547.skip	16
1548.L_udiv64_32_b6:
1549	frcpa.s1	f8,pred=f6,f7;;		// [0]  y0 = 1 / b
1550
1551(pred)	fnma.s1		f9=f7,f8,f1		// [5]  e0 = 1 - b * y0
1552(pred)	fmpy.s1		f10=f6,f8;;		// [5]  q0 = a * y0
1553(pred)	fmpy.s1		f11=f9,f9		// [10] e1 = e0 * e0
1554(pred)	fma.s1		f10=f9,f10,f10;;	// [10] q1 = q0 + e0 * q0
1555(pred)	fma.s1		f8=f9,f8,f8	//;;	// [15] y1 = y0 + e0 * y0
1556(pred)	fma.s1		f9=f11,f10,f10;;	// [15] q2 = q1 + e1 * q1
1557(pred)	fma.s1		f8=f11,f8,f8	//;;	// [20] y2 = y1 + e1 * y1
1558(pred)	fnma.s1		f10=f7,f9,f6;;		// [20] r2 = a - b * q2
1559(pred)	fma.s1		f8=f10,f8,f9;;		// [25] q3 = q2 + r2 * y2
1560
1561	fcvt.fxu.trunc.s1	f8=f8		// [30] q = trunc(q3)
1562	br.ret.sptk.many	b6;;
1563.endp	bn_div_words#
1564#endif
1565