CoCalc -- stwotox.S

GitHub Repository: awilliam/linux-vfio
Path: blob/master/arch/m68k/fpsp040/stwotox.S
¹⁰⁸¹⁷ views
1
|
2
|	stwotox.sa 3.1 12/10/90
3
|
4
|	stwotox  --- 2**X
5
|	stwotoxd --- 2**X for denormalized X
6
|	stentox  --- 10**X
7
|	stentoxd --- 10**X for denormalized X
8
|
9
|	Input: Double-extended number X in location pointed to
10
|		by address register a0.
11
|
12
|	Output: The function values are returned in Fp0.
13
|
14
|	Accuracy and Monotonicity: The returned result is within 2 ulps in
15
|		64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
16
|		result is subsequently rounded to double precision. The
17
|		result is provably monotonic in double precision.
18
|
19
|	Speed: The program stwotox takes approximately 190 cycles and the
20
|		program stentox takes approximately 200 cycles.
21
|
22
|	Algorithm:
23
|
24
|	twotox
25
|	1. If |X| > 16480, go to ExpBig.
26
|
27
|	2. If |X| < 2**(-70), go to ExpSm.
28
|
29
|	3. Decompose X as X = N/64 + r where |r| <= 1/128. Furthermore
30
|		decompose N as
31
|		 N = 64(M + M') + j,  j = 0,1,2,...,63.
32
|
33
|	4. Overwrite r := r * log2. Then
34
|		2**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
35
|		Go to expr to compute that expression.
36
|
37
|	tentox
38
|	1. If |X| > 16480*log_10(2) (base 10 log of 2), go to ExpBig.
39
|
40
|	2. If |X| < 2**(-70), go to ExpSm.
41
|
42
|	3. Set y := X*log_2(10)*64 (base 2 log of 10). Set
43
|		N := round-to-int(y). Decompose N as
44
|		 N = 64(M + M') + j,  j = 0,1,2,...,63.
45
|
46
|	4. Define r as
47
|		r := ((X - N*L1)-N*L2) * L10
48
|		where L1, L2 are the leading and trailing parts of log_10(2)/64
49
|		and L10 is the natural log of 10. Then
50
|		10**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
51
|		Go to expr to compute that expression.
52
|
53
|	expr
54
|	1. Fetch 2**(j/64) from table as Fact1 and Fact2.
55
|
56
|	2. Overwrite Fact1 and Fact2 by
57
|		Fact1 := 2**(M) * Fact1
58
|		Fact2 := 2**(M) * Fact2
59
|		Thus Fact1 + Fact2 = 2**(M) * 2**(j/64).
60
|
61
|	3. Calculate P where 1 + P approximates exp(r):
62
|		P = r + r*r*(A1+r*(A2+...+r*A5)).
63
|
64
|	4. Let AdjFact := 2**(M'). Return
65
|		AdjFact * ( Fact1 + ((Fact1*P) + Fact2) ).
66
|		Exit.
67
|
68
|	ExpBig
69
|	1. Generate overflow by Huge * Huge if X > 0; otherwise, generate
70
|		underflow by Tiny * Tiny.
71
|
72
|	ExpSm
73
|	1. Return 1 + X.
74
|
75

76
|		Copyright (C) Motorola, Inc. 1990
77
|			All Rights Reserved
78
|
79
|       For details on the license for this file, please see the
80
|       file, README, in this same directory.
81

82
|STWOTOX	idnt	2,1 | Motorola 040 Floating Point Software Package
83

84
	|section	8
85

86
#include "fpsp.h"
87

88
BOUNDS1:	.long 0x3FB98000,0x400D80C0 | ... 2^(-70),16480
89
BOUNDS2:	.long 0x3FB98000,0x400B9B07 | ... 2^(-70),16480 LOG2/LOG10
90

91
L2TEN64:	.long 0x406A934F,0x0979A371 | ... 64LOG10/LOG2
92
L10TWO1:	.long 0x3F734413,0x509F8000 | ... LOG2/64LOG10
93

94
L10TWO2:	.long 0xBFCD0000,0xC0219DC1,0xDA994FD2,0x00000000
95

96
LOG10:	.long 0x40000000,0x935D8DDD,0xAAA8AC17,0x00000000
97

98
LOG2:	.long 0x3FFE0000,0xB17217F7,0xD1CF79AC,0x00000000
99

100
EXPA5:	.long 0x3F56C16D,0x6F7BD0B2
101
EXPA4:	.long 0x3F811112,0x302C712C
102
EXPA3:	.long 0x3FA55555,0x55554CC1
103
EXPA2:	.long 0x3FC55555,0x55554A54
104
EXPA1:	.long 0x3FE00000,0x00000000,0x00000000,0x00000000
105

106
HUGE:	.long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
107
TINY:	.long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
108

109
EXPTBL:
110
	.long  0x3FFF0000,0x80000000,0x00000000,0x3F738000
111
	.long  0x3FFF0000,0x8164D1F3,0xBC030773,0x3FBEF7CA
112
	.long  0x3FFF0000,0x82CD8698,0xAC2BA1D7,0x3FBDF8A9
113
	.long  0x3FFF0000,0x843A28C3,0xACDE4046,0x3FBCD7C9
114
	.long  0x3FFF0000,0x85AAC367,0xCC487B15,0xBFBDE8DA
115
	.long  0x3FFF0000,0x871F6196,0x9E8D1010,0x3FBDE85C
116
	.long  0x3FFF0000,0x88980E80,0x92DA8527,0x3FBEBBF1
117
	.long  0x3FFF0000,0x8A14D575,0x496EFD9A,0x3FBB80CA
118
	.long  0x3FFF0000,0x8B95C1E3,0xEA8BD6E7,0xBFBA8373
119
	.long  0x3FFF0000,0x8D1ADF5B,0x7E5BA9E6,0xBFBE9670
120
	.long  0x3FFF0000,0x8EA4398B,0x45CD53C0,0x3FBDB700
121
	.long  0x3FFF0000,0x9031DC43,0x1466B1DC,0x3FBEEEB0
122
	.long  0x3FFF0000,0x91C3D373,0xAB11C336,0x3FBBFD6D
123
	.long  0x3FFF0000,0x935A2B2F,0x13E6E92C,0xBFBDB319
124
	.long  0x3FFF0000,0x94F4EFA8,0xFEF70961,0x3FBDBA2B
125
	.long  0x3FFF0000,0x96942D37,0x20185A00,0x3FBE91D5
126
	.long  0x3FFF0000,0x9837F051,0x8DB8A96F,0x3FBE8D5A
127
	.long  0x3FFF0000,0x99E04593,0x20B7FA65,0xBFBCDE7B
128
	.long  0x3FFF0000,0x9B8D39B9,0xD54E5539,0xBFBEBAAF
129
	.long  0x3FFF0000,0x9D3ED9A7,0x2CFFB751,0xBFBD86DA
130
	.long  0x3FFF0000,0x9EF53260,0x91A111AE,0xBFBEBEDD
131
	.long  0x3FFF0000,0xA0B0510F,0xB9714FC2,0x3FBCC96E
132
	.long  0x3FFF0000,0xA2704303,0x0C496819,0xBFBEC90B
133
	.long  0x3FFF0000,0xA43515AE,0x09E6809E,0x3FBBD1DB
134
	.long  0x3FFF0000,0xA5FED6A9,0xB15138EA,0x3FBCE5EB
135
	.long  0x3FFF0000,0xA7CD93B4,0xE965356A,0xBFBEC274
136
	.long  0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x3FBEA83C
137
	.long  0x3FFF0000,0xAB7A39B5,0xA93ED337,0x3FBECB00
138
	.long  0x3FFF0000,0xAD583EEA,0x42A14AC6,0x3FBE9301
139
	.long  0x3FFF0000,0xAF3B78AD,0x690A4375,0xBFBD8367
140
	.long  0x3FFF0000,0xB123F581,0xD2AC2590,0xBFBEF05F
141
	.long  0x3FFF0000,0xB311C412,0xA9112489,0x3FBDFB3C
142
	.long  0x3FFF0000,0xB504F333,0xF9DE6484,0x3FBEB2FB
143
	.long  0x3FFF0000,0xB6FD91E3,0x28D17791,0x3FBAE2CB
144
	.long  0x3FFF0000,0xB8FBAF47,0x62FB9EE9,0x3FBCDC3C
145
	.long  0x3FFF0000,0xBAFF5AB2,0x133E45FB,0x3FBEE9AA
146
	.long  0x3FFF0000,0xBD08A39F,0x580C36BF,0xBFBEAEFD
147
	.long  0x3FFF0000,0xBF1799B6,0x7A731083,0xBFBCBF51
148
	.long  0x3FFF0000,0xC12C4CCA,0x66709456,0x3FBEF88A
149
	.long  0x3FFF0000,0xC346CCDA,0x24976407,0x3FBD83B2
150
	.long  0x3FFF0000,0xC5672A11,0x5506DADD,0x3FBDF8AB
151
	.long  0x3FFF0000,0xC78D74C8,0xABB9B15D,0xBFBDFB17
152
	.long  0x3FFF0000,0xC9B9BD86,0x6E2F27A3,0xBFBEFE3C
153
	.long  0x3FFF0000,0xCBEC14FE,0xF2727C5D,0xBFBBB6F8
154
	.long  0x3FFF0000,0xCE248C15,0x1F8480E4,0xBFBCEE53
155
	.long  0x3FFF0000,0xD06333DA,0xEF2B2595,0xBFBDA4AE
156
	.long  0x3FFF0000,0xD2A81D91,0xF12AE45A,0x3FBC9124
157
	.long  0x3FFF0000,0xD4F35AAB,0xCFEDFA1F,0x3FBEB243
158
	.long  0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x3FBDE69A
159
	.long  0x3FFF0000,0xD99D15C2,0x78AFD7B6,0xBFB8BC61
160
	.long  0x3FFF0000,0xDBFBB797,0xDAF23755,0x3FBDF610
161
	.long  0x3FFF0000,0xDE60F482,0x5E0E9124,0xBFBD8BE1
162
	.long  0x3FFF0000,0xE0CCDEEC,0x2A94E111,0x3FBACB12
163
	.long  0x3FFF0000,0xE33F8972,0xBE8A5A51,0x3FBB9BFE
164
	.long  0x3FFF0000,0xE5B906E7,0x7C8348A8,0x3FBCF2F4
165
	.long  0x3FFF0000,0xE8396A50,0x3C4BDC68,0x3FBEF22F
166
	.long  0x3FFF0000,0xEAC0C6E7,0xDD24392F,0xBFBDBF4A
167
	.long  0x3FFF0000,0xED4F301E,0xD9942B84,0x3FBEC01A
168
	.long  0x3FFF0000,0xEFE4B99B,0xDCDAF5CB,0x3FBE8CAC
169
	.long  0x3FFF0000,0xF281773C,0x59FFB13A,0xBFBCBB3F
170
	.long  0x3FFF0000,0xF5257D15,0x2486CC2C,0x3FBEF73A
171
	.long  0x3FFF0000,0xF7D0DF73,0x0AD13BB9,0xBFB8B795
172
	.long  0x3FFF0000,0xFA83B2DB,0x722A033A,0x3FBEF84B
173
	.long  0x3FFF0000,0xFD3E0C0C,0xF486C175,0xBFBEF581
174

175
	.set	N,L_SCR1
176

177
	.set	X,FP_SCR1
178
	.set	XDCARE,X+2
179
	.set	XFRAC,X+4
180

181
	.set	ADJFACT,FP_SCR2
182

183
	.set	FACT1,FP_SCR3
184
	.set	FACT1HI,FACT1+4
185
	.set	FACT1LOW,FACT1+8
186

187
	.set	FACT2,FP_SCR4
188
	.set	FACT2HI,FACT2+4
189
	.set	FACT2LOW,FACT2+8
190

191
	| xref	t_unfl
192
	|xref	t_ovfl
193
	|xref	t_frcinx
194

195
	.global	stwotoxd
196
stwotoxd:
197
|--ENTRY POINT FOR 2**(X) FOR DENORMALIZED ARGUMENT
198

199
	fmovel		%d1,%fpcr		| ...set user's rounding mode/precision
200
	fmoves		#0x3F800000,%fp0  | ...RETURN 1 + X
201
	movel		(%a0),%d0
202
	orl		#0x00800001,%d0
203
	fadds		%d0,%fp0
204
	bra		t_frcinx
205

206
	.global	stwotox
207
stwotox:
208
|--ENTRY POINT FOR 2**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
209
	fmovemx	(%a0),%fp0-%fp0	| ...LOAD INPUT, do not set cc's
210

211
	movel		(%a0),%d0
212
	movew		4(%a0),%d0
213
	fmovex		%fp0,X(%a6)
214
	andil		#0x7FFFFFFF,%d0
215

216
	cmpil		#0x3FB98000,%d0		| ...|X| >= 2**(-70)?
217
	bges		TWOOK1
218
	bra		EXPBORS
219

220
TWOOK1:
221
	cmpil		#0x400D80C0,%d0		| ...|X| > 16480?
222
	bles		TWOMAIN
223
	bra		EXPBORS
224

225

226
TWOMAIN:
227
|--USUAL CASE, 2^(-70) <= |X| <= 16480
228

229
	fmovex		%fp0,%fp1
230
	fmuls		#0x42800000,%fp1  | ...64 * X
231

232
	fmovel		%fp1,N(%a6)		| ...N = ROUND-TO-INT(64 X)
233
	movel		%d2,-(%sp)
234
	lea		EXPTBL,%a1	| ...LOAD ADDRESS OF TABLE OF 2^(J/64)
235
	fmovel		N(%a6),%fp1		| ...N --> FLOATING FMT
236
	movel		N(%a6),%d0
237
	movel		%d0,%d2
238
	andil		#0x3F,%d0		| ...D0 IS J
239
	asll		#4,%d0		| ...DISPLACEMENT FOR 2^(J/64)
240
	addal		%d0,%a1		| ...ADDRESS FOR 2^(J/64)
241
	asrl		#6,%d2		| ...d2 IS L, N = 64L + J
242
	movel		%d2,%d0
243
	asrl		#1,%d0		| ...D0 IS M
244
	subl		%d0,%d2		| ...d2 IS M', N = 64(M+M') + J
245
	addil		#0x3FFF,%d2
246
	movew		%d2,ADJFACT(%a6)	| ...ADJFACT IS 2^(M')
247
	movel		(%sp)+,%d2
248
|--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
249
|--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
250
|--ADJFACT = 2^(M').
251
|--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
252

253
	fmuls		#0x3C800000,%fp1  | ...(1/64)*N
254
	movel		(%a1)+,FACT1(%a6)
255
	movel		(%a1)+,FACT1HI(%a6)
256
	movel		(%a1)+,FACT1LOW(%a6)
257
	movew		(%a1)+,FACT2(%a6)
258
	clrw		FACT2+2(%a6)
259

260
	fsubx		%fp1,%fp0		| ...X - (1/64)*INT(64 X)
261

262
	movew		(%a1)+,FACT2HI(%a6)
263
	clrw		FACT2HI+2(%a6)
264
	clrl		FACT2LOW(%a6)
265
	addw		%d0,FACT1(%a6)
266

267
	fmulx		LOG2,%fp0	| ...FP0 IS R
268
	addw		%d0,FACT2(%a6)
269

270
	bra		expr
271

272
EXPBORS:
273
|--FPCR, D0 SAVED
274
	cmpil		#0x3FFF8000,%d0
275
	bgts		EXPBIG
276

277
EXPSM:
278
|--|X| IS SMALL, RETURN 1 + X
279

280
	fmovel		%d1,%FPCR		|restore users exceptions
281
	fadds		#0x3F800000,%fp0  | ...RETURN 1 + X
282

283
	bra		t_frcinx
284

285
EXPBIG:
286
|--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
287
|--REGISTERS SAVE SO FAR ARE FPCR AND  D0
288
	movel		X(%a6),%d0
289
	cmpil		#0,%d0
290
	blts		EXPNEG
291

292
	bclrb		#7,(%a0)		|t_ovfl expects positive value
293
	bra		t_ovfl
294

295
EXPNEG:
296
	bclrb		#7,(%a0)		|t_unfl expects positive value
297
	bra		t_unfl
298

299
	.global	stentoxd
300
stentoxd:
301
|--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT
302

303
	fmovel		%d1,%fpcr		| ...set user's rounding mode/precision
304
	fmoves		#0x3F800000,%fp0  | ...RETURN 1 + X
305
	movel		(%a0),%d0
306
	orl		#0x00800001,%d0
307
	fadds		%d0,%fp0
308
	bra		t_frcinx
309

310
	.global	stentox
311
stentox:
312
|--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
313
	fmovemx	(%a0),%fp0-%fp0	| ...LOAD INPUT, do not set cc's
314

315
	movel		(%a0),%d0
316
	movew		4(%a0),%d0
317
	fmovex		%fp0,X(%a6)
318
	andil		#0x7FFFFFFF,%d0
319

320
	cmpil		#0x3FB98000,%d0		| ...|X| >= 2**(-70)?
321
	bges		TENOK1
322
	bra		EXPBORS
323

324
TENOK1:
325
	cmpil		#0x400B9B07,%d0		| ...|X| <= 16480*log2/log10 ?
326
	bles		TENMAIN
327
	bra		EXPBORS
328

329
TENMAIN:
330
|--USUAL CASE, 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
331

332
	fmovex		%fp0,%fp1
333
	fmuld		L2TEN64,%fp1	| ...X*64*LOG10/LOG2
334

335
	fmovel		%fp1,N(%a6)		| ...N=INT(X*64*LOG10/LOG2)
336
	movel		%d2,-(%sp)
337
	lea		EXPTBL,%a1	| ...LOAD ADDRESS OF TABLE OF 2^(J/64)
338
	fmovel		N(%a6),%fp1		| ...N --> FLOATING FMT
339
	movel		N(%a6),%d0
340
	movel		%d0,%d2
341
	andil		#0x3F,%d0		| ...D0 IS J
342
	asll		#4,%d0		| ...DISPLACEMENT FOR 2^(J/64)
343
	addal		%d0,%a1		| ...ADDRESS FOR 2^(J/64)
344
	asrl		#6,%d2		| ...d2 IS L, N = 64L + J
345
	movel		%d2,%d0
346
	asrl		#1,%d0		| ...D0 IS M
347
	subl		%d0,%d2		| ...d2 IS M', N = 64(M+M') + J
348
	addil		#0x3FFF,%d2
349
	movew		%d2,ADJFACT(%a6)	| ...ADJFACT IS 2^(M')
350
	movel		(%sp)+,%d2
351

352
|--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
353
|--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
354
|--ADJFACT = 2^(M').
355
|--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
356

357
	fmovex		%fp1,%fp2
358

359
	fmuld		L10TWO1,%fp1	| ...N*(LOG2/64LOG10)_LEAD
360
	movel		(%a1)+,FACT1(%a6)
361

362
	fmulx		L10TWO2,%fp2	| ...N*(LOG2/64LOG10)_TRAIL
363

364
	movel		(%a1)+,FACT1HI(%a6)
365
	movel		(%a1)+,FACT1LOW(%a6)
366
	fsubx		%fp1,%fp0		| ...X - N L_LEAD
367
	movew		(%a1)+,FACT2(%a6)
368

369
	fsubx		%fp2,%fp0		| ...X - N L_TRAIL
370

371
	clrw		FACT2+2(%a6)
372
	movew		(%a1)+,FACT2HI(%a6)
373
	clrw		FACT2HI+2(%a6)
374
	clrl		FACT2LOW(%a6)
375

376
	fmulx		LOG10,%fp0	| ...FP0 IS R
377

378
	addw		%d0,FACT1(%a6)
379
	addw		%d0,FACT2(%a6)
380

381
expr:
382
|--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
383
|--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
384
|--FP0 IS R. THE FOLLOWING CODE COMPUTES
385
|--	2**(M'+M) * 2**(J/64) * EXP(R)
386

387
	fmovex		%fp0,%fp1
388
	fmulx		%fp1,%fp1		| ...FP1 IS S = R*R
389

390
	fmoved		EXPA5,%fp2	| ...FP2 IS A5
391
	fmoved		EXPA4,%fp3	| ...FP3 IS A4
392

393
	fmulx		%fp1,%fp2		| ...FP2 IS S*A5
394
	fmulx		%fp1,%fp3		| ...FP3 IS S*A4
395

396
	faddd		EXPA3,%fp2	| ...FP2 IS A3+S*A5
397
	faddd		EXPA2,%fp3	| ...FP3 IS A2+S*A4
398

399
	fmulx		%fp1,%fp2		| ...FP2 IS S*(A3+S*A5)
400
	fmulx		%fp1,%fp3		| ...FP3 IS S*(A2+S*A4)
401

402
	faddd		EXPA1,%fp2	| ...FP2 IS A1+S*(A3+S*A5)
403
	fmulx		%fp0,%fp3		| ...FP3 IS R*S*(A2+S*A4)
404

405
	fmulx		%fp1,%fp2		| ...FP2 IS S*(A1+S*(A3+S*A5))
406
	faddx		%fp3,%fp0		| ...FP0 IS R+R*S*(A2+S*A4)
407

408
	faddx		%fp2,%fp0		| ...FP0 IS EXP(R) - 1
409

410

411
|--FINAL RECONSTRUCTION PROCESS
412
|--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1)  -  (1 OR 0)
413

414
	fmulx		FACT1(%a6),%fp0
415
	faddx		FACT2(%a6),%fp0
416
	faddx		FACT1(%a6),%fp0
417

418
	fmovel		%d1,%FPCR		|restore users exceptions
419
	clrw		ADJFACT+2(%a6)
420
	movel		#0x80000000,ADJFACT+4(%a6)
421
	clrl		ADJFACT+8(%a6)
422
	fmulx		ADJFACT(%a6),%fp0	| ...FINAL ADJUSTMENT
423

424
	bra		t_frcinx
425

426
	|end
427

428
Product

Resources

Company
