Generate more reasonable floating point numbers

32 views
Skip to first unread message

Peter Brady

unread,
May 9, 2016, 11:26:25 AM5/9/16
to sympy
I have some expressions that are being written to a fortran file and sympy is keeping them in a form that involves coefficients that are too large.

For example, here's one of my coefficients:


(6447959512020*c21*r**10 + 89877391198101*c21*r**9 + 482379015493674*c21*r**8 + 1955782912118958*c21*r**7 + 6144750386733894*c21*r**6 + 9734633983924239*c21*r**5 + 3581719348056096*c21*r**4 - 7087617275638458*c21*r**3 - 11332310233710024*c21*r**2 - 6752358207277020*c21*r + 6016737364560*r**10 + 83866633598228*r**9 + 450118807506472*r**8 + 1562757816049524*r**7 + 3409298852122307*r**6 + 2460129818497007*r**5 - 3107226702425082*r**4 - 9732722117772259*r**3 + 12534948043759713*r**2 - 10702947100060800*r + 1022514295602150)/(3387470672547000*r**5 + 21246968829475350*r**4 + 34890947927234100*r**3 - 11310388190004150*r**2 - 62103628996695000*r - 33761791036385100)

The problem is if I write it to a file with a simple transformation of str(expr.n()) I get something troublesome like 

(6447959512020.0*c21*r**10 + 89877391198101.0*c21*r**9 + 482379015493674.0*c21*r**8 + 1.95578291211896e+15*c21*r**7 + 6.14475038673389e+15*c21*r**6 + 9.73463398392424e+15*c21*r**5 + 3.5817193480561e+15*c21*r**4 - 7.08761727563846e+15*c21*r**3 - 1.133231023371e+16*c21*r**2 - 6.75235820727702e+15*c21*r + 6016737364560.0*r**10 + 83866633598228.0*r**9 + 450118807506472.0*r**8 + 1.56275781604952e+15*r**7 + 3.40929885212231e+15*r**6 + 2.46012981849701e+15*r**5 - 3.10722670242508e+15*r**4 - 9.73272211777226e+15*r**3 + 1.25349480437597e+16*r**2 - 1.07029471000608e+16*r + 1.02251429560215e+15)/(3.387470672547e+15*r**5 + 2.12469688294753e+16*r**4 + 3.48909479272341e+16*r**3 - 1.13103881900042e+16*r**2 - 6.2103628996695e+16*r - 3.37617910363851e+16)

I also do some regex processing to add the appropriate '_dp' suffix to all the numbers.

the problem is that c21 and r are O(1) and they are being multiplied by huge numbers and then divided by huge numbers.  I've tried expand and several other options but I can't seem to figure out an automatic way to get sympy to give me more reasonable double precision coefficients.  I can do it "manually" by splitting the expression into the numerator and denominator and then taking the large constant and in the denominator and normalizing the numerator with it before gluing the expression back together but I'm hoping there's a better way.

Here's what I get if do the splitting manually:

(0.342625169147496617*c21*r**10 + 4.77581416328371673*c21*r**9 + 25.632169598556608*c21*r**8 + 103.924419784486855*c21*r**7 + 326.513548464309187*c21*r**6 + 517.269162300637267*c21*r**5 + 190.321790200281699*c21*r**4 - 376.614658232800258*c21*r**3 - 602.164871447902571*c21*r**2 - 358.8*c21*r + 0.319711321605773568*r**10 + 4.45642058838269934*r**9 + 23.9179592041297048*r**8 + 83.0402486340673783*r**7 + 181.159883790403722*r**6 + 130.723897012075825*r**5 - 165.108678569306389*r**4 - 517.167571485358666*r**3 + 666.069426419644693*r**2 - 568.722408026755853*r + 54.3333333333333333)/((r + 1.0)*(180.0*r**4 + 949.0*r**3 + 905.0*r**2 - 1506.0*r - 1794.0))

Is there a better way?


Aaron Meurer

unread,
May 9, 2016, 12:45:33 PM5/9/16
to sy...@googlegroups.com
You could also write a custom StrPrinter subclass that prints floats
how you want. May be more robust than regular expressions.

>
>
> the problem is that c21 and r are O(1) and they are being multiplied by huge
> numbers and then divided by huge numbers. I've tried expand and several
> other options but I can't seem to figure out an automatic way to get sympy
> to give me more reasonable double precision coefficients. I can do it
> "manually" by splitting the expression into the numerator and denominator
> and then taking the large constant and in the denominator and normalizing
> the numerator with it before gluing the expression back together but I'm
> hoping there's a better way.

That's how I would do it. Does dividing by the leading coefficient work?

In [11]: p, q = a.as_numer_denom()

In [12]: print(p/LC(q)/(q/LC(q)))
(0.00190347316193054*c21*r**10 + 0.0265323009071318*c21*r**9 +
0.142400942214203*c21*r**8 + 0.577357887691594*c21*r**7 +
1.81396415813505*c21*r**6 + 2.87371756833687*c21*r**5 +
1.05734327889046*c21*r**4 - 2.09230365684889*c21*r**3 -
3.34536039693279*c21*r**2 - 1.99333333333333*c21*r +
0.00177617400892096*r**10 + 0.0247578921576817*r**9 +
0.132877551134054*r**8 + 0.461334714633707*r**7 +
1.00644379883558*r**6 + 0.726243872289311*r**5 -
0.917270436496146*r**4 - 2.87315317491866*r**3 + 3.70038570233136*r**2
- 3.15956893348198*r + 0.301851851851852)/(1.0*r**5 +
6.27222222222221*r**4 + 10.3*r**3 - 3.3388888888889*r**2 -
18.3333333333333*r - 9.96666666666667)

Aaron Meurer

>
>
> Here's what I get if do the splitting manually:
>
>
> (0.342625169147496617*c21*r**10 + 4.77581416328371673*c21*r**9 +
> 25.632169598556608*c21*r**8 + 103.924419784486855*c21*r**7 +
> 326.513548464309187*c21*r**6 + 517.269162300637267*c21*r**5 +
> 190.321790200281699*c21*r**4 - 376.614658232800258*c21*r**3 -
> 602.164871447902571*c21*r**2 - 358.8*c21*r + 0.319711321605773568*r**10 +
> 4.45642058838269934*r**9 + 23.9179592041297048*r**8 +
> 83.0402486340673783*r**7 + 181.159883790403722*r**6 +
> 130.723897012075825*r**5 - 165.108678569306389*r**4 -
> 517.167571485358666*r**3 + 666.069426419644693*r**2 - 568.722408026755853*r
> + 54.3333333333333333)/((r + 1.0)*(180.0*r**4 + 949.0*r**3 + 905.0*r**2 -
> 1506.0*r - 1794.0))
>
>
> Is there a better way?
>
>
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/ee10ecaf-bd9f-4b21-98e5-68bdf2640a5b%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

Peter Brady

unread,
May 9, 2016, 1:10:49 PM5/9/16
to sy...@googlegroups.com
Thanks Aaron!  That's significantly cleaner than my implementation.

You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/YJoxzw7dU5c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.

To post to this group, send email to sy...@googlegroups.com.
Visit this group at https://groups.google.com/group/sympy.
Reply all
Reply to author
Forward
0 new messages