[Glibc libm] Codeverständniss

27/03/2009 - 10:07 von Markus Wichmann | Report spam
Hi all,

ich habe gerade mal etwas im Quelltext der libm geblàttert und fand dort
folgenden Code für den Sinus. D.h. der Code dort stand in AT&T-Syntax
und hatte ein paar Schwàchen. Ich habe ihn etwass umgeschrieben,
allerding ohne die Semantik zu veràndern. Und ich mag die Intel-Syntax
mehr.

;file sysdeps/i386/fpu/s_sin.S

;x87 Condition Codes map to certain x86 flags, after executing
;fnstsw ax
;sahf
;
;namely:
;C0 - Carry Flag
;C2 - Parity Flag
;C3 - Zero Flag
;C1 would become an unnamed flag and therefore to be testet by bt and so

global sin
section .text
;C style declaration:
;float sin(float x);
sin:
fld dword [esp+4]
fsin
fnstsw ax
sahf
jnp .oor
ret

align 4

.oor:
fldpi
fadd st0
fxch st1
.loop:
fprem1
fntsw ax
sahf
jnp .loop

fstp st1
fsin
ret


OK, den Teil vor dem align verstehe ich noch: Der Sinus wird berechnet.
Wenn der Wert außer Reichweite ist, dann geht es weiter mit .oor (out of
range). Was ich aber nicht verstehe ist das, was danach kommt. OK, bei
.loop sieht der x87-Stack so aus, dass oben x steht und darunter x + pi.
Aber was macht die Schleife? Die Funktionsweise von fprem1 ist mir
nàmlich trotz Doku schleierhaft. Wenn die Rundung des Quotienten nicht
gemacht würde, wàre das Ergebnis dieser Operation in jedem Falle 0, also
kann das echte Ergebnis nur in der Nàhe von 0 liegen. Aber wàre es nicht
sinnvoller dafür eine Art Modulo-Rechnung zu nehmen? Also: man addiere
oder subtrahiere, je nach dem, den Wert 2*pi zur Unbekannten, bis diese
in Reichweite für ein fsin ist.

Vor allem aber die Schleifenbedingung ist mir unklar: das ganze soll
solange wiederholt werden, wie C2 gesetzt ist. Nun steht aber zu fprem1
in der Doku, dass C2 nur dann gelöscht wird, wenn sich die Operanden um
einen Faktor von mehr als 2^64 unterscheiden. Da einer der Operanden pi
ist und bleibt, wird der andere solange reduziert, bis er betragsmàßig
kleiner als (pi/2^64) =~ 1,7e-19 ist? Dessen Sinus kann ich vorhersagen:
Fast 0.

Vielleicht habe ich auch die Schleife vermacht. Der Originalcode ist

.loop: fprem1
fnstsw ax
test eax, 0400h
jnz .loop

Natürlich in AT&T-Syntax und mit numerischen Labels, aber mit diesen
Opcodes.

Kann mir da mal bitte jemand auf die Sprünge helfen?

Danke,
Markus

GUI - ein Hintergrundbild und zwölf XTerms

vim -c "exec \"norm iwHFG#NABGURE#IVZ#UNPXRE\"|%s/#/ /g|norm g??g~~"
 

Lesen sie die antworten

#1 Sebastian Biallas
27/03/2009 - 13:35 | Warnen spam
Markus Wichmann wrote:
Aber wàre es nicht
sinnvoller dafür eine Art Modulo-Rechnung zu nehmen?



Genau das wird doch dort gemacht...

Einfach im Intel-Handbuch mal unter fprem1 nachschlagen:

FPREM1—Partial Remainder (Continued)
Like the FPREM instruction, the FPREM1 computes the remainder through
iterative subtrac-
tion, but can reduce the exponent of ST(0) by no more than 63 in one
execution of the instruc-
tion. If the instruction succeeds in producing a remainder that is less
than one half the modulus,
the operation is complete and the C2 flag in the FPU status word is
cleared. Otherwise, C2 is set,
and the result in ST(0) is called the partial remainder. The exponent of
the partial remainder
will be less than the exponent of the original dividend by at least 32.
Software can re-execute the
instruction (using the partial remainder in ST(0) as the dividend) until
C2 is cleared. (Note that
while executing such a remainder-computation loop, a higher-priority
interrupting routine that
needs the FPU can force a context switch in-between the instructions in
the loop.)
An important use of the FPREM1 instruction is to reduce the arguments of
periodic functions.
When reduction is complete, the instruction stores the three
least-significant bits of the quotient
in the C3, C1, and C0 flags of the FPU status word. This information is
important in argument
reduction for the tangent function (using a modulus of π/4), because it
locates the original angle
in the correct one of eight sectors of the unit circle.

Gruß,
Sebastian

Ähnliche fragen