mersenneforum.org  

Go Back   mersenneforum.org > Math Stuff > Computer Science & Computational Number Theory > PARI/GP

Reply
 
Thread Tools
Old 2020-04-25, 11:29   #34
User133
 
Apr 2020

101002 Posts
Default

Next step would be to make a function of it, so you can call it with starting and stopping values you need. Then, you can move to better "algorithms"
but I don’t understand it at all

Last fiddled with by User133 on 2020-04-25 at 11:35
User133 is offline   Reply With Quote
Old 2020-04-25, 12:01   #35
User133
 
Apr 2020

22×5 Posts
Default

thank you very much
User133 is offline   Reply With Quote
Old 2020-04-25, 12:48   #36
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

25·233 Posts
Default

Quote:
Originally Posted by LaurV View Post
This is for homework thread. It may be moved there later, as it has nothing to do with num theory.
Wrong again, as usual. Do us all a favor. Stop answering questions until you actually
learn some number theory.

Quote:

Anyhow, we are trying to help.
Trying is indeed the correct word. You are not succeeding.

Quote:

2. Do you need to learn pari functions? Or to use pari predefined functions to get the answer? If so, look at help for vecsum() and divisors() functions.
s=0; for(n=3*10^8,4*10^8,z=vecsum(divisors(n)); if(s<z, s=z; print([n,s]))) will give you the answer in 3-5 minutes
The correct way to do this will give the result in seconds and does not use the
divisors function. It does require using the computer in your head along with some
math. It also requires knowing how one computes the sum of divisors
from the prime factorization *without* knowing/computing the non-prime divisors.
Your "brute force/try all possibilities" is simple minded and horribly inefficient.

Last fiddled with by R.D. Silverman on 2020-04-25 at 12:50 Reason: cleanup
R.D. Silverman is offline   Reply With Quote
Old 2020-04-25, 13:25   #37
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

164408 Posts
Default

Quote:
Originally Posted by R.D. Silverman View Post
Wrong again, as usual. Do us all a favor. Stop answering questions until you actually
learn some number theory.



Trying is indeed the correct word. You are not succeeding.



The correct way to do this will give the result in seconds and does not use the
divisors function. It does require using the computer in your head along with some
math. It also requires knowing how one computes the sum of divisors
from the prime factorization *without* knowing/computing the non-prime divisors.
Your "brute force/try all possibilities" is simple minded and horribly inefficient.

Follow on. Give an estimate for the largest prime factor of the
answer without doing any computations at all. You may express your
answer as a function of the input parameters. Justify your answer.

Hint: psi(x) ~ x.

Explain why this question is relevant to the original question.

Last fiddled with by R.D. Silverman on 2020-04-25 at 13:40 Reason: add hint.
R.D. Silverman is offline   Reply With Quote
Old 2020-04-25, 15:35   #38
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

3×2,861 Posts
Default

Quote:
Originally Posted by R.D. Silverman View Post
It also requires knowing how one computes the sum of divisors from the prime factorization *without* knowing/computing the non-prime divisors.
Your "brute force/try all possibilities" is simple minded and horribly inefficient.
Man, don't you think that I worked a bit too much for those aliquot sequences to be unaware of the fact that \sigma(n)=\prod_{i=1}^d{\frac{p_i^{a_i+1}-1}{p_i-1} ? I may not have your math knowledge, but at least see my posts in the past, related to the subject. Do you see what the OP is posting? He's not learning number theory here, but how to set and reset variables in computer... He has to learn not to compute the same thing ever and ever, etc. and you want to teach him academic stuff.

We however agree about the last sentence. This way of computing is "horribly inefficient", regardless of the fact that you use the divisors, like I did, or the prime factors, like you suggested . Because the factorization is the hard part, not combining them in products, or sums...

----
p.s. is matjax broken again? it seems \TeX tag works, but matjax not, or at least, not in previews...

Last fiddled with by LaurV on 2020-04-25 at 15:45
LaurV is offline   Reply With Quote
Old 2020-04-25, 17:51   #39
R.D. Silverman
 
R.D. Silverman's Avatar
 
Nov 2003

11101001000002 Posts
Default

Quote:
Originally Posted by LaurV View Post
Man, don't you think that I worked a bit too much for those aliquot sequences to be unaware of the fact that \sigma(n)=\prod_{i=1}^d{\frac{p_i^{a_i+1}-1}{p_i-1} ? I may not have your math knowledge, but at least see my posts in the past, related to the subject. Do you see what the OP is posting? He's not learning number theory here, but how to set and reset variables in computer...
I can't see any professor giving this computational problem as an exercise in
how to write code or in learning Pari syntax. That could be done by RTFM.

It has to be an exercise in how to derive and code the best method. Code syntax
is too trivial. If it is a weird exercise in Pari syntax then the correct answer is indeed
RTFM.

Quote:




He has to learn not to compute the same thing ever and ever, etc. and you want to teach him academic stuff.
It *IS* "academic stuff".

Quote:
We however agree about the last sentence. This way of computing is "horribly inefficient", regardless of the fact that you use the divisors, like I did, or the prime factors, like you suggested .
The point of such an exercise is to determine that using factorization is the wrong way
to do it. It is in the "Computer Science and CNT" group and not the "programming"
group. It is also learning that this problem requires finding a solution that has
many (necessarily small) prime factors.

Also, the OP said he needed help with a task using GP and not help learning Pari/GP.
Note that there is a subgroup devoted to Pari itself if the point of the exercise were to learn it.

Last fiddled with by R.D. Silverman on 2020-04-25 at 17:56
R.D. Silverman is offline   Reply With Quote
Old 2020-04-25, 18:02   #40
LaurV
Romulan Interpreter
 
LaurV's Avatar
 
Jun 2011
Thailand

3×2,861 Posts
Default

Quote:
Originally Posted by R.D. Silverman View Post
I can't see any professor giving this computational problem as an exercise in
how to write code or in learning Pari syntax. That could be done by RTFM.

It has to be an exercise in how to derive and code the best method. Code syntax
is too trivial.
Agree. We are in violent agreement here. You quoted my post #8, points 1 and 2, and almost all the text, but omitted point 3, which mainly said the same thing as you say:

Quote:
Originally Posted by LaurV View Post
3. Do you need to implement this by yourself? I don't see the teacher who will give a homework that requires using only one line of code, and only 2 functions. When I was in Uni, like >30 years ago, we had to sort a list (in Pascal, or C), and learn different ways to sort it, in the process (bubblesort, flagsort, quicksort, etc). One colleague presented program with one line of code only, he was calling the "sort" function in dBase (at the time) and was very upset when the professor, laughing his lungs out, marked him "failed". He even didn't understand that the requirement was not to sort a stupid, 20 items list, which could be easily sorted by pencil, but to learn the algorithms. True story!

Last fiddled with by LaurV on 2020-04-25 at 18:02
LaurV is offline   Reply With Quote
Old 2020-06-14, 16:25   #41
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

32·61 Posts
Default

There is a program for the computation of Hardy-Littlewood constants for quadratic polynomials in A221712, but it doesn't seem to work:

Code:
/* Auxiliary functions. */
ZetaDN(P, s) = zeta(s) * prod(j = 1, #P, 1 - P[j]^(-s));

LchiN(L, Ebad, s) =
{ my([P, E] = Ebad);
  lfun(L, s) * prod(j = 1, #P, subst(E[j], 'x, P[j]^(-s)));
}

LchiNinit(D, P) =
{ my(Ebad = [], Pbad = []);
  for (j = 1, #P,
    my(p = P[j], s = kronecker(D, p));
    if (s, Ebad = concat(Ebad, 1 - s*'x);
           Pbad = concat(Pbad, p)));
  return ([Pbad, Ebad]);
}

Oddpart(n) = n >> valuation(n,2);

/* The real work; D is fundamental. */
HLW2(D, N) =
{ my(B = getlocalbitprec(), lim, S1, S2, L, P, v, Ebad);

  localbitprec(32); lim = ceil(B*log(2)/log(N/2));
  localbitprec(B + lim + exponent(lim));
  L = lfuninit(D, [1/2, lim, 0]);
  v = vector(lim);
  forfactored(X = 1, lim,
    my([n] = X, S = 0); \\ FIXME: loop over odd divisors?
    fordivfactored(X, Y,
      my([d] = Y);
      if (d % 2, S += moebius(Y) << (n/d)));
    v[n] = S / (2*n);
  );
  P = setunion(factor(abs(D))[,1]~, primes([2, N]));
  Ebad = LchiNinit(D, P);
  S1 = sum(n = 1, lim, v[n] * log(LchiN(L, Ebad, n)));
  S2 = sum(n = 2, lim, (v[n] - if (n%2 == 0, v[n/2]))
                       * log(ZetaDN(P, n)));
  return (S1 + S2);
}

/* Compute the Hardy-Littlewood constant of aX^2+bX+c. */
HardyLittlewood2(A, N = 50) =
{ my(D = poldisc(A), S, P);

  if (poldegree(A) != 2, error("polynomial of degree != 2"));
  my([a, b, c] = Vec(A));
  if (issquare(D) || gcd([2 * a, a + b, c]) > 1, return (0));
  N = max(N, 3);
  /* Take care of the prime p = 2. */
  S = if ((a + b) % 2, 1., 2.);
  /* Take care of odd primes dividing a. */
  P = factor(Oddpart(a))[,1];
  for (j = 1, #P,
    my(p = P[j]);
    S *= if (b % p, (p - 1) / (p - 2), p / (p - 1))
  );
  /* Take care of odd primes dividing the index f. */
  my([D0, f] = coredisc(D, 1));
  P = factor(Oddpart(f))[,1];
  S /= prod(j = 1, #P,
    my(p = P[j]);
    1 - kronecker(D0, p) / (p - 1);
  );
  /* Take care of the primes p <= N. */
  S *= prodeuler(p = 3, N, 1 - kronecker(D0, p) / (p - 1));
  /* Do the real work */
  return (S * exp(-HLW2(D0, N)));
}
Code:
  ***   unexpected character: LchiN(L,Ebad,s)=my([P,E]=Ebad);lfun(L,s)*pro
                                                       ^--------------------
Some help here, please?
mart_r is offline   Reply With Quote
Old 2020-06-14, 17:20   #42
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

328110 Posts
Default

The code translates here. Is you error at runtime? Maybe your should consider upgrading Pari/GP to a later version.

Code:
                   GP/PARI CALCULATOR Version 2.12.0 (alpha)
          amd64 running linux (x86-64/GMP-6.1.2 kernel) 64-bit version
                 compiled: Sep 27 2019, gcc version 9.1.0 (GCC)
                            threading engine: single
                 (readline v7.0 enabled, extended help enabled)
paulunderwood is online now   Reply With Quote
Old 2020-06-14, 22:15   #43
mart_r
 
mart_r's Avatar
 
Dec 2008
you know...around...

32×61 Posts
Default

Okay, I got the newest version and the program can now be read in.

But it is still not working as desired, I'm afraid:
Code:
(23:53) gp > HardyLittlewood2(x^2+x+41)
 ***  at top-level: HardyLittlewood2(x^2+x+41)
 ***                ^--------------------------
 ***  in function HardyLittlewood2: ...,p)/(p-1));return(S*exp(-HLW2(D0,N)))
 ***                                                            ^------------
 ***  in function HLW2: my(B=getlocalbitprec(),lim,S1,S2,L,P,v,EBad);1
 ***                         ^-----------------------------------------
 ***  not a function in function call
Has it something to do with "setting the required precision" (cf. comment in the OEIS entry) - and how?
mart_r is offline   Reply With Quote
Old 2020-06-15, 00:56   #44
paulunderwood
 
paulunderwood's Avatar
 
Sep 2002
Database er0rr

17·193 Posts
Default

On my version:

Code:
? HardyLittlewood2(x^2+x+41)
6.6395463549428433306471137152997759330
What version are you using and does it have the function getlocalbitprec()?

Does ?getlocalbitprec give you an entry?
paulunderwood is online now   Reply With Quote
Reply

Thread Tools


Similar Threads
Thread Thread Starter Forum Replies Last Post
LLL in GP/Pari paul0 Programming 2 2015-11-17 13:04
PARI vs GAP skan Miscellaneous Math 0 2012-12-16 00:13
pari devarajkandadai Programming 21 2012-08-31 18:08
64-bit Pari? CRGreathouse Software 2 2009-03-13 04:22

All times are UTC. The time now is 05:03.

Sat Jul 11 05:03:57 UTC 2020 up 108 days, 2:37, 0 users, load averages: 1.39, 1.33, 1.29

Powered by vBulletin® Version 3.8.11
Copyright ©2000 - 2020, Jelsoft Enterprises Ltd.

This forum has received and complied with 0 (zero) government requests for information.

Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.2 or any later version published by the Free Software Foundation.
A copy of the license is included in the FAQ.