Ask Your Question
1

Support for arbitrarily large numbers?

asked 2018-04-17 17:36:02 +0100

gkanarek gravatar image

Hi, I'm trying to compute the Juggler Sequence for relatively large numbers. Standard python doesn't calculate these correctly; I've discovered that Mathematica does, but I'm more familiar with Python and someone suggested I try out Sage.

I've tried this a couple ways using SageCell, and haven't yet succeeded. Here's the simplest possible implementation:

def juggle(n):
    if n % 2:
        return floor(sqrt(n)**3)
    return floor(sqrt(n))

When I run this on 37, it works perfectly:

x = 37
print x,

while x > 1:
    x = juggle(x)
    print u"\u2192 {}".format(x),

produces

37 → 225 → 3375 → 196069 → 86818724 → 9317 → 899319 → 852846071 → 24906114455136 → 4990602 → 2233 → 105519 → 34276462 → 5854 → 76 → 8 → 2 → 1

which standard python fails (it finds 3374 instead of 3375). However, if I try a hard one, say, x = 48443, it times out after a few iterations:

48443 → 10662193 → 34815273349 → 6496130099313865 → 523578821252958052233532 → 723587466207 → 615512041010804067 → 482897358660562651148793788 → 21974925680433 → 103012783516625098121 → 1045530445028727953685811220915

I can improve this by using logarithms:

def juggle(n):
    if n % 2:
        return floor(exp(1.5*log(n)))
    return floor(sqrt(n))

This implementation is much faster and gets much farther; however, it ends with an error:

TypeError: ECL says: #<a FLOATING-POINT-OVERFLOW>

I'd like to attach the full output as a text file, because it involves some very large numbers, but I don't have the karma yet. The number which causes the code to crash is about 3300 characters long. Here's the traceback, at least:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-a9218cbdcf1b> in <module>()
      8 
      9 while x > Integer(1):
---> 10     x = juggle(x)
     11     print u"\u2192 {}".format(x),

<ipython-input-1-a9218cbdcf1b> in juggle(n)
      1 def juggle(n):
      2     if n % Integer(2):
----> 3         return floor(exp(RealNumber('1.5')*log(n)))
      4     return floor(sqrt(n))
      5 

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/functions/other.pyc in __call__(self, x, **kwds)
    412 
    413         try:
--> 414             return floor(SR(x).full_simplify().canonicalize_radical())
    415         except ValueError:
    416             pass

/home/sc_serv/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.simplify_full (build/cythonized/sage/symbolic/expression.cpp:53709)()
   9698         x = self
   9699         x = x.simplify_factorial()
-> 9700         x = x.simplify_rectform()
   9701         x = x.simplify_trig()
   9702         x = x.simplify_rational()

/home/sc_serv/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.simplify_rectform (build/cythonized/sage/symbolic/expression.cpp:55646)()
   9846 
   9847         """
-> 9848         simplified_expr = self.rectform()
   9849 
   9850         if complexity_measure is None:

/home/sc_serv/sage/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.rectform (build/cythonized/sage/symbolic/expression.cpp:53186)()
   9530 
   9531         """
-> 9532         return self.maxima_methods().rectform()
   9533 
   9534     def unhold(self, exclude=None):

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/symbolic/maxima_wrapper.pyc in __call__(self, *args, **kwds)
     30         """
     31         return super(MaximaFunctionElementWrapper, self).__call__(*args,
---> 32                 **kwds).sage()
     33 
     34 class MaximaWrapper(SageObject):

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, *args, **kwds)
    655 
    656     def __call__(self, *args, **kwds):
--> 657         return self._obj.parent().function_call(self._name, [self._obj] + list(args), kwds)
    658 
    659     def help(self):

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in function_call(self, function, args, kwds)
    576                                        [s.name() for s in args],
    577                                        ['%s=%s'%(key,value.name()) for key, value in kwds.items()])
--> 578         return self.new(s)
    579 
    580     def _function_call_string(self, function, args, kwds):

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in new(self, code)
    345 
    346     def new(self, code):
--> 347         return self(code)
    348 
    349     ###################################################################

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)
    280 
    281         if isinstance(x, string_types):
--> 282             return cls(self, x, name=name)
    283         try:
    284             return self._coerce_from_special_method(x)

/home/sc_serv/sage/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __init__(self, parent, value, is_name, name)
    695                 self._name = parent._create(value, name=name)
    696             except (TypeError, RuntimeError, ValueError) as x:
--> 697                 raise TypeError(x)
    698 
    699     def _latex_(self):

TypeError: ECL says: #<a FLOATING-POINT-OVERFLOW>

Does anyone have any suggestions for me? Thanks in advance!

edit retag flag offensive close merge delete

Comments

1

Welcome to Ask Sage! Thank you for your question!

slelievre gravatar imageslelievre ( 2018-04-17 18:09:26 +0100 )edit
1

The problem is that when x is an integer, sqrt(x) becomes an element in Sage's SymbolicRing.

Hint: always avoid Sage's SymbolicRing if you can!

Stick to the integers!

slelievre gravatar imageslelievre ( 2018-04-17 18:11:55 +0100 )edit

2 Answers

Sort by » oldest newest most voted
2

answered 2018-04-17 18:08:37 +0100

slelievre gravatar image

Use the .isqrt() method to get the floor of the square root of an integer, as an integer.

That seems to work well for the examples suggested in the question.

Define juggle and print_juggle_sequence as follows:

def juggle(n):
    if n % 2:
        return (n^3).isqrt()
    return n.isqrt()

def print_juggle_sequence(x):
    print x 
    while x > 1:
        x = juggle(x)
        print x

Then try it on 37:

sage: print_juggle_sequence(37)
37
225
3375
196069
86818724
9317
899319
852846071
24906114455136
4990602
2233
105519
34276462
5854
76
8
2
1

or on 48443:

sage: print_juggle_sequence(48443)
48443
10662193
34815273349
6496130099313865
523578821252958052233532
723587466207
615512041010804067
482897358660562651148793788
21974925680433
103012783516625098121
1045530445028727953685811220915
...
18392638965492294541243350929469903670518362838411978378223645946128324146
4288664007064705289494967442314405131
8881432808181494071360574181865431001179156074769485216
2980173284925139058877740732
54590963399862
7388569
20083550847
2846169433517467
151841723749529535094491
59167972552621245915902054390059896
243244676308899401
119967888751182910434786797
1314006531652570971499311449697128320012
36249228014573923771
218246929826088259128572354715
101958212868994286046069755788632123549931710
10097435955181606698838
100485998801
31853585940530752
178475729
2384343086560
1544131
1918784550
43803
9167602
3027
166540
408
20
4
2
1
edit flag offensive delete link more

Comments

Fantastic, thank you! That works great. I find the Sage documentation a little hard to navigate, and I hadn't encountered isqrt(). Perfect!

gkanarek gravatar imagegkanarek ( 2018-04-17 18:46:13 +0100 )edit

Once you understand that there are usually more methods than functions, create an object of a given type and explore the available methods with the tab key.

For instance:

sage: a = 5
sage: a.<TAB>

(where <TAB> means pressing the tab key) lets you discover all the methods that can be applied to a Sage integer.

slelievre gravatar imageslelievre ( 2018-04-18 00:51:26 +0100 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-04-17 17:36:02 +0100

Seen: 1,310 times

Last updated: Apr 17 '18