# How can I get Mathematica to produce better Fortran code?

I am new to Mathematica and I am trying to get to it to produce some expressions in Fortran code. However, it seems that Mathematica will output duplicated expressions, that is, expressions the need not be calculated more than once. Is there any good method to define a rule so that Mathematica will recognize duplicated expressions, assign them to a new variable, and then optimize the Fortran code?

P.S:

The following is some code I wrote to test Mathematica's Fortran capabilities. The solution of x, y and z contains many duplicated expressions, for example, (-36*a + 20*a**3 - 216*b)**2.

sol = Solve[{x^2 + y^2 + z^2 == 1, x + y + z == a , x*y*z == b}, {x, y, z}];
xx = x /. sol
yy = y /. sol
zz = z /. sol
Print["Writing Fortran Code . . . : / "];
SetDirectory["F:\\tang\\mathtest"];
strm = OpenWrite["test.f90", FormatType -> FotranForm,  PageWidth -> 70];
(* write subroutine of invisopar*)
WriteString[strm, "subroutine test(x,y,z,a,b)\n"];
WriteString[strm, "implicit none\n"];
WriteString[strm, "real*8::x,y,z,a,b\n"];
nroot = Length[xx];
For[ii = 1, ii <= nroot,
WriteString[strm,
"x = " <> ToString[FortranForm[xx[[ii]]]] <> "\n"]; ii++];

For[ii = 1, ii <= nroot,
WriteString[strm,
"y = " <> ToString[FortranForm[yy[[ii]]]] <> "\n"]; ii++];

For[ii = 1, ii <= nroot,
WriteString[strm,
"z = " <> ToString[FortranForm[zz[[ii]]]] <> "\n"]; ii++];

WriteString[strm, "end subroutine\n"];
Close[strm];
Print["Finished Writing Fortran Code . . . : / "];


When the test.f90 is outputted, the first solution of x is:

x = a/3. + (-6 + 2*a**2)/(3.2*0.6666666666666666*(-36*a + 20*a*3 + Sqrt(4(-6 + 2*a*2)*3 + (-36*a + 20*a**3 - 216*b)**2) - 216*b)**0.3333333333333333) - (-36*a + 20*a*3 + Sqrt(4(-6 + 2*a*2)*3 + (-36*a + 20*a**3 - 216*b)**2) - 216*b)**0.3333333333333333/(6.2*0.3333333333333333)

I would like that the code is as follows:

tmp0 = 216*b
tmp1 = -36*a + 20*a**3
tmp2 = (tmp1 - tmp0 )**2
tmp3 = -6 + 2*a**2
tmp4 = 4*tmp3**3
tmp5 = Sqrt(tmp4 + tmp2 )
tmp6 = tmp1 + tmp5 - tmp0
tmp7 = tmp6**0.3333333333333333
tmp8 = 2**0.6666666666666666
x = a/3. + tmp3 /(3.*tmp8 *tmp7 ) - tmp7/(6.*tmp8)


I noticed that maple has the capability for optimizing Fortran; for example, the following maple code can output optimized Fortran code:

with(codegen, fortran)
A := array(1 .. 2, 1 .. 2, symmetric);
A[1, 1] := log(x); A[1, 2] := 1-log(x); A[2, 2] := 2-log(x);
print(A);
fortran(A, optimized, mode = double)


So I want to know whether Mathematica also can do similar optimization.

• I don't have time to write a full answer, but a couple of things that might help are HornerForm for polynomials and the OptimizeExpression code used by Compile. For the latter, see 1 and 2 and linked questions. Commented Jan 14, 2013 at 7:19
• Hi Simon, thank you very much for your comment. Could you please give me more detailed information about OptimizeExpression? I checked the help and only find the C code example. Thanks, Tang Laoya Commented Jan 14, 2013 at 7:32
• @user58211 I've written up a quick example. Also, can you edit your profile so that it displays your name/alias? People are often more willing to help you if you're more than just a number! Commented Jan 14, 2013 at 9:44
• Did you really write FormatType -> FotranForm in your code? That's not gonna work. Commented Jan 14, 2013 at 12:17

OK, here's a quick go using ExperimentalOptimizeExpression, which is used internally by Compile. I'm not claiming that this is a polished solution to your problem, but hopefully it can be used as a place to start. (Also, forgive any Fortran mistakes, I haven't looked at Fortran for a long time...)

xx = x /. Solve[{x^2 + y^2 + z^2 == 1, x + y + z == a, x*y*z == b}, {x, y, z}] // First;

Quiet[Remove["tmp*"]]  (* otherwise the numbers after 'tmp' keep going up *)
Module[{out},
out = ExperimentalOptimizeExpression[xx, "OptimizationSymbol" -> tmp] /.
HoldPattern[ExperimentalOptimizedExpression[Block[{tempVars__},
CompoundExpression[defs : Set[_, _] ..., expr_]]]] :> Block[{
Set = ToString[#1] <> " = " <> ToString[FortranForm[#2]] &}, List[defs, expr]];
{"subroutine test(x,a,b)", "implicit none", "real*8::x,a,b",
Most@out, "x = " <> ToString@FortranForm@N@Last@out,
"end subroutine"}
] // Flatten // ExportString[#, "Table"] &

subroutine test(x,a,b)
implicit none
real*8::x,a,b
tmp4 = a**2
tmp5 = 2*tmp4
tmp6 = -6 + tmp5
tmp7 = -36*a
tmp8 = a**3
tmp9 = 20*tmp8
tmp12 = -216*b
tmp10 = tmp6**3
tmp11 = 4*tmp10
tmp13 = tmp12 + tmp7 + tmp9
tmp14 = tmp13**2
tmp15 = tmp11 + tmp14
tmp16 = Sqrt(tmp15)
tmp17 = tmp12 + tmp16 + tmp7 + tmp9
x = 0.3333333333333333*a - 0.13228342099734997*tmp17**0.3333333333333333 + \
(0.20998684164914552*tmp6)/tmp17**0.3333333333333333
end subroutine


This can easily be modified to process and output more than just the first solution for x and can also be made to export to a file instead of a string. It's not as compact as hand code would make it, some of the tmp definitions, I would join together, but it's not bad. Also, note that the other solutions for x contain complex numbers and might have to be treated a little more carefully.

• Dear all, Since x have 6 roots, how to let it output all these roots to Fortran code, such as x(1)=... x(2)=... ... Thanks Commented Feb 8, 2014 at 7:16
• Furthermore, I'd like to output the Fortran code to a file, by the For loop, how to do? Commented Feb 8, 2014 at 10:54
• @Simon, how can we get 1.0d0 real*8 value in MMA??? thanks in advance! Commented Oct 29, 2018 at 22:11
• @Simon, your solution will also work with CForm. Is there any way to make these tmp** as Array for Fortran or for C`? Commented Aug 15, 2020 at 7:58

Why not use MathCode?

MathCode makes it possible to develop prototypes in the interactive Mathematica environment which can be automatically translated to fast production code in C++ or Fortran90 and, if necessary, linked to external applications.

• ……Because it's not free? Commented Feb 9, 2014 at 5:12
• @xzczd and Mathematica is free? Commented Jun 22, 2015 at 11:18
• Come on, not every one can spend another $2,495.00 without hesitation. BTW, the discount of MathCode for students is really stingy compared to Mathematica: Mathematica student version only sells for 140 dollars (50.00 dollars in China, at least when I bought it), while MathCode sells for$1,095.00 even in education store! Commented Jun 23, 2015 at 6:48
• if you are looking for free alternatives, then python code can be transformed into fortran using sympy or into c/c++ using many other open-source libs. Commented Jun 25, 2015 at 7:52