Performing arithmetic operations in quadruple precision

General Usage
1

I need to compute simple arithmetic expressions with +,-,*,/,√ in quadruple precision using DecFP where the various coefficients are hard-coded. In order to get an accurate result, all numbers have to be converted to DecFP128 before applying any operators. In the code that would look like

Dec128(5)/Dec128(36)-Dec128(7)*√Dec128(15)/Dec128(30)

which, quite obviously, gives a different result as

Dec128(5/36-7*√15/30)

Is there a more compact way of doing this?

If it weren’t for √ rationals would be a solution, but I need roots…

2

I’m completely unfamiliar with the package but I guess you could do a macro that replaces all integers that it founds with the correspondent DecFP128. For example

function _dec128(x::Expr)
    y = x
    y.args .= _dec128.(y.args)
    return  y
end

_dec128(x::Int) = Dec128(x)

_dec128(x) = x

macro dec128(x)
    return esc(_dec128(x))
end

Now you can just d:

julia> @dec128 5/36-7*√15/30
-7648072252261750509862730377270044E-34

julia> ans == Dec128(5)/Dec128(36)-Dec128(7)*√Dec128(15)/Dec128(30)
true
5 Likes
3

Thanks a lot for the quick reply! I was under the impression that a macro like that would be the solution. I just wasn’t sure how exactly that should look like…

5

Please don’t use convert for this.

6

Okay, i changed the name and got the macro working.
How would i get rid of the eval in the macro?

module ConvertArguments

export @convertarguments

function convertarguments(T, ex::Expr)
   for (i, arg) in enumerate(ex.args)
       if typeof(arg) <: Number
           ex.args[i] = Expr(:call, :convert, T, arg)
       elseif typeof(arg) == Expr
           ex.args[i] = convertarguments(T, arg)
       end
   end
   return :($ex)
end

macro convertarguments(T,ex)
    return convertarguments(T,ex)
end

end #module ConvertArguments

julia> include("ConvertArguments.jl"); using ConvertArguments

julia> @macroexpand @convertarguments BigInt 1+3/4.0*√6
:((ConvertArguments.convert)(ConvertArguments.BigInt, 1) + ((ConvertArguments.convert)(ConvertArguments.BigInt, 3) / (ConvertArguments.convert)(ConvertArguments.BigInt, 4.0)) * √((ConvertArguments.convert)(ConvertArguments.BigInt, 6)))

julia> @convertarguments BigInt 1+3/4.0*√6
2.837117307087383573647963056029418543974460610492502596324519425438220283092974
7

A macro should return an expression. That expression will be eval’ed when the macro is called.

8

So you should pass the symbol T through to the function, and build up an expression that converts things to the type given by the symbol T.

9

I edited my first post to reflect your suggestions. Would that by the ideomatic way to do it?

10

Could you please explain why convert is not a good idea here?

11

Well, it’s by no mean a convertion…

  1. It doesn’t return a type it is converting to
  2. It’ll cause assignment of Expr object to a typed slot to raise a strange error.
  3. It’s probably also a type piracy though that’s much less of a concern…
12

  • This does not handle escape correctly. You must escape all arguments ones and exactly ones.

  • I’m not sure what’s the best way to handle this but you can accidentally handle some Exprs that you don’t want to handle (line numbers for example…)

  • There’s also isa which you should use instead of typeof

  • And :($ex) is a no-op.

1 Like
13

Note that simple binary operators (+,-,*,/) should promote automatically, so in your case, you only need to explicitly use Dec128 for √.

Julia-0.5.1> x1 = Dec128(5)/Dec128(36)-Dec128(7)*√Dec128(15)/Dec128(30)
-7648072252261750509862730377270044E-34

Julia-0.5.1> x2 = Dec128(5)/36-7*√Dec128(15)/30
-7648072252261750509862730377270044E-34

Julia-0.5.1> x1 == x2
true

Also note that DoubleDouble also seems to support the operations you need and could be faster:

Julia-0.5.1> y1 = Double(5)/Double(36)-Double(7)*√Double(15)/Double(30)
Double(-0.7648072252261751, 7.110782587967322e-18)
 - value: -0.76480722522617505098627303772702

Julia-0.5.1> y2 = Double(5)/36-7*√Double(15)/30
Double(-0.7648072252261751, 7.110782587967322e-18)
 - value: -0.76480722522617505098627303772702

Julia-0.5.1> y1 == y2
true

Julia-0.5.1> @btime Dec128(5)/36-7*√Dec128(15)/30
  443.050 ns (0 allocations: 0 bytes)
-7648072252261750509862730377270044E-34

Julia-0.5.1> @btime Double(5)/36-7*√Double(15)/30
  84.492 ns (0 allocations: 0 bytes)
Double(-0.7648072252261751, 7.110782587967322e-18)
 - value: -0.76480722522617505098627303772702
1 Like
14

Also see this: Change default precision - #47 by stevengj

using DecFP
using ChangePrecision

y = Dec128(5)/Dec128(36)-Dec128(7)*√Dec128(15)/Dec128(30)

@changeprecision Dec128 begin
    z = 5/36 - 7*√15/30
end

y == z  # true

EDIT: or more concisely:

@changeprecision Dec128 z = 5/36-7*√15/30
2 Likes