 # Performing arithmetic operations in quadruple precision

#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 `Expr`s 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

``````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
Change default precision
ANN: ChangePrecision.jl