; IEEE Double precision in Blitz2:  written by Andrea Galimberti  26-8-96
;
; You can use this code freely in your source as long as you do not use it
; for commercial purposes, otherwise you have to ask me a permission (my
; address is: via Villoresi, 87   20029 Turbigo (Mi)   Italy).
; You have to write in your work that:
;
;   IEEEDouble is Copyright 1996 by Andrea Galimberti
;
; This work has been modeled on the longreal module for the E language
; written by EA van Breemen (thank you). This is why the double type is
; called "longreal".

; The longreal (double) type is:

NEWTYPE.longreal
  a.l
  b.l
End NEWTYPE

; The following returns true if c$ is a digit

Function isDigit{c$}
  If (Asc(c$)>=48) AND (Asc(c$)<=57)
    Function Return True
  EndIf
  Function Return False
End Function

; This MUST be called to initialize a longreal variable before using it.
; Usage: DEFTYPE.longreal *foo : *foo=dNew{}
; At the end of the program you have to free all variables you initialized,
; e.g., dFree{*foo}.
; Note that you MUST define a longreal variable as a pointer.

Function.l dNew{}
  t.l=AllocVec_(SizeOf.longreal,#MEMF_PUBLIC|#MEMF_CLEAR)
  Function Return t
End Function

; Frees longreal variables initialized with dNew{}.

Statement dFree{*x.longreal}
  FreeVec_ *x
End Statement

; Converts an integer (long) to a double: i  ->  lr

Statement dFloat{i.l,*lr.longreal}
  w.l=IEEEDPFlt_(i)
  PutReg d0,*lr\a
  PutReg d1,*lr\b
End Statement

; Converts a double to an integer (long)
; Returns the integer value

Function.l dFix{*lr.longreal}
  w.l=IEEEDPFix_(*lr\a,*lr\b)
  Function Return w
End Function

; Tests a double against the value 0.0
; Returns:  1  ->  x>0
;           0  ->  x=0
;          -1  ->  x<0

Function dTst{*x.longreal}
  flag.b=IEEEDPTst_(*x\a,*x\b)
  Function Return flag
End Function

; Compares two doubles
; Returns:  1  ->  x>y
;           0  ->  x=y
;          -1  ->  x<y

Function dCompare{*x.longreal,*y.longreal}
  flag.b=IEEEDPCmp_(*x\a,*x\b,*y\a,*y\b)
  Function Return flag
End Function

; Adds: x+y  ->  x

Statement dAdd{*x.longreal,*y.longreal}
  w.l=IEEEDPAdd_(*x\a,*x\b,*y\a,*y\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Subtracts: x-y  ->  x

Statement dSub{*x.longreal,*y.longreal}
  w.l=IEEEDPSub_(*x\a,*x\b,*y\a,*y\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Multiplies: x*y  ->  x

Statement dMul{*x.longreal,*y.longreal}
  w.l=IEEEDPMul_(*x\a,*x\b,*y\a,*y\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Division: x/y  ->  x

Statement dDiv{*x.longreal,*y.longreal}
  w.l=IEEEDPDiv_(*x\a,*x\b,*y\a,*y\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Nearest integer less than or equal to x  ->  x

Statement dRound{*x.longreal}
  w.l=IEEEDPFloor_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Nearest integer greater than or equal to x  ->  x

Statement dRoundUp{*x.longreal}
  w.l=IEEEDPCeil_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; -1*x  ->  x

Statement dNeg{*x.longreal}
  w.l=IEEEDPNeg_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Absolute value of x  ->  x

Statement dAbs{*x.longreal}
  w.l=IEEEDPAbs_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Copy y  ->  x

Statement dCopy{*x.longreal,*y.longreal}
  *x\a=*y\a
  *x\b=*y\b
End Statement

; Power: x^y  ->  x

Statement dPow{*x.longreal,*y.longreal}
  w.l=IEEEDPPow_(*y\a,*y\b,*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Converts IEEE Single to IEEE Double: x  ->  y

Statement dDouble{x.l,*y.longreal}
  w.l=IEEEDPFieee_(x)
  PutReg d0,*y\a
  PutReg d1,*y\b
End Statement

; Converts IEEE Double to IEEE Single
; Returns: single value

Function.l dSingle{*x.longreal}
  w.l=IEEEDPTieee_(*x\a,*x\b)
  Function Return w
End Function

; Square root of x  ->  x

Statement dSqrt{*x.longreal}
  w.l=IEEEDPSqrt_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Pi value  ->  x

Statement dPi{*x.longreal}
 *x\a=$400921FB             ;Dirty but quick 8-)
 *x\b=$54442D18
End Statement

; Converts degrees to radians: x  ->  x

Statement dRad{*x.longreal}
  DEFTYPE.longreal *s,*t

  *s=dNew{} : *t=dNew{}

  dPi{*t}
  dFloat{180,*s}

  dDiv{*t,*s}
  dMul{*t,*x}
  *x\a=*t\a
  *x\b=*t\b

  dFree{*s} : dFree{*t}
End Statement

; Sine of x  ->  x

Statement dSin{*x.longreal}
  w.l=IEEEDPSin_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Cosine of x  ->  x

Statement dCos{*x.longreal}
  w.l=IEEEDPCos_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Tangent of x  ->  x

Statement dTan{*x.longreal}
  w.l=IEEEDPTan_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; ArcSine of x  ->  x

Statement dASin{*x.longreal}
  w.l=IEEEDPAsin_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; ArcCosine of x  ->  x

Statement dACos{*x.longreal}
  w.l=IEEEDPAcos_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; ArcTan of x  ->  x

Statement dATan{*x.longreal}
  w.l=IEEEDPAtan_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Hyperbolic Sine of x  ->  x

Statement dSinh{*x.longreal}
  w.l=IEEEDPSinh_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Hyperbolic Cosine of x  ->  x

Statement dCosh{*x.longreal}
  w.l=IEEEDPCosh_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Hyperbolic Tangent of x  ->  x

Statement dTanh{*x.longreal}
  w.l=IEEEDPTanh_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Exponential: e^x  ->  x

Statement dExp{*x.longreal}
  w.l=IEEEDPExp_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Natural logaritm: ln x  ->  x

Statement dLn{*x.longreal}
  w.l=IEEEDPLog_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement

; Logaritm in base 10: log x  ->  x

Statement dLog{*x.longreal}
  w.l=IEEEDPLog10_(*x\a,*x\b)
  PutReg d0,*x\a
  PutReg d1,*x\b
End Statement


;*************************************************************
;* Converts an ascii representation to a longreal            *
;*                                                           *
;* Accepts the exponential notation: 1E23 or 23.34e8         *
;*                                                           *
;* Params:  s$  =  string with the ascii representation      *
;*          x   =  converted longreal                        *
;*************************************************************

Statement a2d{s$,*x.longreal}
  DEFTYPE.longreal *divider,*fraction,*ten,*tmp,*longexp
  DEFTYPE.l i,expo,expsign,sign

  If s$=""
    *x\a=0 : *x\b=0
    Statement Return
  EndIf

  *divider=dNew{} : *fraction=dNew{} : *ten=dNew{}
  *tmp=dNew{} : *longexp=dNew{}

  tmpbuffer$=""

  dFloat{0,*x}
  dFloat{10,*ten}
  i=1
  sign=1
  If Mid$(s$,i,1)="-"
    sign=-1
    i+1
  Else
    If Mid$(s$,i,1)="+" Then i+1
  EndIf

  While isDigit{Mid$(s$,i,1)}
    dFloat{Asc(Mid$(s$,i,1))-48,*tmp}
    dMul{*x,*ten}
    dAdd{*x,*tmp}
    i+1
  Wend


  If (Mid$(s$,i,1)=".") AND isDigit{Mid$(s$,i+1,1)}
    i+1
    dFloat{1,*divider}
    dFloat{0,*fraction}
    While isDigit{Mid$(s$,i,1)}
      dMul{*fraction,*ten}
      dFloat{Asc(Mid$(s$,i,1))-48,*tmp}
      dAdd{*fraction,*tmp}
      dMul{*divider,*ten}
      i+1
    Wend
    dDiv{*fraction,*divider}
    dAdd{*x,*fraction}
  EndIf
  dFloat{sign,*tmp}
  dMul{*x,*tmp}

  If (Mid$(s$,i,1)="E") OR (Mid$(s$,i,1)="e")
    i+1
    If Mid$(s$,i,1)="-"
      expsign=-1
      i+1
    Else
      expsign=1
      If Mid$(s$,i,1)="+" Then i+1
    EndIf
    expo=0
    While isDigit{Mid$(s$,i,1)}
      expo=expo*10+Asc(Mid$(s$,i,1))-48
      i+1
    Wend
    dFloat{expo*expsign,*longexp}
    dPow{*ten,*longexp}
    dMul{*x,*ten}
  EndIf

  dFree{*divider} : dFree{*fraction} : dFree{*ten}
  dFree{*tmp} : dFree{*longexp}
End Statement


;*******************************************************
;* Converts a longreal x to ascii with num digits      *
;* Only for fraction numbers (nonnegative)             *
;*                                                     *
;* Params:   x    =  longreal to be converted          *
;*           num  =  number of digits                  *
;* Returns:  string with the ascci representation      *
;*******************************************************

Function$ dFormat{*x.longreal,num.l}
  DEFTYPE.longreal *c,*d,*e
  pos.l=0

  *c=dNew{} : *d=dNew{} : *e=dNew{}

  v.l=dFix{*x}
  s$+Str$(v)+"."
  dCopy{*c,*x}
  For k.l=1 To num
    dCopy{*d,*c}
    dRound{*d}
    dSub{*c,*d}
    dFloat{10,*d}
    dMul{*c,*d}
    f.l=dFix{*c}+48
    s$+Chr$(f)
  Next k

  c.w=Len(s$)
  While Mid$(s$,c,1)="0"
    c-1
  Wend
  If Mid$(s$,c,1)="." Then c-1
  If c<Len(s$) Then s$=Left$(s$,c)

  dFree{*c} : dFree{*d} : dFree{*e}
  Function Return s$
End Function


;***********************************************************************
;* Converts a longreal x to ascii with num digits and using the        *
;* exponential notation if power>maxpow.                               *
;* Also for 'large' numbers                                            *
;*                                                                     *
;* Params:   x       =  longreal to be converted                       *
;*           num     =  number of digits                               *
;*           maxpow  =  power above which exponential notation is used *
;* Returns:  string with ascii representation                          *
;***********************************************************************

Function$ dLFormat{*x.longreal,num.l,maxpow.l}
  DEFTYPE.longreal *a,*one,*ten
  sign.l=1

  *a=dNew{} : *one=dNew{} : *ten=dNew{}

  dFloat{10,*ten}
  dFloat{1,*one}
  dCopy{*a,*x}
  power.l=0
  If dTst{*a}=0
    s$=dFormat{*a,num}
    dFree{*a} : dFree{*one} : dFree{*ten}
    Function Return s$
  EndIf
  If (dTst{*a}=-1)
    sign=-1
    dNeg{*a} : dNeg{*x}
  EndIf
  If dCompare{*a,*one}=-1
    While dCompare{*a,*one}=-1
      dMul{*a,*ten}
      power-1
    Wend
  Else
    While dCompare{*a,*ten}=1
      dDiv{*a,*ten}
      power+1
    Wend
  EndIf
  If Abs(power)<=maxpow
    dCopy{*a,*x}
  EndIf
  s$=dFormat{*a,num}
  If sign=-1
    s$="-"+s$
  EndIf
  If Abs(power)>maxpow
    s$+"E"+Str$(power)
  EndIf

  dFree{*a} : dFree{*one} : dFree{*ten}
  Function Return s$
End Function