StrE(FPnumber, sig_digits)

[akj]

I sometimes have a need to output numbers in scientific format such as 3.1415926535E-24 which StrD() and StrF() cannot appear to do. The code below fills this gap:

Code: Select all

Procedure.s StrE(fp.d, length=16)
; Convert FP number to fixed length scientific format
; The length is the number of significant mantissa digits

; The fundamental concept for a given floating point number fp:
; fp = sign*mantissa[radix 2]* 2^exponent[radix 2]
;    = sign*mantissa[radix10]*10^exponent[radix10]
; Where the first format is how the FP number is stored internally
;   and the second format is how it is shown as a string
; The purpose of StrE(fp) is to convert between these formats

; The hard bit is accurately calculating 10^exponent[radix10]
; This is done via the assembler instruction F2XM1 which
;   calculates 2^(x)-1 where x = -1.0 .. +1.0
;   and the relationship: 10^x = 2^(x*log2(10))

Protected mnt.d, exp ; mantissa, unbiased exponent [radix 10]
Protected fp$ ; The result
Protected test.d ; Test whether F2XM1 argument is in range

If length<0 Or length>20: length = 20: EndIf ; Sanity check
; Special case
If fp=0.0
  ProcedureReturn "+"+StrD(0.0, length-1)+" E+0"
EndIf
; Extract sign
If fp<0.0: fp = -fp: fp$ = "-": Else: fp$ = "+": EndIf
; Calculate radix 10 mantissa and exponent
EnableASM
  FNINIT ; Round to nearest, all exceptions masked, 64-bit precision
  FLD qword [p.v_fp]
  FXTRACT ; Exponent (radix 2) -> ST1, mantissa -> ST0
  FXCH ; Exponent (radix 2) -> ST0, mantissa -> ST1
  FLD ST0 ; Push ST0 (duplicate exp2)
  FLDLG2 ; Push Log10(2) = 0.301
  FMULP ; Approx exponent (radix 10) -> ST0
  FRNDINT ; Round exponent to integer
          ;   Must round (not truncate) for F2XM1 to work
  FIST exp ; Approx exponent (radix 10), may be wrong by 1
  ; Calculate mantissa correction = 2^(exp2-exp*log2(10))
  FLDL2T ; Push Log2(10) = 3.322
  FMULP ; exp2*log2(10)
  FSUBP ; exp2-exp*log2(10)  Hopefully in the range -1.0..+1.0
  FST test ; Magnitude upper limit is 0.5*Log2(10)=1.661
  If Abs(test)<1 ; If already in range -1.0..+1.0
    F2XM1 ; (2^x)-1
    FLD1 ; +1
    FADDP ; Mantissa correction
  Else ; Use 2^x = Square_of(2^(x/2))
    FLD1 ; 1
    FLD ST0 ; Duplicate
    FADDP ; 2
    FDIVP ; x/2
    F2XM1 ; (2^(x/2))-1
    FLD1 ; +1
    FADDP ; [Square root of mantissa correction]
    FLD ST0 ; Duplicate
    FMULP ; Square to give mantissa correction
  EndIf
  FMULP ; Mantissa
  FSTP mnt
DisableASM
; Normalise mantissa to be >=1 and <10 (i.e. 1 digit before DP)
; In theory at most one adjustment will ever be needed
While mnt<1.0: mnt * 10.0: exp - 1: Wend
While mnt>=10.0: mnt / 10.0: exp + 1: Wend
; Build result
fp$ + StrD(mnt, length-1)+" E"
If exp>=0: fp$ + "+": EndIf
ProcedureReturn fp$ + Str(exp)
EndProcedure

; Examples
Debug StrE(1/81, 9) ; +1.2345...E-2
Debug StrE(Pow(2, -52)) ; +2.2204...E-16
Debug StrE(-#PI, 20) ; -3.1415...E+0

I suspect that the FNINIT instruction should be replaced by something less drastic and that the If..Else..Endif block within the assembler code should be replaced by conditional jumps in pure assembler, but I don't know what would be best. Can anyone advise me?

[Psychophanta]

Sometimes ASM is not a worth for certain things, and this could be one.
Not sure about your doubts, but i used this

Code: Select all

Procedure$ StrE(fp.d,dec.b=6)
  Protected E.b
  If dec<0:dec=0:ElseIf dec>20:dec=20:EndIf;<- Sanity check
  While Abs(fp)>1E1
    fp*1E-1:E+1
  Wend
  While Abs(fp)<1E-1
    fp*1E1:E-1
  Wend
  If E:ProcedureReturn StrD(fp,dec)+"E"+Str(E)
  Else:ProcedureReturn StrD(fp,dec)
  EndIf
EndProcedure

And this one is probably better and with results more similar to your ones (i mean no 0. ....), and no loops. Is the one i use currently:

Code: Select all

Procedure$ StrE(fp.d,dec.b=6)
  Protected E.b
  If dec<0:dec=0:ElseIf dec>20:dec=20:EndIf;<- Sanity check
  E=Log10(fp)
  If E
    fp*Pow(10,-E)
    ProcedureReturn StrD(fp,dec)+"E"+Str(E)
  EndIf
  ProcedureReturn StrD(fp,dec)
EndProcedure

[dioxin]

akj,

I suspect that the FNINIT instruction should be replaced by something less drastic

You could save and restore the full FPU state so it's not messed up by your code but this is a bit slow:

Code: Select all

'save the full FPU state
!sub esp,108     'make room for 108 bytes of FPU status
!fsave [esp]

..do what FPU stuff you like here..

'restore the full FPU state to what it was before you started
!frstor [esp]         'restore FPU state
!add esp,108          'clear up stack

'or you could just set the control word to what you need and then restore it when finished:

Code: Select all

!push &h1f7f0000       'the FP control word needed is in the first 4 digits
!fstcw [esp]           'save the current FP control word
!fldcw [esp+2]         'set the FP control word to what you need
..
your other ASM code goes here
..
!fldcw [esp]           'restore the original control word 
!add esp,4             'clean up CPU stack

For your IF THEN ELSE ENDIF block you could use something like this:

Code: Select all

'if abs(test)<1
!fabs    'test is on the top of the stack. Get its ABS value
!fld1    'load the number 1
!fcomip  'compare ST0 and ST1 and pop the stack (remove the 1). Set CPU flags on the result
!jna short label1  'skip the IF code block
'your IF code block goes here
..
..

!jmp short label2  'skip the ELSE code block
label1:
'your ELSE code block goes here
..
..

label2: 

[akj]

@Psychophanta :

I like your algorithms as they are very straightforward.

However, I wonder whether they could result in loss of accuracy from either the repeated multiplications or the Pow() function, the latter having been reported as having problems in this area. For example:
http://www.purebasic.fr/english/viewtopic.php?t=34238

In my assembler method, I was doing my best to retain full accuracy.
I think I will spend some time comparing the results returned by the three routines to see whether any errors are creeping in.


@dioxin:

Thank you for your feedback which is most welcome.

I am inexperienced in using assembler under Windows, although I have programmed in Z80 assembler (Amstrad CPC) in the past.

I shall certainly implement your suggestions regarding the FPU control word and the conditional jumps. I have always struggled with the latter, in trying to decide which of the jump instructions best fits what I'm trying to achieve. It would have taken me ages to realise that JNA was what I needed, so thanks for that.

Is there any documentation that describes PureBasic's implementation of assembler? For example, where does it document:
* Having to use ST0 rather then ST(0)
* The meanings of dword, tword, qword, etc
* When it is necessary to use p.v_Variable rather than the simple variable name.
Does such documentation exist in the public domain?

[Trond]

Is there any documentation that describes PureBasic's implementation of assembler? For example, where does it document:
* Having to use ST0 rather then ST(0)
* The meanings of dword, tword, qword, etc
* When it is necessary to use p.v_Variable rather than the simple variable name.
Does such documentation exist in the public domain?

Have you read the help file? Press F1 and look up asm in the index.

[akj]

@Trond:

I'm sorry, I see what you mean. I now realise I need to read the Flat assembler documentation, in particular http://flatassembler.net/docs.php?article=manual