Beefy Boxes and Bandwidth Generously Provided by pair Networks
We don't bite newbies here... much
 
PerlMonks  

[OT] Function to do x**y for |y|<1 w/o ** operator.

by swampyankee (Parson)
on Feb 19, 2013 at 18:55 UTC ( #1019644=perlquestion: print w/ replies, xml ) Need Help??
swampyankee has asked for the wisdom of the Perl Monks concerning the following question:

In another forum (one that's science, not software, oriented), there was a poster asking for programming advice. He's trying to port some code to a PLC, the development environment of which can't do exponentiation.

Alas, he needs to do non-integer exponents. So, my question for the Mathematical Monks is does anybody know of an exponentiation function? Because of the way exponens work, the only concern is fractional powers, i.e. those with absolute values less than one. I'm digging through my copy of Knuth and Numerical Recipes, but since real languages can do $x = $y ** $z; with some restrictions, I've never looked into that deeply into the numerical libraries.


Information about American English usage here and here. Floating point issues? Please read this before posting. — emc

Comment on [OT] Function to do x**y for |y|<1 w/o ** operator.
Re: [OT] Function to do x**y for |y|<1 w/o ** operator.
by BrowserUk (Pope) on Feb 19, 2013 at 19:06 UTC

    Does he have log and exp available?

    print exp( log( 100 ) *.5 );; 10 print exp( log( 1000 ) *.3333333333333333333333 );; 10

    Ie.

    sub pow{ exp( log( $_[0] ) * $_[1] ) };; printf "%2d : %20.17f\n", $_, pow( 10**$_, 1/$_ ) for 1 .. 20;; 1 : 10.00000000000000200 2 : 10.00000000000000200 3 : 9.99999999999999820 4 : 10.00000000000000200 5 : 10.00000000000000200 6 : 9.99999999999999820 7 : 9.99999999999999820 8 : 10.00000000000000200 9 : 9.99999999999999820 10 : 10.00000000000000200 11 : 10.00000000000000200 12 : 9.99999999999999820 13 : 10.00000000000000200 14 : 9.99999999999999820 15 : 9.99999999999999820 16 : 10.00000000000000200 17 : 10.00000000000000200 18 : 9.99999999999999820 19 : 10.00000000000000200 20 : 10.00000000000000200 printf "%2d : %20.17f\n", $_, pow( 12345**$_, 1/$_ ) for 1 .. 20;; 1 : 12345.00000000000500000 2 : 12345.00000000000500000 3 : 12345.00000000000500000 4 : 12345.00000000000500000 5 : 12345.00000000000500000 6 : 12345.00000000000500000 7 : 12345.00000000000500000 8 : 12345.00000000000500000 9 : 12345.00000000000500000 10 : 12345.00000000000500000 11 : 12345.00000000000500000 12 : 12345.00000000000500000 13 : 12345.00000000000500000 14 : 12345.00000000000500000 15 : 12344.99999999998400000 16 : 12345.00000000000500000 17 : 12344.99999999998400000 18 : 12345.00000000000500000 19 : 12344.99999999998400000 20 : 12345.00000000000500000

    With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
    Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
    "Science is about questioning the status quo. Questioning authority".
    In the absence of evidence, opinion is indistinguishable from prejudice.

      He has log, but not exp.


      Information about American English usage here and here. Floating point issues? Please read this before posting. — emc

        If you don't need fast or accurate:

        #! perl -slw use strict; sub myExp{ my $e = 1 + ( shift() / 10000000 ); my $r = $e; $r *= $e for 2 .. 10000000; $r; } printf "%30.20g : %30.20g\n", exp( $_ ), myExp( $_ ) for map log( 10**$_ ), -10 .. 10; __END__ C:\test>exp 9.9999999999999965e-011 : 9.9997349036884958e-011 1.0000000000000007e-009 : 9.9997852777366009e-010 9.999999999999982e-009 : 9.9998303391048947e-009 9.9999999999999943e-008 : 9.9998700988200641e-008 1.0000000000000004e-006 : 9.9999045679358269e-007 9.9999999999999974e-006 : 9.9999337241854084e-006 0.00010000000000000009 : 9.9999575897369183e-005 0.0010000000000000002 : 0.00099999761423453593 0.010000000000000005 : 0.0099999893930861859 0.10000000000000002 : 0.099999973530420908 1 : 1 10.000000000000002 : 9.9999973561653093 100.00000000000004 : 99.999893882330838 999.99999999999977 : 999.99761406302275 10000.000000000009 : 9999.9575911775646 100000.00000000001 : 99999.337176962814 999999.99999999953 : 999990.45646420098 10000000.000000006 : 9999870.1098560225 100000000.00000018 : 99998303.311974689 999999999.99999928 : 999978527.31280696 10000000000.000004 : 9999734913.5273266

        With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.

        Update:Changed constants for slightly greater accuracy for the sake of another couple of iterations.

        Here's a somewhat more accurate, and much faster version:

        #! perl -slw use strict; sub myExp{ my $e = 1 + ( shift() / 4294967296 ); $e *= $e for 1 .. 32; $e; } printf "%30.20g : %30.20g\n", exp( $_ ), myExp( $_ ) for map log( 10**$_ ), -10 .. 10; __END__ C:\test>exp 9.9999999999999965e-011 : 1.0000000285332257e-010 1.0000000000000007e-009 : 1.0000000829163549e-009 9.999999999999982e-009 : 1.0000002130576865e-008 9.9999999999999943e-008 : 9.9999976329917627e-008 1.0000000000000004e-006 : 9.9999969112813618e-007 9.9999999999999974e-006 : 9.9999977582604716e-006 0.00010000000000000009 : 9.999998681099902e-005 0.0010000000000000002 : 0.00099999984556405617 0.010000000000000005 : 0.0099999993405499293 0.10000000000000002 : 0.099999996702749588 1 : 1 10.000000000000002 : 9.9999952742988558 100.00000000000004 : 100.00000084944983 999.99999999999977 : 999.99959325720954 10000.000000000009 : 10000.000169889967 100000.00000000001 : 99999.966690555739 999999.99999999953 : 1000000.1401405634 10000000.000000006 : 9999996.1589527112 100000000.00000018 : 100000003.39779937 999999999.99999928 : 999999662.81900096 10000000000.000004 : 10000002874.419859

        With the rise and rise of 'Social' network sites: 'Computers are making people easier to use everyday'
        Examine what is said, not who speaks -- Silence betokens consent -- Love the truth but pardon error.
        "Science is about questioning the status quo. Questioning authority".
        In the absence of evidence, opinion is indistinguishable from prejudice.
        He has log, but not exp.

        Perhaps there's also something helpful in Math::BigFloat's _pow() sub:
        # Taylor: | u u^2 u^3 | # x ** y = 1 + | --- + --- + ----- + ... | # |_ 1 1*2 1*2*3 _|
        where (I gather)
        u = y * ln x
        Cheers,
        Rob
          He has log, but not exp.

        Seriously. What kind of brain-dead platform has one but not the other. That's truly ridiculous.

        Alex / talexb / Toronto

        "Groklaw is the open-source mentality applied to legal research" ~ Linus Torvalds

Re: [OT] Function to do x**y for |y|<1 w/o ** operator.
by moritz (Cardinal) on Feb 19, 2013 at 21:43 UTC
Re: [OT] Function to do x**y for |y|<1 w/o ** operator.
by salva (Abbot) on Feb 20, 2013 at 09:18 UTC
    If you have log available, you can use Newton method to easyly find the solution to the equation f(z) = log(z) - y * log x = 0 :
    #!/usr/bin/perl use strict; use warnings; my ($X, $Y, $n) = @ARGV; $n ||= 20; my $Z = $X ** $Y; # f(z) = log z - y * log x = 0; # f'(z) = f1(z) = 1/z my $ylogx = $Y * log $X; my $z = 1; for my $i (1..$n) { my $fz = log($z) - $ylogx; my $f1z = 1 / $z; $z = $z - $fz / $f1z; my $e = abs($z - $Z); printf "%3d: z=%8f, e=%8f\n", $i, $z, $e; }
    This method does not handle correctly values of y < 0 but this is not a problem because z(y) = 1 / z(-y).

    Also, note that even if the algorithm converges in very few iterations, it uses log that is usually an expensive operation and so, other methods requiring more iterations may be actually faster.

Re: [OT] Function to do x**y for |y|<1 w/o ** operator.
by salva (Abbot) on Feb 20, 2013 at 11:37 UTC
    Another solution based on the fact that a floating point number is actually a sum of powers of 2.

    For 0 <= y <= 0, it holds that y = b0*20 + b1*2-1 + b2*2-2 + b3*2-3 + ...

    And so, xy = xb0*20 + b1*2-1 + b2*2-2 + b3*2-3 + ... = xb0*20 * xb1*2-1 * xb2*2-2 * xb3*2-3 * ... = Πi, bi≠0 x2-i

    Which in Perl becomes...

    #!/usr/bin/perl use strict; use warnings; my ($X, $Y, $n) = @ARGV; $n ||= 20; my $Z = $X ** $Y; my ($neg, $y) = ($Y < 0 ? (1, -$Y) : (0, $Y)); $y > 1 and die "Y is out of range"; my $bit = 1; my $x = $X; my $z = 1; while ($y) { if ($y >= $bit) { $y -= $bit; $z *= $x; } $bit /= 2; $x = sqrt($x); } $z = 1/$z if $neg; my $e = abs($z - $X ** $Y); print "z=$z, e=$e\n";

    sqrt is also an expensive operation. If you have to calculate xy for different values of y while x stays constant, you may be able to speed up the process creating a table with the values of x2-i.

    Alternatively, you can write a function to calculate ex (using a table with the values of e2i), and then calculate xy as elog(x)*y.

    update: there are also several implementations of exp(x) freely available, for instance, the one in OpenBSD is here.

Log In?
Username:
Password:

What's my password?
Create A New User
Node Status?
node history
Node Type: perlquestion [id://1019644]
Approved by BrowserUk
help
Chatterbox?
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others cooling their heels in the Monastery: (8)
As of 2014-12-26 04:02 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    Is guessing a good strategy for surviving in the IT business?





    Results (165 votes), past polls