Bzoj 1101: [Poi2007]zap

Description

FGD is cracking a cipher, and he needs to answer a number of similar questions: for a given integer, a, B, and D, how many positive integers are to X, Y, X<=a,y<=b, and gcd (x, y) =d. As a classmate of FGD, FGD hopes to get your help.

Data range: Multiple groups of queries, 1<=n<= 50000, each group of 1<=d<=a,b<=50000solution

Test instructions: 1<=x<=a,1<=y<=b, satisfies the (x, y) logarithm of gcd (x, y) = d

In order to learn Möbius inversion from the Huang to find the template problem

The most immediate idea is of course

For I 1 to a

For J 1 to B

if (Gcd==d) count++;

We write it in the form of summation.

Make

There is a formula in Möbius inversion

This formula can be arbitrarily pushed with a two-item theorem.

Observe the above two formulas

Found that the two is exactly the perfect fit when and only when the number of Gcd==d is added to the limit conditions

So we consider replacing the GCD (I,J) with a binary, and changing the original from N to gcd (I,J)

At this point, we get a new formula.

So we change the way we enumerate

Enumerate gcd possible factor d at the outermost layer, and then take a multiple of two d from a ', B '

The number of such schemes is ⌊a '/d⌋⌊b '/d⌋, and the nature of the conversion enumeration is the same as the above

Then get

Then μ also has the method of linear sieve, can learn

The problem is almost solved, the complexity is about N*a, N is the number of data groups

But looking at the data 50000*50000 is not going to make it.

Think again, for ⌊a '/d⌋, must be a number of consecutive fixed value, that is, for a period of time d,⌊a '/D⌋ value will not change, ⌊b '/D⌋ the same

You can then record the μ prefix and then block it to ask for

But ⌊a '/d⌋ theoretically have √a ' value, ⌊b '/D⌋ Likewise, ⌊a '/d⌋⌊b '/D⌋ value is not √a ' *b ', complexity doesn't seem to be optimized

In fact, because the divisor d of both is the same at the same time, so at a certain time in D ⌊a '/D⌋ only one ⌊b '/d⌋ value only √a ' +√b '

So the complexity becomes n*√a '

Problem solving

1#include <map>2#include <cmath>3#include <ctime>4#include <queue>5#include <stack>6#include <cstdio>7#include <climits>8#include <iomanip>9#include <cstring>Ten#include <cstdlib> One#include <iostream> A#include <algorithm> -  - #defineMAXP 50000 the #defineMAXN 50000+5 - #defineSet (a B) memset (A, (b), sizeof (a)) - #defineFR (i,a,b) for (ll i= (a), _end_= (b); i<=_end_;i++) - #defineRF (I,b,a) for (ll i= (a), _end_= (b); i>=_end_;i--) + #defineFe (I,A,B) for (int i=first[(b)],_end_= (a); i!=_end_;i=s[i].next) - #defineFEC (I,A,B) for (int &i=cur[(b)],_end_= (a); i!=_end_;i=s[i].next) +  A using namespacestd; at  -typedefLong Longll; -  - intprime[maxn],pri[maxn],tot=0; - intSUM[MAXN],MIU[MAXN]; - intans; in intn,m; -  to voidRead () + { - #ifndef Online_judge theFreopen ("1101.in","R", stdin); *Freopen ("1101.out","W", stdout); $ #endifPanax Notoginseng   //cin >> N; -scanf"%d",&n); the } +  A voidWrite () the {} +  - voidprint () $ { $   //cout << ans << endl; -printf"%d\n", ans); - } the  - voidMobius ()Wuyi { themiu[1]=1; -Fr (I,2, MAXP) { Wu     if(!prime[i]) pri[++tot]=i,miu[i]=-1; -     intj=1; About      while(J<=tot && pri[j]*i<=Maxp) { $prime[pri[j]*i]=1; -       if(i%pri[j]==0 ){ -miu[pri[j]*i]=0; -      　　  Break; A       } +miu[pri[j]*i]=-Miu[i]; theJ + +; -     } $   } theFr (I,1, MAXP) thesum[i]=sum[i-1]+Miu[i]; the } the  - intCalintXinty) in { the   intres=0, pos,i=1; the   if(x>y) Swap (x, y); About    while(i<=x) { thePos=min (x/(x/i), y/(y/i)); theres+= (sum[pos]-sum[i-1]) * (x/i) * (y/i); thei=pos+1; +   } -   returnRes; the }Bayi  the voidWork () the { - Mobius (); -Fr (I,1, N) { the     inta,b,d; the     //cin >> a >> b >> D; thescanf"%d%d%d",&a,&b,&d); theAns=cal (a/d,b/d); - print (); the   } the } the 94 intMain () the { the read (); the Work ();98 write (); About   return 0; -}

Bzoj 1101: [Poi2007]zap